- Email: [email protected]
Accuracy of the wave equation in predicting arterial pulse propagation
Mathematical and Computer Modelling 57 (2013) 460–468
Contents lists available at SciVerse ScienceDirect
Mathematical and Computer Modelling journal homepage: www.elsevier.com/locate/mcm
Accuracy of the wave equation in predicting arterial pulse propagation Ahmed Al-Jumaily a,∗ , Andrew Lowe b a
Institute of Biomedical Technologies, Auckland University of Technology, Private Bag 92006, Auckland, New Zealand
b
Pulsecor Limited, Auckland, New Zealand
article
info
Article history: Received 2 June 2011 Received in revised form 15 December 2011 Accepted 13 June 2012 Keywords: Blood pressure Aorta Mathematical model Arterial pulses Geometric tapering Arterial stiffening
abstract This paper investigates the effect of the various terms in the one-dimensional acoustic wave equation on the pulse characteristics within the aorta. To mimic the physiological nature of the systemic arteries, the aorta is modelled as an elastic conical tube. The frequency spectrum is used to study the effect of different terms in this equation on the pressure ratio between the aortic root at the heart exit and the iliac bifurcation. For validation, the effective reflection distance calculated using this model is within 7% error of clinical observations, suggesting that the model is able to mimic physiological pressure propagation with a high degree of accuracy and can therefore be used to generate and test hypotheses. This work demonstrates that: (i) tapering in the aorta lumen radius causes supplementary amplification of the pressure pulses in the system and increases the propagation velocity; however, tapering in the aorta wall thickness generates opposite effects, (ii) increasing the wall stiffness causes a change in the natural frequency of the system and increases the propagation velocity, and (iii) inclusion of either the advective momentum correction term or the viscosity term insignificantly affects the pressure ratio. The last observation suggests that the flow pattern does not influence the pressure propagation characteristics. © 2012 Elsevier Ltd. All rights reserved.
1. Introduction Blood flow induces body forces and stresses in the arterial walls due to complex fluid–structure interactions. These forces and stresses play an important role in the onset and progression of many acquired and congenital cardiovascular diseases such as atherosclerosis and aneurysms. Atherosclerosis involves the accumulation of plaque in the intima of the arterial wall, which reduces arterial lumen and increases local arterial stiffness. There is substantial evidence on the localisation of these plaque deposits at sites with haemodynamic conditions commonly characterised by low wall shear stress [1,2]. Aneurysms, on the other hand, involve the degradation of local arterial wall tissues resulting in lowering of local arterial stiffness and enlargement of local vessel cross-sections. If, in extreme cases, the wall stresses induced due to the transient fluid–structure interactions exceed the strength limit of the dilated artery wall, they cause vessel rupture leading to death from internal haemorrhage which has been reported to be in more than 80% of the cases [3]. Aneurysms are common in locations with secondary flow and flow recirculation even in normal resting conditions [4–6]. Considering the correlations between the various haemodynamic conditions and the onset and progression of different cardiovascular diseases, there is worldwide consensus on the need for enhancements in the current understanding of cardiovascular mechanopathobiology [7]. Since blood flow and pressure are dependent on a number of factors in addition to physical phenomena, experimental studies demand rigorous screening, which make them expensive and time consuming. A preferred mode of investigation is by computer simulation using mathematical models.
∗
Corresponding author. Tel.: +64 99219777; fax: +64 99219973. E-mail address: [email protected] (A. Al-Jumaily).
0895-7177/$ – see front matter © 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.mcm.2012.06.023
A. Al-Jumaily, A. Lowe / Mathematical and Computer Modelling 57 (2013) 460–468
461
Fig. 1. The axisymmetrical cylindrical 1-D domain; section S (x, t ) is defined at a location x at any time t with a cross-sectional area A(x, t ).
Quateroni et al. [8] and Taylor et al. [9] provide excellent reviews of current modelling techniques. Principles of conservation of mass, momentum and energy form the theoretical bases of all mathematical models that have been developed for the study of the physical aspects of blood flow. Van de Vosse and Stergiopulos [10] presents an excellent review which focuses on the physical and mathematical modelling of pulse wave propagation based on a general fluid dynamic approach. However, an analytical solution of the full form of these equations has not been developed and only solutions for special cases are available today. Proper application of numerical techniques for solving these equations demands the familiarity with assumptions and approximations that can be made to achieve a reasonable solution. This paper develops a mathematical model capable of representing pressure propagation in human arteries and uses this model to investigate the significance of different physical terms within this model on pressure propagation in the aorta. Although the development of 1-D mathematical models for studying pressure propagation in arteries and the effect of some phenomena have been previously discussed in the literature [11–19], to the best of our knowledge the application of a structured sensitivity analysis is presented for the first time in this work. Further, previous publications assume the aorta to be straight and either linearly or exponentially tapered which is not anatomically representative. For a better shape regeneration, bicubic splines are introduced for both wall thickness and lumen radius and this explains the higher accuracy of the model predictions of complex physiological parameters associated with arterial wave propagation. 2. Model development The pulse propagation in the arteries can be simulated as the propagation of pressure pulses in an incompressible, viscous or inviscid fluid (blood) flowing in a deformable tube (artery) [9,19]. To develop a reasonable simulation, the following assumptions are made based on clinical observations: the arterial pulses have a significantly longer wavelength, in the order of metres [20], when compared to the millimetre-scale arterial dimensions [21]. Since the aorta expands and contracts around 10% (±5% around the mean) with each beat [22], the amplitude of these pulses can be assumed to be very small. Shear rates in the proximal arteries are large enough to permit the assumption that the blood is Newtonian [23]. Velocity measurements in the proximal arteries have confirmed that the axial velocity is considerably larger than the radial velocity and thus 1-D equations can be used. Further, the artery wall is assumed thin, homogenous, elastic and impermeable, whereby a zero-slip condition may be applied to these equations. Under these assumptions, applying the one dimensional conservation of mass and momentum equations leads to the following governing equations of blood flow,
∂A ∂Q + =0 ∂t ∂x 2 ∂Q ∂ Q A ∂p + α + + kR u = 0 ∂t ∂x A ρ ∂x
(1) (2)
where A (x, t ) is the cross-sectional area, Q (x, t ) is the volume flow rate across the section, u (x, t ) = Q (x, t )/A (x, t ) is the axial velocity of the blood averaged across the section, ρ is the density of blood which is considered as a constant in this investigation, α is the advective momentum correction coefficient, p (x, t ) is the mean pressure at the section and kR is the loss of momentum per unit length due to viscosity of blood (ν). For a Poiseuille flow profile, α = 1.33 and kR = 8π ν while for the clinically measured flow profile, a value of α = 1.1 and kR = 21π ν has been suggested [11,24–26]. α = 1 assumes a completely flat profile, but it is often used as it permits significant mathematical simplification [27] (see Fig. 1). From physical theories, it is known that as the wave pulses move through the arterial tree, they cause local distension of the arterial wall. A relationship linking this distension and the pressure at a cross-section can be followed assuming the aorta as a thin walled pressurised cylinder where the stresses induced remain constant through the thin arterial wall.
462
A. Al-Jumaily, A. Lowe / Mathematical and Computer Modelling 57 (2013) 460–468
Neglecting inertial and viscoelasticity effects [24], the area can be written as A = Ao + ϕ(pnet )
(3)
where Ao is the resting cross-sectional area, ϕ (pnet ) is the area increment caused by the increment in pressure, pnet and
ϕ(pnet ) =
2Ao r Ew tw
· 1pnet =
2Ao r Ew t w
· ∆(p − pext ) = Ao c · ∆(p − pext )
(4)
where r is the internal radius of the artery, pext is the extramural pressure, Ew is the modulus of elasticity of the artery wall, tw is the thickness of the artery wall and c is the compliance of the cross-section. Since the extramural pressure, acting on the artery remains reasonably unchanged as the aorta is not rigidly held in position, pext ≈ constant, and hence the expansion and contraction of the artery is in response to internal pressure change p. Since the primary interest in this study is in pressure propagation, it can be assumed that the cross-sectional area, the pressure and the velocity of the flow perturbates around a steady state value (e.g. A = A0 + Ap ). The subscript 0 defines a steady state operating-point (mean) value and the subscript p defines a small perturbation in the system variables. In vivo measurements have shown that the blood flow (and hence all other coupled variables) oscillates periodically around a mean value [28] with small pulses compared to the steady state values. Previous studies [29] have confirmed that the propagation of the arterial waves can be simulated using first order differential equations. Although this system is nonlinear in nature, perturbation theories are used to simplify the differential equations. This procedure has been widely used in the literature and it has proven to generate the standard wave equation which is an effective tool in understanding wave propagations in the arterial system. Adding nonlinear complexity may improve the mathematical formulation but definitely will not have significant additions to understanding those waves. To obtain the final governing equations, substitute Eqs. (1) and (2) into Eqs. (3) and (4), implement perturbation technique for each of the variables around a steady state value, neglect high order terms, Laplace transform and after some mathematical manipulation obtain A0 csP + (cu0 P + Ucp0 + U )
(1 + cp0 )sU + u0 (1 + cp0 )
dA0 dx
dU dx
+ cu0 A0
+
dP dx
+ A0 (1 + cp0 )
(1 + cp0 ) dP p
dx
+
kR A0
dU dx
=0
U =0
(5) (6)
where U and P are the Laplace transforms of the instantaneous velocity and pressure oscillations respectively. The spatial derivatives are approximated using forward finite differences and matrix partition techniques are applied to reduce the equations to a [A] · [B] = [C ] form as shown below kp (n, s)Pn − [uo An c ]Pn−1 + [(2An − An−1 )(1 + po c )]Un − [An (1 + po c )]Un−1 = 0 kv (n, s)Un − [uo (1 + po c )]Un−1 +
1 + po c
ρ
Pn −
1 + po c
ρ
(7)
P n −1 = 0
(8)
where s is the Laplace variable, h is the step size, the subscript n represents the variable at the nth step and the subscript n − 1 represents the variable at the n − 1th step and kv (n, s) = h(1 + po c )s + uo (1 + po c ) + kp (n, s) = An hcs + cuo (2An − An−1 ).
hkR An
(9) (10)
Eqs. (7) and (8) are sufficient to represent fluid flow through a linearly tapered element. As these equations are homogeneous, the results drawn from them are continuous in one-dimensional Gaussian space. Replacing s by iω in these equations transfers the solution from the s domain to the frequency domain, where i is the imaginary unit and ω is the angular frequency. It is worth mentioning here that each one of the elements in the finite difference formulation is very small compared to the overall tapered aorta and is assumed to be of a uniform shape. Thus, the aorta tapering is reduced to small step changes and is not needed to be included in the above equations when applied to each element. The system is closed by specifying a relationship between pressure and flow at any one location along the aorta. The equations are then solved by implementing a modified Gauss elimination (L-U decomposition) strategy to give transfer functions (TFs) between pressure and flow at different locations along the aorta. 3. Simulation specifications This section is aimed at accurately describing the different parameters used in the simulation process to ensure that the results of this paper are reproducible.
A. Al-Jumaily, A. Lowe / Mathematical and Computer Modelling 57 (2013) 460–468
463
Table 1 Geometric domain description. Source: Data adapted from Westerhof et al. [21]. Artery segment
Length (cm)
Internal radius (cm)
Wall thickness (cm)
Young’s modulus (×105 N/m2 )
Ascending aorta (1) Ascending aorta (2) Aortic arch (1) Aortic arch (2) Thoracic aorta (1) Thoracic aorta (2) Thoracic aorta (3) Abdominal aorta (1) Abdominal aorta (2) Abdominal aorta (3) Common iliac artery (1)
2 2 2 3.9 5.2 5.2 5.2 5.3 5.3 5.3 5.8
1.470 1.440 1.120 1.070 0.999 0.675 0.645 0.610 0.580 0.548 0.368
0.164 0.161 0.132 0.127 0.120 0.090 0.087 0.084 0.082 0.078 0.063
4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0
Fig. 2. Invasively measured ascending aortic pressure after digitisation.
3.1. Geometric domain description Table 1 shows geometric properties at discrete locations along the aorta [22]. Although the modulus of elasticity of the physiological system is not uniform, the data presented in this table gives an acceptable approximation for linear modelling purpose. These data were interpolated using cubic splines to generate the geometric domain of the model. 3.2. Frequency spectrum of interest Within the experimental context of this investigation, the pressure fluctuation at the ascending aorta was invasively measured during catheterisation and was extracted in a beat-to-beat manner. The recorded pressure was digitised at 100 Hz and the resultant waveform is shown in Fig. 2. This beat was analysed to check its frequency contents. It was observed that a significant portion (94.1%) of the wave energy was contained in frequencies below 15.6 Hz. Previously reported values indicate that the range of 0–15 Hz is physiologically relevant. Thus, the present data agrees well with previous readings. 3.3. Impedance boundary condition Although branching plays an important role in wave reflections, this paper focuses on the significance of each term for a single tube. Previous simulations for analysing blood pressure and flow use arbitrary boundary conditions. For example, flow boundary conditions [12,17,18], pressure boundary conditions [30,31], modified forms of Windkessels as boundary conditions [13] and use of the characteristic impedance boundary conditions [15,32] have been used. The use of the impedance boundary condition has the advantage that the model results are generalised; however, the use of pressure or flow specifications are better suited for patient-specific models. Olufsen et al. have proposed the use of the fractal approach to formulate a ‘‘structured tree’’ boundary condition [17,18,33]. However, the assumed similarity laws in the derivation of this boundary condition are difficult to be physiologically realised. In light of the above-mentioned facts, invasively measured impedance is used as a boundary condition as it gives the best representation of the downstream vasculature. At the end of the aorta, the invasively measured impedance can be written as Z = f (ω).
(11)
Mills et al. measured the impedance at the iliac artery up to 13 Hz [34]. Based on the Windkessel theory, the impedance at high frequencies approaches a constant value and in this research, the impedance at frequencies at and above 13 Hz is assumed to remain constant. Without mention of the precise location of Mills’ measurement, it is further assumed that the measurement was made at the start of the iliac artery. The impedance at the end of the aorta can be calculated by applying the parallel impedance analogy from electrical circuit analysis. Each one of the finite difference elements needs a starting condition and end condition. The last element in the simulation uses the values of the impedance given above.
464
A. Al-Jumaily, A. Lowe / Mathematical and Computer Modelling 57 (2013) 460–468
Fig. 3. Map showing the simulated pressure waveforms (mm Hg) at different locations along the aorta when the model is used to simulate the pressure propagation in the geometric domain specified in Table 1. A is an arc tracing the timing of peak pressure, which may be used as a reference to calculate PWV. B is a region of low pressure along the aorta and C is a region with post-systole pressure fluctuation.
4. Model validation The invasively measured pressure waveform shown in Fig. 2 is broken down into its harmonics and represented as a Fourier series by, yx=0 = Ao +
n [Ai cos ωi t − Bi sin ωi t ].
(12)
i=1
Previous research [32] has shown that the first 5 harmonics contain over 90% of the energy of the pressure waveform. The pressure wave is assumed to contain 10 harmonics so as to contain more than 99% of the energy. Using the pressure transfer functions (TFx ) for different locations along the aorta and input Fourier series, the blood pressure at different locations can be calculated as p(x, t ) = A0 +
∞ [|TF|x An sin(nωt + ]TFx ) + Bn cos(nωt + ]TFx )]
(13)
n=1
where p (x, t ) is the pressure wave at any distance x from the aortic root, An and Bn are the coefficients of the nth harmonic of the pressure and TFx is the ratio of pressure at distance x from the aortic root to that at the aortic root, |TF| is the derived normalised amplitude ratio and |TF| is the associated phase lag. Fig. 3 shows a temporal and spatial map of the pressure predicted along the aorta. To validate the model, comparisons are made between the model outputs and available clinical data. The basic parameters for comparison are the pulse wave velocity (PWV), blood pressure trends and the effective reflecting distance (ERD). The model-predicted PWVs at different locations along the aorta are within the physiologically measured range of 4.6–6.5 m/s [35] and they progressively increase as the wave travels to the periphery (refer to Fig. 3, arc A) in agreement with clinical observations. From the pressure waveforms it is observed that the wave steepens as it propagates away from the heart and that a region of minimum pressure and consequently minimum wall shear stress (for a cylindrical thin-wall vessel, the wall shear stress is a function of the pressure) develops in the thoracic aorta (refer to Fig. 3, ellipse B). Careful simulation has shown that increasing the aortic stiffness decreases the minimum pressure and expands the region of low pressure. Previous research has shown that such sites, with low wall shear stress, are prone to plaque deposition [1,2]. Consequently, this demonstrates that the risk for the development of atherosclerosis considerably increases with increasing aortic stiffness (which normally takes place with ageing). A considerably different pressure pattern is noticed in the abdominal aorta where a distinct diastolic oscillation is observed (refer Fig. 3, ellipse C) which becomes more pronounced with the increase of aortic stiffness. This diastolic oscillation highlights the presence of the reflected wave coming from the periphery that is coupled to a negative flow wave. This negative flow wave slows down local fluid particles and may cause instantaneous flow reversal. Locations demonstrating these characteristics are commonly associated with vascular remodelling leading to the development of aneurysms, which are more commonly reported in the abdominal aorta [4–6]. At this stage, it may be difficult to identify in-plane locations of plaque deposition and aneurysm formation using 1-D models; however, it is possible to link the axial locations of arterial remodelling with the predicted pressure fields.
A. Al-Jumaily, A. Lowe / Mathematical and Computer Modelling 57 (2013) 460–468
465
Fig. 4. Effect of aortic radius and wall thickness tapering on the pressure propagation characteristics.
ERD is calculated based on the reflected wave transit time (RWTT) principle using the formula ERD = PWV × RWTT/2 which was proposed by Murgo et al. [36]. Using the available numerical data for the aorta [20], the model predicts ERD to be 47.5 cm from the heart while previous mathematical models predict ERD to be 39.8 cm away from the heart [32]. The present result correlates very well with previously reported values based on invasive data by Mills et al. [34] who estimated ERD to be 50.9 cm away from the heart. RWTT analysis carried out on the Framingham Cohort showed a mean ERD of 50 cm in women and 59 cm in men [37], to which the present model prediction compares favourably. The lack of precise data required to reconstruct the arterial model from the Framingham Cohort restricts the investigation to this level; however it is believed that if patient-specific data is available, the model can be used to predict the precise location of ERD. Therefore, the model prediction matches major clinical observations and hence it is able to simulate the propagation of pressure waves in the human aorta with a high degree of accuracy. 5. Pressure propagation under different conditions Considering each condition individually, the significance of the different physical phenomena on the pressure propagation in the aorta is studied by investigating the change in the pressure TF between the aortic root and the iliac bifurcation. 5.1. Geometric tapering Fig. 4 shows the effect of tapering of the artery internal radius and wall thickness on the TF for three different cases: (i) in a straight tube of mean internal radius and wall thickness, 1.009 cm and 0.121 cm respectively, (ii) a tube of constant internal radius with linearly tapering wall thickness and (iii) a tube with linearly tapering internal radius and wall thickness. The degrees of tapering in internal radius and wall thickness respectively were calculated based on the physiological data given in Table 1. Fig. 4 shows that tapering in the internal radius causes the frequency spectrum of the TF to shift towards higher frequencies (indicating an increase of resonant or natural frequency of the aorta). There is also an amplification of the pulse amplitude and a reduction in the phase delay, indicating an increase in PWV. However, tapering in the wall thickness has the effect of reducing the pulse amplitude, increasing the phase delay indicating a decrease in PWV and increasing the natural frequency of the system. The difference between the cases is observable for frequencies above 4 Hz, which are within the physiological range. These changes may be explained by the fact that tapering changes the compliance of the system following Eq. (4). It is also observed that tapering in the internal radius influences the pressure propagation characteristics of the aorta to a much higher degree than the tapering of wall thickness. This may be due to the fact that the internal radius tapers down at 1.2758° while the wall thickness tapers down at 0.1190° over the length of the aorta. 5.2. Aortic stiffness The aortic stiffness under the normal condition is assumed as 4 × 105 N/m2 [20]. Fig. 5 shows the change in TF caused by doubling the aortic stiffness. It is seen that increasing the aortic stiffness shifts the natural frequency to higher values which is in agreement to the quarter wavelength formula [38] while a reduced phase delay suggests the increase of PWV. These changes are observed because the system compliance is dependent on the aortic stiffness according to Eq. (4).
466
A. Al-Jumaily, A. Lowe / Mathematical and Computer Modelling 57 (2013) 460–468
Fig. 5. Effect of change of arterial stiffness.
Fig. 6. Effect of change of the momentum correction factor, α .
5.3. Changing the assumed flow profile The two model parameters used to specify the flow pattern in the aorta are the advective momentum correction factor
α and the momentum loss factor kR from Eq. (2).
Fig. 6 shows the cases when the advective acceleration term is included or disregarded in the simulation. It is indicated that the advective acceleration term itself is very small compared to the other terms in Eq. (2) and does not affect the pressure propagation in the aorta. Minor changes in the momentum correction factor, i.e. Poiseuille to the experimentally measured flow profile with 1α < 1, should therefore have smaller effects on the propagation characteristics of the aorta than complete omission, which is studied here. The viscosity of blood is 3.2 × 10−3 Pa-s [39] and is incorporated into the simulation through the term kR . Fig. 7 shows the change in the propagation characteristics of the aorta by omitting kR . The graph shows a small change in the amplitude of the TF peaks with no change in the overall trend. There is no observable change in the propagation velocity of the pulses confirming that the viscosity of the blood insignificantly influences pressure propagation in the elastic arteries. Therefore, these observations indicate that the changes in α and kR have an insignificant effect on the TF and hence the assumed flow profiles should not change the pressure propagation characteristics. Future studies can assume Poiseuillian flow as it leads to considerable mathematical simplifications. 6. Conclusions In this paper, blood flow in the aorta is simulated using a 1-D fluid flow model to investigate the effect of various physical terms on the pulse wave propagation characteristics. Clinical parameters predicted by the simulation are close to invasive
A. Al-Jumaily, A. Lowe / Mathematical and Computer Modelling 57 (2013) 460–468
467
Fig. 7. Effect of change of the viscous correction term, kR .
measurements. The model has shown that the changes in the geometry and material properties strongly influence the pressure propagation characteristics while the flow profile has an insignificant effect. Acknowledgement We would like to acknowledge IBTec former employee Ashis Mookerjee for his computational contribution to this article. References [1] C.G. Caro, J. Fritz-Gerald, R. Schroter, Atheroma and arterial wall shear: observation, correlation and proposal for a shear dependent mass transfer mechanism for atherogenesis, Proceedings of the Royal Society of London B. Biological Science 177 (1971) 109–159. [2] C.A. Taylor, C.P. Cheng, L.A. Espinosa, B.T. Tang, D. Parker, R.J. Herfkens, In vivo quantification of blood flow and wall shear stress in the human abdominal aorta during lower limb exercise, Annals of Biomedical Engineering 30 (3) (2002) 402–408. [3] D. Katz, J. Cronenwett, The cost effectiveness of early surgery verses watchful waiting in the management of small abdominal aortic aneurysms, Journal of Vascular Surgery 19 (1994) 980–990. [4] C.M. Scotti, A.D. Shkolnik, S.C. Muluk, E.A. Finol, Fluid–structure interaction in abdominal aortic aneurysms: effects of asymmetry and wall thickness, Biomedical Engineering Online 4 (1) (2005) 64. [5] J.C. Lasheras, The biomechanics of arterial aneurysms, Annual Review of Fluid Mechanics 39 (2007) 293–319. [6] R. Peattie, T. Riehle, E. Bluth, Pulsatile flow in fusiform models of abdominal aortic aneurysms: flow fields, velocity patterns and flow-induced wall stresses, Journal of Biomechanical Engineering 126 (4) (2004) 438–446. [7] M. O’Rourke, J. Hashimoto, Mechanical factors in arterial ageing-a clinical perspective, Journal of the American College of Cardiology 50 (1) (2007) 1–13. [8] A. Quarteroni, M. Tuveri, A. Veneziani, Computational vascular fluid dynamics: problems, models and methods, Computing and Visualization in Science 2 (2000) 163–197. [9] C.A. Taylor, M. Draney, Experimental and computational methods in cardiovascular fluid dynamics, Annual Review of Fluid Mechanics 36 (2004) 197–231. [10] F.N. van de Vosse, N. Stergiopulos, Pulse wave propagation in the arterial tree, Annual Review of Fluid Mechanics 43 (2011) 467–499. [11] P. Reuderink, H. Hoogstraten, P. Sipkema, B. Hillen, N. Westerhof, Linear and nonlinear one-dimensional models of pulse wave transmission at high Wormsley numbers, Journal of Biomechanics 22 (8–9) (1989) 819–827. [12] S.J. Sherwin, L. Formaggia, J. Peiro, V. Franke, Computational modelling of 1D blood flow with variable mechanical properties and its application to the simulation of wave propagation in the human arterial system, International Journal for Numerical Methods in Fluids 43 (2001) 673–700. [13] P. Nardinocchi, G. Pontrelli, L. Teresi, A one-dimensional model for blood flow in prestressed vessels, European Journal of Mechanics A/Solids 24 (2005) 23–33. [14] J.R. LaCourse, G. Monahakrishnan, K. Sivaprasad, Simulations of arterial pressure pulses using a transmission line model, Journal of Biomechanics 19 (9) (1986) 771–780. [15] J. Wang, K. Parker, Wave propagation in a model of the arterial circulation, Journal of Biomechanics 37 (2004) 457–470. [16] K. Matthys, J. Alastruey, J. Peiro, A. Khir, P. Segers, P. Verdonck, et al., Pulse wave propagation in a model human arterial network: assessment of 1-D numerical simulations against in vitro measurements, Journal of Biomechanics 40 (15) (2007) 3476–3486. [17] M.S. Olufsen, Structured tree outflow conditions for blood flow in large systemic arteries, American Journal of Physiology 276 (1999) H257–H268. [18] M.S. Olufsen, C. Peskin, W. Kim, E. Pederson, A. Nadim, J. Larsen, Numerical simulation and experimental validation of blood flow in arteries with structured tree outflow conditions, Annals of Biomedical Engineering 28 (11) (2000) 1281–1299. [19] J.B. Grotberg, O.E. Jensen, Biofluid mechanics in flexible tubes, Annual Review of Fluid Mechanics 36 (2004) 121–147. [20] P. Segers, S. Carlier, A. Pasquet, S.I. Rabben, L.R. Hellevik, E. Remme, et al., Individualising the aorto-radial pressure transfer: feasibility of a model based approach, American Journal of Physiology 279 (2000) H542–H549. [21] N. Westerhof, F. Bosman, C. de Vries, A. Noodergraaf, Analog studies of the human systemic arterial tree, Journal of Biomechanics 2 (1969) 121–143. [22] P. Boutouyrie, S. Laurent, A. Benetos, X.J. Girerd, A.P. Hoeks, M. Safar, Opposing effects of ageing on distal and proximal large arteries in hypertensives, Journal of Hypertension Supplements 10 (1992) S87–S91. [23] F.J.H. Gijsen, F.N. van de Vosse, J.D. Jassen, The influence of the non-Newtonian properties of blood on the flow in large arteries: steady flow in a carotid bifurcation model, Journal of Biomechanics 32 (1999) 601–608. [24] S. Balar, T. Rogge, D. Young, Computer simulation of blood flow in the human arm, Journal of Biomechanics 22 (6–7) (1989) 691–697. [25] L. Formaggia, D. Lamponi, A. Quarteroni, One-dimensional models for blood flow in arteries, Journal of Engineering Mathematics 47 (2003) 251–276.
468
A. Al-Jumaily, A. Lowe / Mathematical and Computer Modelling 57 (2013) 460–468
[26] E.M. Pedersen, W.Y. Kim, N.H. Staalsen, J.M. Hasenkam, H. Nygaard, P.K. Paulsen, Development of velocity profiles and retrograde flow in the porcine abdominal aorta under different haemodynamic conditions, Scandinavian Cardiovascular Journal 33 (1999) 206–214. [27] L. Formaggia, J. Gerbeau, F. Nobile, A. Quarteroni, On the coupling of 3D and 1D Navier–Stokes equations for flow problems in compliant vessels, Computer Methods in Applied Mechanics and Engineering 191 (2001) 561–582. [28] J. Moore, D. Ku, Pulsatile velocity measurements in a model of the human abdominal aorta under simulated exercise and postprandial conditions, Journal of Biomechanical Engineering 116 (1) (1994) 107–111. [29] J. Womersley, The mathematical analysis of the arterial circulation in a state of oscillatory motion, Wright Air Development Centre, 1957, Report No. WADC-TR56-614. [30] T. Mabotuwana, L. Cheng, A. Pullan, A model of blood flow in the mesenteric arterial system, Biomedical Engineering Online 6 (1) (2007) 17. [31] M.S. Olfusen, A. Nadim, L.A. Lipsitz, Dynamics of cerebral blood flow regulation explained using a lumped parameter model, American Journal of Physiology 282 (2002) R611–R622. [32] A.P. Avolio, Multi-branched model of the human arterial system, Medical and Biological Engineering and Computing 18 (6) (1980) 709–718. [33] B. Steele, M.S. Olufsen, C. Taylor, Fractal network model for simulating abdominal and lower extremity blood flow during resting and exercise conditions, Computer Methods and Programs in Biomedicine 10 (1) (2007) 39–51. [34] C.J. Mills, I.T. Gabe, J.H. Gault, D.T. Mason, J. Ross, E. Braunwald, et al., Pressure-flow relationships and vascular impedance in man, Cardiovascular Research 4 (1970) 405–417. [35] A. Avolio, S. Chen, R. Wang, C. Zhang, M. Li, M. O’Rourke, Effects of ageing on changing arterial compliance and left ventricular load in a northern Chinese urban community, Circulation 68 (1983) 50–58. [36] J.P. Murgo, N. Westerhof, J.P. Giolma, S.A. Altobelli, Aortic input impedance in normal man: relationship to pressure wave forms, Circulation 62 (1980) 105–116. [37] G.F. Mitchell, H. Parise, E.J. Benjamin, M.G. Larson, M.J. Keyes, J.A. Vita, et al., Changes in arterial stiffness and wave reflection with advancing age in healthy men and women: the framingham heart study, Hypertension 43 (2004) 1239–1245. [38] P. Segers, P. Verdonck, Role of tapering in aortic wave reflection: hydraulic and mathematical model study, Journal of Biomechanics 33 (2000) 299–306. [39] A.R. Pries, D. Neuhaus, P. Gaehtgens, Blood viscosity in tube flow: dependance on diameter and hematocrit, American Journal of Physiology 263 (1992) H1770–H1778.
Contents lists available at SciVerse ScienceDirect
Mathematical and Computer Modelling journal homepage: www.elsevier.com/locate/mcm
Accuracy of the wave equation in predicting arterial pulse propagation Ahmed Al-Jumaily a,∗ , Andrew Lowe b a
Institute of Biomedical Technologies, Auckland University of Technology, Private Bag 92006, Auckland, New Zealand
b
Pulsecor Limited, Auckland, New Zealand
article
info
Article history: Received 2 June 2011 Received in revised form 15 December 2011 Accepted 13 June 2012 Keywords: Blood pressure Aorta Mathematical model Arterial pulses Geometric tapering Arterial stiffening
abstract This paper investigates the effect of the various terms in the one-dimensional acoustic wave equation on the pulse characteristics within the aorta. To mimic the physiological nature of the systemic arteries, the aorta is modelled as an elastic conical tube. The frequency spectrum is used to study the effect of different terms in this equation on the pressure ratio between the aortic root at the heart exit and the iliac bifurcation. For validation, the effective reflection distance calculated using this model is within 7% error of clinical observations, suggesting that the model is able to mimic physiological pressure propagation with a high degree of accuracy and can therefore be used to generate and test hypotheses. This work demonstrates that: (i) tapering in the aorta lumen radius causes supplementary amplification of the pressure pulses in the system and increases the propagation velocity; however, tapering in the aorta wall thickness generates opposite effects, (ii) increasing the wall stiffness causes a change in the natural frequency of the system and increases the propagation velocity, and (iii) inclusion of either the advective momentum correction term or the viscosity term insignificantly affects the pressure ratio. The last observation suggests that the flow pattern does not influence the pressure propagation characteristics. © 2012 Elsevier Ltd. All rights reserved.
1. Introduction Blood flow induces body forces and stresses in the arterial walls due to complex fluid–structure interactions. These forces and stresses play an important role in the onset and progression of many acquired and congenital cardiovascular diseases such as atherosclerosis and aneurysms. Atherosclerosis involves the accumulation of plaque in the intima of the arterial wall, which reduces arterial lumen and increases local arterial stiffness. There is substantial evidence on the localisation of these plaque deposits at sites with haemodynamic conditions commonly characterised by low wall shear stress [1,2]. Aneurysms, on the other hand, involve the degradation of local arterial wall tissues resulting in lowering of local arterial stiffness and enlargement of local vessel cross-sections. If, in extreme cases, the wall stresses induced due to the transient fluid–structure interactions exceed the strength limit of the dilated artery wall, they cause vessel rupture leading to death from internal haemorrhage which has been reported to be in more than 80% of the cases [3]. Aneurysms are common in locations with secondary flow and flow recirculation even in normal resting conditions [4–6]. Considering the correlations between the various haemodynamic conditions and the onset and progression of different cardiovascular diseases, there is worldwide consensus on the need for enhancements in the current understanding of cardiovascular mechanopathobiology [7]. Since blood flow and pressure are dependent on a number of factors in addition to physical phenomena, experimental studies demand rigorous screening, which make them expensive and time consuming. A preferred mode of investigation is by computer simulation using mathematical models.
∗
Corresponding author. Tel.: +64 99219777; fax: +64 99219973. E-mail address: [email protected] (A. Al-Jumaily).
0895-7177/$ – see front matter © 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.mcm.2012.06.023
A. Al-Jumaily, A. Lowe / Mathematical and Computer Modelling 57 (2013) 460–468
461
Fig. 1. The axisymmetrical cylindrical 1-D domain; section S (x, t ) is defined at a location x at any time t with a cross-sectional area A(x, t ).
Quateroni et al. [8] and Taylor et al. [9] provide excellent reviews of current modelling techniques. Principles of conservation of mass, momentum and energy form the theoretical bases of all mathematical models that have been developed for the study of the physical aspects of blood flow. Van de Vosse and Stergiopulos [10] presents an excellent review which focuses on the physical and mathematical modelling of pulse wave propagation based on a general fluid dynamic approach. However, an analytical solution of the full form of these equations has not been developed and only solutions for special cases are available today. Proper application of numerical techniques for solving these equations demands the familiarity with assumptions and approximations that can be made to achieve a reasonable solution. This paper develops a mathematical model capable of representing pressure propagation in human arteries and uses this model to investigate the significance of different physical terms within this model on pressure propagation in the aorta. Although the development of 1-D mathematical models for studying pressure propagation in arteries and the effect of some phenomena have been previously discussed in the literature [11–19], to the best of our knowledge the application of a structured sensitivity analysis is presented for the first time in this work. Further, previous publications assume the aorta to be straight and either linearly or exponentially tapered which is not anatomically representative. For a better shape regeneration, bicubic splines are introduced for both wall thickness and lumen radius and this explains the higher accuracy of the model predictions of complex physiological parameters associated with arterial wave propagation. 2. Model development The pulse propagation in the arteries can be simulated as the propagation of pressure pulses in an incompressible, viscous or inviscid fluid (blood) flowing in a deformable tube (artery) [9,19]. To develop a reasonable simulation, the following assumptions are made based on clinical observations: the arterial pulses have a significantly longer wavelength, in the order of metres [20], when compared to the millimetre-scale arterial dimensions [21]. Since the aorta expands and contracts around 10% (±5% around the mean) with each beat [22], the amplitude of these pulses can be assumed to be very small. Shear rates in the proximal arteries are large enough to permit the assumption that the blood is Newtonian [23]. Velocity measurements in the proximal arteries have confirmed that the axial velocity is considerably larger than the radial velocity and thus 1-D equations can be used. Further, the artery wall is assumed thin, homogenous, elastic and impermeable, whereby a zero-slip condition may be applied to these equations. Under these assumptions, applying the one dimensional conservation of mass and momentum equations leads to the following governing equations of blood flow,
∂A ∂Q + =0 ∂t ∂x 2 ∂Q ∂ Q A ∂p + α + + kR u = 0 ∂t ∂x A ρ ∂x
(1) (2)
where A (x, t ) is the cross-sectional area, Q (x, t ) is the volume flow rate across the section, u (x, t ) = Q (x, t )/A (x, t ) is the axial velocity of the blood averaged across the section, ρ is the density of blood which is considered as a constant in this investigation, α is the advective momentum correction coefficient, p (x, t ) is the mean pressure at the section and kR is the loss of momentum per unit length due to viscosity of blood (ν). For a Poiseuille flow profile, α = 1.33 and kR = 8π ν while for the clinically measured flow profile, a value of α = 1.1 and kR = 21π ν has been suggested [11,24–26]. α = 1 assumes a completely flat profile, but it is often used as it permits significant mathematical simplification [27] (see Fig. 1). From physical theories, it is known that as the wave pulses move through the arterial tree, they cause local distension of the arterial wall. A relationship linking this distension and the pressure at a cross-section can be followed assuming the aorta as a thin walled pressurised cylinder where the stresses induced remain constant through the thin arterial wall.
462
A. Al-Jumaily, A. Lowe / Mathematical and Computer Modelling 57 (2013) 460–468
Neglecting inertial and viscoelasticity effects [24], the area can be written as A = Ao + ϕ(pnet )
(3)
where Ao is the resting cross-sectional area, ϕ (pnet ) is the area increment caused by the increment in pressure, pnet and
ϕ(pnet ) =
2Ao r Ew tw
· 1pnet =
2Ao r Ew t w
· ∆(p − pext ) = Ao c · ∆(p − pext )
(4)
where r is the internal radius of the artery, pext is the extramural pressure, Ew is the modulus of elasticity of the artery wall, tw is the thickness of the artery wall and c is the compliance of the cross-section. Since the extramural pressure, acting on the artery remains reasonably unchanged as the aorta is not rigidly held in position, pext ≈ constant, and hence the expansion and contraction of the artery is in response to internal pressure change p. Since the primary interest in this study is in pressure propagation, it can be assumed that the cross-sectional area, the pressure and the velocity of the flow perturbates around a steady state value (e.g. A = A0 + Ap ). The subscript 0 defines a steady state operating-point (mean) value and the subscript p defines a small perturbation in the system variables. In vivo measurements have shown that the blood flow (and hence all other coupled variables) oscillates periodically around a mean value [28] with small pulses compared to the steady state values. Previous studies [29] have confirmed that the propagation of the arterial waves can be simulated using first order differential equations. Although this system is nonlinear in nature, perturbation theories are used to simplify the differential equations. This procedure has been widely used in the literature and it has proven to generate the standard wave equation which is an effective tool in understanding wave propagations in the arterial system. Adding nonlinear complexity may improve the mathematical formulation but definitely will not have significant additions to understanding those waves. To obtain the final governing equations, substitute Eqs. (1) and (2) into Eqs. (3) and (4), implement perturbation technique for each of the variables around a steady state value, neglect high order terms, Laplace transform and after some mathematical manipulation obtain A0 csP + (cu0 P + Ucp0 + U )
(1 + cp0 )sU + u0 (1 + cp0 )
dA0 dx
dU dx
+ cu0 A0
+
dP dx
+ A0 (1 + cp0 )
(1 + cp0 ) dP p
dx
+
kR A0
dU dx
=0
U =0
(5) (6)
where U and P are the Laplace transforms of the instantaneous velocity and pressure oscillations respectively. The spatial derivatives are approximated using forward finite differences and matrix partition techniques are applied to reduce the equations to a [A] · [B] = [C ] form as shown below kp (n, s)Pn − [uo An c ]Pn−1 + [(2An − An−1 )(1 + po c )]Un − [An (1 + po c )]Un−1 = 0 kv (n, s)Un − [uo (1 + po c )]Un−1 +
1 + po c
ρ
Pn −
1 + po c
ρ
(7)
P n −1 = 0
(8)
where s is the Laplace variable, h is the step size, the subscript n represents the variable at the nth step and the subscript n − 1 represents the variable at the n − 1th step and kv (n, s) = h(1 + po c )s + uo (1 + po c ) + kp (n, s) = An hcs + cuo (2An − An−1 ).
hkR An
(9) (10)
Eqs. (7) and (8) are sufficient to represent fluid flow through a linearly tapered element. As these equations are homogeneous, the results drawn from them are continuous in one-dimensional Gaussian space. Replacing s by iω in these equations transfers the solution from the s domain to the frequency domain, where i is the imaginary unit and ω is the angular frequency. It is worth mentioning here that each one of the elements in the finite difference formulation is very small compared to the overall tapered aorta and is assumed to be of a uniform shape. Thus, the aorta tapering is reduced to small step changes and is not needed to be included in the above equations when applied to each element. The system is closed by specifying a relationship between pressure and flow at any one location along the aorta. The equations are then solved by implementing a modified Gauss elimination (L-U decomposition) strategy to give transfer functions (TFs) between pressure and flow at different locations along the aorta. 3. Simulation specifications This section is aimed at accurately describing the different parameters used in the simulation process to ensure that the results of this paper are reproducible.
A. Al-Jumaily, A. Lowe / Mathematical and Computer Modelling 57 (2013) 460–468
463
Table 1 Geometric domain description. Source: Data adapted from Westerhof et al. [21]. Artery segment
Length (cm)
Internal radius (cm)
Wall thickness (cm)
Young’s modulus (×105 N/m2 )
Ascending aorta (1) Ascending aorta (2) Aortic arch (1) Aortic arch (2) Thoracic aorta (1) Thoracic aorta (2) Thoracic aorta (3) Abdominal aorta (1) Abdominal aorta (2) Abdominal aorta (3) Common iliac artery (1)
2 2 2 3.9 5.2 5.2 5.2 5.3 5.3 5.3 5.8
1.470 1.440 1.120 1.070 0.999 0.675 0.645 0.610 0.580 0.548 0.368
0.164 0.161 0.132 0.127 0.120 0.090 0.087 0.084 0.082 0.078 0.063
4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0
Fig. 2. Invasively measured ascending aortic pressure after digitisation.
3.1. Geometric domain description Table 1 shows geometric properties at discrete locations along the aorta [22]. Although the modulus of elasticity of the physiological system is not uniform, the data presented in this table gives an acceptable approximation for linear modelling purpose. These data were interpolated using cubic splines to generate the geometric domain of the model. 3.2. Frequency spectrum of interest Within the experimental context of this investigation, the pressure fluctuation at the ascending aorta was invasively measured during catheterisation and was extracted in a beat-to-beat manner. The recorded pressure was digitised at 100 Hz and the resultant waveform is shown in Fig. 2. This beat was analysed to check its frequency contents. It was observed that a significant portion (94.1%) of the wave energy was contained in frequencies below 15.6 Hz. Previously reported values indicate that the range of 0–15 Hz is physiologically relevant. Thus, the present data agrees well with previous readings. 3.3. Impedance boundary condition Although branching plays an important role in wave reflections, this paper focuses on the significance of each term for a single tube. Previous simulations for analysing blood pressure and flow use arbitrary boundary conditions. For example, flow boundary conditions [12,17,18], pressure boundary conditions [30,31], modified forms of Windkessels as boundary conditions [13] and use of the characteristic impedance boundary conditions [15,32] have been used. The use of the impedance boundary condition has the advantage that the model results are generalised; however, the use of pressure or flow specifications are better suited for patient-specific models. Olufsen et al. have proposed the use of the fractal approach to formulate a ‘‘structured tree’’ boundary condition [17,18,33]. However, the assumed similarity laws in the derivation of this boundary condition are difficult to be physiologically realised. In light of the above-mentioned facts, invasively measured impedance is used as a boundary condition as it gives the best representation of the downstream vasculature. At the end of the aorta, the invasively measured impedance can be written as Z = f (ω).
(11)
Mills et al. measured the impedance at the iliac artery up to 13 Hz [34]. Based on the Windkessel theory, the impedance at high frequencies approaches a constant value and in this research, the impedance at frequencies at and above 13 Hz is assumed to remain constant. Without mention of the precise location of Mills’ measurement, it is further assumed that the measurement was made at the start of the iliac artery. The impedance at the end of the aorta can be calculated by applying the parallel impedance analogy from electrical circuit analysis. Each one of the finite difference elements needs a starting condition and end condition. The last element in the simulation uses the values of the impedance given above.
464
A. Al-Jumaily, A. Lowe / Mathematical and Computer Modelling 57 (2013) 460–468
Fig. 3. Map showing the simulated pressure waveforms (mm Hg) at different locations along the aorta when the model is used to simulate the pressure propagation in the geometric domain specified in Table 1. A is an arc tracing the timing of peak pressure, which may be used as a reference to calculate PWV. B is a region of low pressure along the aorta and C is a region with post-systole pressure fluctuation.
4. Model validation The invasively measured pressure waveform shown in Fig. 2 is broken down into its harmonics and represented as a Fourier series by, yx=0 = Ao +
n [Ai cos ωi t − Bi sin ωi t ].
(12)
i=1
Previous research [32] has shown that the first 5 harmonics contain over 90% of the energy of the pressure waveform. The pressure wave is assumed to contain 10 harmonics so as to contain more than 99% of the energy. Using the pressure transfer functions (TFx ) for different locations along the aorta and input Fourier series, the blood pressure at different locations can be calculated as p(x, t ) = A0 +
∞ [|TF|x An sin(nωt + ]TFx ) + Bn cos(nωt + ]TFx )]
(13)
n=1
where p (x, t ) is the pressure wave at any distance x from the aortic root, An and Bn are the coefficients of the nth harmonic of the pressure and TFx is the ratio of pressure at distance x from the aortic root to that at the aortic root, |TF| is the derived normalised amplitude ratio and |TF| is the associated phase lag. Fig. 3 shows a temporal and spatial map of the pressure predicted along the aorta. To validate the model, comparisons are made between the model outputs and available clinical data. The basic parameters for comparison are the pulse wave velocity (PWV), blood pressure trends and the effective reflecting distance (ERD). The model-predicted PWVs at different locations along the aorta are within the physiologically measured range of 4.6–6.5 m/s [35] and they progressively increase as the wave travels to the periphery (refer to Fig. 3, arc A) in agreement with clinical observations. From the pressure waveforms it is observed that the wave steepens as it propagates away from the heart and that a region of minimum pressure and consequently minimum wall shear stress (for a cylindrical thin-wall vessel, the wall shear stress is a function of the pressure) develops in the thoracic aorta (refer to Fig. 3, ellipse B). Careful simulation has shown that increasing the aortic stiffness decreases the minimum pressure and expands the region of low pressure. Previous research has shown that such sites, with low wall shear stress, are prone to plaque deposition [1,2]. Consequently, this demonstrates that the risk for the development of atherosclerosis considerably increases with increasing aortic stiffness (which normally takes place with ageing). A considerably different pressure pattern is noticed in the abdominal aorta where a distinct diastolic oscillation is observed (refer Fig. 3, ellipse C) which becomes more pronounced with the increase of aortic stiffness. This diastolic oscillation highlights the presence of the reflected wave coming from the periphery that is coupled to a negative flow wave. This negative flow wave slows down local fluid particles and may cause instantaneous flow reversal. Locations demonstrating these characteristics are commonly associated with vascular remodelling leading to the development of aneurysms, which are more commonly reported in the abdominal aorta [4–6]. At this stage, it may be difficult to identify in-plane locations of plaque deposition and aneurysm formation using 1-D models; however, it is possible to link the axial locations of arterial remodelling with the predicted pressure fields.
A. Al-Jumaily, A. Lowe / Mathematical and Computer Modelling 57 (2013) 460–468
465
Fig. 4. Effect of aortic radius and wall thickness tapering on the pressure propagation characteristics.
ERD is calculated based on the reflected wave transit time (RWTT) principle using the formula ERD = PWV × RWTT/2 which was proposed by Murgo et al. [36]. Using the available numerical data for the aorta [20], the model predicts ERD to be 47.5 cm from the heart while previous mathematical models predict ERD to be 39.8 cm away from the heart [32]. The present result correlates very well with previously reported values based on invasive data by Mills et al. [34] who estimated ERD to be 50.9 cm away from the heart. RWTT analysis carried out on the Framingham Cohort showed a mean ERD of 50 cm in women and 59 cm in men [37], to which the present model prediction compares favourably. The lack of precise data required to reconstruct the arterial model from the Framingham Cohort restricts the investigation to this level; however it is believed that if patient-specific data is available, the model can be used to predict the precise location of ERD. Therefore, the model prediction matches major clinical observations and hence it is able to simulate the propagation of pressure waves in the human aorta with a high degree of accuracy. 5. Pressure propagation under different conditions Considering each condition individually, the significance of the different physical phenomena on the pressure propagation in the aorta is studied by investigating the change in the pressure TF between the aortic root and the iliac bifurcation. 5.1. Geometric tapering Fig. 4 shows the effect of tapering of the artery internal radius and wall thickness on the TF for three different cases: (i) in a straight tube of mean internal radius and wall thickness, 1.009 cm and 0.121 cm respectively, (ii) a tube of constant internal radius with linearly tapering wall thickness and (iii) a tube with linearly tapering internal radius and wall thickness. The degrees of tapering in internal radius and wall thickness respectively were calculated based on the physiological data given in Table 1. Fig. 4 shows that tapering in the internal radius causes the frequency spectrum of the TF to shift towards higher frequencies (indicating an increase of resonant or natural frequency of the aorta). There is also an amplification of the pulse amplitude and a reduction in the phase delay, indicating an increase in PWV. However, tapering in the wall thickness has the effect of reducing the pulse amplitude, increasing the phase delay indicating a decrease in PWV and increasing the natural frequency of the system. The difference between the cases is observable for frequencies above 4 Hz, which are within the physiological range. These changes may be explained by the fact that tapering changes the compliance of the system following Eq. (4). It is also observed that tapering in the internal radius influences the pressure propagation characteristics of the aorta to a much higher degree than the tapering of wall thickness. This may be due to the fact that the internal radius tapers down at 1.2758° while the wall thickness tapers down at 0.1190° over the length of the aorta. 5.2. Aortic stiffness The aortic stiffness under the normal condition is assumed as 4 × 105 N/m2 [20]. Fig. 5 shows the change in TF caused by doubling the aortic stiffness. It is seen that increasing the aortic stiffness shifts the natural frequency to higher values which is in agreement to the quarter wavelength formula [38] while a reduced phase delay suggests the increase of PWV. These changes are observed because the system compliance is dependent on the aortic stiffness according to Eq. (4).
466
A. Al-Jumaily, A. Lowe / Mathematical and Computer Modelling 57 (2013) 460–468
Fig. 5. Effect of change of arterial stiffness.
Fig. 6. Effect of change of the momentum correction factor, α .
5.3. Changing the assumed flow profile The two model parameters used to specify the flow pattern in the aorta are the advective momentum correction factor
α and the momentum loss factor kR from Eq. (2).
Fig. 6 shows the cases when the advective acceleration term is included or disregarded in the simulation. It is indicated that the advective acceleration term itself is very small compared to the other terms in Eq. (2) and does not affect the pressure propagation in the aorta. Minor changes in the momentum correction factor, i.e. Poiseuille to the experimentally measured flow profile with 1α < 1, should therefore have smaller effects on the propagation characteristics of the aorta than complete omission, which is studied here. The viscosity of blood is 3.2 × 10−3 Pa-s [39] and is incorporated into the simulation through the term kR . Fig. 7 shows the change in the propagation characteristics of the aorta by omitting kR . The graph shows a small change in the amplitude of the TF peaks with no change in the overall trend. There is no observable change in the propagation velocity of the pulses confirming that the viscosity of the blood insignificantly influences pressure propagation in the elastic arteries. Therefore, these observations indicate that the changes in α and kR have an insignificant effect on the TF and hence the assumed flow profiles should not change the pressure propagation characteristics. Future studies can assume Poiseuillian flow as it leads to considerable mathematical simplifications. 6. Conclusions In this paper, blood flow in the aorta is simulated using a 1-D fluid flow model to investigate the effect of various physical terms on the pulse wave propagation characteristics. Clinical parameters predicted by the simulation are close to invasive
A. Al-Jumaily, A. Lowe / Mathematical and Computer Modelling 57 (2013) 460–468
467
Fig. 7. Effect of change of the viscous correction term, kR .
measurements. The model has shown that the changes in the geometry and material properties strongly influence the pressure propagation characteristics while the flow profile has an insignificant effect. Acknowledgement We would like to acknowledge IBTec former employee Ashis Mookerjee for his computational contribution to this article. References [1] C.G. Caro, J. Fritz-Gerald, R. Schroter, Atheroma and arterial wall shear: observation, correlation and proposal for a shear dependent mass transfer mechanism for atherogenesis, Proceedings of the Royal Society of London B. Biological Science 177 (1971) 109–159. [2] C.A. Taylor, C.P. Cheng, L.A. Espinosa, B.T. Tang, D. Parker, R.J. Herfkens, In vivo quantification of blood flow and wall shear stress in the human abdominal aorta during lower limb exercise, Annals of Biomedical Engineering 30 (3) (2002) 402–408. [3] D. Katz, J. Cronenwett, The cost effectiveness of early surgery verses watchful waiting in the management of small abdominal aortic aneurysms, Journal of Vascular Surgery 19 (1994) 980–990. [4] C.M. Scotti, A.D. Shkolnik, S.C. Muluk, E.A. Finol, Fluid–structure interaction in abdominal aortic aneurysms: effects of asymmetry and wall thickness, Biomedical Engineering Online 4 (1) (2005) 64. [5] J.C. Lasheras, The biomechanics of arterial aneurysms, Annual Review of Fluid Mechanics 39 (2007) 293–319. [6] R. Peattie, T. Riehle, E. Bluth, Pulsatile flow in fusiform models of abdominal aortic aneurysms: flow fields, velocity patterns and flow-induced wall stresses, Journal of Biomechanical Engineering 126 (4) (2004) 438–446. [7] M. O’Rourke, J. Hashimoto, Mechanical factors in arterial ageing-a clinical perspective, Journal of the American College of Cardiology 50 (1) (2007) 1–13. [8] A. Quarteroni, M. Tuveri, A. Veneziani, Computational vascular fluid dynamics: problems, models and methods, Computing and Visualization in Science 2 (2000) 163–197. [9] C.A. Taylor, M. Draney, Experimental and computational methods in cardiovascular fluid dynamics, Annual Review of Fluid Mechanics 36 (2004) 197–231. [10] F.N. van de Vosse, N. Stergiopulos, Pulse wave propagation in the arterial tree, Annual Review of Fluid Mechanics 43 (2011) 467–499. [11] P. Reuderink, H. Hoogstraten, P. Sipkema, B. Hillen, N. Westerhof, Linear and nonlinear one-dimensional models of pulse wave transmission at high Wormsley numbers, Journal of Biomechanics 22 (8–9) (1989) 819–827. [12] S.J. Sherwin, L. Formaggia, J. Peiro, V. Franke, Computational modelling of 1D blood flow with variable mechanical properties and its application to the simulation of wave propagation in the human arterial system, International Journal for Numerical Methods in Fluids 43 (2001) 673–700. [13] P. Nardinocchi, G. Pontrelli, L. Teresi, A one-dimensional model for blood flow in prestressed vessels, European Journal of Mechanics A/Solids 24 (2005) 23–33. [14] J.R. LaCourse, G. Monahakrishnan, K. Sivaprasad, Simulations of arterial pressure pulses using a transmission line model, Journal of Biomechanics 19 (9) (1986) 771–780. [15] J. Wang, K. Parker, Wave propagation in a model of the arterial circulation, Journal of Biomechanics 37 (2004) 457–470. [16] K. Matthys, J. Alastruey, J. Peiro, A. Khir, P. Segers, P. Verdonck, et al., Pulse wave propagation in a model human arterial network: assessment of 1-D numerical simulations against in vitro measurements, Journal of Biomechanics 40 (15) (2007) 3476–3486. [17] M.S. Olufsen, Structured tree outflow conditions for blood flow in large systemic arteries, American Journal of Physiology 276 (1999) H257–H268. [18] M.S. Olufsen, C. Peskin, W. Kim, E. Pederson, A. Nadim, J. Larsen, Numerical simulation and experimental validation of blood flow in arteries with structured tree outflow conditions, Annals of Biomedical Engineering 28 (11) (2000) 1281–1299. [19] J.B. Grotberg, O.E. Jensen, Biofluid mechanics in flexible tubes, Annual Review of Fluid Mechanics 36 (2004) 121–147. [20] P. Segers, S. Carlier, A. Pasquet, S.I. Rabben, L.R. Hellevik, E. Remme, et al., Individualising the aorto-radial pressure transfer: feasibility of a model based approach, American Journal of Physiology 279 (2000) H542–H549. [21] N. Westerhof, F. Bosman, C. de Vries, A. Noodergraaf, Analog studies of the human systemic arterial tree, Journal of Biomechanics 2 (1969) 121–143. [22] P. Boutouyrie, S. Laurent, A. Benetos, X.J. Girerd, A.P. Hoeks, M. Safar, Opposing effects of ageing on distal and proximal large arteries in hypertensives, Journal of Hypertension Supplements 10 (1992) S87–S91. [23] F.J.H. Gijsen, F.N. van de Vosse, J.D. Jassen, The influence of the non-Newtonian properties of blood on the flow in large arteries: steady flow in a carotid bifurcation model, Journal of Biomechanics 32 (1999) 601–608. [24] S. Balar, T. Rogge, D. Young, Computer simulation of blood flow in the human arm, Journal of Biomechanics 22 (6–7) (1989) 691–697. [25] L. Formaggia, D. Lamponi, A. Quarteroni, One-dimensional models for blood flow in arteries, Journal of Engineering Mathematics 47 (2003) 251–276.
468
A. Al-Jumaily, A. Lowe / Mathematical and Computer Modelling 57 (2013) 460–468
[26] E.M. Pedersen, W.Y. Kim, N.H. Staalsen, J.M. Hasenkam, H. Nygaard, P.K. Paulsen, Development of velocity profiles and retrograde flow in the porcine abdominal aorta under different haemodynamic conditions, Scandinavian Cardiovascular Journal 33 (1999) 206–214. [27] L. Formaggia, J. Gerbeau, F. Nobile, A. Quarteroni, On the coupling of 3D and 1D Navier–Stokes equations for flow problems in compliant vessels, Computer Methods in Applied Mechanics and Engineering 191 (2001) 561–582. [28] J. Moore, D. Ku, Pulsatile velocity measurements in a model of the human abdominal aorta under simulated exercise and postprandial conditions, Journal of Biomechanical Engineering 116 (1) (1994) 107–111. [29] J. Womersley, The mathematical analysis of the arterial circulation in a state of oscillatory motion, Wright Air Development Centre, 1957, Report No. WADC-TR56-614. [30] T. Mabotuwana, L. Cheng, A. Pullan, A model of blood flow in the mesenteric arterial system, Biomedical Engineering Online 6 (1) (2007) 17. [31] M.S. Olfusen, A. Nadim, L.A. Lipsitz, Dynamics of cerebral blood flow regulation explained using a lumped parameter model, American Journal of Physiology 282 (2002) R611–R622. [32] A.P. Avolio, Multi-branched model of the human arterial system, Medical and Biological Engineering and Computing 18 (6) (1980) 709–718. [33] B. Steele, M.S. Olufsen, C. Taylor, Fractal network model for simulating abdominal and lower extremity blood flow during resting and exercise conditions, Computer Methods and Programs in Biomedicine 10 (1) (2007) 39–51. [34] C.J. Mills, I.T. Gabe, J.H. Gault, D.T. Mason, J. Ross, E. Braunwald, et al., Pressure-flow relationships and vascular impedance in man, Cardiovascular Research 4 (1970) 405–417. [35] A. Avolio, S. Chen, R. Wang, C. Zhang, M. Li, M. O’Rourke, Effects of ageing on changing arterial compliance and left ventricular load in a northern Chinese urban community, Circulation 68 (1983) 50–58. [36] J.P. Murgo, N. Westerhof, J.P. Giolma, S.A. Altobelli, Aortic input impedance in normal man: relationship to pressure wave forms, Circulation 62 (1980) 105–116. [37] G.F. Mitchell, H. Parise, E.J. Benjamin, M.G. Larson, M.J. Keyes, J.A. Vita, et al., Changes in arterial stiffness and wave reflection with advancing age in healthy men and women: the framingham heart study, Hypertension 43 (2004) 1239–1245. [38] P. Segers, P. Verdonck, Role of tapering in aortic wave reflection: hydraulic and mathematical model study, Journal of Biomechanics 33 (2000) 299–306. [39] A.R. Pries, D. Neuhaus, P. Gaehtgens, Blood viscosity in tube flow: dependance on diameter and hematocrit, American Journal of Physiology 263 (1992) H1770–H1778.