Analogue simulation of pharyngeal airflow response to Twin Block treatment in growing patients with Class II1 and mandibular retrognathia

The flow dynamics of respiratory airflow is the basic factor that influences the ventilation function of the upper airway. This research aimed to investigate the pharyngeal flow field characteristics after Twin Block (TB) treatment in growing patients with Class II1 and mandibular retrognathia by computation fluid dynamics (CFD) simulation. Cone beam computed tomography (CBCT) scans of patients who have completed TB treatment (n = 30) and about to accept TB treatment (n = 30) were reconstructed. After CFD simulation, correlations between the pharyngeal pressure drop and morphological parameters were further analyzed. During inspiration, we found that the pressure minimum occurred in the hypopharynx, while the maximum pressure drop and velocity was located in the oropharynx. After TB treatment, the oropharynx and hypopharynx showed significant differences in airflow features, and the most obvious change was observed in the oropharynx. A significant correlation was discovered between the change amount of oropharyngeal pressure drop and volume (r = 0.694, p = 0.001), mean cross-sectional area (r = 0.859, p = 0.000), and ratio of the minimum and mean cross-sectional area (r = 0.898, p = 0.000) of the oropharynx. Our research suggested that the pharyngeal airflow characteristics response positively to mandibular advancement with the enlargement in volume, cross-sectional area and more uniform oropharyngeal area distribution.

Many studies evaluated the pharyngeal flow dynamics response followed by the morphological changes in adult patients with OSAS [11][12][13] . They reported that the upper airway resistance is inversely correlated with the pharyngeal dimension. W. Vos. et al. measured the anatomical properties and simulated the airway flow of the upper airway in OSAS patients, and confirmed that there exists a relationship between the morphological parameters of the upper airway and apnea hypopnea index (AHI) in these patients 14 . However, there have been no studies related to pharyngeal flow field characteristics due to morphological changes after mandibular advancement in growing patients. Since the upper airway of children is anatomically different from that of the adults 15 , simulation of the upper airway response to mandibular advancement in growing patients with Class II 1 and mandibular retrognathia needs to be established.
Computational Fluid Dynamics (CFD) is an appropriate method in analogue simulating the flow field characteristics of the upper airway. It could accurately calculate the functional parameters such as pressure drop and flow resistance before and after treatment based on the patient specific geometrical models of the upper airway derived from 3D imaging technology [16][17][18] . Thereby provide a better understanding in the effect of the morphological and respiratory airflow changes on the pharyngeal function. Therefore, the objective of this study is to evaluate the pharyngeal airflow response to TB treatment in growing patients with class II 1 and mandibular retrognathia by means of CFD analogue simulation.

Methods
The study was critically reviewed and approved by the Ethical Committee of Shandong University Dental School. All the written informed consents were received from parents or legal guardians of 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.
For the subjects involved in our research are growing patients, their personal growth development have a negligible impact on the upper airway during treatment period 19 . For this reason, a control group is needed. Body mass index (BMI) was calculated for each patient to quantify the development status. Thirty growing patients (13 boys and 17 girls, mean age 11.57 ± 0.94 years, BMI 17.86 ± 1.25) who were diagnosed as class II 1 malocclusion with mandibular retrognathia and have finished TB treatment were randomly selected as TB group. The treatment period of TB group was 13.67 ± 1.51 months. The control group (13 boys and 17 girls, mean age 11.72 ± 0.86 years, BMI 18.14 ± 0.97) included the patients with the same diagnosis and matched well with the patients of TB group through age, gender, BMI and growth pattern. By this method, factors that influence the upper airway anatomy and morphology other than appliance therapy could be excluded. All the patients had no other potential airway abnormalities.
Image acquisition and 3D model reconstruction. All the CBCT scans were performed using the CBCT scanner (KaVo Dental GmbH, Bismarckring, Germany) on each wake patient with the Frankfort horizontal plane parallel to the floor. During scanning, patients were required to close their mouth with maximum intercuspation and the tongue was in contact with anterior palate. The scanning area was from the basis cranii level to the third cervical vertebra level (scan time: 8.9 seconds; slice thickness: 0.4 mm, 120 kV, 5 mA). The CBCT data were saved as DICOM (Digital Imaging and Communications in Medicine) format. The pre-and post-treatment CBCT scans of TB group were collected as T1 and T2 data, while the scans of control group before function therapy were collected as control data.
All 3D models of the upper airway in TB and control groups were reconstructed using MIMICS 16.0 (Materialism's Interactive Medical Image Control System) software. The upper airway was highlighted by setting the threshold between − 1024 Hounsfield Units (HU) and − 480 HU. The paranasal sinuses were erased manually for there were almost no airflow in the cavity. After image segmenting and region growing, the patient-specific upper airway 3D geometry was obtained.
The upper airway model consists of nasal cavity and pharynx. After identifying the basis cranii, the PNS point (posterior nasal spine), the superior border of the epiglottis and the C 3 point (the most anterior inferior point of the third cervical vertebra) in the midsagittal plane, the pharynx is further divided into three parts: the nasopharynx, oropharynx, and hypopharynx by the corresponding cross-sectional slices (Fig. 1). All 3D models were exported as stereolithography (STL) files.

CFD mesh generation.
After 3D reconstruction, the upper airway geometry of STL format was loaded into ANSYS ICEM CFD (ANSYS 12.0) for model trimming and mesh generation. The cross-sections of nostrils were selected as the inlet of the upper airway, while the outlet was located at the cross-section in C 3 level. The whole pharyngeal cavity was defined as body of the geometry. Then, an unstructured tetrahedral volume mesh was generated in the 3D upper airway model to form a computational domain of fluid field (Fig. 2). The maximum grid edge length was set as 1.0 mm while the minimum was 0.5 mm, which resulted in a computational grid of pre-, post-treatment and the control models typically consisting of 8.81 ± 0.09 billion, 12.04 ± 0.10 billion and 9.63 ± 0.08 billion computational cells, respectively. The maximum scale factor of the grids was 0.8.
Airflow simulation. The 3D mesh was exported and read into a Reynolds Average Navier Stokes CFD solver (Fluent 14.0, Fluent Inc.) for numerical simulation. The air was assumed to be a newtonian and homogeneous fluid. Considering the Mach number of airflow during eupnea is smaller than 0.3, the flow field within the upper airway was incompressible. Pressure boundary conditions are applied for the upper airway constriction have no impact on the static pressure difference between the lungs and atmosphere. A standard atmospheric pressure (gauge pressure was 0 Pa) was set for the inlet and in the outlet, a pressure of − 20 Pa was applied. The wall surface of the geometry was assumed to be stationary and non-slip. The effects of the temperature, humidity and vibrissae factors on the fluid were neglected in the simulation, while the gravity factor was included.
Scientific RepoRts | 6:26012 | DOI: 10.1038/srep26012 A maximum inlet volume flow rate 166 ml/s (10 L/min) was set for the inhaled airflow in an awakened state 18 . The maximum and minimum values of pressure and velocity were artificial limited to unified the reference legends. Laminar approach was selected to model the pharyngeal airflow. The calculation residual was set as 10 −3 and the iteration numbers were 1000 steps. Second-order discretization schemes were used and the pressure-velocity coupling was solved through the SIMPLE scheme.
Outcome parameters. The upper airway morphology was quantified by volume (V) and mean cross-sectional area (S mean ) of each section. The minimum cross-sectional area (S min ) could be defined and S min /S mean was calculated to determine the pharyngeal morphological uniform 11 . The pharyngeal shape was determined by ratio of the largest lateral (LR) and anteroposterior (AP) dimensions (LR/AP) of each crosssectional slice.
After simulation of the airflow, all the cross sections from nasopharynx to hypopharynx were selected along Z-axis every 1 mm for parameters calculation. The maximum pressure (P max ), minimum pressure (P min ) and maximum velocity (v max ) of the upper airway were extracted from the outcome of the CFD analysis. The airflow pressure drop (Δ P) is computed by Δ P = P max − P min , while the pharyngeal resistance (R) is obtained through the formula: R = Q/Δ P (Q was pharyngeal inlet volume flow rate). As Q was constant in our study, Δ P could directly represent the change of R in the following data statistics. All the parameters were measured three times by the same researcher, and the average value was applied to the study. Statistical analysis. All the measurements were repeated for each randomly selected patient, with a 1-week interval, to assess intra-rater reliability through paired t-test. The method error (ME) was calculated as: = ∑d n ME /2 2 (d is deviation between the two measurements; n for the number of paired double measurements) 20 .
The final measuring values were loaded into SPSS 17.0 for statistical analysis. The mean value and standard deviation (SD) were calculated for each variable. The statistical differences between the flow dynamics parameters of the upper airway in TB group were analyzed by paired-t test, while the T2 and control data were compared by the independent t-test. By this way, sections with the most significant changes in flow field characteristics of the upper airway due to Twin Block treatment could be found out. The relationship between the pharyngeal pressure drop and morphological parameters (V, S mean , S min /S mean , LP/AP) were further assessed by Pearson correlation coefficient.

Results
The age distribution and BMI showed no statistical differences between the TB and control group (p = 0.914, 0.852). There were no significant differences between the first and repeated measurements (p = 0.873), which proved high intra-rater reliability. The maximum ME of pharyngeal volume and area measurements were 0.21 mm 3 and 7.85 mm 2 , respectively. In the numerical simulation of pharyngeal airflow, the maximum ME of pressure measurements was 1.65 Pa, while for velocity measurements was 0.59 m/s. Table 1 gives an overview of the measured and computed hydrodynamic parameters in the numerical simulation. The inspiratory pressure falls continuously with the respiratory airflow move downwards from the nasal cavity in all subjects. The pressure minimum occurs in the hypopharynx while airflow travels around the outlet of the pharyngeal model. The maximum pressure drop and velocity, however, were located in the oropharynx region where the minimum cross-section area lies in all of the subjects modeled.
The airflow characteristics of the oropharynx and hypopharynx changed obviously through comparing the T2 data of TB group and the control data, while the nasopharyngeal flow field shows no statistical differences after mandibular advanced treatment. The most significant change of airflow patterns was still observed in the oropharynx. Figure 3a-c shows the changes in pressure profile of the airflow with the variation of the pharyngeal morphology. The pressure gradient in the narrowest region of oropharynx was significantly lower in the T2 data as the oropharyngeal pressure drop decreased from 67.43 ± 15.12 Pa to 40.18 ± 8.26 Pa after TB treatment (Fig. 3d-f). From Fig. 4a significantly reduced number of local eddy was also observed in the oropharynx of T2 data as the velocity of airflow slipped from 11.24 ± 1.92 m/s to 7.14 ± 1.92 m/s after mandibular repositioning. In addition, we noticed that as the inhaled airflow accelerates through the oropharynx and into the narrowest region, the most negative pressure and the maximum velocity were generally located along the anterior wall of the upper airway. A significant negative correlation was observed in the change amount of oropharyngeal pressure drop and V, S mean , and S min /S mean of the oropharynx ( Table 2).

Discussion
In this study, we explored the biomechanical reaction of obstructed upper airway using patient-specific 3D CFD models derived from CBCT scans before and after mandibular advancement. It is shown that the flow field inside the pharynx especially in oropharynx response positively to TB treatment in growing patients with class II 1 and mandibular retrognathia. In addition, a significant correlation was observed between the upper airway patency and V, S mean and S min /S mean of the oropharynx, suggesting that the anatomic morphology may predict the functional properties of the upper airway in clinical treatment.
During inspiration, as increase of the kinetic energy and dissipative energy loss in the respiratory airflow upstream, the pharyngeal luminal pressure decreases generally from 0 Pa in nostrils and reaches subatmospheric values in the oropharynx and hypopharynx region 21 . The pressure gradient was relatively large in the oropharynx for the maximum pressure drop that reflects airway resistance and quantitatively indicates the OSAS severity was observed here 18 . whereas, the pressure fluctuation in the narrow region of oropharynx was negligible and the pressure drop became much lower after mandibular advancement compared to the pre-treatment and the control group. This may result in decrease in pharyngeal compliant and thereby an improvement in pharyngeal collapse prevention capability 22 . The expansion in the oropharyngeal space could help reduce the kinetic energy of airflow and maintain a higher luminal pressure that permit a freer passage of air 21 .
Given the inhaled volume flow rate, the respiratory airflow velocity reached a peak value in the narrow region of the oropharynx. This is consistent with previous studies about adult OSAS patients, showing that the maximum velocity of airflow occurred in the velopharynx, the most occlusive area of the upper airway 23,24 . The restricted region of the oropharynx could directly induce a jet flow with significant increase in the pressure drop as well as the airflow velocity, as a linear correlation has been proved between the flow rate and the square root of the maximum pressure drop in several flow sensitive studies 18,25 .
The wall shear stresses, which depended on the pharyngeal airflow, was the crucial factor controlling the defense mechanisms of the upper airway 26,27 . For patients with narrowed airway, the excessive wall shear stresses   Table 2. Correlations between oropharyngeal pressure drop and morphological parameters.
generated by the high speed airflow along the pharyngeal wall may adversely affect the functioning of epithelial cells 28 . Children with insufficient neuromuscular response that maintains airway patency may result in airway collapse and even induce OSAS 21,29 . In our study, the airflow became more smoothly after TB treatment may indicate a significantly lower level of wall shear stresses that enhance epithelial barrier function and thus draw a protection on the mucous membrane of the pharyngeal wall 27,30 . V, S mean and S min /S mean are discovered as good parameters to predict the changes in the upper airway resistance, which confirms the results of previous researches 11,31-33 . Holsbeke. et al. 11 pointed out that S min and S min /S mean play a more dominant role in the decrease of pharyngeal resistance than other morphological parameters, as S min is more sensitive to the rotation of the tongue that caused by mandibular repositioning. Whereas, an individual difference in the location and size of S min in eupnoea was reported in the research of Ye. et al. 34 . They reckoned that the measurement errors between individuals in different measuring planes could be effectively reduced by calculation of S min /S mean . In the present study, S min /S mean shows a strong link with pharyngeal patency, indicating that the upper airway with a rather uniform area distribution react more positively than pre-treated airways which are more concave and obstructed.
The nasopharynx was the only segment that showed no statistical differences in flow field features after CFD simulation. A relatively constant anatomical morphology in nasopharynx due to TB treatment was reported in our published research 35 , which further confirmed the interaction between the pharyngeal shape and ventilation function of the upper airway.
In our research, the numerical simulation was conducted with several limiting but reasonable assumptions, including laminar model, steady-state airflow and neglecting the influence of temperature and humidity. In the particle imaging velocimetry (PIV) experiment of Hörschler. et al. 36 , the respiratory airflow state was proved to be laminar in quiet inspiratory and expiratory. The steady-state airflow inside the upper airway was verified through steady and unsteady computation in a comparative research between growing OSAS patients and normal children 33 . The temperature and humidity of airflow were considered as negligible factors for no significant effects were discovered on the pharyngeal flow field in a scaled up model of nasopharynx 37 . In spite of the limitations, these assumptions make the complicated flow geometry tractable since the flow field inside the pharynx is time varying and unknown and the flow rate measurements were not available for the patients involved.

Conclusion
The pharyngeal airflow pressure distribution and resistance response positively after mandibular advancement with the enlargement in volume, cross-sectional area and more uniform area distribution, and ultimately result in an improvement in ventilation function of the upper airway in growing patients with class II 1 and mandibular retrognathia.