Aerodynamic Simulation of Small Airway Resistance: A New Imaging Biomarker for Chronic Obstructive Pulmonary Disease

Introduction

Chronic obstructive pulmonary disease (COPD) is a leading cause of death around the world, characterized by airflow obstruction primarily caused by the disease of the small airways, airways less than 2 mm in internal diameter.^1,2 Small airway disease occurs early in the disease process, prior to emphysema.^3–5 Accurate assessment of small airway structure and function is crucial for subclinical and early COPD warnings.

The pulmonary function test is currently the gold standard for COPD diagnosis and can even identify individuals at high risk of developing COPD.^6–8 However, its value in assessing small airway obstruction remains limited.⁹ Airway obstruction can be reflected by an increase in airway resistance, clinically measured by whole-body plethysmography. However, this method does not exclusively measure the small airway resistance, but also the resistance to flow in the larger airways, which vary greatly among individuals.

There is a lack of consensus on the most effective indices suggestive of early-stage small airway dysfunctions. While the parametric response map (PRM) has offered some new insight,^10,11 its widespread implementation in clinical practice has been limited due to its resource-demanding nature and greater radiation exposure due to the paired inspiratory and expiratory scannings.

Recent advances in image-based computational fluid dynamics (CFD) have shed light on the precise quantification of airway aerodynamics, which has greatly enhanced the understanding of human upper airway physiology.¹² The CFD technique has been applied in various studies on the physiological and pathological conditions of the respiratory tract, including the studies on the characterizations of the aerodynamic behaviours of particle transportation and deposition,^13–15 and the aerodynamic feature characterization of COPD and other respiratory diseases.^16–22 Some of these studies also validated the CFD calculations with various methods, including 3D-printed in vitro models^20,22 and the combination of radioactive tracer and SPECT/CT technique.¹⁶ Nevertheless, few existing studies attempted to quantify the small airway resistance which would contribute directly to the diagnosis and management of COPD patients.

To better characterize the early changes in the small airways resistance (R_SA), we aimed to develop a novel strategy to decouple R_SA from the total airway resistance obtained from body plethysmography, by aerodynamic simulation of the large airways reconstructed from respiratory CT images incorporating PFT results as the boundary conditions. We implemented the developed strategy on a total of 24 patients classified by the PFT results as normal, high-risk, and COPD, and compared the derived R_SA against pulmonary function parameters, including the ratio of the first-second forced expiratory volume (FEV₁) and forced vital capacity (FVC), ie, FEV₁/FVC, FEV₁ measured value as a percentage of the predicted value (FEV₁% predicted), the measured value of MEF50 or 75 (maximal expiratory flow after 50 or 75% of the FVC has not been exhaled) and the MMEF75/25 (maximal mid expiratory flow) as a percentage of the predicted value (MEF50 or 75% predicted, MMEF75/25% predicted).

Materials and Methods

This retrospective study was approved by the institutional review committee of Changzheng Hospital (No.2020SL028), and the patient’s informed consent was waived due to its retrospective nature. Patient data is confidential and compliant with the Declaration of Helsinki.

Case Selection

From August 2020 to December 2020, a total of 24 subjects who attended our outpatient clinic and underwent PFT and chest CT scans were consecutively included (http://www.chictr.org.cn, ChiCTR2000035283). Inclusion criteria were as follows: patients underwent PFT and inspiratory phase CT with a maximal period of 4 weeks between the two examinations. Exclusion criteria were: 1) underlying lung diseases such as severe pulmonary interstitial fibrosis, lung cancer, and massive pulmonary infection; 2) pleural effusion; 3) a chest surgery history; 4) marked respiratory motion or metal artifact on CT images; 5) CT images having a minimal slice thickness greater than 1 mm; 6) an active pulmonary disease other than COPD; 7) an acute respiratory disease episode (ie, the development of new or increased respiratory symptoms) within a month before enrolment.

According to the GOLD (Global Initiative for Chronic Obstructive Lung Disease) recommendations, the subjects were classified into a non-COPD group (FEV₁/FVC ≥70%) and a COPD group (FEV₁/FVC<70% after inhalation of bronchodilators). Within the non-COPD group, subjects with a FEV₁% predicted value<95% were deemed at high risk of developing COPD,²³ and subjects with a FEV₁% predicted value≥95% were identified as normal.

CT Scanning and Image Evaluation

Non-contrast enhanced chest CT scanning was performed at the end of full inspiration using a 256-layer CT scanner (Brilliance iCT, Philips Medical Systems, The Netherlands). All patients underwent breath-hold training before CT scanning. The following CT scanning parameters were used: 120kV, 200mAs, 128*0.625mm collimation, scan pitch 0.915, slice thickness 1 mm, layer pitch 1mm, 500-msec rotation time.

The images were analyzed by a thoracic radiologist with 5 years’ experience who was blinded to the clinical information and PFT results of the patients.

Pulmonary Function Test

The PFT was performed 30 minutes after inhalation of the bronchodilator Ventolin 200ug using a spirometer (Jaeger, Würzburg, Germany) by a well-trained respiratory technologist. We recorded PFT results including FEV₁, FEV₁/FVC, FEV₁% predicted, MEF50% predicted, MEF75% predicted and MMEF 75/25% predicted. Besides, body plethysmography was carried out in a pressure box (Jaeger, Würzburg, Germany), in which the subjects were located. When the pressure of the box stayed stabilized, the subjects were instructed to insert a mouthpiece and hold it over the nose with a clip, and then they breathed quietly. The airway valve was closed during subjects’ normal exhalation and the subjects were instructed to gasp gently into the occlusion. The parameter of body plethysmography, airway resistance (R_BP) at the expiration phase, was recorded.

Aerodynamic Simulations

We reconstructed the conducting airway model from the end-inspiratory CT images with a 3D slicer and generated tetrahedral computational grids with ICEM-CFD (ANSYS, Canonsburg, USA). The maximal sizes for the surface and volume elements, which were designed to meet the stability of CFD calculations, were set as 0.05 and 0.1 mm, respectively. The maximal size must be smaller than the minimum diameter of the conducting airway in order to accurately calculate the flow velocity and pressure inside the airway. The ANSYS-CFX solver was used to resolve the flow physiology within the first second of forced expiratory. The outlet boundary condition at the tracheal opening was specified as the standard atmosphere (ie, 101,325 Pa), whereas the inlet boundary condition at the distal end of each bronchiole was derived from the FEV₁ following a splitting method at each bifurcation under all distal ends were reached:

(1)

where Q_i is the air flowrate at a daughter bronchiole with a diameter D_i, and n is the total number of daughter bronchioles at an airway bifurcation. A steady-state CFD simulation was performed with the density and viscosity of air assumed to be 1.15 kg/m³ and 1.87×10⁻⁵ kg/m·s.

Derivation of the Small Airway Resistance

A convergence state of the simulation was reached when the residuals of the normalised velocity, pressure, and continuity were no greater than 10⁻⁵. We calculated the effective resistance of the reconstructed airway model based on the topology by analogy with an electric circuit, where the resistance of a bronchiole segment between two neighbour bifurcations was calculated as

(2)

where ∆P represents the pressure drop across the neighbouring bifurcations and Q_i is the average airflow within the bronchiole segment. We assumed the difference between the total resistance measured by body plethysmography (R_BP) and the effective resistance of the reconstructed airway model (level I to VI) calculated by CFD (R_CFD) to be the resistance of the small airways (R_SA):

(3)

Statistical Analyses

The two-tailed Student’s T-test was used to compare the aerodynamic parameters across different subject groups. A variance was used to measure differences in demographics among the three groups. The relationships between the derived R_SA and pulmonary function parameters were evaluated by the Pearson correlation coefficient or Spearman correlation coefficient, depending on whether the data was normal distribution or not. Throughout the study, p < 0.05 was considered statistically significant. All statistical analyses were done with the MedCalc^® Statistical Software Version 22.007 (MedCalc Software Ltd, Ostend, Belgium; https://www.medcalc.org; 2023).

Results

Basic Demographic Data and Pulmonary Function Tests

Of the 24 subjects enrolled in the study, there were 10 normal subjects, 6 high-risk subjects and 8 COPD patients. The detailed demographics of the included subjects are shown in Table 1. There were no significant differences among the three groups in gender (the control group: men 20%; the high-risk group: men 50%; the COPD group: men 63%), age (the control group: mean 51.1 years; the high-risk group: mean 58.8 years; the COPD group: mean 60.3 years), and body mass index (the control group: mean 23.9 kg/m²; the high-risk group: mean 27.6 kg/m²; the COPD group: mean 25.1 kg/m²) (P>0.05). The PFT parameters, including FEV₁ (the control group: mean 2.7 L; the high-risk group: mean 2.6 L; the COPD group: mean 2.1 L), FEV₁% predicted (the control group: mean 106.2%; the high-risk group: mean 88.9%; the COPD group: mean 82.7%) and FEV₁/FVC (the control group: mean 84%; the high-risk group: mean 82.2%; the COPD group: mean 65.4%), were significantly different among the three groups (P=0.04, P<0.001, P<0.001).

Table 1 Patient Demographics

Airway Aerodynamic Characteristics

We compared the total airway resistance measured by body plethysmography (R_BP), resistance of the reconstructed airway model estimated by CFD simulation (R_CFD), and the derived resistance of the small airways (R_SA) among the three groups and found that the R_BP (mean 0.45 kPa s/L), R_CFD(mean 0.34 kPa s/L), R_SA(mean 0.13 kPa s/L) and R_CFD/R_SA as a percentage of total airway resistance (R_CFD%/R_SA%, mean 72.8%/27.2%) in the normal control group were significantly lower than those (R_BP: mean 0.72 kPa s/L; R_CFD: mean 0.21 kPa s/L; R_SA: mean 0.52 kPa s/L; R_CFD%: mean 27.63%; R_SA%: mean 72.38%) in the COPD group(p=0.012; p=0.022; p<0.001; p<0.001; p<0.001). It was shown that the R_SA of the high-risk group had a significantly greater absolute value (mean 0.31 kPa s/L) and proportion (mean 54%) (p=0.002, p<0.001) than the control group. The R_SA and R_SA% in the high-risk group were significantly lower than those in the COPD group (p=0.020; p<0.001) (Figure 1, Table 2).

Table 2 Total Airway Resistance Measured by Body Plethysmography (RBP), Resistance of the Reconstructed Airway Model Estimated by CFD Simulation (RCFD), and the Derived Resistance of the Small Airways (RSA)

Figure 1 Hydrodynamic simulations performed in the airway model. Visualisation of the airflow streamlines (Top), and pressure (Bottom) distributions of three representative cases respectively for the normal control (Left), patients at a high risk of developing COPD (Middle) and diagnosed with COPD (Right).

Correlation of Airway Aerodynamics with Parameters Obtained from PFTs

As illustrated in Figure 2, R_SA was significantly negatively correlated with MEF75% (r = −0.71, p < 0.001) and moderately negatively correlated with FEV₁/FVC (r = −0.70, p < 0.001), FEV₁% (r = −0.60, p = 0.002), MEF50% (r = −0.64, p = 0.001) and MMEF75/25% (r = −0.64, p = 0.001). R_SA% was moderately negatively correlated with FEV₁/FVC (r = −0.67, p < 0.001), FEV₁% (r = −0.57, p = 0.004), MEF50% (r = −0.64, p = 0.001), MEF75% (r = −0.60, p = 0.002) and MMEF75/25% (r = −0.64, p = 0.001).

Figure 2 Correlations of the estimated small-airway resistance (RSA) and the small airway resistance as a percentage of the total airway resistance (RSA%) with pulmonary function parameters (FEV1/FVC ratio, FEV1% predicted, MEF50% predicted, MEF75% predicted and MMEF75/25% predicted).

In females, R_SA was moderately negatively correlated with FEV₁% (r = −0.65, p = 0.011) and significantly negatively correlated with MEF50% (r = −0.71, p = 0.005), MEF75% (r = −0.72, p = 0.003) and MMEF75/25% (r = −0.72, p = 0.003). R_SA% was moderately negatively correlated with MEF50% (r = −0.55, p = 0.042) and MMEF75/25% (r = −0.60, p = 0.023). (Figure 3)

Figure 3 Correlations of the estimated small-airway resistance (RSA) and the small airway resistance as a percentage of the total airway resistance (RSA%) with pulmonary function parameters (FEV1/FVC ratio, FEV1% predicted, MEF50% predicted, MEF75% predicted and MMEF 75/25% predicted) in females.

In males, R_SA was moderately negatively correlated with MEF75% (r = −0.65, p = 0.041). R_SA% was significantly negatively correlated with FEV₁/FVC (r = −0.79, p = 0.006), MEF50% (r = −0.76, p = 0.011), MEF75% (r = −0.76, p = 0.011) and MMEF75/25% (r = −0.78, p = 0.008) and moderately negatively correlated with FEV₁% (r = −0.69, p = 0.028). (Figure 4)

Figure 4 Correlations of the estimated small-airway resistance (RSA) and the small airway resistance as a percentage of the total airway resistance (RSA%) with pulmonary function parameters (FEV1/FVC ratio, FEV1% predicted, MEF50% predicted, MEF75% predicted and MMEF 75/25% predicted) in males.

Discussion

In the present study, we have reconstructed subject-specific airway models for ten normal controls, six high-risk subjects of developing COPD, and eight COPD patients. We performed CFD for each patient to resolve the resistance of flow passing through the airway model, incorporating the FEV₁. This then allows the effective resistance of the small airways distal to the generated models to be quantified, by comparing the derived resistance with the total airflow resistance measured by PFT. Results suggested that R_SA and the R_SA% obtained from CFD can already identify high-risk subjects when total airway resistance from pulmonary function tests has not yet changed. R_SA and the R_SA% were in correlation with pulmonary function parameters, including FEV₁/FVC ratio, FEV₁% predicted, MEF50% predicted, MEF75% predicted and MMEF 75/25% predicted.

Although the R_SA was calculated based upon the total flow resistance measured by body plethysmography, it is worth noting that these two are not linearly correlated and the decoupling of R_SA requires the pressure loss in the larger airway system to be characterised in a patient-specific manner. Likewise, FEV₁ does not immediately determine the resistance to flow imposed by the reconstructed larger airway structure, although it has been employed as the boundary condition to set up the aerodynamic simulation. Changes in pressure loss are more likely to be affected by the structure of a large airway reconstructed from CT images. In this regard, combining CFD with patient-specific chest CT imaging and PFT enables the estimation of R_SA, which is otherwise not quantifiable by any existing non-invasive approach.

Patients with COPD have chronic inflammation in the lungs, which leads to bronchial remodeling. It has been reported that small airway disease occurs prior to emphysema in COPD.^3–5 COPD patients and high-risk smokers who have not yet met the diagnostic criteria for COPD could have small airway dysfunction.^24,25 Early detection of small airway abnormalities is imperative for the early diagnosis of COPD.

To quantify small airway abnormalities, we explored a new approach, combining fluid dynamics based on CT images and lung function, to calculate the airway resistance. Our study successfully reconstructed airway models in which the airflow velocity of high-risk individuals and COPD patients decreased in varying degrees, which is consistent with clinical findings. Based on this, we found that it is feasible to use the R_SA to reflect airway abnormalities since it tends to increase as COPD progresses, and it is significantly different among normal controls, high-risk groups, and COPD patients. Besides, the small airway resistance is in moderate or significant correlation with lung function parameters, including MEF50% predicted, MEF75% predicted and MMEF 75/25% predicted, which reflect small airway obstruction. It suggests that small airway resistance derived from airway CFD can be used as a new imaging biomarker to participate in the early warning of COPD.

Compared to previous studies on the CFD simulation of aerodynamics in patients with respiratory diseases, our study is the first to apply the CFD analysis to decouple the small airway resistance from the total airway resistance. To our knowledge, there were no earlier attempts to measure the resistance of small airways in COPD patients in this way. Besides, due to the high cost and time-consuming nature of digital reconstruction of human airway geometry from the medical imaging data, most studies using this methodology included very few subjects (typically <10 subjects);^16–20 the number of subjects included in our experiment is second only to the study published by Van de Moortele et al.²² The expansion of data scale can better reveal the relationship of aerodynamic parameters and the development of COPD, and also lead to more robust conclusions.

Our study had several limitations. First, from a clinical perspective, this study still had a small sample size due to its single-center and retrospective nature, although the included subject number is higher than in many similar studies. Future studies would benefit from testing the reported findings in a larger group of patients with a greater diversity of COPD severity. Second, the aerodynamic simulation was performed in a steady-state fashion, without having considered the changes in airway flow resistance over the entire respiratory cycle. This is partly due to the potential high computational cost of transient CFD simulations, as well as the associated non-negligible change in the airway structure, reconstruction of which is still beyond the capacity of current CT technologies. Third, despite correlations between R_SA and PFT were found, the superiority of R_SA in predicting patients at high risk of developing COPD remains inconclusive, because the high-risk population of COPD included in this study was defined only by pulmonary function tests, without clinical information, such as smoking history. The data volume should be expanded to include high-risk populations defined by diverse criteria or the patient’s clinical follow-up to further validate the value of small airway resistance in identifying the COPD high-risk population. Fourth, this method requires patients to undergo spirometry, plethysmography, and CT scans, which limits its use outside of a hospital setting. Finally, the included subjects did not undergo paired respiratory phase chest CT scanning, therefore PRM parameters could not be obtained. Furthermore, since this is a retrospective study, it is difficult to collect information on the patients’ clinical symptoms during the examinations, such as cough and dyspnea. As a result, there is a lack of analysis of correlations between airway resistance parameters and PRM parameters as well as clinical symptoms. Further prospective studies will be conducted to analyze the correlations to confirm the clinical application value of small airway resistance.

Conclusion

In conclusion, airway CFD combined with PFT is a valuable method for estimating small airway resistance. The derived R_SA helps to differentiate high-risk populations from normal populations and COPD patients, suggesting its potential in the early diagnosis of COPD. Further validation of our method on an independent, large follow-up data set is required to investigate accuracy.