Airway Resistance and Respiratory Distress in Laryngeal Cancer: A Computational Fluid Dynamics Study

Obstructive upper airway pathologies are a great clinical challenge for the airway surgeon. Protection against acute obstruction is critical, but avoidance of unnecessary tracheostomy must also be considered. Decision‐making regarding airway, although supported by some objective findings, is largely guided by subjective experience and training. This investigation aims to study the relationship between clinical respiratory distress and objective measures of airway resistance in laryngeal cancer as determined by computational fluid dynamic (CFD) and morphometric analysis.


BACKGROUND
Obstructive upper airway pathologies are a great clinical challenge for the airway surgeon.Protection of the respiratory tract against the possibility of acute obstruction is always considered paramount. 1However, the introduction of a surgical airway like a tracheostomy brings with it a significant negative impact on patient quality of life. 2,35][6] An additional option of endoscopic tumor debulking exists as well for select neoplastic cases, but its effects are temporary and it does not always preclude the need for tracheostomy before or during definitive cancer treatment. 7,8irway analysis is therefore an essential part of the evaluation for a new laryngeal cancer diagnosis.Although some cases present in acute airway distress requires emergent management, 7,9,10 the majority of cases can proceed with a standard work-up, including history and physical with flexible fiber-optic laryngoscopy, computed tomography (CT) AE magnetic resonance imaging (MRI) of the primary tumor and neck, and in-office or operative biopsy. 11Proper airway analysis is a complex task.It requires integration of clinical observations like dyspnea and stridor, measured data like airway caliber, blood gas testing, and radiologic findings.The goal for airway management is to make decisions based on objective data.Currently, however, much of this decision-making is based on subjective experience and training.
Recent computational advances have introduced more sophisticated techniques for modeling the human upper aerodigestive tract.The pathology most extensively studied in this manner is obstructive sleep apnea, where 3D shape analysis, [12][13][14] computational fluid dynamics (CFDs), [15][16][17] and fluid-structure interaction 18,19 focused on the pharynx have shown an ability to predict disease severity and response to treatment.][24] With respect to the larynx, some previous work has experimentally validated methods of approximating flow through a healthy glottis using calculated coefficients. 25,26ther laryngeal analysis has been able to demonstrate a correlation between regions of higher fluid pressure and regions prone to developing cancer. 27Computational models are another experimentally validated approach, 28 and indeed, CFD analysis lends itself well to the larynx and trachea given their relatively rigid boundaries and complex geometry, particularly when pathology is present.From this, a series of clinically oriented studies have evaluated fluid flow with CFDs in subglottic stenosis [29][30][31][32] and vocal cord immobility. 33,34Again, CFD analysis has shown its applicability to clinical cases both for surgical planning and in its correlation to clinical findings in small sample sizes.
There are no reported studies to date that evaluate fluid flow around laryngeal tumors.In addition, there is overall a very small number of samples that have been used to correlate CFD model results against clinical findings and outcomes in laryngotracheal disease.This study investigates the relationship between clinical respiratory distress in laryngeal cancer and airflow resistance determined by CFD analysis.It will additionally study the correlation between morphometric measurements, such as length and degree of stenosis, and clinical outcomes.Our hypothesis is that the computational fluid model will produce flow parameters that are significantly correlated to patient outcomes, demonstrating its effectiveness as a potential risk prediction tool.

Patient Selection
This was designed as a retrospective case-control study of patients newly diagnosed with laryngeal cancer at a single tertiary care institution in Montreal, Canada.Patients were diagnosed by one of five academic Otolaryngologists.Included cases were adults with a biopsy-proven laryngeal cancer of any subsite who required urgent surgical management or hospitalization for acute airway obstruction.Cases were excluded if they had no available CT neck or larynx prior to definitive oncologic treatment, if a stabilizing airway surgery performed prior to CT scan significantly affecting obstruction geometry (i.e., tumor debulking), or if the patient had a separate concurrent upper airway pathology.Controls were obtained based on the same criteria as above except with no surgical or inpatient care performed prior to definitive oncologic treatment.
The study protocol was evaluated and approved by the research ethics board of the McGill University Health Centre (2021-7178).
A total of 20 cases and 20 controls were obtained from January 2014 to December 2020, baseline characteristics for which can be seen in Table I.Controls were matched to cases based on the tumor staging, split into either early (T1 or T2 according to the American Joint Committee on Cancer (AJCC) 8th edition 35 ) or locally advanced (T3 or T4) disease.After categorizing based on matching data, controls were determined at random.

Airway Reconstruction
All CT scans were obtained using a GE Revolution scanner (GE HealthCare, Chicago, IL, USA).The highest resolution acquisition available for each patient was used for further processing, with either 1 mm (neck protocol) or 0.5 mm (larynx protocol) voxel size.De-identified CT scans were processed in Dragonfly v2021.1 (ORS, Montreal, Canada) for the purposes of airway reconstruction.Airway anatomy was segmented from the nasal/oral inlet down to trachea at the level of the sternal notch.Nasal, oral, or both inlets were incorporated into the reconstruction based on the patency of the oral cavity and velopharynx, respectively, and if neither were patent (i.e., tracheostomy) then an artificial velopharyngeal airway was approximated manually.Thresholding was used to differentiate air from surrounding tissues, followed by region of interest analysis to isolate the main structure from background and non-contiguous structures.If a tracheostomy or nasogastric tube were present, a watershed transform applied on a Sobel-filtered image identified the relevant segments, 36 which were incorporated into the airway.For the purposes of the fluid simulation, the walls of the airway and the outlet were converted into 1 mm thick volumes (i.e., a conduit and an end piece) using dilate and subtract operations.The airway, airway walls, and outlet were then converted into surface meshes, and a Laplacian smoothing operation 36 was finally used to remove any voxel edge artifacts.

Computational Fluid Dynamics
Simulations were performed in PowerFLOW v6-2019 (Dassault Systèmes, Vélizy-Villacoublay, France), which is based on the lattice Boltzmann method for modeling fluid flow.Steadystate flow was simulated using a constant inspiratory flow rate _ Q of 15 L/min, which is the one of the most commonly used flow rates in CFD literature for approximating adult breathing at rest. 31,37Global parameters, including viscosity, specific heat capacity, and molecular weight, were all set to standard values of air at standard ambient temperature and pressure (SATP).PowerFLOW features a renormalization group (RNG) form of the k À ε equations as its turbulence model, and for these simulations the default turbulence intensity was set to 1% and the turbulence length scale 2 mm.The imported airway, airway boundaries, and outlet surface meshes were placed in a bounding box filled with air.The following boundary conditions were used: 1. constant mass flow at the tracheal outlet of 0.0003 kg/s (15 L/ min equiv.at SATP); 2. constant static pressure inlet anterior to the nostrils at standard pressure (101 kPA); 3. stationary walls with no-slip condition (zero velocity at interface).
A convergence study was used to determine optimal discretization voxel size.For this, a series of simulations on a normal airway were performed starting with a voxel length of 1 mm and decreasing by one-half until the steady-state total pressure in the mid-trachea converged (<5% change in subsequent simulations).This was repeated on a critically narrowed airway with a thickness of 1 mm to ensure minimal discretization error in the stenotic segment.The discretization was subdivided into variable resolution regions, with the stenotic segment the highest, followed by the rest of the airway at half resolution, and the background air at one sixteenth resolution.
To verify the model, a series of five normal airways were first processed prior to the cases and controls.These patients had no upper airway pathology noted in their chart or on their CT scans.The results of these simulations were analyzed in the same fashion as the others, calculating the pressure drop ΔP across the larynx (from mid-pharynx to mid-trachea), maximum velocity at the narrowest portion (in this case the glottis), and the calculated resistance R ¼ ΔP= _ Q.Given the relatively large volume of computation required for each simulation, a parallelized solution was implemented and executed on Calcul Québec (www.calculquebec.ca),now a part of the Compute Canada (www.computecanada.ca)supercomputer networks.Computations were parallelized using up to four nodes (ranging from 32 to 40 cores each) of the supercomputer CPU (Intel Gold 6148 Skylake, 2.4GHz).This step decreased the total computation time to convergence for each simulation from approximately 1-5 days to 3-12 h and allowed multiple simulations to be evolved in parallel.

Morphometric Analysis
As a secondary analysis, airway geometry was compared to clinical and simulation outcomes by using several radiologic measurements.All morphometrics were collected using the Dragonfly platform.Measurements were made on axial sections from the standard image acquisitions.Although some previous literature has reported morphometrics based on angled planes, 29,31 this study aimed to model the relationship of what a clinician may measure during their analysis on the standard image planes.
First, on the axial sections, the minimum cross-sectional area of the airway was determined, A min .On the same slice, the approximate minimum thickness of the airway was also measured, d min , using the short axis if the shape was ovoid or rectangular.This represents the airway size as reported by some clinicians (e.g., a 2 mm airway).Inferior to this, the crosssectional area of the mid-trachea, A mt , was also recorded from axial slices.The stenosis-tracheal ratio, a quantity previously described in subglottic stenosis CFD literature, 29 was then approximated from these quantities, RA st ¼ A min =A mt : Finally, the length of the stenotic segment was determined as an approximated measurement on parasagittal slices, L s .

Statistical Analysis
All statistical analysis was performed in R v4.0.2 (R Foundation for Statistical Computing, Vienna, Austria).The primary analysis aimed to evaluate the relationship between airway resistance and the clinical endpoint of acute airway obstruction requiring surgical or inpatient management.This was performed first by univariable logistic regression, followed by multivariable regression combining resistance measurements with patient characteristics (age, obesity, COPD).In secondary analysis, a comparison of morphometric measurements (A min , t min , RA st , L s ) to clinical airway outcomes was performed using linear regression.This was followed by performing a statistical regression of the geometric measurements to the CFD results, using thickness to the negative fourth power according to Poiseuille's law. 38And finally, the correlation between morphometrics and resistance calculations was evaluated, with the intention of verifying basic geometric measurements as a predictor of fluid dynamics.

RESULTS
The initial convergence study was performed on a normal airway not included in study data, with minimum airway thickness at the glottis of 15 mm.The results demonstrated minimal change in the inlet to midtracheal pressure drop from a voxel length of 250 μm.After completion of segmentation and morphometric analysis over the study data set, a second convergence study was performed on a case with a critically narrowed larynx (d min = 1.96 mm).This was to ensure minimal discretization error when a high degree of stenosis was present.Convergence was obtained for an isometric voxel length of 125 μm at the stenosis.The Reynolds number computed during this convergence simulation was 7611, which is well above the threshold considered for fully turbulent flow (Re > 4000).This demonstrated the further need for turbulence modeling in the CFD analysis, and subsequent simulation of cases demonstrated Reynolds numbers in the transition or fully turbulent ranges (from 3390 to 10018).
Simulations performed on the five normal airways yielded an average pressure drop of 5.9 Pa (range 0.4-14.9Pa), average resistance of 0.0236 Pa s/mL (range 0.0016-0.0596Pa s/mL), and a max flow velocity of 2.76 m/s (range 0.69-5.21m/s).For the pathologic cases, a representative sample of CFD simulations from nostril to trachea are illustrated in Figure 1, ranging from minimal to severe obstruction.In mild cases, most of the airway resistance was generated by the nasal cavity, with the glottis only contributing a few pascals of pressure loss.With significant laryngeal narrowing, however, pressure drop rapidly increased relative to decreasing airway caliber.This is emphasized between Figure 1C,D, where a decrease in a fraction of a millimeter resulted in an order of magnitude larger pressure drop.
Measurements from the CFD analysis included pressure drop and maximum air velocity, reported in Table II.Resistance was calculated according to R ¼ ΔP= _ Q. Due to the quadratic nature of aerodynamic resistance in duct systems, CFD measurements spanned several orders of magnitude (range: 0.008-328.1Pa s/mL), which contributed to the high variance of collected data.Despite this, both resistance and maximum velocity were significantly correlated to the primary outcome through logistic regression.Morphometric analysis similarly showed statistically significant correlation for measurements evaluating the airway cross-section.Length of stenosis, however, was not significant.
The primary statistical endpoint, airway obstruction as a function of calculated airway resistance, was combined into a multivariable regression model with relevant patient characteristics (Table III).The resistance odds ratio continued to show a statistically significant relationship in multivariable analysis.Age and existing diagnosis of COPD both showed a trend as expected of increased odds with less favorable status, however they did not show statistical significance.Obesity (defined as body mass index (BMI) > 30 kg/m 2 ) had a relationship opposing expectation, but again without statistical significance.
In secondary analysis, the relationship between CFD-obtained airway resistance and morphometric measurements was evaluated (Fig. 2).A power law to the negative fourth degree, y = ax À4 + b, was statistically fit to verify the efficacy of simple thickness measurements as a predictor of resistance.The coefficients of the fit resulted as a = 253.09and b = 7.15, the adjusted R 2 value 0.832, and overall, it achieved statistical significance ( p < 0.001).This is in keeping with the Hagen-Poiseuille equation for flow through a conduit.It also illustrated a thickness measurement threshold of approx.2.5 mm-below which resistances climb exponentially.In Figure 3, resistance distributions in each laryngeal subsite are shown.No statistically significant difference was found between groups, but the glottic distribution had the largest proportion of lesions with clinically significant resistances.

DISCUSSION
Major advances in computer processing technology over the recent decades has revolutionized our ability to solve complex problems in many fields including engineering, math, and basic sciences.In clinical medicine, however, the use of computational methods remains infrequent.This is a particularly underutilized tool in the field of Otolaryngology, where many functional outcomes such as hearing, balance, respiration, and phonation are based on measurable physical quantities for which accurate modeling could influence patient care.Some of the currently used computational methods in the head and neck include radiation dosimetry and planning, which are required as part of the current standard of care in radiation oncology, and CFD analysis for nasal cavity airflow, which has recently become one of the clinical tests available for quantifying and localizing nasal obstruction. 39omputational models of laryngeal flows have been studied as well, but thus far sample sizes have been small, and it have not been validated for use in adults.This study aimed to further evaluate CFD models for airflow through a pathologic human larynx, with a focus on airway obstruction in laryngeal cancer.
The results demonstrated a significant relationship between calculated airway resistance and acute airway obstruction.It confirmed this result on multivariable regression, and additionally showed a significant relationship between maximum air velocity and obstruction.This suggests that CFD analysis could be an effective model for estimating resistance and predicting airway distress in a pathologic larynx.The clinical factors included in the multivariable analysis (COPD, obesity, and age) did not reach significance, although this study was not powered to properly assess these variables.Increasing sample size could be used to assess the suspected correlations.Subgroup analysis did not find a significant difference of the laryngeal subsites in terms of calculated resistance or clinical outcome, although it did suggest a possible trend toward higher resistances in the glottis.This would be consistent with airway anatomy, where the glottis has the narrowest cross-section between nasopharynx and carina and would intuitively be most vulnerable to obstruction.Morphometric analysis also yielded multiple relevant relationships with the primary outcome, where all measurements based on the minimum airway crosssection demonstrated significance.The length of stenosis measurement was not significant, although this may be partly because stenosis length is difficult to approximate in cases of mild obstruction, where subtle superficial tissue changes can be difficult to approximate.Morphometric data were also shown to be significantly related to resistance calculations on a regression analysis, and when observing this relationship (demonstrated in Fig. 2), a threshold airway thickness less than 2.5 mm appears to be associated with a rapid rise in resistance, likely to lead to obstruction.
Previous studies evaluating CFDs in the normal human airway from pharynx to trachea have reported pressure drops of 3.83 Pa, 32 4 Pa, 18 and 12.89 Pa, 40 and resistances of 0.0191 Pa s/mL 18 and 0.02 Pa s/mL. 32To assess the validity of this model, five normal airways were processed and analyzed.Although there is notable variability stemming from the range of normal airway geometry, the average pressure drop and resistance of 5.9 Pa and 0.0236 Pa s/mL, respectively, fall well within these ranges.To help further assess the validity of the model, an inverse power law to the fourth power (i.e., 1/ d 4 ) was fit using the stenotic thickness to resistance relationship.This demonstrated statistical significance (p < 0.001), showing that the model is in keeping with the  commonly accepted Hagen-Poiseuille equation for flow through a conduit, Δp ¼ 8 μL Q=πr 4 (where pressure drop, and by extension resistance, is related to conduit radius the negative fourth power).CFD analysis, although very powerful in its ability to model fluid flow through complex geometry, requires intensive computational resources and long computation times.This study used a supercomputer network, a resource not generally available in a clinical setting, and computations times were still in the order of hours (from days on a single CPU).Given the urgency of treatment decisions in airway distress, the model presented in its current form would not be fast enough to be clinically useful in most cases.Further work will be needed to develop this into a valuable clinical tool, and there are multiple forms that it could take.One option could be to develop a large data set that represents the spectrum of laryngotracheal stenosis in its various forms and have a registration algorithm that matches the current patient's airway to a similar, previously analyzed airway.This could provide resistance and risk estimates with a muchreduced computation time.Another would involve a machine learning algorithm, where a large dataset of patients presenting with laryngeal cancer prior to airway distress would be analyzed to develop a risk assessment tool, predicting future tracheostomy need.
There were several limitations in the current study.For patient selection, the number of patients who met criteria over the studied period was quite limited, which prevented the matching of patients beyond a simple binary early versus advanced tumor status.To improve the limited sample size, patients who underwent tracheostomy prior to CT (medically unstable) were included if the operative team did not significantly alter airway geometry (i.e., tumor debulking).Some sampling bias is still likely present from cases that were excluded for this reason.Data collection was performed retrospectively, and therefore was susceptible to its inherent biases, and the primary outcome was a product of each surgeon's decision-making, which is only partly based on objective data and is likely to have some inter-rater variability.Other important clinical variables such as subjective dyspnea, laryngoscopy findings, and vital signs were not included, which likely limited the efficacy of the statistical model.Given that  laryngoscopy is generally considered more clinically relevant than imaging in airway analysis, future work should consider its implementation.Subjective dyspnea questionnaires and numerical BMI would likely also contribute significantly.For the image segmentation, volumes were generated from a single image capture with the patient in supine position, which may be distorted from baseline if the patient was not following protocol (i.e., phonating, swallowing).For the computational fluid analysis, rigid boundary conditions may have affected results in the oropharynx and supraglottis given their softer tissue characteristics, although they should still apply well in advanced tumor disease as laryngeal cancer is often fixed and rigid.Regarding CFD parameters, all CFD simulations were computed using a single flow rate representing normal breathing at rest.This may have led to some error with cases of moderate to severe stenosis, as some of the pressure drops calculated would not be feasibly generated by human lungs (>100 cm H 2 O).Adjusting for this though would be a challenge, as patients in obstruction naturally adjust their tidal volume and flow rate to optimize gas exchange -adjustments that are difficult to predict.Given that the primary calculation, resistance, is naturally adjusted for flow rate, however, this error is likely minimal.

CONCLUSIONS
CFD analysis applied to the human larynx in the setting of laryngeal cancer demonstrated a significant relationship between calculated airflow resistance and an acute airway obstruction.Morphometric measurements were also shown to have a significant relationship both to the primary outcome and to the calculated resistance, suggesting that a simplified geometric statistical model could be used to help predict resistance and outcome.One measurement in particular, minimum airway thickness less than 2.5 mm, was shown to be correlated with exponential increases in airflow resistance.Future work would validate these findings and attempt to build a predictive model of airway resistance that combines clinical parameters with calculated CFD and/or morphometric measurements.This would provide further objective data to help clinicians make informed airway decisions.

Fig. 1 .
Fig. 1.Axial CT images and computational fluid dynamic pressure distributions from nostril to trachea.(A) Early stage glottic lesion and (B) advanced transglottic lesion, managed as outpatient.(C) and (D) severe and very severe obstructions at the level of the glottis/subglottis, requiring urgent (within 48 h) and emergent (immediate) airway surgery respectively.

Fig. 2 .
Fig. 2. Scatter plot of resistance data as a function of stenosis thickness, with overlying curves statistically fit to a negative exponential function to the fourth degree.

Fig. 3 .
Fig. 3. Distribution of resistances (converted to log base 10 for visualization) in each laryngeal subsite.

TABLE I .
Baseline Characteristics.

TABLE II .
Study Measurements, Statistically Correlated to the Primary Outcome by Simple Logistic Regression.A min = minimum cross-sectional area; d min = minimum airway thickness; L st = length of stenosis; R = resistance; RA st = ratio of minimum to mid-tracheal cross-sectional area; V max = maximum air velocity.

TABLE III .
Univariable and Multivariable Analysis, Evaluating the Primary Measurement-Resistance as Obtained from CFD Analysis-Along With Patient Risk Factors as They Relate to Airway Obstruction.