REDUCED MODEL FOR FORCED EXPIRATION AND ANALYSIS OF ITS SENSITIVITY

REDUCED MODEL FOR FORCED EXPIRATION AND ANALYSIS OF ITS SENSITIVITY Janusz Mroczka and Adam G. Polak Chair of Electronic and Photonic Metrology, Wroclaw University of Technology ul. B. Prusa 53/55, 50-317 Wroclaw, Poland

Abstract: Pathological changes in the lung modify the shape of the flow-volume curve registered during forced expiration. Computational models allowing simulation of the test results are too complex for the estimation of their parameters. In this study a complex model was reduced by introduction of the functions scaling airway properties. Influence of the individual parameters of the reduced model on the flow-volume curve was evaluated by means of the sensitivity analysis. The conclusion is that the parameters of the scaling functions and elastic properties of lung tissue affect the measured data most significantly and that the descending part of the curve should be used to assess them. Copyright © 2006 IFAC Keywords: Biomedical systems, Modelling, Model reduction, Sensitivity analysis, Inverse problem.

1. INTRODUCTION Forced expiration is the most common test of the lung function. It has been shown that the registered maximal expiratory flow-volume (MEFV) curve is effort-independent and simultaneously sensitive to respiratory disorders (Hyatt, et al., 1958). The observed connection between the respiratory state and the shape of the MEFV curve encourage to elaborate a method for the quantitative evaluation of lung parameters using the flow-volume data. Solving this so-called inverse problem requires application of a mathematical model that couples the respiratory mechanical properties with the MEFV curve. Three useful approaches to model and simulate the forced expiration, incorporating morphological data and main physiological phenomena, can be found in the literature (Lambert, et al., 1982; Elad, et al., 1988; Polak, and Lutchen, 2003). They are, however, too complex for the estimation of the respiratory system parameters. The early study by Lambert (1984) exposed these difficulties and finally failed in solving the problem. Recently, Lambert and co-workers were

successful predicting individual flow-volume curves with the use of their model (Lambert, et al., 2004; Lambert, and Beck, 2004). It has been, however, proved that the set of parameters determining maximal expiratory flow is much bigger than the maximal airway areas adjusted in those studies. The aim of this study was to transform the chosen complex computational model for the forced expiration into a model characterized by a reduced number of free parameters and then to analyse its sensitivity. These activities are the first step towards solving, if possible, the inverse problem in the forced expiration. Such an inverse model would have been very useful clinically (e.g. in case of asthmatic or COPD patients) giving the possibility to retrieve airway and/or lung tissue properties from the MEFV curve. The paper is organised as follows. At the beginning (section 2) the complex model chosen for further analysis is briefly described. Then a method of its reduction, engaging sigmoid functions that rescale

airway mechanical properties distributed along the bronchial tree, is proposed. Introduction of the scaling functions results in the reduced number of the free model parameters. Next, the sensitivity of the reduced model is analysed. Section 3 presents results of the model reduction and sensitivity analysis, and conclusions exposing the most influential model parameters as well as other insights following the study are drawn in section 4. 2. METHODS 2.1 Complex computational model A complex model with a symmetrically bifurcating airway structure has been chosen since it describes airway morphology correctly, allows independent adjustment of their mechanical properties and simultaneously is computationally efficient (comparing to other approaches). This model has been discussed extensively elsewhere (Lambert et al., 1982; Polak, 1998) and will be quoted briefly here. Its bronchial tree structure, following the geometry proposed by Weibel (1963), consists of 24 symmetrical generations, each consisting of identical airways arranged in parallel, with the trachea categorised as generation 0 (Fig. 1). The mechanical properties of the airways can be specified independently for each generation. First of all, they include parameters describing the dependence of the airway lumen area A on transmural pressure Ptm , i.e.the tube law; Lambert et al., 1982): ⎧⎪ Amα 0 (1 − Ptm P1 )− n , Ptm ≤ 0, (1) A(Ptm ) = ⎨ −n ⎪⎩ Am 1 − (1 − α 0 )(1 − Ptm P2 ) , Ptm > 0, 1

[

2

]

where Am is the maximal lumen area, α0 is the normalized area (i.e. A/Am) at zero transmural pressure, and n1 and n2 are shape-changing scalars. P1 and P2 are pressure asymptotes given by: P1 = n1α 0 α 0′ ,

(2)

P2 = n2 (α 0 − 1) α 0′ ,

where α'0 denotes the slope of the A-Ptm curve (i.e. dA/dPtm) at Ptm = 0 . Additionally, intrapleural bronchi lengths alter with the lung volume (VL): 1

⎛ V + VL ⎞ 3 ⎟⎟ , l (V L ) = l1.5 ⎜⎜ t ⎝ Vt + 1.5 ⎠

(3)

where l1.5 is the length of an airway at lung volume of 1.5 dm3, and Vt is lung tissue volume. Weibel has proposed exponential equations to represent the airway diameter and length dependences on the generation number (Weibel, 1963).

Q Ptm l

A

g=0 (trachea) g=1

g=2

Fig. 1. Schematic representation of the first 3 generations of the bronchial tree: Q is expiratory flow, g is the generation number, A is the airway lumen area, l is the airway length, and Ptm is transmural pressure. Another model feature is the use of a nonlinear characteristics for the lung recoil (Bogaard, et al., 1995; see Fig. 2):

⎧VL − V0 , VL ≤ Vtr , ⎪ C E ⎪ Pst (VL ) = ⎨ ⎪Vm − Vtr ln⎛⎜ Vm − Vtr ⎜V −V ⎪ C st L ⎝ m ⎩

⎞ Vtr − V0 ⎟⎟ + , VL > Vtr , C st ⎠ (4)

where Pst is static lung recoil pressure, VL is lung volume, Vm and V0 are maximal and minimal lung volumes, Vtr is transition volume, and Cst is lung compliance at zero recoil pressure. During simulations, pressure drops along compliant airways are calculated by numerical integration of the formula for the pressure gradient (including both wave-speed flow limitation and viscous pressure losses) at every lung volume analysed (Lambert et al., 1982). Since the bronchial tree is fully symmetrical, the calculated pressure drop is simultaneously the pressure loss along a generation that the given airway belongs to. Integration is interrupted at each junction between airway generations for the evaluation of convective acceleration of gas according to the Bernoulli law. The calculated total pressure drop across the bronchial tree is then summarized with the pressure loss in the upper airways which depends on the scaling coefficients Ru and r (Polak, 1988), and equated to the driving pressure produced by the lung elastic recoil and expiratory muscle. The driving pressure (Pd) is the function of expiration time and lung volume (Polak, 1988): ⎛ V − RV ⎞ Pd (t ,VL ) = Pm (1 − e −t τ )⎜ L ⎟, ⎝ VC ⎠

(5)

where Pm is maximal expiratory pressure, t stands for time of expiration, τ is a time constant of expiratory muscle, RV denotes residual volume, and VC is vital capacity.

A total amount of 159 parameters is used in the complex model: 144 describe properties of 24 airway generations, 11 parameters for the lung recoil characteristics and other features of the respiratory system (including VC), and 4 physical or semiphysical constants. Using known forced vital capacity (FVC) instead of VC, 154 of the parameters should be regarded as the free ones. The Lambert model for the forced expiration has been proven to be a valuable tool in the analysis of the respiratory system for the last two decades. Its computational abilities have been enhanced (Polak, 1998) by applying a method of succeeding approximations when finding maximal airflows in the quasi-static conditions at the consecutive lung volumes in the whole range of vital capacity. 2.2 Reduction of the complex model Proposed by Weibel equations expressing the relationship between airway dimensions and the generation number follow the fact that the mechanical properties of neighbouring airways are not independent. In the complex model, airway elastic properties, described by Eqs 1-2, include five parameters specified for each generation. This huge number of parameters can be reduced by the use of scaling functions that would capture changes in distribution of the airway mechanical properties along the bronchial tree (Habib, et al, 1994; Polak, and Mroczka, 1998). Such changes may result from possible pathology or intersubject variability. Let us denote the baseline values of the airway parameters generally as θ* and the rescaled ones as θ, then the following rescaling, depending on the generation number (g), can be proposed: Am (g ) =

2 Am* (g ) , 1 + exp( pa1 g + pa 2 )

(6)

α 0 (g ) =

1 + exp(0.33 g − 2.4) * α 0 (g ) , 1 + exp( pz1 g + pz 2 )

(7)

α 0′ (g ) =

2 α m′* (g ) , 1 + exp( pc1 g + pc 2 )

(8)

n1 (g ) = ( pn11 g + pn12 ) ⋅ n1* (g ) ,

(9)

n2 (g ) = ( pn 21 g + pn 22 )n2* (g ) ,

(10)

where pa1, pa2, pz1, pz2, pc1, pc2, pn11, pn12, pn21, pn22 are the scaling parameters. Two features, Am and α'0, are rescaled by sigmoid functions with two degrees of freedom and values between 0 and 2 (Eqs 6 and 8). Since α0 cannot exceed 1 (it is a relative parameter), the form of the numerator in Eq. (7) has been chosen to satisfy this relationship for all generations.

Table 1: Parameter values of the reduced model for the normal and diseased lung. H was presumed as equal to 176 cm, and VC as equal to 5.5 dm3 (normal lungs) and 4 dm3 (diseased lungs) Symbol Am(0) (mm2) α0(0) α'0(0) (kPa-1) n1(0) n2(0) pa1 pa2 pz1 pz2 pc1 pc2 pn11 pn12 pn21 pn22 RV (dm3) ∆V0 (dm3) ∆Vtr (dm3) ∆Vm (dm3) Vt (dm3) Cst (dm3kPa-1) Pm (kPa) τ (s) Ru r

Parameter value Normal lung Diseased lung 237 237 0.882 0.882 0.11 0.11 0.5 0.5 10 10 0.0 0.0 0.0 0.0 0.33 0.36 -2.4 -2.4 0.0 0.04 0.0 1.1 0.0 0.0 1.0 1.0 0.0 0.0 1.0 1.0 1.5 4.0 0.0 0.0 2.5 1.0 5.8 4.1 0.9 0.9 3.5 6.0 24 24 0.2 0.2 0.11 0.11 1.68 1.68

To include, however, their intersubject variability, the airway lengths at lung volume of 1.5 dm3 (l1.5) are related to the known patient’s height (H) (Habib, et al, 1994; Polak, and Mroczka, 1998), assuming that the baseline data derive from a man of 176 cm: l1.5 (g ) =

H * l1.5 (g ) . 176

(11)

The effect of variation in residual volume (RV), seen as a parallel shift of the recoil pressure – lung volume relationship, on the MEFV curve is small (Lambert, 1984), and it has been investigated in this study by introducing the following reparameterisation into the V0 = RV + ∆V0 , Vtr = RV + ∆Vtr and model: Vm = RV + ∆Vm , where ∆V0, ∆Vtr, and ∆Vm are new “shift” parameters applied instead of V0, Vtr, and Vm used originally. All the parameters of the reduced model are listed in Table 1. 2.3 Sensitivity analysis The sensitivity analysis enables to figure out if the model is suitable for estimation of any of its parameters, which of them, what accuracy may be excepted and which part of data should be used to

∂y , ∂θ

(12)

where y denotes a vector of output data (expiratory flow in our case) and θ is a vector of the model parameters. Since the parameters may differ in the magnitude by several orders, so may the sensitivity vectors (i.e. the columns in matrix X). It is better to use the normalized sensitivity XN , X N = X ⋅ diag(diag(X T X )) , −1 2

(13)

to observe the connection of the individual parameters with specific fractions of the output data (Thomaseth, and Cobelli, 1999; operator diag transforms a matrix into a vector including the matrix diagonal elements or a vector into a diagonal matrix). Another advantage following the sensitivity analysis is the possibility to evaluate correlations between the sensitivity vectors. High correlation indicates that the parameters influence the measured model output in a very similar way and their values cannot be estimated precisely from noisy data. The sensitivity of the reduced model derived in section 2.2 was analysed with the intention to assess the model abilities in solving the inverse problem and to plan further investigations. The sensitivity vectors could not be derived analytically, so they were computed numerically. 2.4 Values of the parameters The complex model was simulated using the published baseline data (Lambert, 1984; Polak, 1998). Being aware of the dependence of the nonlinear model analysis on the considered point in the parameter space, the investigations were performed in two cases: for the normal lung with baseline values of the parameters and for the diseased lung characterized by hyperinflation, loss of lung elastic recoil and airway obstruction, the symptoms found in asthmatic and emphysematic patients (Gelb, and Zamel, 2000; Baldi, et al., 2001). The relevant lung static recoil characteristics are shown in Fig. 2. Taking into account the quantitative effects of airway constriction (Gelb, and Zamel, 2000; Morlion, and Polak, 2005), the obstruction was simulated by adequate alternation of the scaling parameters pz1, pz2, pc1, and pc2. This modification reduced α0 twice in generation 23 (leaving it unchanged in the generation 1), and α'0 was decreased from two-fold in the generation 1 up to four-fold in the generation 23. The resulting MEFV curves are presented in Fig. 3.

In the reduced model, the values of the intrapleural airway properties (i.e. generations 1 to 23) are calculated with the use of the functions rescaling the baseline values, so the model possesses 25 free parameters: 5 of the trachea (generation 0), 10 of the rescaling functions, 2 of upper airway resistance, 4 of the lung recoil characteristics, and 4 others. For the parameter values of the normal lung (Table 1), the complex and reduced model outputs overlap (Fig. 3, solid line), since the scaling functions equal unity for all generations. The outputs of the models would overlap also in case of the diseased lung, if the airway parameters of the complex model were the same as baseline values multiplied by the rescaling functions with parameter values given in Table 1 (the diseased lung column). The impact of the individual parameters on the flowvolume data is shown in Figs 4 and 5. Parameters describing the extrapleural airways and expiratory muscle determine only the ascending part of the MEFV curve (see Fig. 3), both for the normal and diseased lung (Fig. 4). On the contrary, parameters of the intrapleural airways and lung recoil of the normal as well as diseased lung affect expiratory flow in the whole range of VC (Fig. 5). 8 7

Lung volume (dm3)

X=

3. RESULTS

6 5 4 3

Normal lungs Diseased lungs

2 1

0

0.5

1

1.5

2

Static recoil pressure (kPa)

2.5

3

Fig. 2. Static recoil pressure – lung volume curves of the normal and diseased lung used in the study. 9 8

Normal lungs Diseased lungs

7

Flow (dm3s-1)

this end. The primary point of the analysis consists in the determination of the output sensitivity (X) to the parameters:

6 5 4 3 2 1 0

0

1

2

3

4

Expired volume (dm3)

5

6

Fig. 3. Maximum expiration flow-volume curves of the normal and diseased lung.

(A)

0.5

0.15 0.1

0.3

Normalized sensitivity

Normalized sensitivity

0.4

0.2 0.1 0 -0.1 -0.2 -0.3

(A)

0.2

0.05 0 -0.05 -0.1 -0.15 -0.2 -0.25

0

1

2

3

4

5

Expired volume (dm3)

-0.3

6

0

1

2

3

(B)

0.5

4

5

Expired volume (dm3)

6

(B) 0.25

0.4

Normalized sensitivity

Normalized sensitivity

0.2 0.3 0.2 0.1 0 -0.1 -0.2

0.15 0.1 0.05 0 -0.05 -0.1 -0.15 -0.2

-0.3 -0.4

-0.25 0

0.5

1

1.5

2

2.5

Expired volume (dm3)

3

3.5

4

Fig. 4. Normalized sensitivities to parameters (XN) characterizing extrapleural airways and expiratory muscle of normal (A) and diseased lungs (B). The vertical dashed line indicates the volume related to peak expiratory flow (PEF). Some of the sensitivity vectors of the reduced model appeared strongly correlated (absolute values of the relevant correlation coefficients were greater than 0.99), and the number of the correlated pairs was bigger in case of the diseased lungs. 4. DISCUSSION AND CONCLUSIONS Application of the functions scaling airway properties according to their generation number goes back as far as the first Weibel’s models of the bronchial tree geometry (Weibel, 1963). Since then such an approach has been practised many times demonstrating that the properties of the neighbouring generations are related to each other. In the present work, the proposed functions rescale airway baseline properties rather than follow strictly the former approaches. Their form implies that the reduced model will efficiently mimic a serial distribution of pathological changes in the bronchial tree. Sigmoid shape of the functions rescaling Am and α'0, used also by others (Habib, et al., 1994), enables capturing nonlinear decrease or increase of these parameters along the generations. Additionally, more complicated form of the function scaling α0 allows imitating its non-monotonic variations.

0

0.5

1

1.5

2

2.5

Expired volume (dm3)

3

3.5

4

Fig. 5. Normalized sensitivities to parameters (XN) characterizing intrapleural airways and lung recoil of normal (A) and diseased lungs (B). Recent experimental results show that airway obstruction during metacholine challenge may have just such serial character (Lambert, and Beck, 2004). Nevertheless, the case of selective constriction of individual airway generations, investigated by Lambert and Beck (2004), cannot be well described by this model. The sensitivity analysis (Fig. 4) has revealed that the descending part of the MEFV curve is insensitive to the parameters describing the extrathoracic airways and expiratory muscle. As shown by others, this portion of flow-volume data is set by the flow limiting mechanism and the resulting expiratory flow does not depend on the properties of the airways situated downstream to the site of flow limitation (Dawson, and Elliott, 1977). Simultaneously, this part of the curve is effort-independent (Hyatt, et al., 1958) and the parameters describing muscle action cannot influence its shape. The conclusion is that the data following peak expiratory flow (PEF) should be used to estimate the parameters of the intrapleural airways and the lung recoil characteristics. This is in line with the practice of earlier trials (Lambert, 1984; Lambert, and Beck, 2004). The proposed form of the reduced model implies that all processes modifying the MEFV curve shape will be seen as variations in the free parameters during identification of the model, and then interpreted as a

change of lung recoil or as serially distributed alternations of airway mechanics. Prior diagnosis and knowledge of occurring pathological processes will predestine the use of the model in some diseases or, on the contrary, will devalue it in other disorders. In the latter case, however, it is still possible to propose another parameterisation of the reduced model, e.g. by the choice of other rescaling functions that will fit the pathogenesis better, and then to repeat the procedure applied in this study. Linear methods have been used to investigate the nonlinear model, so the results are correct only in the neighbourhood of the point in the parameter space, at which the model was analysed. It has been shown, however, that such linear approximations may yield quite accurate description (Yuan, et al., 1998), justifying this approach. Achieving the complete results would require the analysis to be done in the whole parameter space consisting of an infinitive number of points. To generalise the outcomes, we have performed investigations at two representative states of the respiratory system (the normal and diseased lung) with hope that common conclusions would be also true for other possible conditions. Analysis of the model sensitivity at other characteristic points unquestionably needs additional studies. On the other hand, the very precise quantitative outcomes seem to be unnecessary at this stage of investigations. Summarizing, this study has shown that the model for the forced expiration can be reduced to 25 free parameters describing intrapleural airway mechanics and the lung static recoil. Only these parameters influence the descending part of the MEFV curve and it should be chosen for the model identification. High correlations between some of the sensitivity vectors indicate that the reduced model cannot be properly identified and this reveals the need for further investigations. One of possible approaches is the selection of the most influential and less-correlated parameters for the estimation (Polak, 2001). This should be the next step towards elaboration of the inverse model for the forced expiration. REFERENCES Baldi, S., M. Miniati, C.R. Bellina, L. Battolla, G. Catapano, E. Begliomini, D. Giustini and C. Giuntini (2001). Relationship between extent of pulmonary emphysema by high-resolution computed tomography and lung elastic recoil in patients with chronic obstructive pulmonary disease. Am. J. Respir. Crit. Care Med., 164, 585-589. Bogaard, J.M., S.E. Overbeek, A.F.M. Verbraak, C.Vons, H.T.M. Folgering, Th.W. van der Mark, C.M. Roos, P.J. Sterk, and the Dutch CNSLD study group (1995). Pressure-volume analysis of

the lung with an exponential and linearexponential model in asthma and COPD. Eur. Respir. J., 8, 1525-1531. Dawson, S.D. and E.A. Elliott (1977). Wave-speed limitation on expiratory flow - a unifying concept. J. Appl. Physiol.: Respirat. Environ. Exercise Physiol., 43, 498-515. Elad, D., R.D. Kamm and A.H. Shapiro (1988). Mathematical simulation of forced expiration. J. Appl. Physiol., 65, 14-25. Gelb, A.F. and N. Zamel (2000). Unsuspected pseudophysiologic emphysema in chronic persistent asthma. Am. J. Respir. Crit. Care Med., 162, 1778-1782. Habib, R.H., R.B. Chalker, B. Suki and A.C. Jackson (1994). Airway geometry and wall mechanical properties estimated from subglottal input impedance in humans. J. Appl. Physiol., 77, 441451. Hyatt, R.E., D.P. Schilder and D.L. Fry (1958). Relationship between maximum expiratory flow and degree of lung inflation. J. Appl. Physiol., 13, 331-336. Lambert, R.K. (1984). Sensitivity and specificity of the computational model for maximal expiratory flow. J. Appl. Physiol.: Respirat. Environ. Exercise Physiol., 57, 958-970. Lambert, R.K. and K.C. Beck (2004). Airway area distribution from the forced expiration maneuver. J. Appl. Physiol., 97, 570-578. Lambert, R.K., R.G. Castile and R.S. Tepper (2004). Model of forced expiratory flows and airway geometry in infants. J. Appl. Physiol., 96, 688692. Lambert, R.K., T.A. Wilson, R.E. Hyatt and J.R. Rodarte (1982). A computational model for expiratory flow. J. Appl. Physiol.: Respirat. Environ. Exercise Physiol., 52, 44-56. Morlion, B. and A.G. Polak (2005). Simulation of lung function evolution after heart-lung transplantation using a numerical model. IEEE Trans. Biomed. Eng., 52, 1180-1187. Polak, A.G. (1998). A forward model for maximum expiration. Comput. Biol. Med., 28, 613-625. Polak, A.G. (2001). Indirect measurements: combining parameter selection with ridge regression. Meas. Sci. Technol., 12, 278-287. Polak, A.G. and K.R. Lutchen (2003). Computational model for forced expiration from asymmetric normal lungs. Ann. Biomed. Eng., 31, 891-907. Polak, A.G. and J. Mroczka (1998). A metrological model for maximum expiration. Measurement, 23, 265-270. Thomaseth, K. and C. Cobelli (1999). Generalized sensitivity functions in physiological system identification. Ann. Biomed. Eng., 27, 607-616. Yuan, H., B. Suki and K.R. Lutchen (1998). Sensitivity analysis for evaluating nonlinear models of lung mechanics. Ann. Biomed. Eng., 26, 230-241. Weibel, E.R. (1963). Morphometry of the Human Lung. Springer, Berlin.