Airflow in the respiratory system is complicated as it goes through various regions with different geometries and mechanical properties. Three-dimensional (3-D) simulations are typically limited to local areas of the system because of their high computational cost. On the other hand, the one-dimensional (1-D) equations of flow in compliant tubes offer a good compromise between accuracy and computational cost when a global assessment of airflow in the system is required. The aim of the current study is to apply the 1-D formulation in space and time variables to study the propagation of a pulse wave in human airways; first in a simple system composed of just one bifurcation, trachea-main bronchi, according to the symmetrical Weibel model. Then extending the system to include a further generation, the bronchi branches. Pulse waveforms carry information about the functionality and morphology of the respiratory system and the 1-D modelling, in terms of space and time variables, represents an innovative approach for respiratory response interpretation. 1-D modelling in space-time variables has been extensively applied to simulate blood pressure and flow in the cardiovascular system. This work represents the first attempt to apply this formulation to study pulse waveforms in the human bronchial tree

Alastruey, J

Clavica, F

Khir, AW

Sherwin, SJ

English

Brunel University Research Archive

    Abstract— Airflow in the respiratory system is complicated as it goes through various regions with different geometries and mechanical properties. Three-dimensional (3-D) simulations are typically limited to local areas of the system because of their high computational cost. On the other hand, the one-dimensional (1-D) equations of flow in compliant tubes offer a good compromise between accuracy and computational cost when a global assessment of airflow in the system is required. The aim of the current study is to apply the 1-D formulation in space and time variables to study the propagation of a pulse wave in human airways; first in a simple system composed of just one bifurcation, trachea-main bronchi, according to the symmetrical Weibel model. Then extending the system to include a further generation, the bronchi branches. Pulse waveforms carry information about the functionality and morphology of the respiratory system and the 1-D modelling, in terms of space and time variables, represents an innovative approach for respiratory response interpretation. 1-D modelling in space-time variables has been extensively applied to simulate blood pressure and flow in the cardiovascular system. This work represents the first attempt to apply this formulation to study pulse waveforms in the human bronchial tree.   I. INTRODUCTION  nformation on the morphology and functionality of the cardiovascular and respiratory systems can  be derived by analyzing pulse waveforms in arteries [1, 2] and airways, respectively. Understanding the underlying mechanisms of pulse wave propagation in normal conditions and the impact of disease and anatomical variation on the patterns of propagation is therefore relevant to improve prevention, diagnosis and treatment of disease. The one-dimensional (1-D) equations of flow in compliant vessels offer a good  Manuscript received  April 23, 2009 F. Clavica is with Brunel Institute for Bioengineering, Brunel University,  Uxbridge UB8 3PH, UK  (Francesco.Clavica@brunel.ac.uk) J. Alastruey is with the Department of Bioengineering, Imperial College London, South Kensington Campus, SW7 2AZ London, UK  (jordi.alastruey-arimon@imperial.ac.uk)   S. J. Sherwin is with the Department of Aeronautics, Imperial College London, South Kensington Campus, SW7 2AZ London, UK . (s.sherwin@imperial.ac.uk)  A. W. Khir is with Brunel Institute for Bioengineering, Brunel University,  Uxbridge UB8 3PH, UK  (ashraf.khir@brunel.ac.uk)  compromise between accuracy and computational cost when a global assessment of the system is required; 3-D simulations are typically limited to local areas of the system because of their high computational cost [2]. The 1-D model in space-time variables provides information regarding the distance and the magnitude of potential blockages in the airways respectively derived from wave timing and amplitude. The aim of this paper is to apply the 1-D formulation to study air wave propagation in the first and second generation of the human airways using different boundary conditions at the terminal segments.   II. METHODOLOGY A.1 Governing equations Conservation of mass and momentum applied to a 1-D impermeable and elastic tubular control volume of Newtonian incompressible fluid leads to the system of hyperbolic partial differential equations [1]:                      0)( =∂∂+∂∂xAUtA ,                                             (1)                     AfxpxUUtUρρ=∂∂+∂∂+∂∂ 1 ,                                 (2) where x is the axial coordinate of the tube, t is the time, A(x,t) is the cross-sectional area of the tube, U(x,t) is the average axial velocity of the fluid, p(x,t) is the average internal pressure over the cross-section, ρ is the density of the air ( 3/204.1 mkg ) , and f(x,t) is the friction force per unit length. Equations (1) and (2) can be completed with the pressure area relationship [1, 2]:                       )( eext AApp −+= β ,                                 (3) where:                eeAEhx )1()( 2νpiβ−=                                       (4) is a parameter associated with the mechanical properties of the tube wall. It depends on the wall thickness he(x) and the sectional area Ae(x) at the equilibrium state (p,U)=(pext,0). E is the Young’s modulus, pext is the constant external pressure, and ν=0.5 is the Poisson’s ratio (the vessel wall is considered to be incompressible). Equations (1) to (3) are solved using a discontinuous Galerkin scheme with a spectral/hp spatial discretisation, a second order Adams-Bashforth time-integration scheme, and the initial conditions One-Dimensional Modelling of Pulse Wave Propagation in Human Airway Bifurcations in Space–Time Variables  Francesco Clavica, Student Member, IEEE EMBS, Jordi Alastruey, Spencer J. Sherwin,              Ashraf W. Khir I 548231st Annual International Conference of the IEEE EMBSMinneapolis, Minnesota, USA, September 2-6, 2009978-1-4244-3296-7/09/$25.00 ©2009 IEEEAuthorized licensed use limited to: Brunel University. Downloaded on January 28, 2010 at 10:09 from IEEE Xplore.  Restrictions apply.   (p,U)=(0,0) everywhere in the system. The parameter β is related to the speed of pulse wave propagation c(x,t) through:                                   2/122AAceρβ=                              (5) A detailed description of the numerical solution is given in earlier work [1].  A.2  Linear analysis of wave reflections   Considering an incident wave that encounters a reflection site a component of the wave will be reflected and another transmitted. The reflection coefficient (R) can be defined as: R = (δp − ∆p) /∆p , where δp-∆p is the change of pressure across the reflected wave and ∆p is the change of pressure in the incident wave [1].  Considering a bifurcation where the parent vessel is 0 and the daughter vessels are 1 and 2 (Fig. 1), R takes the form:                     )/()/()/()/()/()/(221100221100cAcAcAcAcAcAR++−−=                        (6) Since pressure is constant at the bifurcation, the transmission coefficient can be calculated as: T=1+R. The terminal reflection coefficient (Rt), in a system with a single resistance (Rs) coupled to the outflow of a 1-D model terminal segment, can be expressed:                                   00ZRZRRsst+−=                                      (7) Where Z0 is the characteristic impedance of the terminal segment (eAcZ /0 ρ= ); 11 <<− tR .  Rt =-1 corresponds to an open-end condition, Rt =0 non-reflective condition (Rs=Z) and Rt =1 corresponds to a total blockage of the terminal segment.  B. Geometrical Model  In the human respiratory system the elementary airways branch in an irregular dichotomy way. However, in this study we follow the model of Weibel [3], which assumes symmetry at each bifurcation, and thus allows for an easier interpretation of our preliminary results. The parameters we used to build the model are: lengths (l) taken from Weibel [3], volume corrected diameter (DFRC) (applied to Weibel model to obtain diameter at functional residual capacity FRC) according to [7], wall thickness (h) from Montaudon et al. [4], and Young’s modulus (E) defined according to the content of cartilage and soft tissue of each generation [5-7].   TABLE I Lengths (l), Diameter at FRC (DFRC), wall thickness (h), Young’s modulus (E) and wave speed (c) for the considered Weibel generations (Gen.)  Gen. l(mm) DFRC(mm) h(mm) E (MPa) c(m/s) 0 120 17.83 1.4 2.9 506.5 1 47.6 12.1 1.3 2.9 592.3 2 19 8.05 1.3 1.4 514.9 Figure 1 and 2 show the bronchial tree models under consideration; in Figure 1 the one generation model (model 1) has Rt =0 and Rt =1 as terminal reflection coefficients for the two terminal segments (main bronchi). Rt=0 Rt=1021 G0 Fig. 1. Model 1: model of the trachea and main bronchi (first generation) with Rt =1 in one of the two main bronchi. G describes a bifurcation with subscript that denotes the number of the mother tube.  Figure 2 shows the model for a two generations system respectively with one (model 2A) and two (model 2B) blockages (Rt=1) in the terminal segments. The remaining segments are set with Rt =0.  Rt=0 Rt=0 Rt=0 Rt=1023 4 5 61A)G0G1 G2Rt=1 Rt=0 Rt=0 Rt=1023 4 5 61B)G0G1 G2 Fig. 2. (A) Model 2A: model of the trachea, first and second generation; Rt=1 in one of the four terminal segments. (B) Model 2B:  model of the first and second generation, Rt =1 in two of the four terminal segments. G describes a bifurcation with subscript that denotes the number of the corresponding mother tube.  We enforced a flow pulse (Q) of 0.1 l/s at the inlet of the trachea. This delta wave was approximated with the Gaussian function:                      ))01.0(*10(1 27.9*10)( −−−= tetQ                                    (8) The time of 0.01s corresponds to the time of the initial Gaussian pulse reaching its peak. We impose that, after the pulse, the inlet behaves as a reflective boundary with Rt =1. Calculations were made at the middle point of each segment. In studying the waves running through the bifurcation we refer to the mother tube as 0 and the daughters as tubes 1 and 2 according to Figure 1. In Figure 2  the mother tube is 0, the two segments of the first generation are 1 and 2 while segments of the second generation are identified with 3-4-5-6. Point G describes a bifurcation with subscript that denotes the number of the corresponding mother tube. The waves travelling forward and backward will be denoted by ‘+’ and ‘-‘ respectively. Therefore, a wave that runs in the path 0-0 means that it runs forward from the inlet of the mother tube, is reflected at the bifurcation G0  then runs backward to the inlet of the mother tube. The path 01-1-0 means that the 5483Authorized licensed use limited to: Brunel University. Downloaded on January 28, 2010 at 10:09 from IEEE Xplore.  Restrictions apply.   wave runs from the inlet of the mother tube towards the bifurcation G0, is transmitted into the daughter tube 1, is reflected backwards to the bifurcation G0 and finally runs backwards to the inlet of the mother tube [8]. III. RESULTS  Figure 3, 4 and 5 show the simulated pressure waves  at middle point of the mother tube respectively according to model 1 (Fig. 1), model 2A (Fig. 2A) and model 2B (Fig. 2B). The waves shown in the figures are associated with the corresponding paths. For example, along the path 0, the time (s) of arrival of incident wave at middle point of the mother tube is given by:  t(0)=0.01s+ )*2/( 00 cl ; comparing this theoretical time with the time shown in the figures it is possible to associate ‘Wave A’ with the path 0; same method is applied to determine the time of arrival for all the other waves according to the corresponding path and wave speed. Table II and IV provide the paths, times and the simulated amplitude for the waves associated to model 1 (Figure 1, Table II), model 2A and 2B (Figure 2A, B, Table IV). ABCDEFGHIL MTime (s)Pressure (Pa)0.01 0.0102 0.0104 0.0106 0.0108 0.011 0.0112-100-50050100150200250Pressure (Pa)Pressure (Pa) Fig. 3. Pressure waves at the middle point of segment 0 (trachea) according to model 1 (Fig. 1)  TABLE II Wave paths with theoretical (t=l/c) arrival time at the middle point of the mother tube (0). (-) denotes a backward wave. The numbering of path refers to Fig. 1. Computational amplitudes refer to peaks in Fig. 3, for theoretical amplitudes see text for explanation.   IV. DISCUSSION AND CONCLUSIONS  The expression of flow pulse (Eq. 8) generates an incident wave of approximately 0.06 ms duration. Due to the high wave speed in air (Table I), this short pulse is necessary to assure that the arrival of reflected waves appears after the end of the incident wave.  From Figure 3 and Table II, the ‘wave B’ (path 0-0) corresponds to the wave reflected at the bifurcation between the mother tube and the two segments of the first generation (point G0); it is possible to compare its computational amplitude (Table II) with the theoretical amplitude (IB) calculated by:                     5.28*0== forwardGAB RII Pa                    (9) In which IA is the amplitude of incident Wave (‘Wave A’) and RG0forward is the reflection coefficient in the forward direction at the bifurcation point G0, according to Table III. Using the same method, for the amplitude of wave C (path 02-2-0, Figure 3) we have:     )1(**)1(*00 backwardGtforwardGAC RRRII ++=         (10) Rt is the terminal reflection coefficient that in model 1 (Fig. 1) is equal to 0 and 1 for the two terminal segments. Being the ‘Wave C’ (path 02-2-0), the wave associated to the reflection from the terminal segment, the knowledge of the time and wave amplitude provides information regarding the distance and the magnitude of the blockage. Theoretical amplitudes for each path of model 1 are shown in Table II.   TABLE III Reflection coefficients (R) evaluated in forward and backward directions at the bifurcation points G0, G1, G2 (Fig. 1-2) according to Eq. 6.  R  Forward Backward G0 0.118 -0.559 G1 or G2 -0.009 -0.495  Applying the same concepts, it is possible to consider and compare a more complicated system composed by two generations in different conditions of peripheral blockage according to model 2A (one blockage, Rt=1, Figure 2A) and model 2B (two blockages, Figure 2B) to study the effects on wave amplitude (Table IV). Using Table IV it is possible to observe that ‘Wave D’ in model 2A is composed by two different sub-waves (path 0-00 and 026-6-2-0). In model 2B the same ‘Wave D’ derives from the combination of three waves (path 0-00, 026-6-2-0 and 013-3-1-0) that arrive simultaneously at the middle point of the mother tube. The difference in amplitude for this Wave, in model 2A and model 2B (Table IV), agrees with the theoretical expectations, using the same principles described in Eq. 9 and 10. Wave E (wave path 026-66-6-2-0 in model 2A and wave path 026-66-6-2-0 plus 013-33-3-1-0 in model 2B, Table IV) is the expansion wave generated by the re-reflection in the backward direction at point G2 (model 2A) and points G1, G2 (model 2B), R=-0.495 from Table III. In model 2B the amplitude of this wave appears to double compared to model 2A (Table IV), for the symmetry of the system according to theoretical expectations. Wave Wave path Time (s) Theoretical Amplitude (Pa)Computational  Amplitude (Pa)A 0 0.010118 240.32 240.32B 0-0 0.010355 28.55 28.32C 02-2-0 0.010516 118.46 115.63D 0-00 0.010592 28.55 28.98E 02-22-2-0 0.010676 -66.26 -64.91F 02-2-00 0.010752 118.46 115.21G 02-22-22-2-0 0.010837 37.06 38.13H 02-22-2-00 0.010913 -66.26 -63.890-002-2-0 0.010989 14.0702-22-22-22-2-0 0.010998 -20.7302-2-00-0 0.010989 14.07 L 02-22-22-2-00 0.011074 37.06 38.1602-2-002-2-0 0.011151 58.3902-22-22-22-22-2-0 0.011158 11.5902-22-2-00-0 0.011151 -7.8614.83IM 49.995484Authorized licensed use limited to: Brunel University. Downloaded on January 28, 2010 at 10:09 from IEEE Xplore.  Restrictions apply.   Pressure (Pa)BCDEFGHI LATime (s)0.01 0.0102 0.0104 0.0106 0.0108 0.011 0.0112-100050100150200250-50Pressure (Pa) Fig. 4. Pressure waves at the middle point of segment 0 (trachea) according to model 2A (Fig. 2A). ATime (s)Pressure (Pa)BCDEFGHI L0.01 0.0102 0.0104 0.0106 0.0108 0.011 0.0112-100-50050100150200250Pressure (Pa) Fig. 5. Pressure waves at the middle point of segment 0 (trachea) according to model 2B (Fig. 2B)  The 1-D modelling, applied to human airways bifurcations, has been shown to catch the theoretical timing and amplitude of waves even when the fluid is not blood. These preliminary results are encouraging and indicate that once the 1-D model is extended to study more bifurcating generations of airways, it would be a useful tool to predict and to better understand physiological and pathological respiratory conditions and their effects on velocity and pressure in each segment.   A. Limitations  In these investigations we tested a limited number of bifurcations for the preliminary results. We used only total reflective (Rt=1) and non-reflective (Rt=0) reflection coefficients however we acknowledge that terminal reflections may be different. The inflow pulse used in this study was chosen to be very short in time which we acknowledge to be much shorter than a normal breathing cycle.    TABLE IV Wave paths with theoretical (t=l/c) arrival time to the middle point of the mother tube (0). (-) denotes a backard wave. The numbering of path refers to fig. 2A-B. Wave amplitude refers to computational results for model 2A (figure 2A) and model 2B (figure 2B).  Wave Wave path Time (s)Computational Amplitude (Pa)                  model 2AComputational Amplitude (Pa)  model 2B                  A 0 0.01011 240.32 240.32B 0-0 0.01035 28.32 28.32 02-2-0 0.01051601-1-0 0.010516  0-00 0.010592026-6-2-0 0.010589013-3-1-0 0.010589026-66-6-2-0 0.01066 -28.78013-33-3-1-0 0.01066026-66-66-6-2-0 0.01073 14.18013-33-33-3-1-0 0.010730-00-0 0.010829026-6-2-00 0.010826026-6-226-6-2-0 0.010824026-66-66-66-6-2-0 0.010811013-3-1-00 0.010826013-3-113-3-1-0 0.010824013-33-33-33-3-1-0 0.010811026-6-213-3-1-0 0.010824013-3-126-6-2-0 0.010824026-66-6-2-00 0.010905026-66-66-66-66-6-2-0 0.010884026-66-6-226-6-2-0 0.010898013-33-3-1-00 0.010905013-33-33-33-33-3-1-0 0.010884013-33-3-113-3-1-0 0.010898026-66-66-6-2-00 0.010974026-66-66-66-66-66-2-0 0.01095013-33-33-3-1-00 0.010974013-33-33-33-33-33-1-0 0.01095026-6-2-00-0 0.011063026-6-226-6-2-00 0.0110609013-3-1-00-0 0.011063012-3-113-3-1-00 0.011060943.6385.66I3.05L 21.37-48.23H-12.925.7622.57C -3.14G 107.601143.511E -57.5128.72FD-2.64 REFERENCES [1] S. J. Sherwin,. V. Franke, J. Peirό, and K. H. Parker, “One-dimensional modelling of a vascular network in space-time variables” Journal of Engineering Mathematics vol.47, pp.217-250., 2003 [2] K.S Matthys, J. Alastruey, J. Peirό, A. W. Khir, P. Segers, P. R. Verdonck, K.H. Parker, and S.J. Sherwin “Pulse wave propagation in a model human arterial network: Assessment of 1-D numerical simulations against in vitro measurements” Journal of Biomechanics vol.40, pp.3476-3486, 2007 [3] Weibel, E.R. “Morphometry of the human lung”. Berlin, Springer-Verlag, 1963.  [4] M. Montaudon, P. Desbarats, P. Berger, G. de Dietrich, R.  Marthan, and F. Laurent, “Assessment of bronchial wall thickness and lumen diameter in human adults using multi-detector computed tomography: comparison with theoretical models” J. Anat. vol. 211, pp. 579- 588, 2007. [5] R.H. Habib, R.B. Chalker, B. Suki and A.C. Jackson, “Airway geometry and wall mechanical properties estimated from subglottal input impedance in humans” J. Appl. Physiol., vol.77 (1), pp.441-451, 1994.  [6] R. K. Lambert, E. M. Baile, R. Moreno, J. Bert and P.D. Pare’, “A method for estimating the Young’s Modulus of complete tracheal cartilage rings” J. Appl. Physiol.  vol. 70 (3), pp. 1152-1159, 1991. [7] B.R. Wiggs, R. Moreno, J.C. Hogg, C. Hillian., and P.D. Pare’, “A model of the mechanics. of airway narrowing” J. Appl. Physiol. vol. 69 (3), pp. 849-860, 1990. [8] A. W. Khir, K.H. Parker, “Measurements of wave speed and reflected waves in elastic tubes and bifurcations”  Journal of Biomechanics vol. 35,  pp 775-783, 2002.   5485Authorized licensed use limited to: Brunel University. Downloaded on January 28, 2010 at 10:09 from IEEE Xplore.  Restrictions apply. 

One-dimensional modelling of pulse wave propagation in human airway bifurcations in space-time variables

http://bura.brunel.ac.uk/bitstream/2438/9615/1/Fulltext.pdf

One-dimensional modelling of pulse wave propagation in human airway bifurcations in space-time variables

Abstract

Similar works

Full text

Available Versions

Brunel University Research Archive