Estimating Airway Resistance from Forced Expiration in Spirometry

Featured Application: Authors are encouraged to provide a concise description of the speciﬁc application or a potential application of the work. This section is not mandatory. Abstract: Spirometry is the gold standard to detect airﬂow limitation, but it does not measure airway resistance, which is one of the physiological factors behind airﬂow limitation. In this study, we describe the dynamics of forced expiration in spirometry using a deﬂating balloon and using this model. We propose a methodology to estimate ζ (zeta), a dimensionless and e ﬀ ort-independent parameter quantifying airway resistance. In N = 462 (65 ± 8 years), we showed that ζ is signiﬁcantly ( p < 0.0001) greater in COPD (2.59 ± 0.99) than healthy smokers (1.64 ± 0.18), it increased signiﬁcantly ( p < 0.0001) with the severity of airﬂow limitation and it correlated signiﬁcantly ( p < 0.0001) with airway resistance (r = 0.55) and speciﬁc conductance (r = − 0.60) obtained from body-plethysmography. ζ also showed signiﬁcant associations ( p < 0.001) with di ﬀ usion capacity (r = − 0.64), air-trapping (r = 0.68), and CT densitometry of emphysema (r = 0.40 against % below − 950 HU and r = − 0.34 against 15th percentile HU). Moreover, simulation studies demonstrated that an increase in ζ resulted in lower airﬂows from baseline. Therefore, we conclude that ζ quantiﬁes airway resistance from forced expiration in spirometry—a method that is more abundantly available in primary care than traditional but expensive methods of measuring airway resistance such as body-plethysmography and forced oscillation technique.


Introduction
Airflow limitation or obstruction is the characteristic feature of obstructive lung diseases such as chronic obstructive pulmonary disease (COPD) and asthma. In COPD, it arises out of the dual effects of permanent parenchymal destruction (emphysema) and small airway dysfunction due to prolonged exposure of the lungs to cigarette smoke [1]. While in asthma, it arises due to narrowing of airways from excessive inflammation and mucus production as a result of bronchial hyper-responsiveness to allergens [2].
In clinical practice, spirometry is the gold standard to measure airflow limitation [3]. During a spirometry test, a subject performs a forced expiration using a spirometer that measures the flow and volume of exhaled air. Several indices of diagnostic importance, such as forced expiratory volume in 1 second (FEV1), forced vital capacity (FVC), and Tiffeneau's index (FEV1/FVC) are obtained from spirometry [4]. Two widely accepted criteria for the presence of airflow limitation are a FEV1/FVC ratio below a fixed cutoff of 0.7 [3] or below the lower limits of normal (LLN, defined as the 5th percentile of a normally distributed set of values of FEV1/FVC for a population of non-smoking, normal individuals) [5]. The plot of expiratory flow against volume measurements is termed as maximal expiratory flow-volume curve (MEFVC), and it is also known to be associated with different pathological states [6].
Although spirometry is used to indicate airflow limitation, the latter itself is an end-result of many factors.
One of these factors is an increased airway resistance, and therefore, measuring it may directly provide additional insight into the occurrence of airflow limitation. Traditionally, body-plethysmography has been used to determine airway resistance (Raw) as a ratio of driving alveolar pressure to airflow. Other methods of measuring airway resistance include forced oscillation technique (FOT), which estimates respiratory system resistance (R RS ) and reactance (X RS ), and the interrupter technique which measures interrupter resistance (R int ) at a single moment of tidal breathing [7].
All of the above-mentioned techniques measure airway resistance during tidal breathing. Nevertheless, factors affecting airway resistance during a forced expiration such as loss of lung elasticity [8] and dynamic airway collapse cannot be underestimated [9], and often compose the main reason for airflow limitation and clinical symptoms [10]. Measures of airway resistance obtained during forced expiration will provide better insight into these mechanisms, which are potentially amendable for intervention. Currently, there exists no method in the literature to calculate dynamic airway resistance from the forced expiratory manoeuvre. The consensus is that spirometry cannot be used to estimate airway resistance since it does not measure alveolar pressure [7,11,12]. In the past, researchers have proposed parameters reflective of airway resistance [13] that only describes the shape of MEFVC, without any physiological basis.
In this study, we describe the dynamics of forced expiration using a model of a deflating balloon, and thereby propose a method to estimate a parameter representative of airway resistance from spirometry. We hypothesized that this resistance parameter would be significantly higher in COPD compared to healthy smokers, increase with the severity of airflow limitation in COPD patients [14] and show a significant correlation with Raw and its flow standardized inverse, specific airway conductance (sGaw). We investigated the relationship of this parameter with other pulmonary function tests (PFT) parameters such as diffusion capacity for carbon monoxide (DLCO) [15] and the ratio of residual volume (RV) to total lung capacity (TLC), an index that quantifies air-trapping in COPD subjects [16]. We also associated this parameter with computed tomography (CT)-based densitometric quantification of emphysema, a condition of severe airflow limitation in COPD when lung parenchyma is permanently destroyed [17]. Finally, we studied how changes in the resistance parameters affect the shape of the MEFVC.

Study Population
We used data of 462 COPD individuals and healthy smokers from the Leuven COPD cohort. All subjects were Caucasian with an age range between 50 and 90 years, and with a smoking history of at least 15 pack-years. Individuals with a suspicion or diagnosis of asthma, with exacerbations due to COPD within the last 6 weeks of enrolment or with a diagnosis of other respiratory diseases were excluded. These individuals performed complete PFT including post-bronchodilator spirometry, body-plethysmography, and diffusion capacity at the time of enrolment, while a CT scan was carried out within one year of enrolment. COPD was defined based on a post-bronchodilator ratio of FEV1/FVC < 70 % [3], while the severity airflow limitation was based on FEV1 expressed as a percentage of healthy reference value (FEV1 %pred) as described in the Global Initiative of Chronic Obstructive Lung Diseases (GOLD) [14]. The study design of the Leuven COPD cohort can be found at www.clinicaltrials.gov (NCT00858520).

Spirometry
Subjects performed spirometry using a MasterScreen Pneumo spirometer (available from Vyaire Medical Inc., Illinois, USA). Post-bronchodilator spirometry manoeuvres, the ones with the highest sum of FEV1 and FVC among three acceptable and repeatable manoeuvres, were used for the study [4].

Whole-body Plethysmography
Whole-body plethysmography was carried out using MasterScreen Body (available from Vyaire Medical Inc., Illinois, USA). Airway resistance (Raw), defined as a ratio of the difference of alveolar pressure and mouth pressure (or total driving pressure) to flow rate, was determined by a surrogate ratio of specific resistance (sRaw, the slope of shift volume to flow rate) to functional residual capacity (FRC) under tidal breathing conditions. Specific conductance (sGaw) was determined as a ratio of inverse of Raw to FRC. Finally, lung volumes including residual volume (RV) and total lung capacity (TLC) were also measured, and their ratio was defined as an index for air-trapping [1].

Diffusion Capacity
Diffusing capacity (DLCO) was measured by the single-breath carbon monoxide gas transfer method and corrected for alveolar ventilation but not hemoglobin concentration [18].

Computed Tomography Densiometry
CT scans were obtained in a routine setting using different multi-detector row scanners with different acquisition parameters. Patients were examined in supine position and the scans were conducted at end-inspiration. The severity of emphysema was calculated using the percentage of total voxels with an X-ray attenuation value below -950 HU, and using the attenuation value of the 15th percentile along a histogram of voxels [17].

Deflating Balloon Model of Forced Expiration
During forced expiration, a person empties his lungs forcefully from total lung capacity (TLC) until a residual volume (RV) is reached [4]. As such, one can draw an analogy with an elastic balloon that deflates from a volume of TLC until RV (Figure 1) or from a volume of FVC until 0 (since FVC is the difference between TLC and RV), under the influence of a driving pressure. At a given time t, we denote the volume of the balloon as x, and the velocity of deflation as .
x. If we ignore the effect of intra-thoracic gas compression [19], then the volume of air coming out of the balloon is FVC-x and the magnitude of airflow is .
x, which we assume to be laminar. Furthermore, we considered a coordinate system such that all vectors in the direction of deflation are negative. Model of a balloon that deflates from total lung capacity (or a volume of FVC) at time time (t) = 0 until residual volume (or a volume of 0) at t = t n with a peak airflow resulting at t = t 1 . The deflation process is driven by alveolar pressure, which arises due to a combination of expiratory muscle effort and elastic recoil of the lungs; while it is opposed by a resistive pressure, which depends on the width of the outlet.
First, we write the equation of motion of the deflating balloon as: where I ..
x is the inertial pressure, P alv is the alveolar pressure, and P res is the resistive pressure. Equation (1) implies that the inertia of a deflating balloon is the result of P alv , which is a driving pressure and P res , which opposes deflation [9,20]. The initial conditions of Equation (1) at t = 0 are: x(0) = 0 (2) and final conditions at t = t n (time at the end of forced expiration) are: At the beginning of the forced expiration, flow increases from zero to peak expiratory flow (PEF) within a very short interval t 1 (0 < t 1 << t n ). While the airflows in this phase are greatly influenced by expiratory muscle effort [21], the airflows during t 1 until t n , generally after 50% of expired volume are known to be effort-independent [22]. To avoid effort-dependent biased estimates, we drop the former phase and consider only the later phase (t 1 < t < t n ) by modifying the initial conditions as: where ∆v is the volume of air expired between 0 and t 1 , and can be written in terms of measured airflow F(t) as: In Equation (1), the driving pressure P alv is a sum total of pleural pressure (P pl ), which arises purely out of expiratory muscle effort, and elastic recoil pressure (P st ), which is a restoring pressure generated by the inherent elasticity in the lungs [9]. According to the data produced by Agostani et al. [23], P alv decreases non-linearly with lung volume from a peak until zero during forced expiration. However, we make a simplifying assumption by representing P alv as a linear function of x with a coefficient k. Thus, The resistive pressure P res in Equation (1) is determined by the width of the outlet of the balloon ( Figure 1a). P res decreases as the width of the outlet increases. We can represent it as a viscous damping pressure that dissipates the energy supplied to the system by P alv , and denote it as a proportion of deflating velocity .
x with a damping coefficient c: x.
Putting Equations (6) and (7) in Equation (1) and rearranging the terms, we arrive at a familiar second-order linear differential equation describing free vibrations of a single-degree-of-freedom mass-spring-damper system: ..
is a dimensionless coefficient representative of resistance in the balloon while ω 2 = k I is a coefficient of driving pressure per unit inertance I. It is interesting to note that ζ turns out to be an effort-independent parameter unlike ω. It quantifies the resistance of the small airways in the lungs. In the following section, we describe a method to estimate ζ and ω.

Estimation of ζ
Lung volumes decrease from TLC to RV during FE, unless noisy artifacts such as cough are present which we neglect in this study. So, we can consider the system in Equation (8) as an overdamped system, which implies ζ > 1 [24]. In an overdamped system, the solution to Equation (8) is as follows: with There are four unknown terms, C 1 , C 2, ζ, and ω in solution 9, which can theoretically be determined using the two initial conditions (Equation (4)) and the two final conditions (Equation (3)). Applying the initial conditions (Equation (4)), we can rewrite C 1 and C 2 in terms of ζ and ω as follows (see Appendix A): We now use an optimization strategy to obtain the unknown parameters, ζ and ω. There are two reasons for this step; the first is that final conditions (Equation (3)) could have anyway been solved numerically and the second, the estimate of the parameter ζ may not be robust if only a single data point at time t n is considered. Therefore, we consider the flow and volume of the entire manoeuver from t 1 until t n and express our cost function as follows: where V(t) and F(t) are measured flow and volume of air during forced expiration. Then, we formulate the optimization as a minimization problem as follows: We can solve the above optimization problem for ζ and ω using any suitable method. In our study, we used differential evolution as implemented in Python by the original authors [25]. We defined a search space of [1,5] and (0, 5] for ζ and ω, respectively. We wrote a Python script to automate the calculation of ζ and ω from flow-volume measurements of all subjects [26].

Statistical Analysis
We compared ζ between COPD and healthy smokers using a one-sided t-test. We studied ζ across different GOLD stages using one-way ANOVA followed by a post-hoc test using Tukey's honest significant difference test. We explored correlation of ζ against FEV1, FVC, FEV1/FVC, Raw, sGaw, DLCO, RV/TLC, and CT densitometry using Pearson's correlation coefficient. The significance level was set at 0.01 and all analyses were carried out in R Studio software [27]. Values were expressed as mean ± standard deviation (SD), unless specified otherwise. We used one-way ANOVA and logistic regression tests to assess group differences for continuous and categorical baseline variables.

Study Population
The baseline characteristics of the cohort are shown in Table 1. This cohort comprised of 63 GOLD I, 111 GOLD II, 92 GOLD III, 61 GOLD IV, and 135 healthy smokers. In this cohort, FEV1/FVC decreased with severity of airflow limitation from healthy (74 ± 4%) to GOLD I (64 ± 4%) until GOLD IV (30 ± 6%). In addition, Raw increased with severity of airflow limitation while other PFT parameters such as DLCO and air-trapping index (RV/TLC) also showed expected patterns. The extent of emphysema as quantified in CT by % below −950 HU and 15th percentile HU increased considerably with airflow limitation.

Model Fit
The fit of the deflating balloon model was very good as evident by a low mean squared error (MSE) for both volume (0.01 ± 0.01) (R 2 = 0.99 ± 0.01)) and flow (0.02 ± 0.02) (R 2 = 0.98 ± 0.03). Figure 2 shows the fit of model flow-volume against actual measurements for a healthy and a COPD subject.

Raw and sGaw
ζ and Raw were positively correlated (r = 0.55, p < 0.0001), and both increased with severity of airflow limitation, thus confirming our hypothesis. Furthermore, it was negatively correlated with sGaw (r = −0.60, p < 0.0001), which was indeed a very logical observation as sGaw is inversely proportional to Raw. The trends of ζ, Raw, and sGaw against severity of airflow limitation can be seen in Figure 3.

CT Densitometry
We considered a cohort of COPD subjects only (N = 327), since emphysema is considered a trait of COPD only. We observed that ζ was significantly correlated with 15th percentile HU (r = −0.34, p < 0.0001) and with % below -950 HU (r = 0.40, p < 0.0001), although the magnitude of correlations was weak. Interestingly, the known parameters of airway resistance, Raw and sGaw, did not correlate with any of the densitometry indices (p > 0.05).

ζ simulation Studies
We studied how the shape of MEFVC changes when ζ changes from a baseline value. For a given manoeuver, we simulated MEFVCs by changing the baseline ζ by 20% while keeping ω as constant. Figure 4a,b shows the simulation of model flow and volume in the case of a healthy and a COPD individual for the baseline models shown in Figure 2a,b, respectively. As expected, an increase in ζ resulted in decreased airflows with a lower simulated FEV1 and FVC, while a decrease in ζ led to increased airflows with a higher simulated FEV1 and a slightly elevated FVC. The slight increase in FVC is seen as the model slightly underestimates the observed FVC. Furthermore, we also observed an increase in the concavity of the MEFVC as ζ increased in both the simulation studies.

Discussion
In this study, we described the dynamics of forced expiration in spirometry by using a model of a deflating balloon. With this model, we proposed a methodology to estimate a parameter ζ representative of airway resistance in the lungs. We proved our hypothesis by demonstrating that ζ was significantly (p < 0.0001) greater in COPD than in healthy smokers. It increased significantly with the severity of airflow limitation in COPD and it significantly correlated with Raw and sGaw, which are known indices of airway resistance. ζ also showed significant correlations with DLCO and air-trapping (RV/TLC), and with CT densitometric indices of emphysema in COPD patients. Moreover, simulation studies showed that an increase in ζ resulted in decreased airflows, and therefore a reduced FEV1 and FVC, which further validates our claim that ζ quantifies airway resistance.
The parameter ζ is not dependent on expiratory effort, but rather, it stems from the narrowing of the small airways (<2 mm diameter). Narrower airways result in increased airway resistance that is reflected by an elevated ζ. In COPD, narrowing of the small airways during forced expiration is a result of a multitude of factors such as loss of elasticity in lung parenchyma leading to dynamic.
Airway collapse [9]; and airway abnormalities such as cellular inflammation [28], smooth muscle contraction [29], mucosal thickening [30], and fibrotic retraction [31]. Finally, a higher value of ζ also indicates the severity of airflow limitation, which itself is a functional consequence of a complex interaction between airway abnormalities and destruction of parenchyma elasticity.
Our method is the first to estimate airway resistance from forced expiration in spirometry. The use of a second-order linear differential equation in our method further validates the existence of second-order transfer functions to describe the dynamics of forced expiration [32]. Nevertheless, the application of second-order differential equations to describe respiratory mechanics is not a new concept [33]. In fact, the estimation of airway resistance during mechanical ventilation [34], body-plethysmography, and the interrupter technique is made possible by a simplification of these equations under tidal breathing conditions [7]. However, our method has an important distinction in that it does not require a measurement of alveolar pressure. Since spirometry does not measure alveolar pressure, we make a reasonable assumption on the variation of the latter (Equation (6)) that nonetheless, provides a good estimate of airway resistance. We believe that this is an important development in the field of respiratory mechanics because until now, the consensus has been that measuring alveolar pressure is indispensable for the calculation of airway resistance [7,35,36].
The clinical usefulness of ζ can be illustrated by the fact that it quantifies an important physiological phenomenon using spirometry alone, which is a simple but widely used pulmonary function test.
In fact, our results show that ζ is superior to its counterparts, Raw and sGaw, in their association with CT densitometry, which is considered as the most direct method for assessing the severity of emphysema [37]. It is also better correlated with DLCO, which is an index for alveolar destruction and loss of capillary bed in COPD [38]. Compared to the interrupter technique, which estimates airway resistance using measurements obtained at a single instant of time under tidal breathing [7], we achieve a more robust value of airway resistance as we consider the entire manoeuver of forced expiration. This motivates the application of our methodology, especially, in primary care where spirometers are more abundantly available than bulky and expensive equipment such as a body-box or FOT machine. A major limitation of our methodology is that it does not shed any light on expiratory muscle effort and elastic recoil pressure, which are the components of alveolar pressure that drive forced expiration. The term ω proves to be of little use as it masks these underlying components. An interesting result would have been to estimate elastic recoil coefficient and then, to correlate it with elasticity loss in emphysema. However, such a step would have required pleural pressure measurements using an invasive and impractical esophageal balloon [39]. Another drawback is that ζ may not be sensitive to upper airway resistance. During forced expiration, the emptying of the upper airways occurs in the beginning, a phase that was neglected in our model (Equation (4)). It is also worth mentioning that validation of ζ by comparing it against Raw may not be completely appropriate. The latter quantifies total resistance during tidal breathing when resistance of the larger airways and viscoelastic properties of lung parenchyma account dominate over smaller airways resistance [40,41]. A better validation of ζ would have been to compare it against the difference of FOT resistance at 5 Hz and 20 Hz (R 5 -R 20 ), metric that reflects small airway resistance [42]. Finally, our study was focused on airflow limitation associated with COPD only. A more inclusive study is required to evaluate the efficacy of ζ by considering other obstructive cohorts such as asthma, asthma-COPD overlap (ACO), and obstructive asthma.

Conclusions
In this study, we describe the dynamics of forced expiration in spirometry using a deflating balloon and using this model, we proposed a methodology to estimate ζ, a dimensionless and effort-independent parameter quantifying airway resistance in the lungs. ζ correlated significantly with Raw and sGaw, which are traditional indices of airway resistance from body-plethysmography. It increases with the severity of airflow limitation as expected from an airway resistance, and in some cases, even showing better associations with DLCO and CT densitometry of emphysema than Raw and sGaw. Moreover, simulation studies showed that an increase in ζ resulted in decreased airflows, and therefore a reduced FEV1 and FVC. Thus, we conclude that ζ estimated from forced expiration in spirometry reflects airway resistance in the lungs. In primary care, ζ may be clinically useful, as spirometry is more widely available than traditional but expensive methods to measure airway resistance such as body-plethysmography or FOT. In the future, we aim to validate ζ in a more inclusive cohort of obstructive airway diseases with longitudinal follow-up, as well as study the association of ζ against airway resistance parameters obtained from the interrupter technique and FOT.

Patents
Based on the present study, we have filed a patent titled "Methods and apparatus for determining airflow limitation" (PCT/EP2019/059903) on 11 May 2019. similarly, multiplying (A1) by s 1 throughout, subtracting (A1) from it and rearranging for C 2 , we get 11b