1 Introduction Pressure and flow pulse waveforms carry information on the functionality of the cardiovascular system and the morphology of the arterial network, which can be valuable for the diagnosis and treatment of disease. Modelling allows us to study the effect on pulse waveforms of the physical properties of the system, such as arterial geometry and distensibility, cardiac output, heart rate and peripheral impedance to flow, and analyse the mechanisms underlying clinically relevant changes (Stergiopulos et al., 1992; Wang and Parker, 2004; Mynard and Nithiarasu, 2008; Liang et al., 2009). In particular, it allows us to test clinical diagnostic techniques that are based on pulse wave analysis and disentangle their underlying mechanisms (Karamanoglu et al., 1994; Trachet et al., 2010; Alastruey et al., 2006; Alastruey, 2011). These studies can be extremely challenging in vivo for technical and physiological reasons; e.g. some vessels are inaccessible, manipulation of the properties of interest can be dangerous or can elicit reflex compensation, and several parameters of interest are not directly measurable. We have previously shown the ability of the nonlinear 1-D equations of pulse wave propagation in elastic vessels to capture the main features of pressure and flow waveforms measured in well-defined experimental 1:1 replicas of the larger conduit arteries in the human systemic circulation (Segers et al., 1998; Segers and Verdonck, 2000; Alastruey, 2006; Matthys et al., 2007). Matthys et al. (2007) reported average relative root-mean-square errors between numerical and experimental waveforms smaller than 3.5% for pressure and 19% for the flow rate at 70 locations in the tapered silicone network sketched in Fig. 1. Much of these errors arose from relatively high-frequency oscillations in the peripheral vessels predicted in the numerical model but not seen in the experimental measurements. Wall viscosity was suggested to play an important role in damping oscillations in the numerical pulse waveforms, but this hypothesis could not be tested since a purely elastic tube law was used in the 1-D model. Recent works have shown the benefits of modelling the arterial wall using a visco-elastic law. Bessems et al. (2008) tested a visco-elastic 1-D formulation against in vitro data measured in single tapered polyurethane vessels, showing that wall visco-elasticity is necessary to accurately predict the propagation and attenuation of pressure and flow waves. Reymond et al. (2009) compared the predictions of a visco-elastic 1-D distributed model of the human arterial tree with average pressure and flow waveforms measured at several arterial locations in a group of young subjects. Qualitative comparison of these waveforms and quantitative comparison of systolic, diastolic and mean pressure and flow indicated that the visco-elastic effects may be significant, especially in peripheral branches. They argued that energy losses and damping effects due to wall visco-elasticity are of the same order of magnitude as wall friction in large and medium size vessels. However, the parameters of their model were not specific to each subject. Several works (Armentano et al., 1995a,b; Craiem et al., 2005; Čanić et al., 2006) have used Voigt-type visco-elastic models to simulate the visco-elastic behaviour of human and animal arteries, and improve the agreement between measured and simulated pulse waveforms. These models feature hysteresis (different stress–strain relationship for loading and unloading) and creep (continuous extension at constant load). The aim of this work is to study the role of visco-elasticity in explaining the discrepancies of the results presented in Matthys et al. (2007). We will use a nonlinear, time-domain 1-D formulation that simulates the arterial wall as a Voigt-type visco-elastic material. Except for the shape of the velocity profile, the parameters required by the numerical algorithm will be directly measured in the in vitro setup and no data fitting will be involved. We will, therefore, provide an assessment of nonlinear visco-elastic 1-D modelling by comparison against in vitro measurements in a well-defined arterial network with tapered vessels. The notation and SI units of the quantities used in this paper are given in Table 1. 2 Methods 2.1 Visco-elastic 1-D formulation Conservation of mass and momentum applied to a 1-D impermeable and deformable tubular control volume of incompressible and Newtonian fluid flowing with a constant axisymmetric velocity profile yields (Sherwin et al., 2003) (1) ∂ A ∂ t + ∂ ( AU ) ∂ x = 0 , ∂ U ∂ t + U ∂ U ∂ x + 1 ρ ∂ P ∂ x = f ρ A , where x is the axial coordinate along the vessel, t the time, A(x, t) A ( x , t ) the cross-sectional area of the lumen, U(x, t) U ( x , t ) the average axial velocity, P(x, t) P ( x , t ) the average internal pressure over the cross-section, and f = − 22μ π U f = − 22 μ π U (Alastruey, 2006) the friction force per unit length. The constant density and viscosity of the 65–35% water–glycerol mixture used to mimic blood were ρ = 1050 kg m − 3 ρ = 1050 kg m − 3 and μ = 2.5 mPa s μ = 2.5 mPa s , respectively. We have modelled the silicone used to make the experimental arteries as a Voigt-type material; i.e. its tensile stresses τ and strains ϵ are related through (Fung, 1993) (2) τ = E ϵ + φ d ϵ d t , where E is the Young's modulus and φ the viscosity of silicone. Assuming the in vitro arterial walls to be thin, isotropic, homogeneous and incompressible, and to deform axisymmetrically, each circular cross-section independently of the others, τ can be related to the internal pressure P and ϵ to the change of radius of the vessel wall η = r − r 0 η = r − r 0 using Laplace's law (Olufsen et al., 2000) (3) τ = r ( P − P ext ) h , ϵ = 1 1 − σ 2 η r 0 , where r = A / π , h(x) is the wall thickness, σ = 0.5 σ = 0.5 the Poisson's ratio and r 0(x) r 0 ( x ) the luminal radius at P = P ext P = P ext ; P ext being the constant external pressure (Matthys et al., 2007). Combining Eqs. (2) and (3) and assuming that 1 / r 1 / r can be approximated by 1 / r 0 1 / r 0 yields (4) P = P ext + 4 Eh 3 ( r 0 ) 2 η + 4 φ h 3 ( r 0 ) 2 ∂ η ∂ t . Expressing η as a function of A = π r 2 A = π r 2 and A 0 = π(r 0)2 A 0 = π ( r 0 ) 2 , yields the visco-elastic tube law (5) P = P ext + β A 0 ( A − A 0 ) + Γ A 0 A ∂ A ∂ t , β = 4 3 π Eh , Γ = 2 3 π φ h . The parameters β(x) β ( x ) and Γ(x) Γ ( x ) are related to the elastic and visco-elastic properties of the arterial wall, respectively. It is important to remark that Eq. (5) can also be derived from a generalised string model (Quarteroni et al., 2000; Formaggia et al., 2003). 2.2 Numerical solution Eqs. (1) and (5) were solved in the 37 arterial segments of the hydraulic model in Fig. 1 using a discontinuous Galerkin scheme with a spectral/hp spatial discretisation. This is a convenient scheme for high-order discretisation of convection-dominated flows (Cockburn and Shu, 1998), such as arterial flows. It allows us to propagate waves of different frequencies without suffering from excessive dispersion and diffusion errors. Eqs. (1) and (5) can be written in the following conservative form: ∂ U ∂ t + ∂ F ∂ x = S U , with U = A U , S U = 0 f ρ A , F = F e + F v = AU U 2 2 + P ext + P e ρ + 0 − Γ ρ A 0 A ∂ ( AU ) ∂ x , and P e = ( β / A 0 ) ( A − A 0 ) being the elastic component of pressure. The flux F was separated into an elastic (F e ) ( F e ) and a viscous (F v ) ( F v ) term, and the mass conservation ∂ A / ∂ t = − ∂ (AU) / ∂ x ∂ A / ∂ t = − ∂ ( AU ) / ∂ x was applied to change the time derivative to a spatial derivative in the viscous term of Eq. (5). The discrete form of this conservative law in a domain Ω discretised into a mesh of N el elemental non-overlapping regions Ω e = [ x e l , x e u ] , such that x e u = x e + 1 l for e = 1, …, N el − 1 e = 1 , … , N el − 1 , and ⋃ e = 1 N el Ω e = Ω is given by (Karniadakis and Sherwin, 2005; Alastruey, 2006) (6) ∑ e = 1 N el ∂ U δ ∂ t , ψ δ Ω e + ∂ F ( U δ ) ∂ x , ψ δ Ω e + [ ψ δ · { F u − F ( U δ ) } ] x e l x e u = ∑ e = 1 N el ( S U δ , ψ δ ) Ω e , for all ψ δ ψ δ in V δ V δ , where (u, v) Ω = ∫ Ω uv d x ( u , v ) Ω = ∫ Ω uv d x is the standard L 2(Ω) L 2 ( Ω ) inner product, U δ U δ and ψ δ ψ δ denote the approximation of U and test functions ψ , respectively, in the finite space V δ V δ of piecewise polynomial vector functions (they may be discontinuous across inter-element boundaries), and F u = F e u + F v u is the approximation of the flux at the interface. The term F e u was treated through the solution of a Riemann problem as described by Alastruey (2006). The term F v u requires a different treatment. Various ways of dealing with this term were analysed by Zienkiewicz et al. (2003) and the references therein. Here, F v u at the inter-element boundaries was approximated as F v u | x e u = F v u | x e + 1 l = 1 2 ( F v | x e u + F v | x e + 1 l ) , e = 1 , … , N el − 1 , with F v u | x 1 l = F v | x 1 l at the inlet of the domain and F v u | x N el u = F v | x N el u at the outlet, so that F v u − F v ( U δ ) = 0 at both boundaries. The expansion bases were selected to be a polynomial space of order 𝒫 and the solution was expanded on each region Ω e Ω e in terms of Legendre polynomials L p (ξ) L p ( ξ ) ; i.e. (7) U δ | Ω e ( x e ( ξ ) , t ) = ∑ p = 0 P L p ( ξ ) U ^ e p , where U ^ e p ( t ) are the expansion coefficients. Legendre polynomials are particularly convenient because the basis is orthogonal with respect to the L 2(Ω e ) L 2 ( Ω e ) inner product. Following standard finite element techniques, the following elemental affine mapping was introduced: x e ( ξ ) = x e l ( 1 − ξ ) 2 + x e u ( 1 + ξ ) 2 , with ξ in the reference element Ω st = { − 1 ≤ ξ ≤ 1} Ω st = { − 1 ≤ ξ ≤ 1 } . The choice of discontinuous discrete solution and test functions allows us to decouple the problem on each element, the only link coming through the boundary fluxes. Substitution of (7) into (6) and letting ψ δ | Ω e = U δ | Ω e ψ δ | Ω e = U δ | Ω e , yields 2𝒫 2 P equations to be solved for each Ω e Ω e , e = 1, …, N el e = 1 , … , N el , (8) d U ^ i , e p d t = F ( U δ | Ω e ) , p = 0 , … , P , i = 1 , 2 , where U ^ i , e p , i = 1 , 2 are the two components of U ^ e p , F ( U δ | Ω e ) = − ∂ F i ∂ x , L p Ω e − 1 J e [ L p [ F i u − F i ( U δ | Ω e ) ] ] x e l x e u + ( S i ( U e δ ) , L p ) Ω e , and J e = 1 2 ( x e u − x e l ) is the Jacobian of the elemental mapping from Ω st Ω st . For every Ω e Ω e , U δ | Ω e U δ | Ω e was evaluated at Gauss–Lobatto–Legendre quadrature points of order 𝒬 to evaluate the integrals (∂ F i / ∂ x, L p ) Ω e ( ∂ F i / ∂ x , L p ) Ω e and (S i , L p ) Ω e ( S i , L p ) Ω e . All spatial derivatives were calculated using collocation differentiation at the quadrature points (Karniadakis and Sherwin, 2005). Arterial segments were divided into non-overlapping elements with a 2 cm length, when possible, and a polynomial and quadrature order of 3. Elements or segments shorter than 1.5 cm were given a polynomial and quadrature order of 2. Zero pressures and velocities were considered as initial conditions. The discretisation in time was performed by a second-order Adams–Bashforth scheme. The time step was Δ t = 100 μ s Δ t = 100 μ s for the elastic case and Δ t = 20 μ s Δ t = 20 μ s for the visco-elastic case, since the convection Δ t Δ t scales like the square of the polynomial order 𝒫, whereas the diffusion Δ t Δ t scales like twice the square of 𝒫 (Karniadakis and Sherwin, 2005). Although Δ t Δ t had to be reduced in the visco-elastic model, it is still practical to use an explicit scheme because of the small 𝒫 used and the fact that wave speeds are much larger than flow velocities. Visco-elasticity was neglected at the boundary conditions of the network and the junctions, which were implemented as described in Matthys et al. (2007, Section 2.2.3) for both the purely elastic and visco-elastic models. 2.3 Numerical parameters from in vitro data Table 2 shows the geometry, elasticity and boundary conditions required by the elastic and visco-elastic formulations, which were measured in the silicone network without any parameter fitting, as we detailed in Matthys et al. (2007). Most of the arterial segments were linearly tapered and their wall thickness decreased toward distal locations. Their Young's modulus was approximately constant (1.2 MPa), since they were made of the same silicone material. We assumed P ext = 0 P ext = 0 since the experimental setup was not submerged in water. Here we estimated φ (which was not measured in Matthys et al., 2007) from an uniaxial extension test in a sample of the same silicone material used in the experimental setup. Fig. 2 (top) shows the extension ψ and measured load F used in the test, which were normalised by their corresponding maximum values F max and ψ max ψ max . The test was carried out using a tensometer and 1 kN load cell (Instron 5542, High Wycombe, UK), loading and unloading at a rate of 20 mm min − 1 min − 1 to a maximum extension of 1 mm. Three consecutive loading and unloading cycles were followed by 30 s at zero extension, loading to maximum extension which was held for 60 s before the final unloading to zero extension. A cubic spline was fitted to F(t) and ψ(t) ψ ( t ) , separating loading and unloading. Stresses τ(t) τ ( t ) were calculated as the ratio of F and the cross-sectional area of the sample, and strains ϵ(t) ϵ ( t ) as the ratio of ψ(t) ψ ( t ) and the initial sample length. Fig. 2 (bottom) shows the F ψ-loop F ψ -loop for the first loading–unloading cycle (only the first cycle is shown for clarity; the following cycles were similar with no evidence of conditioning). We calculated φ as (9) φ = ∫ ϵ ini ϵ end τ d ϵ ∫ t ini t end d ϵ d t 2 d t , which follows from integrating both sides of Eq. (2) over one load–unload cycle and using ∫ ϵ ini ϵ end E ϵ d ϵ = 0 . The limits of the integrals are the initial (t ini , ϵ ini ϵ ini ) and final (t end , ϵ end ϵ end ) time and strain in a load–unload cycle, with ϵ ini = ϵ end ϵ ini = ϵ end . We obtained φ = 3.0±0.3 kPa s φ = 3.0 ± 0.3 kPa s . 2.4 Error calculations To assess the accuracy of the numerical predictions of pressure P and volume flow rate Q=AU, the following relative errors were calculated at each measuring site, E P = 1 N ∑ i = 1 N P i exp − P i P i exp 2 , E Q = 1 N ∑ i = 1 N ( Q i exp − Q i ) 2 max ( Q exp ) , (10) E P ¯ = ( mean ( P exp ) − mean ( P ) ) 2 mean ( P exp ) , E Q ¯ = ( mean ( Q exp ) − mean ( Q ) ) 2 mean ( Q exp ) , where N is the total number of samples in one cardiac cycle, the subscript i indicates the sampling point (the sampling rate was 1 kHz), the superscript exp indicates experimental data, the operator max(·) max ( · ) calculates the maximum value within i = 1, …, N i = 1 , … , N , and mean(·) mean ( · ) calculates the mean value within i = 1, …, N i = 1 , … , N . In the frequency domain, we calculated (11) E P ^ = 1 M ∑ j = 1 M ( P ^ j exp − P ^ j ) 2 P ^ 1 exp , E Q ^ = 1 M ∑ j = 1 M ( Q ^ j exp − Q ^ j ) 2 Q ^ 1 exp , where M is the total number of harmonics and P ^ j and Q ^ j are the amplitudes of the pressure and flow jth harmonic, respectively. 3 Results Fig. 3 compares the ensemble averaged experimental waveforms at four locations with the numerical predictions obtained using both the elastic and the visco-elastic 1-D formulations. As in the purely elastic study, we matched the onset of the experimental and numerical systolic ejections. These comparisons are representative of the results obtained at other measurement sites and show the ability of visco-elasticity to reduce the high-frequency oscillations produced by the elastic formulation, especially in peripheral vessels. Wall visco-elasticity reduces relative errors for pressure and the flow rate in the time domain for all the 70 measurement sites. Table 3 shows these errors for both the (original) purely elastic and the (new) visco-elastic formulations in the midpoint of the arterial segments, and Table 4 shows average relative errors over all 70 measuring sites, and over all the sites in the aorta and the vessels of the first, second and third generation of bifurcations. For both formulations, relative errors are larger for the flow rate than for pressure. If visco-elasticity is modelled, the averaged pressure ℰ P E P and flow ℰ Q E Q errors over all measuring sites are significantly decreased (from 3.0% to 2.5% and from 15.7% to 10.8%, respectively). Variations in φ of ±0.3 kPa s ± 0.3 kPa s (the standard deviation of the measured φ) leads to non-significant changes in the previous errors of less than 0.05% for pressure and 0.2% for the flow rate. Changes in the errors for mean pressures E P ¯ and flow rates E Q ¯ are not significant if visco-elasticity is modelled (Table 4). In the frequency domain, the visco-elastic formulation produces greater decreases in the relative errors of the amplitudes from the 5th to the 20th pressure E P ^ and flow E Q ^ harmonics than for the first five harmonics (Tables 3 and 4). Fig. 4 compares the amplitudes of the first 20 experimental flow harmonics in the midpoint of a vessel of the first generation of bifurcations, with the numerical predictions obtained using the elastic and the visco-elastic formulations. The latter captures these amplitudes better, especially with the increasing frequency. Similar results were obtained at the other locations studied. The A–P curves obtained using the visco-elastic model present an area of hysteresis that is much smaller than the area obtained from in vivo measurements (Armentano et al., 1995a,b; Valdez-Jasso et al., in press), as Fig. 5 shows in an aortic location. The numerical formulation predicts an increase in this area with the increasing wall viscosity φ (Fig. 5). 4 Discussion We have shown that simulation of wall visco-elasticity decreases the relative root-mean-square errors of numerical pressure and flow predictions in the aorta and all the generations of bifurcations of the silicone model, in both the frequency and time domains. According to our frequency–domain analysis, the damping effect due to visco-elasticity is more significant from the 5th flow harmonic (Table 4). This leads to a decrease in the high-frequency oscillations of the purely elastic model and supports the hypothesis suggested in Matthys et al. (2007) that visco-elasticity has a more significant effect on damping oscillations predicted by the purely elastic formulation (with higher amplitudes than those observed in vivo) than energy losses at junctions and peripheral compliances in the overflow reservoirs. The damping effect in the visco-elastic model is greater in distal vessels, which is in agreement with the increase in the viscous term of the tube law (5) with the decreasing cross-sectional area toward distal locations. The minor reductions in the relative errors of mean pressure and flow rate when visco-elasticity was modelled suggest that the effect of visco-elasticity on mean pressures and flow rates is secondary. The relative errors of flow rate predictions were larger than corresponding pressure predictions, which is consistent with the larger uncertainty of the experimental flow measurements compared with the experimental pressure measurements discussed in Matthys et al. (2007). Additional uncertainties were introduced by changes in the boundary conditions of the experimental setup every time the pressure and flow transducers were repositioned. Indeed, significantly different errors were obtained between measured and predicted mean flow rates E Q ¯ at several locations within the same vessel (e.g. Segments 31–33 in Table 3). The uncertainties in the measurement of the visco-elastic modulus were relatively small and did not introduce large errors in the numerical results. However, the uniaxial extension tests were carried out at a slower speed (about 6 s for a load–unload cycle) than the experimental inflow imposed by the pump (70 beats per minute). This difference may be responsible for some of the errors observed in the visco-elastic pressure and flow predictions. Another limitation of this study is the use of a linear model to describe the dynamics of the vessel walls (Eq. (2)) that does not account for the stress relaxation observed in the uniaxial tests (Fig. 2, top). Models accounting for stress relaxation (Bessems et al., 2008; Devault et al., 2008; Valdez-Jasso et al., 2009), and the nonlinear behaviour of the wall (Reymond et al., 2009; Valdez-Jasso et al., in press) have been proposed, but their applicability to patient-specific 1-D modelling is more challenging, since more parameters must be estimated than in the Voigt model. We believe that a linear constitutive model is a reasonable approach, since wall deformations are small under normal physiological conditions. According to the 1-D visco-elastic results, changes in the diameter of the ascending aorta over the cardiac cycle are less than 6.5% of its mean diameter. These changes decrease toward distal locations. Although the areas of hysteresis produced by the visco-elastic formulation (Fig. 5) are smaller than those observed in vivo, they are in agreement with the corresponding areas measured experimentally (Fig. 2, bottom). We could not assess the accuracy of the diameter waveforms predicted by the 1-D model, since we did not measure them experimentally. However, the area of hysteresis increases if a greater wall visco-elasticity φ = 24 kPa s φ = 24 kPa s is used. This is a more physiological value according to data in Armentano et al. (1995b) and Gariepy et al. (1993) for the carotid and femoral arteries of normotensive men. Using an anatomically realistic model of the arterial network, we have shown that wall viscosities about an order of magnitude smaller than those measured in the human play a significant role in shaping the high-frequency components of the pressure and flow waveforms, especially in the periphery. This result supports the qualitative observations by Reymond et al. (2009) in the human and our previous conclusion that energy losses at bifurcations have a secondary effect on the blood flow in large arteries compared to visco-elasticity (Matthys et al., 2007). We have provided further evidence to support the conclusions given in the previous works (Segers et al., 1998; Segers and Verdonck, 2000; Olufsen et al., 2000; Alastruey, 2006; Matthys et al., 2007; Steele et al., 2007; Bessems et al., 2008; Reymond et al., 2009) about the ability of the 1-D formulation to capture the main features of pressure and flow wave propagation in large arteries (and hence the use of this mathematical model to simulate clinically relevant problems), if accurate measurements of the parameters of the model are provided. The 1-D formulation is computationally inexpensive (it takes the order of minutes to solve one cardiac cycle for a whole body model) and, hence, offers a good balance between accuracy and computational cost when a global assessment of pressure wave propagation in the cardiovascular system is required. Conflict of interest statement There are no conflicts of interest between the authors of this paper and other external researchers or organisations that could have inappropriately influenced this work.
