Computational fluid dynamics simulation of the upper airway response to large incisor retraction in adult class I bimaxillary protrusion patients

The changes of the upper airway after large retraction of the incisors in adult class I bimaxillary protrusion patients were assessed mainly focused on the anatomic variation and ignored the functional changes. This study aimed to investigate the changes of the upper airway in adult class I bimaxillary protrusion patients after extraction treatment using the functional images based on computational fluid dynamics (CFD). CFD was implemented after 3D reconstruction based on the CBCT of 30 patients who have completed extraction treatment. After treatment, pressure drop in the minimum area, oropharynx, and hypopharynx increased significantly. The minimum pressure and the maximum velocity mainly located in the hypopharynx in pre-treatment while they mostly occured in the oropharynx after treatment. Statistically significant correlation between pressure drop and anatomic parameters, pressure drop and treatment outcomes was found. No statistical significance changes in pressure drop and volume of nasopharynx was found. This study suggested that the risk of pharyngeal collapsing become higher after extraction treatment with maximum anchorage in bimaxillary protrusion adult patients. Those adverse changes should be taken into consideration especially for high-risk patients to avoid undesired weakening of the respiratory function in clinical treatment.

drops in the pharynx were more closely associated with the presence and severity of OSAS than anatomical endpoints, which may be a better tool to predict the risk of the upper airway collapsing 10 .
Computational Fluid Dynamics (CFD) is a validated method to accurately compute aerodynamic flow characteristics of the upper airway [11][12][13][14][15] , which allows describing regional flow and pressure profiles in the upper airway. Literature had reported that a CFD model is more sensitive or accurate than anatomical parameters alone in evaluating the effects of airway restriction in OSAS 10 . Therefore, CFD simulation together with 3D upper airway model makes it possible to provide a more comprehensive understanding of the variety in pharynx after extraction treatment.
The aim of this study is to evaluate the pharyngeal airflow response to large incisor retraction in class I adult bimaxillary protrusion patients using CFD simulation based on the patient-specific CBCT dataset.

Subjects. The study was reviewed and approved by School of Stomatology Shandong University Research
Ethic Board (protocol number 20160103). All the written informed consents were received from the patients, and the study was conducted according to the tenets of the Declaration of Helsinki for research involving human subjects. The methods were carried out in accordance with the approved guidelines of SCIENTIFIC REPORTS. The sample of this study consisted 30 (19 females, 11 males, age 25.87 ± 0.78 years, BMI 20.56 ± 1.48 kg/m 2 ) class I bimaxillary dentoalveolar protrusion patients, who were treated by extraction of four first premolars using maximum anchorage. The exclusion criteria included: (1) patients have any other severe craniofacial anomalies except class I bimaxillary dentoalveolar protrusion, (2) age under 18 years old, (3) a history of serious diseases, or acquiring serious diseases during treament, such as serious cardiovascular disease, respiratory system disease, serious craniofacial deformity, severe nasal diseases and relative surgical history, (4) weight increased apparently into obese (BMI > 30 kg/m 2 ) after treatment. The patients were informed consent, including the potential risk associated with CT radiation and miniscrew methodologies. Oriental pre-adjusted appliance KOSAKA slot brackets (OPA-K, Tomy; Fukushima-ken, Japan) were used, and miniscrews were implanted as anchorage for the retraction and intrusion of the teeth. A force of 150 g per side of elastic chains was applied from the miniscrew to the upper crimpable hook to retract and intrude the upper anterior tooth. The patients were examined at 1-month intervals until the retraction of the upper and lower anterior teeth completed.

Reconstruction of the model and mesh generation.
The three-dimensional models used in the calculation were developed from pre-and post-treatment CBCT data from patients selected. All CBCT scans were implemented with each patient awake in the position of the Frankfort horizontal plane parallel to the floor using the same CBCT scanner (KaVo Dental GmbH, Bismarckring, Germany. scan time: 8.9 seconds; slice thickness: 0.4 mm, 120 kV, 5 mA) by the same operator. Patients were guided to close their mouths with the maximum intercuspation and the tongues touching their hard palates during the scanning. The pre-and post-treatment CBCT scans of patients were collected as T1 and T2 data. Subsequently, the dataset was exported in DICOM (Digital Imaging and Communications in Medicine) file format, and then was read into MIMICS16.0 (Materialism's Interactive Medical Image Control System) software. After that, the models were reorientated in three planes: in the coronal view, the most inferior points on the infraorbital margin (orbitale) lie on the same horizontal plane; in the sagittal view, the models were reorientated to make the Frankfort plane horizontal; in the axial view, the models were reorientated to make the line through the crista galli and the midpoint on the anterior margin of foramen magnum (basion) vertical 16 . The upper airway of interest was segmented by setting the threshold between − 1024 Hounsfield Units (HU) and − 480 HU, and the 3D anatomically accurate patient-specific models were reconstructed. An appropriate smoothing algorithm was used to transform the 3D model into a smooth one without the loss of the patient-specific characters in the shape of the upper airway (Fig. 1). After that, the Unstructured tetrahedral volume mesh was generated in ANSYS ICEM CFD (ANSYS 16.0) using hybrid mesh scheme (Fig. 2). A grid convergence analysis was performed by repeating the solution with five different element size meshes (Grid 1, Grid 2, Grid 3, Grid 4 and Grid5) to establish grid independence solutions. Similar element size was used in the pre-reatment and post-treatment models for reliable comparisons of the results. Changes in the average pressure and average velocity in a selected plane were used as convergence criteria. when the changes in those criteria were less than 1%, an acceptable level of grid-independence was achieved. As shown in Fig. 3, changes of the average pressure and average velocity between Grid 3 and Grid 4 were 0.32%, 0.93% for pretreatment model, and 0.38%, 0.02% for post-treatment model. It was suggested that changes in those criteria become practically negligible if Grid3 is used. Thus, all the data presented in the current work were from a simulation with Grid3, resulting in a computational grid of pre-and post-treatment typically consisting of 1.10 ± 0.05 million and 1.01 ± 0.09 million elements, respectively. The maximum skewness of the grids was 0.87. Then the 3D models were imported into the Reynolds Average Navier Stokes CFD solver (Fluent 16.0, Fluent Inc.).

Computational fluid dynamics. Flow simulations were performed using the Reynolds Average Navier
Stokes CFD solver (Fluent 16.0, Fluent Inc.). Inspiratory upper airway flow was modeled, as inspiration is associated with negative pressures causing airway collapse in OSAS, greater airflow pressure and velocity gradients. Steady laminar flow field was applied in this study. Previous studies showed that fully developed turbulence would not be reached until a flow rate of 30 L/min 17 . The flow rate in the current study was less than 24 L/min, therefore  the airflow was presumed as incompressible and laminar. Wootton et al. 11 have validated steady model results with unsteady flow simulation, and the pressure distribution in the steady model at peak flow was very similar to the unsteady model, the most changes in pressure was less than 2 Pa. To reduce the computational resources needed, steady state condition was used in this study. The magnitude of the peak flow rate was calculated by the formula V = Q/A. Q was the volume flow rate of the model, and A was the area of the bottom of the hypopharynx. In this study, the respiratory tidal volume was 600 mL, and the respiratory cycle was T = 4 s, an approximate sinusoidal variation in the airflow rate with time was observed in the rhinomanometer test, and the ratio of the inspiratory time and expiratory time was 1:1 (Fig. 4), , ω = π /2 is the angular frequency, Q tidal volum = 600 ml. For inspiration process, an ambient static pressure condition was set at the nostril (the flow outlet), the negative value of velocity at the peak flow rate was set at the bottom of the hypopharynx (the flow inlet), the wall surface of the geometry was assumed to be non-slip.
The effects of some factors such as the temperature, humidity and vibrissae on the fluid were neglected in the simulation, while the gravity was included. The calculation residual was set as 10 −4 and the iteration numbers were 1000 steps when all residuals fell below 0.2%. Second-order discretization schemes were used with pressure-velocity coupling solved through the SIMPLE scheme.
Outcome parameters. Two planes: a line passing from anterior nasal spine to posterior nasal spine and a line passing from the anterosuperior edge of the fourth cervical vertebra (C4) to menton, which have been proved to be reliable and reproducible landmarks for defining upper airway margins 16 were uesd to divided the whole pharynx into three parts (nasopharynx, oropharynx, and hypopharynx). In this study, volumes of the whole pharynx (V) and each part of the pharynx: nasopharynx (V na ), oropharynx (V or ), and hypopharynx (V hy ) was calculated. The minimum cross-sectional area (A min ) and its location was defined. Parameters above were used to assess the size changes of the pharynx. In addition, the mean cross-sectional area (A mean ) was computed by ratio of the pharyngeal volume and the length of the pharynx (V/L). A min /A mean was also obtained to explore the uniform of the cross-sectional changes. Upper central incisor retraction amount at edge in the horizontal direction (ICE) were achieved through registration of the pre-and post-treatment models, and the process for registration was presented in details in the work of Chen Y et al. 7 .
Once the simulation was completed, all of the cross sections from choana to the bottom of hypopharynx were selected per millimeter along Z-axis for parameters calculation. Area-weighted average pressure and velocity in each cross-section was extracted to evaluate the location of the minimum area-weighted average pressure and the maximum velocity in the pharynx. The area-averaged pressure drop from choana to the minimum area (Δ P max ) was defined to quantify the effect of the pharyngeal airway narrowing on the flow field. Airflow pressure drops of nasopharynx (Δ P na ), oropharynx (Δ P or ), and laryngopharynx (Δ P hy ) were counted by Δ P = P max − P min . The resistance of the airway (R) was achieved by the formula: R = Q/Δ P (Q was pharyngeal inlet volume flow rate). Because the flow rate was constant in our study, the Δ P reflected the changes of the resistance for each pharyngeal part.
All the parameters were measured twice by the same researcher, and the average value was applied to the study.
Statistical analysis. All statistical analyses were performed using the SPSS 17.0 software package. For each variable themean and the standard deviation (SD) were calculated. The differences between the flow dynamics and anatomical parameters of the upper airway between T1 and T2 were analyzed by paired-t test, and the relationships between the ΔP and ICE were assessed using Pearson correlation coefficient. To assess the measurement errors, all the variables were measured twice with the manual method by the same operator with a 2-week interval. The error of method was calculated as follows: = Σd n ME /2 2 (where d indicates deviations between the two measurements, and n indicates number of paired objects). Intraexaminer reliability was calculated by intraclass correlation coefficient (ICC) for the measurements obtained by the examiner at both times. For all the analyses, P < 0.05 was considered significant. Stepdown Bonferroni multiple testing was used to correct the P-values.

Results
No significant difference between the first and repeated measurements were observed, which proved the reliability between the first and second measurements (ICC from0.92 to 0.97). The maximum ME of pharyngeal volume and area measurements were 0.15 mm 3 and 0.31 mm 2 , respectively, while the maximum ME of pressure measurements was 0.47 Pa in the numerical simulation of pharyngeal airflow, and 0.24 mm of the upper central incisor retraction amount at edge in the horizontal direction.
A statistically significant decrease between the pre and post-treatment was observed in the volume and the value of the A min and A min /A mean (Table 1) of the oropharynx and hypopharynx, while the difference in the nasopharynx between pre and post-treatment was not significant. The location of the minimum of the cross-section moved upward from the hypopharynx to the oropharynx in the post-treatment patients compared to the pre-treatment patients. Upper central incisor retraction amount at edge in the horizontal direction (ICE) was 6.84 ± 1.68 mm after the extraction treatment. Figures 5 and 6 showed the pressure and velocity distribution in upper airway for a representative patient. In the pre-treatment subjects, the pressure had no abrupt changes from the choana to the hypopharynx, and the minimum pressure and maximum velocity occurred in the hypopharynx. After extraction treatment, the pressure decreased from the choana to the hypopharynx, and the minimum pressure and maximum velocity occurred in the oropharynx. In addition, the mainstream of airflow with high velocity flowed along the posterior wall of the pharynx, while the flow velocity was lower in the anterior wall. The local swirls occurred mostly in the places with abrupt changes in morphology. The value of the Δ P max , Or-Δ P, Hy-Δ P increased significantly after the extraction treatment ( Table 1). The pressure drop in nasopharynx was statistically insignificant.
As shown in Table 2, a significant negative correlation was found between Or-Δ P, Hy-Δ P and the relevant volume. There was a significant positive correlation between Or-Δ P, Hy-Δ P, Δ P max and A min , A min /A mean . In Table 3, there was a significant positive correlation between Or-Δ P, Hy-Δ P, Δ P max and ICE.

Discussion
The present study was designed to assess the airflow characteristics variation of the upper airway in adult class I bimaxillary protrusion patients using functional imaging based on CFD after extraction treatment. The pressure drop from choana to the minimum cross-section, and the resistance in oropharynx and hypopharynx increased significantly in post-treatment subjects compared to the pre-treatment subjects. These changes were correlative with the relevant anatomical changes and the treatment outcomes.
Previous studies have reported a close relationship between pharyngeal narrowing and OSAHS [18][19][20][21] . Δ P max and A min have been proved to be positively correlated with the severity of the OSAS 10,20,21 , which are parameters to quantify clinically significant flow-and pressure-related effects of anatomical narrowing. These two variables may serve as a guide for clinical treatment. In this study, the correlation between the A min and Δ P max was confirmed. After extraction treatment, the pharyngeal minimum cross section was diminished significantly, which induced additional negative internal pressure and higher airflow velocity due to the "Bernoulli effect" near the smallest cross-sectional area of the upper airway. Bigger pressure gradient between the internal pressure and the external pressure from surrounding tissues made the upper airway a more collapsible vessel. Due to the reduction of the   smallest cross-sectional area, pharyngeal luminal pressure decreased more due to an increase in kinetic energy and a dissipative energy loss upstream during inspiration. Consequently, the Δ P max raised significantly, which may suggest higher probability of pharyngeal obstruction. The correlation between Δ P max and A min /A mean was stronger than that between A min and Δ P max , suggesting that the relative narrowing of the upper airway had a stronger effect on the collapsibility than the absolute narrowing 15 .
In pre-treatment subjects, the airflow passed through the pharynx smoothly without sharp changes. The pressure reduced uniformly from nasopharynx to hypopharynx, and the minimum pressure existed in the hypopharynx. After the extraction treatment, the minimum pressure moved up along the pharynx. This occurred mostly in the root of tongue, next in the back of the soft palate where the pharyngeal obstruction occurred commonly in OSAS patients [22][23][24] . The resistance of both the oropharynx and hypopharynx raised, which was significantly correlated with the relevant volumetric decrease and the amount of the retraction of the incisors. The retraction of the incisor diminished the length of the oral and enforced the tongue backward, therefore the upper airway diminished in size, especially in the base of the tongue, next in the back of the soft palate 9 . These findings were concordant with Chen Y et al.' studies 6,7,9 , but not with Al Maaitah et al.' fingings 8,24,25 , possibly due to the difference in the selected subjects and growth.
It has been demonstrated that the fast changing period of the upper airway is under 20 years old 26 . Since all of the patients in the present study were older than 20, the growth potential impact might be ignored.
The limitation of this study is that the CFD model assumed the wall of the upper airway as hard wall, while the upper airway was surrounded by soft tissues, such as tongue and soft palate in the anterior wall. Therefore, a 3D fluid-structure interaction simulation of the upper airway, which considers the interaction between the airflow and around soft tissues, is needed to assess the functional changes of the upper airway [27][28][29] . Another limitation is that the OSAS occurs during sleeping, when the upper airway dilator muscles andrespiratory pump muscles are less activated thus cannot adjust the shape and posture of the soft tissues. However, the CBCT has been taken while the patient is awake, and the size and morphology of the upper airway may be different between asleep and awake 30,31 . In the future, a FSI model of the sleeping patient should be established to assess the functional changes of the upper airway after the extraction treatment.

Conclusion
The oropharynx and hypopharynx increase in flow resistance and the risk of collapsing of the pharynx become higher after the extraction treatment with the maximum anchorage in bimaxillary protrusion adult patients. Those adverse changes in upper airway owing to the extraction should be taken into consideration to avoid weakening of the respiratory function.  Table 3. The correlation between pressure drop in minimum cross section, oropharynx, hypopharynx and the amount of incisor retraction.