- Email: [email protected]
Innovative model to simulate exhalation phase in human respiratory system
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 4 ( 2 0 1 1 ) 300–305
journal homepage: www.intl.elsevierhealth.com/journals/cmpb
Innovative model to simulate exhalation phase in human respiratory system Tommaso Sbrana a,∗ , Alberto Landi b , Giosuè Angelo Catapano c a b c
Centro Interdipartimentale di Ricerca E. Piaggio, Faculty of Engineering, University of Pisa, Italy Department of Energy and Systems Engineering, University of Pisa, Italy Tuscany Foundation ‘Giuseppe Monasterio’ CNR, Pisa, Italy
a r t i c l e
i n f o
a b s t r a c t
Article history:
In this paper, we present a mathematical model, which mimics the bronchial resistances
Received 13 November 2010
of human’s lung in an expiratory act. The model is implemented in Matlab. The inputs that
Received in revised form
are used in this model derive from spirometry test. This model is able to study a physiologic
22 January 2011
condition, a pathologic one and the patient’s follow up after drug treatment. We split our
Accepted 4 February 2011
study into two parts. The first one focuses the analysis on the gas fluido dynamic inside of the respiratory pathways. The second part takes care of the pressure equilibrium in the
Keywords:
exchange zone. We use the outputs that derive from the second subsystem to solve the
Model
Bernoulli’s equation of the first part. The model was validated with data provided from
Lung
“Clinical Physiology Institute” of CNR and G. Monasterio Foundation of Pisa.
Respiratory system
© 2011 Elsevier Ireland Ltd. All rights reserved.
COPD Exhalation phase
1.
Introduction
Mathematical models in medicine are often used to evaluate the information provided by diagnostic tests. There are several advantages in this approach, like the possibility to test a hypothesis, to mimic a pathologic condition and/or a drug action avoiding side effects in patients. From the point of view of a model builder two different states have to be usually considered: basal or pathological conditions. The key problem is: how to move from a model emulating average conditions of the physiologic system in a basal state to a model able to represent the system in a pathological-state. Solving this problem is relevant from a diagnostic and prognostic point of view,
because it allows the physician to prove his diagnosis using the model, reducing delays in obtaining a correct diagnosis and avoiding invasive additional tests. In recent years many mathematical models have been proposed in the literature for describing human physiological systems. One of the main problems of many models is their inadequacy for helping medical diagnostic techniques and new therapeutic protocols: sometimes they need too many unknown parameters derived from an accurate description of the physiological system, so that they cannot be easily identified. On the other hand simplified models can be used for explaining and teaching how a biological system works [1], but they are usually too approximate and simple to help physicians. Many attempts have been reported in the literature to
Abbreviations: Pi , pressure “i”; Zi , height “i” from earth; hf , parameter relative to weight losses; L, length of the tube; D, diameter of cross section surface of the tube; v, medium rate of velocity of the air flow along the tube; f, frequency of the respiratory act; Vt , tidal volume; Vd , death space volume. ∗ Corresponding author. Tel.: +39 334 901 760. E-mail address: [email protected] (T. Sbrana). 0169-2607/$ – see front matter © 2011 Elsevier Ireland Ltd. All rights reserved. doi:10.1016/j.cmpb.2011.02.005
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 4 ( 2 0 1 1 ) 300–305
model and simulate human respiratory physiology. Although it is rather simplistic to classify these models, because of the complexity of respiratory physiology, different approaches can be distinguished, e.g., electrical models, mechanical models, morphometric models and mathematical models. Electrical models [2] describe the system using a Maxwell or a Firestone analogy in order to transform the physiological system into equivalent electrical circuits. The problem involved with this approach is that the transformed system is often too simplified and that a circuital study does not match clinical outputs that well correlate with medical findings. Mechanical models transform the physiological system into an equivalent mechanical system, based on pistons, gears and springs [3]. This approach has led to models able to obtain outputs well-correlated with doctor expectancy, but unable to mimic pathological conditions for their structural rigidity that makes unfeasible the emulation of cases different from the basic one. Morphometric models need inputs usually unknown and difficult to obtain without invasive tests [4]. Mathematical models can describe physiological system by using a datadriven (or black-box) approach [1], useful for fitting curves to input–output data when our understanding is limited and inadequate to produce a detailed model, or a system-based (or white-box) approach, formulated on the basis of comprehensive models able to adequately describe the underlying physiology. In the complex case of respiratory physiology a combination of two approaches (or gray-box), e.g. [5], is the optimal solution for joining the explicit structure of parts of model based on physiology and available data. In this paper we follow a reference model, i.e. a white-box approach, oriented to correlate spirometric data with the patient condition. It has to be highlighted that white-box approach needs some expedient to extrapolate all required parameters, as it will be explained in Section 2. Model presents study of bronchial resistances of human lung in an expiratory act and is composed by static and linear equations, which are implemented in Matlab M files, with the aim to obtain a simple method to correlate spirometric data with the patient condition.
2.
Materials and methods
The human respiratory system is divided into two parts: conductive and respiratory zone [6]. This division can be done, considering the exchange of oxygen between bronchial tubes and other tissues. It is interesting to notice that gas passes the endothelial membrane only in respiratory zone, while in conductive zone there is no gas exchange. In this hypothesis the model is split into two parts: one mimics the conductive part, and the other is deputed to the respiratory zone. The two parts are linked each other in a junction zone called “A” in Fig. 1. Some comments regarding this approach must be considered. Since the aim is to build a simple model able to correlate spirometric data with the patient condition, the choice to split into two parts the respiratory system was due to the attempt of maintaining essential functions of the system, reducing as much as possible the number of required anatomical data. Furthermore, this approach avoids a detailed – difficult to extract – description of the lung structure and of bronchial tubes. In our opinion a more detailed description of the system may
301
Fig. 1 – Design of the system. The two parts are depicted in different colors (blue for the conductive zone, violet for the respiratory one) and are divided by junction zone A. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)
lead to a model complexity that is excessively high, compared to the available data.
2.1.
Materials and methods: conductive zone
The first part of the model is considered as a rigid smoothwalled tube with a circular cross section surface, where the air can flow. The assumption of rigid wall derives from the absence of a pulsating movement of the wall, like that one of the arteries. The assumption of smooth surface derives from the presence of mucus along the bronchial tubes [6,7]. In physiologic systems the bronchial diameters become smaller and smaller from the mouth to the deepest regions of the bronchial tree and at the same time the total surface grows up. These phenomena are considered in the literature with an overestimation of the tube diameter [8] (an average value can be 2.2 cm2 ) and the total length of the conductive zone is estimated in the literature [8] as L = 30 cm. The circuit is studied by a classical fluid-dynamic approach [9]. Several studies have shown that the air flux in the conductive zone is turbulent [6] and the physiological system can be modeled using a fluid-dynamic analysis in a turbulence state. Therefore the Bernoulli’s expression is introduced: v2 v2 P2 P1 + g · Z1 + 1 = + g · Z2 + 2 + hf 2 2
(1)
Pi is the pressure “i”, Zi is the height “i” from earth, hf is the parameter relative to weight losses, vi is the velocity “i”, g is the gravity acceleration. The subscript 1 indicates values data correlated to the mouth, the subscript 2 considers values linked to the respiratory zone. The desired output of first sub-model is P1 . The height values are estimated for each case: the mouth is positioned 15 cm under the scalp. Z2 is correlated with the length of the tube, so that it can be assumed as a known data. The load losses (hf ) are composed by two elements: one due to friction (hf1 ) and the second to branches along the bronchial tree (hf2 ).
302
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 4 ( 2 0 1 1 ) 300–305
In (2) we express the relation between the two parameters hf1 and hf2 : hf = hf 1 + hf 2
(2)
From the Moody’s diagram it is possible to obtain Funning’s factor (f), for estimating the term hf1 [9]: hf 1 = 2 · f ·
L · v2 D
(3)
L is the length of the tube, D is the diameter of cross section surface of the tube, v is the medium rate of velocity of the air flow along the tube. The load losses due to the branches in a circuit in turbulence flow conditions [9] is:
hf 2 = ε ·
v2 2
(4)
where ε is the proportional coefficient typical of turbulence flows. The ε parameter can be estimated if the angle between two pathways is assumed of 60◦ . In physiology, this value represents an angles average between two pathways, as it is well described in [7]. This value implies ε equal to 0.3 [10]. Therefore, in order to obtain P1 , unknown parameters of the Bernoulli’s expression are the velocities v1 , v2 , and the pressure called P2 . The last one is linked to the junction region between the respiratory zone and the conductive one. The air flow is defined by: V˙ = v
(5)
The output flow at mouth is estimated by spirometry measure of forced expiratory volume in 1 second (FEV1 ). In order to obtain velocity at mouth (v1 ), it can be used:
v1 =
V˙ s
(6)
where S is the cross section surface constant for previous hypothesis. The spirometry test does not permit to obtain data linked to the airflow between the junction region. There is an empirical relation [11] which links the tidal volume, which is estimated via spirometry, and the patient’s anatomical death space, in order to extrapolate medium air flow rate in respira˙ tory zone (Va): ˙ = f · (Vt − Vd ) Va
(7)
f is the frequency of the respiratory act, Vt is the tidal volume, Vd is the death space volume. Frequency used in simulation is 1 Hz, because the aim of this model is to study a single expiratory act. Once flux is found, (6) can be used, in order to obtain v2 . P2 is the last unknown parameter for computing P1 from (1) and can be estimated from respiratory zone model.
Fig. 2 – Second part of the model, that mimics respiratory zone.
2.2.
Materials and methods: respiratory zone
The system reported in Fig. 2 is expressed through balance of pressures at equilibrium. Respiratory zone is structured as one bubble. Inside this structure there are an alveolar pressure which tries to expand the sac and an elastic pressure which is linked to the viscoelastic properties of the tissue. The last one tries to minimize the variation of the alveolar volume. Out of the bubble there is a pleural pressure, which describes the action of the outside environment. It pushes the wall and it limits the expansion of the bubble. During the expiratory phase pleural pressure acts with the elastic pressure in order to throw out stored air volume in inspiratory-phase. Alveolar pressure counteracts this movement, in order to avoid the collapse of the system [3].
Palveolar − Ppleuric = Pelastic Pelastic = PEEP + V/C
(8)
Note that elastic pressure term is composed by two parts: one is time invariant and it is calculated at t = 0 and the other one is estimated by a ratio between volume V and compliance C. At t = 0, we consider a steady state after a previous expiratory phase. The elastic pressure estimated at t = 0 is also called PEEP. In non-pathologic subjects this pressure is assumed to be zero [12]. The relationship between volume and pressure is expressed by a compliance term, representing tissue visco-elasticity. In the literature, compliance is described as a time-variable or a constant parameter, depending on the physiological aspects that have to be described. We use a static (constant) compliance, because model describes only one expiratory act and we are not interested in dynamic changes in a long period. In the literature the static compliance is equal to 0.002 dm3 /Pa [13]. The term V represents difference between tidal volume (Vt ) and death space volume (Vd ).
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 4 ( 2 0 1 1 ) 300–305
303
Table 1 – Input data for model validation. N
Height [cm]
FEV1 [l/s]
Vt [ml]
Vd [ml]
1 2 3 4 5 6 7 8 9 10
162 170 173 173 177 160 168 170 158 180
3.55 3.92 4.66 4.14 4.83 3.18 4.54 3.86 3.15 3.94
400 440 470 370 570 338 607 690 473 643
60 54 90 75 155 88 106 35 89 250
If these two volumes are known, system (8) can be solved and alveolar pressure can be derived. This is the output of second part of the model. The value assigned to pleural parameter is obtained from the literature and it is typically equal to 500 Pa. On the other hand it is impossible to measure pleural pressure in a patient, because a too invasive test would be required. The alveolar pressure is used as an input of the first part of the model to solve the Bernoulli’s equation.
3.
Overview of the model
The overview of the model is shown in Fig. 3, where inputs and outputs are highlighted. P = R · V˙
(9)
The term R is the resistance to be estimated. The airflow at mouth is extracted by spyrometric test and it is called FEV1. In conclusion (10) is the resistance of bronchial tubes: R=
P ˙ 2 FEV1
(10)
where P is the difference between alveolar pressure and pressure relative to mouth.
4.
Data analysis
The model has been validated comparing the outputs of the model with gold standards that can be found in the literature. In particular, the estimated value for the mouth pressure P1 is compared with a range of values acceptable for a physiologic respiratory system (typically they are varying from 2.5 Pa to 5 Pa [13]). The “Clinical Physiology Institute” of Pisa provided the inputs of the model: FEV1, tidal volume and death space volume for 10 different subjects. Spirometry test previously ascertained that all subjects were healthy. Table 1 summarizes the input data for model validation. The model estimates alveolar pressure, mouth pressure and resistance of a bronchial resistance to air flow in exhalation phase as it is summarized in Table 2. It can be noticed that P1 values obtained from simulation phase fall in acceptable range typical of literature for mouth pressure [13]. Moreover it can be argued that accuracy of model outputs is 90%. This percentage means that 9 subjects on 10 show P1 values that are enclosed in the acceptable
Fig. 3 – Summary of outputs from the respiratory zone, used in conductive zone for obtaining P1 and the input of the model. When pressures P1 and P2 are obtained, the model is able to estimate the resistance of bronchial tubes to air flow. The model studies the system in turbulence state, and (9) summarizes the relationship between pressures and flow.
range and only one is very closed to it (5.1 Pa vs. the limit of 5 Pa). In the literature the gold standard values for resistances of different levels of bronchial tree to airflow are well-known [6]. Thus we compared resistance value that was estimated for each subject of Table 2 with a curve that represents bronchial resistance in relation with system level (Fig. 4).
304
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 4 ( 2 0 1 1 ) 300–305
Table 2 – Outputs estimated from the model corresponding to the input data of Table 1. N
P2 [Pa]
P1 [Pa]
R [cm H2 O/(l × s)2 ]
1 2 3 4 5 6 7 8 9 10
670 693 690 647.5 707.5 625 750.5 827.5 692 718
4.53 4.13 4.21 4.87 3.83 5.1 2.75 0.08 4.14 3.59
0.528 0.448 0.316 0.375 0.302 0.613 0.363 0.555 0.693 0.458
Fig. 4 – Average value of resistance measured at each level of bronchial tree. The red point represents the value estimated from model and his level in the system. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)
From Fig. 4 we obtain an interesting result: the model studies the resistance of duct with medium cross section surface. This is very important, as it will be discussed in the next section.
5.
Discussion and conclusion
The model is implemented in a M-file in Matlab 7.0 environment. The equations describing the model are time-invariant, leading to algebraic simple computations. The input data are collected in an Excel file and passed to Matlab environment. The estimation of output is therefore very fast and it is organized in a table, which allows a simple and fast information for the users. The model is innovative with respect to the previous literature, in terms of estimated outputs. It was validated with 10 different healthy subject’s data provided from “Clinical Physiology Institute” of Pisa. It can be noticed that P1 values, which we obtained from simulation phase, fall in acceptable range that the literature shows for mouth pressure (2.5–5.0 Pa) [13]. Moreover it can be argued that model accuracy is 90%, because 9 subjects on 10 show values that are enclosed in the range and the last one is very closed to it (5.1 Pa). Moreover, resistances estimated by the model are correlated with medium caliber ducts. This is a very important fact to go into depth. Currently the common diagnostic tests through spirometry and body plethysmography to measure the airway resistance (RAW), are unable to give data linked to a particular level of bronchial tree. The outputs of these tests are concerned with whole respiratory system; in particular they are influenced in height percentage by big caliber ducts. The airway resistances “RAW” are influenced in 80% by ducts of large cross section surface [14], as shown in Fig. 5. The red bar indicates tree level sections that influence in higher percentage the output of test. As it can be argued from Fig. 5, the proposed model provides an output which is influenced in high percentage by resistance of medium cross section ducts.
Fig. 5 – “A” represents how the output of modern diagnostic methods is influenced by difference levels of tree (red part: high influence, blue one: low influence). “B” shows that the output of the model proposed is influenced only by the ducts of medium cross section surface (red part). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 4 ( 2 0 1 1 ) 300–305
305
Chronic obstructive pulmonary disease (COPD) is a pathology that causes airflow limitation as a consequence of airways obstruction. Nowadays functional diagnosis is based on spirometry tests and RAW tests [15]. These tests are gold
mated if post therapy airways resistances show a lower value than before therapy. This implies that drug therapy is effective for the subject.
standard as they are non invasive, standardized, reproducible and objective. Raw tests estimates correlation between air flow, measured at mouth and output pressure. This method permits evaluation of viscous mechanic properties, in normal condition of expiration. A limitation of these tests is correlated with the functional standpoint of COPD, which is characterized by a progressive airflow obstruction [15]. It is well known that airways obstruction is not homogeneous, but it is maximum in medium-small caliber ducts [15]. For this reason, RAW test is not the best method to make a specific diagnosis, because the ducts influence the output of this test [15], as shown in Fig. 5. Moreover, a pathology like the COPD interests in early-stage the duct of little cross section surface, then it spreads to other bronchial tree sections. When this type of disease is detected by diagnostic methods currently in use, outputs of tests are influenced in high percentage by big section ducts and only when pathology is widespread to all ducts of respiratory system it is detected from diagnostic tests [15]. The trend in research is to develop methods that allow the diagnosis of COPD, before all bronchial tree is affected, in order to anticipate treatment and improve its efficacy. Following the trend, the proposed model allows the study of a particular section of bronchial tree. This model will be able to detect the pathology when it interests the medium cross section surface duct, so diagnosis can be more precise and specific. This expected result requires a validation of model with pathological subject’s data to be statistically consistent and it represents a future step of the research. Preliminary results obtained in cooperation with the Clinical Physiology Institute of CNR and G. Monasterio Foundation in Pisa are strongly encouraging and confirm that proposed model is effective in the detection of the pathology at an early-stage in patients. To our knowledge, other models that can be found in the literature are unable to study pathologic conditions, because of their complex structure and the need of data collected from too invasive tests. The set of equations used is time independent. This fact, although hides dynamics, allows the clinician to use the model, following steady-state slow variations, for studying in a simple way the follow up of the patient. Future work will be devoted to monitor the patient making new tests after drug therapy and using the results as inputs of the model. The post-therapy outputs will be compared with the pre-therapy outputs, e.g., taking into account the airways resistances, we could estimate the ratio between the airways resistances calculated before therapy and those evaluated after drug administration. A general condition improvement will be esti-
Conflict of interest The authors declare that they have no conflict of interest.
references
[1] Claudio, Cobelli e Ewart Carson “Introduction to Modelling in Physiology and Medicine”, Academic Press, 2008. [2] M.C.K. Khoo, Physiological Control Systems: Analysis, Simulation and Estimation, Prentice Hall of India Pvt. Ltd., New Delhi, 2001. [3] J.G. Chase, T. Yuta, K.J. Mulligan, G.M. Shaw, B. Horn, A novel mechanical lung model of pulmonary diseases to assist with teaching and training, Biomed Central, 20 August 2006. [4] K. Sohn, W.J. Warwick, Y.W. Lee, J. Lee, J.E. Holte, Investigation of non-uniform air flow signal oscillation during high frequency chest compression, Biomedical Engineering, on line 2005. http://www.biomedicalengeneering-online.com/content/4/1/31. [5] S. Mesic, R. Babuska, H.C. Hoogsteden, A.F.M. Verbraak, Computer-controlled mechanical simulation of the artificially ventilated human respiratory system, IEEE Transactions on Biomedical Engineering 50 (June) (2003). [6] J.B. West, Respiratory Physiology: The Essentials, Lippincott-Raven, 1979. [7] H.S. Martin, Physiological Basis of Respiratory Disease, BC Decker Inc Hamilton, 2005. [8] J. Chretien, “Pneumologia” annual report, Masson Italia Editori, 1978. [9] S. Middleman, An Introduction to Fluid Mechanics: Principles of Analysis and Design, Wiley, 1998. [10] Internet website: http://www.caleffi.it http://www.caleffi.it/Caleffi:hydronic solutions. [11] Internet website: http://www.cfww.org/ipg-cf/article/195/, “Physiotherapy in the treatment of cystic fibrosis: from infant to adult” International Physiotherapy Group for Cystic Fibrosis, 2009. [12] K.A. Cope, M.T. Watson, W.M. Foster, S.S. Sehnert, T.H. Risby, Effect of ventilation on the collection of exhaled breath collection, American Physiological Society (2003). [13] W.S. Fowler, Lung function studies. II. The respiratory dead space, J. Appl. Physiol. 96 (2004) 1371–1379. [14] Internet website: http://en.wikipedia.org/wiki/Body plethysmography. [15] T.J. Wilt, D. Niewoehner, C. Kim, et al., Use of spirometry for case finding, diagnosis, and managment of chronic obstructive polmunary disease (COPD), Evidence Report/Technology Assessment No. 121 AHRQ Publication No. 05-E017-2. Rockville, MD. Agency for Healthcare Research and Quality, September, 2005.
journal homepage: www.intl.elsevierhealth.com/journals/cmpb
Innovative model to simulate exhalation phase in human respiratory system Tommaso Sbrana a,∗ , Alberto Landi b , Giosuè Angelo Catapano c a b c
Centro Interdipartimentale di Ricerca E. Piaggio, Faculty of Engineering, University of Pisa, Italy Department of Energy and Systems Engineering, University of Pisa, Italy Tuscany Foundation ‘Giuseppe Monasterio’ CNR, Pisa, Italy
a r t i c l e
i n f o
a b s t r a c t
Article history:
In this paper, we present a mathematical model, which mimics the bronchial resistances
Received 13 November 2010
of human’s lung in an expiratory act. The model is implemented in Matlab. The inputs that
Received in revised form
are used in this model derive from spirometry test. This model is able to study a physiologic
22 January 2011
condition, a pathologic one and the patient’s follow up after drug treatment. We split our
Accepted 4 February 2011
study into two parts. The first one focuses the analysis on the gas fluido dynamic inside of the respiratory pathways. The second part takes care of the pressure equilibrium in the
Keywords:
exchange zone. We use the outputs that derive from the second subsystem to solve the
Model
Bernoulli’s equation of the first part. The model was validated with data provided from
Lung
“Clinical Physiology Institute” of CNR and G. Monasterio Foundation of Pisa.
Respiratory system
© 2011 Elsevier Ireland Ltd. All rights reserved.
COPD Exhalation phase
1.
Introduction
Mathematical models in medicine are often used to evaluate the information provided by diagnostic tests. There are several advantages in this approach, like the possibility to test a hypothesis, to mimic a pathologic condition and/or a drug action avoiding side effects in patients. From the point of view of a model builder two different states have to be usually considered: basal or pathological conditions. The key problem is: how to move from a model emulating average conditions of the physiologic system in a basal state to a model able to represent the system in a pathological-state. Solving this problem is relevant from a diagnostic and prognostic point of view,
because it allows the physician to prove his diagnosis using the model, reducing delays in obtaining a correct diagnosis and avoiding invasive additional tests. In recent years many mathematical models have been proposed in the literature for describing human physiological systems. One of the main problems of many models is their inadequacy for helping medical diagnostic techniques and new therapeutic protocols: sometimes they need too many unknown parameters derived from an accurate description of the physiological system, so that they cannot be easily identified. On the other hand simplified models can be used for explaining and teaching how a biological system works [1], but they are usually too approximate and simple to help physicians. Many attempts have been reported in the literature to
Abbreviations: Pi , pressure “i”; Zi , height “i” from earth; hf , parameter relative to weight losses; L, length of the tube; D, diameter of cross section surface of the tube; v, medium rate of velocity of the air flow along the tube; f, frequency of the respiratory act; Vt , tidal volume; Vd , death space volume. ∗ Corresponding author. Tel.: +39 334 901 760. E-mail address: [email protected] (T. Sbrana). 0169-2607/$ – see front matter © 2011 Elsevier Ireland Ltd. All rights reserved. doi:10.1016/j.cmpb.2011.02.005
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 4 ( 2 0 1 1 ) 300–305
model and simulate human respiratory physiology. Although it is rather simplistic to classify these models, because of the complexity of respiratory physiology, different approaches can be distinguished, e.g., electrical models, mechanical models, morphometric models and mathematical models. Electrical models [2] describe the system using a Maxwell or a Firestone analogy in order to transform the physiological system into equivalent electrical circuits. The problem involved with this approach is that the transformed system is often too simplified and that a circuital study does not match clinical outputs that well correlate with medical findings. Mechanical models transform the physiological system into an equivalent mechanical system, based on pistons, gears and springs [3]. This approach has led to models able to obtain outputs well-correlated with doctor expectancy, but unable to mimic pathological conditions for their structural rigidity that makes unfeasible the emulation of cases different from the basic one. Morphometric models need inputs usually unknown and difficult to obtain without invasive tests [4]. Mathematical models can describe physiological system by using a datadriven (or black-box) approach [1], useful for fitting curves to input–output data when our understanding is limited and inadequate to produce a detailed model, or a system-based (or white-box) approach, formulated on the basis of comprehensive models able to adequately describe the underlying physiology. In the complex case of respiratory physiology a combination of two approaches (or gray-box), e.g. [5], is the optimal solution for joining the explicit structure of parts of model based on physiology and available data. In this paper we follow a reference model, i.e. a white-box approach, oriented to correlate spirometric data with the patient condition. It has to be highlighted that white-box approach needs some expedient to extrapolate all required parameters, as it will be explained in Section 2. Model presents study of bronchial resistances of human lung in an expiratory act and is composed by static and linear equations, which are implemented in Matlab M files, with the aim to obtain a simple method to correlate spirometric data with the patient condition.
2.
Materials and methods
The human respiratory system is divided into two parts: conductive and respiratory zone [6]. This division can be done, considering the exchange of oxygen between bronchial tubes and other tissues. It is interesting to notice that gas passes the endothelial membrane only in respiratory zone, while in conductive zone there is no gas exchange. In this hypothesis the model is split into two parts: one mimics the conductive part, and the other is deputed to the respiratory zone. The two parts are linked each other in a junction zone called “A” in Fig. 1. Some comments regarding this approach must be considered. Since the aim is to build a simple model able to correlate spirometric data with the patient condition, the choice to split into two parts the respiratory system was due to the attempt of maintaining essential functions of the system, reducing as much as possible the number of required anatomical data. Furthermore, this approach avoids a detailed – difficult to extract – description of the lung structure and of bronchial tubes. In our opinion a more detailed description of the system may
301
Fig. 1 – Design of the system. The two parts are depicted in different colors (blue for the conductive zone, violet for the respiratory one) and are divided by junction zone A. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)
lead to a model complexity that is excessively high, compared to the available data.
2.1.
Materials and methods: conductive zone
The first part of the model is considered as a rigid smoothwalled tube with a circular cross section surface, where the air can flow. The assumption of rigid wall derives from the absence of a pulsating movement of the wall, like that one of the arteries. The assumption of smooth surface derives from the presence of mucus along the bronchial tubes [6,7]. In physiologic systems the bronchial diameters become smaller and smaller from the mouth to the deepest regions of the bronchial tree and at the same time the total surface grows up. These phenomena are considered in the literature with an overestimation of the tube diameter [8] (an average value can be 2.2 cm2 ) and the total length of the conductive zone is estimated in the literature [8] as L = 30 cm. The circuit is studied by a classical fluid-dynamic approach [9]. Several studies have shown that the air flux in the conductive zone is turbulent [6] and the physiological system can be modeled using a fluid-dynamic analysis in a turbulence state. Therefore the Bernoulli’s expression is introduced: v2 v2 P2 P1 + g · Z1 + 1 = + g · Z2 + 2 + hf 2 2
(1)
Pi is the pressure “i”, Zi is the height “i” from earth, hf is the parameter relative to weight losses, vi is the velocity “i”, g is the gravity acceleration. The subscript 1 indicates values data correlated to the mouth, the subscript 2 considers values linked to the respiratory zone. The desired output of first sub-model is P1 . The height values are estimated for each case: the mouth is positioned 15 cm under the scalp. Z2 is correlated with the length of the tube, so that it can be assumed as a known data. The load losses (hf ) are composed by two elements: one due to friction (hf1 ) and the second to branches along the bronchial tree (hf2 ).
302
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 4 ( 2 0 1 1 ) 300–305
In (2) we express the relation between the two parameters hf1 and hf2 : hf = hf 1 + hf 2
(2)
From the Moody’s diagram it is possible to obtain Funning’s factor (f), for estimating the term hf1 [9]: hf 1 = 2 · f ·
L · v2 D
(3)
L is the length of the tube, D is the diameter of cross section surface of the tube, v is the medium rate of velocity of the air flow along the tube. The load losses due to the branches in a circuit in turbulence flow conditions [9] is:
hf 2 = ε ·
v2 2
(4)
where ε is the proportional coefficient typical of turbulence flows. The ε parameter can be estimated if the angle between two pathways is assumed of 60◦ . In physiology, this value represents an angles average between two pathways, as it is well described in [7]. This value implies ε equal to 0.3 [10]. Therefore, in order to obtain P1 , unknown parameters of the Bernoulli’s expression are the velocities v1 , v2 , and the pressure called P2 . The last one is linked to the junction region between the respiratory zone and the conductive one. The air flow is defined by: V˙ = v
(5)
The output flow at mouth is estimated by spirometry measure of forced expiratory volume in 1 second (FEV1 ). In order to obtain velocity at mouth (v1 ), it can be used:
v1 =
V˙ s
(6)
where S is the cross section surface constant for previous hypothesis. The spirometry test does not permit to obtain data linked to the airflow between the junction region. There is an empirical relation [11] which links the tidal volume, which is estimated via spirometry, and the patient’s anatomical death space, in order to extrapolate medium air flow rate in respira˙ tory zone (Va): ˙ = f · (Vt − Vd ) Va
(7)
f is the frequency of the respiratory act, Vt is the tidal volume, Vd is the death space volume. Frequency used in simulation is 1 Hz, because the aim of this model is to study a single expiratory act. Once flux is found, (6) can be used, in order to obtain v2 . P2 is the last unknown parameter for computing P1 from (1) and can be estimated from respiratory zone model.
Fig. 2 – Second part of the model, that mimics respiratory zone.
2.2.
Materials and methods: respiratory zone
The system reported in Fig. 2 is expressed through balance of pressures at equilibrium. Respiratory zone is structured as one bubble. Inside this structure there are an alveolar pressure which tries to expand the sac and an elastic pressure which is linked to the viscoelastic properties of the tissue. The last one tries to minimize the variation of the alveolar volume. Out of the bubble there is a pleural pressure, which describes the action of the outside environment. It pushes the wall and it limits the expansion of the bubble. During the expiratory phase pleural pressure acts with the elastic pressure in order to throw out stored air volume in inspiratory-phase. Alveolar pressure counteracts this movement, in order to avoid the collapse of the system [3].
Palveolar − Ppleuric = Pelastic Pelastic = PEEP + V/C
(8)
Note that elastic pressure term is composed by two parts: one is time invariant and it is calculated at t = 0 and the other one is estimated by a ratio between volume V and compliance C. At t = 0, we consider a steady state after a previous expiratory phase. The elastic pressure estimated at t = 0 is also called PEEP. In non-pathologic subjects this pressure is assumed to be zero [12]. The relationship between volume and pressure is expressed by a compliance term, representing tissue visco-elasticity. In the literature, compliance is described as a time-variable or a constant parameter, depending on the physiological aspects that have to be described. We use a static (constant) compliance, because model describes only one expiratory act and we are not interested in dynamic changes in a long period. In the literature the static compliance is equal to 0.002 dm3 /Pa [13]. The term V represents difference between tidal volume (Vt ) and death space volume (Vd ).
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 4 ( 2 0 1 1 ) 300–305
303
Table 1 – Input data for model validation. N
Height [cm]
FEV1 [l/s]
Vt [ml]
Vd [ml]
1 2 3 4 5 6 7 8 9 10
162 170 173 173 177 160 168 170 158 180
3.55 3.92 4.66 4.14 4.83 3.18 4.54 3.86 3.15 3.94
400 440 470 370 570 338 607 690 473 643
60 54 90 75 155 88 106 35 89 250
If these two volumes are known, system (8) can be solved and alveolar pressure can be derived. This is the output of second part of the model. The value assigned to pleural parameter is obtained from the literature and it is typically equal to 500 Pa. On the other hand it is impossible to measure pleural pressure in a patient, because a too invasive test would be required. The alveolar pressure is used as an input of the first part of the model to solve the Bernoulli’s equation.
3.
Overview of the model
The overview of the model is shown in Fig. 3, where inputs and outputs are highlighted. P = R · V˙
(9)
The term R is the resistance to be estimated. The airflow at mouth is extracted by spyrometric test and it is called FEV1. In conclusion (10) is the resistance of bronchial tubes: R=
P ˙ 2 FEV1
(10)
where P is the difference between alveolar pressure and pressure relative to mouth.
4.
Data analysis
The model has been validated comparing the outputs of the model with gold standards that can be found in the literature. In particular, the estimated value for the mouth pressure P1 is compared with a range of values acceptable for a physiologic respiratory system (typically they are varying from 2.5 Pa to 5 Pa [13]). The “Clinical Physiology Institute” of Pisa provided the inputs of the model: FEV1, tidal volume and death space volume for 10 different subjects. Spirometry test previously ascertained that all subjects were healthy. Table 1 summarizes the input data for model validation. The model estimates alveolar pressure, mouth pressure and resistance of a bronchial resistance to air flow in exhalation phase as it is summarized in Table 2. It can be noticed that P1 values obtained from simulation phase fall in acceptable range typical of literature for mouth pressure [13]. Moreover it can be argued that accuracy of model outputs is 90%. This percentage means that 9 subjects on 10 show P1 values that are enclosed in the acceptable
Fig. 3 – Summary of outputs from the respiratory zone, used in conductive zone for obtaining P1 and the input of the model. When pressures P1 and P2 are obtained, the model is able to estimate the resistance of bronchial tubes to air flow. The model studies the system in turbulence state, and (9) summarizes the relationship between pressures and flow.
range and only one is very closed to it (5.1 Pa vs. the limit of 5 Pa). In the literature the gold standard values for resistances of different levels of bronchial tree to airflow are well-known [6]. Thus we compared resistance value that was estimated for each subject of Table 2 with a curve that represents bronchial resistance in relation with system level (Fig. 4).
304
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 4 ( 2 0 1 1 ) 300–305
Table 2 – Outputs estimated from the model corresponding to the input data of Table 1. N
P2 [Pa]
P1 [Pa]
R [cm H2 O/(l × s)2 ]
1 2 3 4 5 6 7 8 9 10
670 693 690 647.5 707.5 625 750.5 827.5 692 718
4.53 4.13 4.21 4.87 3.83 5.1 2.75 0.08 4.14 3.59
0.528 0.448 0.316 0.375 0.302 0.613 0.363 0.555 0.693 0.458
Fig. 4 – Average value of resistance measured at each level of bronchial tree. The red point represents the value estimated from model and his level in the system. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)
From Fig. 4 we obtain an interesting result: the model studies the resistance of duct with medium cross section surface. This is very important, as it will be discussed in the next section.
5.
Discussion and conclusion
The model is implemented in a M-file in Matlab 7.0 environment. The equations describing the model are time-invariant, leading to algebraic simple computations. The input data are collected in an Excel file and passed to Matlab environment. The estimation of output is therefore very fast and it is organized in a table, which allows a simple and fast information for the users. The model is innovative with respect to the previous literature, in terms of estimated outputs. It was validated with 10 different healthy subject’s data provided from “Clinical Physiology Institute” of Pisa. It can be noticed that P1 values, which we obtained from simulation phase, fall in acceptable range that the literature shows for mouth pressure (2.5–5.0 Pa) [13]. Moreover it can be argued that model accuracy is 90%, because 9 subjects on 10 show values that are enclosed in the range and the last one is very closed to it (5.1 Pa). Moreover, resistances estimated by the model are correlated with medium caliber ducts. This is a very important fact to go into depth. Currently the common diagnostic tests through spirometry and body plethysmography to measure the airway resistance (RAW), are unable to give data linked to a particular level of bronchial tree. The outputs of these tests are concerned with whole respiratory system; in particular they are influenced in height percentage by big caliber ducts. The airway resistances “RAW” are influenced in 80% by ducts of large cross section surface [14], as shown in Fig. 5. The red bar indicates tree level sections that influence in higher percentage the output of test. As it can be argued from Fig. 5, the proposed model provides an output which is influenced in high percentage by resistance of medium cross section ducts.
Fig. 5 – “A” represents how the output of modern diagnostic methods is influenced by difference levels of tree (red part: high influence, blue one: low influence). “B” shows that the output of the model proposed is influenced only by the ducts of medium cross section surface (red part). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 4 ( 2 0 1 1 ) 300–305
305
Chronic obstructive pulmonary disease (COPD) is a pathology that causes airflow limitation as a consequence of airways obstruction. Nowadays functional diagnosis is based on spirometry tests and RAW tests [15]. These tests are gold
mated if post therapy airways resistances show a lower value than before therapy. This implies that drug therapy is effective for the subject.
standard as they are non invasive, standardized, reproducible and objective. Raw tests estimates correlation between air flow, measured at mouth and output pressure. This method permits evaluation of viscous mechanic properties, in normal condition of expiration. A limitation of these tests is correlated with the functional standpoint of COPD, which is characterized by a progressive airflow obstruction [15]. It is well known that airways obstruction is not homogeneous, but it is maximum in medium-small caliber ducts [15]. For this reason, RAW test is not the best method to make a specific diagnosis, because the ducts influence the output of this test [15], as shown in Fig. 5. Moreover, a pathology like the COPD interests in early-stage the duct of little cross section surface, then it spreads to other bronchial tree sections. When this type of disease is detected by diagnostic methods currently in use, outputs of tests are influenced in high percentage by big section ducts and only when pathology is widespread to all ducts of respiratory system it is detected from diagnostic tests [15]. The trend in research is to develop methods that allow the diagnosis of COPD, before all bronchial tree is affected, in order to anticipate treatment and improve its efficacy. Following the trend, the proposed model allows the study of a particular section of bronchial tree. This model will be able to detect the pathology when it interests the medium cross section surface duct, so diagnosis can be more precise and specific. This expected result requires a validation of model with pathological subject’s data to be statistically consistent and it represents a future step of the research. Preliminary results obtained in cooperation with the Clinical Physiology Institute of CNR and G. Monasterio Foundation in Pisa are strongly encouraging and confirm that proposed model is effective in the detection of the pathology at an early-stage in patients. To our knowledge, other models that can be found in the literature are unable to study pathologic conditions, because of their complex structure and the need of data collected from too invasive tests. The set of equations used is time independent. This fact, although hides dynamics, allows the clinician to use the model, following steady-state slow variations, for studying in a simple way the follow up of the patient. Future work will be devoted to monitor the patient making new tests after drug therapy and using the results as inputs of the model. The post-therapy outputs will be compared with the pre-therapy outputs, e.g., taking into account the airways resistances, we could estimate the ratio between the airways resistances calculated before therapy and those evaluated after drug administration. A general condition improvement will be esti-
Conflict of interest The authors declare that they have no conflict of interest.
references
[1] Claudio, Cobelli e Ewart Carson “Introduction to Modelling in Physiology and Medicine”, Academic Press, 2008. [2] M.C.K. Khoo, Physiological Control Systems: Analysis, Simulation and Estimation, Prentice Hall of India Pvt. Ltd., New Delhi, 2001. [3] J.G. Chase, T. Yuta, K.J. Mulligan, G.M. Shaw, B. Horn, A novel mechanical lung model of pulmonary diseases to assist with teaching and training, Biomed Central, 20 August 2006. [4] K. Sohn, W.J. Warwick, Y.W. Lee, J. Lee, J.E. Holte, Investigation of non-uniform air flow signal oscillation during high frequency chest compression, Biomedical Engineering, on line 2005. http://www.biomedicalengeneering-online.com/content/4/1/31. [5] S. Mesic, R. Babuska, H.C. Hoogsteden, A.F.M. Verbraak, Computer-controlled mechanical simulation of the artificially ventilated human respiratory system, IEEE Transactions on Biomedical Engineering 50 (June) (2003). [6] J.B. West, Respiratory Physiology: The Essentials, Lippincott-Raven, 1979. [7] H.S. Martin, Physiological Basis of Respiratory Disease, BC Decker Inc Hamilton, 2005. [8] J. Chretien, “Pneumologia” annual report, Masson Italia Editori, 1978. [9] S. Middleman, An Introduction to Fluid Mechanics: Principles of Analysis and Design, Wiley, 1998. [10] Internet website: http://www.caleffi.it http://www.caleffi.it/Caleffi:hydronic solutions. [11] Internet website: http://www.cfww.org/ipg-cf/article/195/, “Physiotherapy in the treatment of cystic fibrosis: from infant to adult” International Physiotherapy Group for Cystic Fibrosis, 2009. [12] K.A. Cope, M.T. Watson, W.M. Foster, S.S. Sehnert, T.H. Risby, Effect of ventilation on the collection of exhaled breath collection, American Physiological Society (2003). [13] W.S. Fowler, Lung function studies. II. The respiratory dead space, J. Appl. Physiol. 96 (2004) 1371–1379. [14] Internet website: http://en.wikipedia.org/wiki/Body plethysmography. [15] T.J. Wilt, D. Niewoehner, C. Kim, et al., Use of spirometry for case finding, diagnosis, and managment of chronic obstructive polmunary disease (COPD), Evidence Report/Technology Assessment No. 121 AHRQ Publication No. 05-E017-2. Rockville, MD. Agency for Healthcare Research and Quality, September, 2005.