Novel wave intensity analysis of arterial pulse wave propagation accounting for peripheral reflections

We present a novel analysis of arterial pulse wave propagation that combines traditional wave intensity analysis with identification of Windkessel pressures to account for the effect on the pressure waveform of peripheral wave reflections. Using haemodynamic data measured in vivo in the rabbit or generated numerically in models of human compliant vessels, we show that traditional wave intensity analysis identifies the timing, direction and magnitude of the predominant waves that shape aortic pressure and flow waveforms in systole, but fails to identify the effect of peripheral reflections. These reflections persist for several cardiac cycles and make up most of the pressure waveform, especially in diastole and early systole. Ignoring peripheral reflections leads to an erroneous indication of a reflection-free period in early systole and additional error in the estimates of (i) pulse wave velocity at the ascending aorta given by the PU–loop method (9.5% error) and (ii) transit time to a dominant reflection site calculated from the wave intensity profile (27% error). These errors decreased to 1.3% and 10%, respectively, when accounting for peripheral reflections. Using our new analysis, we investigate the effect of vessel compliance and peripheral resistance on wave intensity, peripheral reflections and reflections originating in previous cardiac cycles.


INTRODUCTION
Blood pressure and flow waveforms in systemic arteries carry valuable information for the diagnosis and treatment of cardiovascular disease and play a significant role in clinical conditions such as hypertension. The waveforms result from a complex ventricular-vascular interaction involving cardiac contraction, impedance of large and medium-sized distensible arteries and resistance of smaller arteries and arterioles. Blood behaves as an incompressible fluid in arteries, which distend to accommodate the sudden increase in blood volume delivered by cardiac contraction. When elastic energy stored in the distended arterial walls is released, arteries contract. The regular expansion and contraction of arteries (the pulse) that follows cardiac contraction propagates in the form of Figure 1. Blood pressure (P ) waveform measured in vivo along the human (left, modified from [1]) and rabbit (right) aortas. (a) Measurements were made every 10 cm down the aorta starting approximately 5 cm from the aortic valve. Each waveform is an ensemble average of continuous pressure measurements over 1 min using the peak of the R-wave of the electrocardiogram as the reference time. The circles indicate the time of the diastolic pressure (P d ) after which pressure increases due to left ventricular ejection of blood. The slope of the dotted lines connecting the circles indicates the pulse wave speed (6.9 m s 1 ) with which the pressure wavefront at the start of cardiac contraction propagates down the aorta. Systole is the phase of the cardiac cycle when the heart muscle contracts, and diastole is the phase when the heart muscle relaxes. (b) Measurements were made every 1 cm starting approximately 1 cm from the aortic valve, as described in Section 2.1. Each waveform is an ensemble average of ten continuous cardiac cycles using the foot of the flow waveform in the aortic root as the reference time.
pulse waves. These produce continuous changes in blood pressure and flow that can be studied as pressure and flow wavefronts (infinitesimal changes in pressure and flow) ‡ running forwards and backwards (away from and towards the heart, respectively), with backward wavefronts originating from reflected forward wavefronts at sites of vascular impedance mismatch. Figure 1 shows typical blood pressure waveforms measured in vivo along human (left) and rabbit (right) aortas, from the root to the aorto-iliac bifurcation, under normal conditions. The slope of the line joining the feet of these waveforms shows clearly that the pressure wavefront originated at the start of cardiac contraction propagates away from the heart; the measured space-averaged speed is 6.9 m s 1 in the human and 6.1 m s 1 in the rabbit. Thus, during a typical cardiac cycle, which takes about 1 s in the human and 0.25 s in the rabbit, a pulse wave has sufficient time to travel from the heart to the arterial vasculature and back multiple times.
Several studies have used wave intensity analysis (WIA) to investigate the role of wave reflections in shaping in vivo pressure and flow waveforms in systemic arteries [2][3][4][5][6][7][8][9][10][11], including the coronary circulation [12][13][14]. Given simultaneous measurements of blood pressure and flow velocity with time at an arbitrary location in the arterial network, we can calculate the local pulse wave velocity (PWV) and apply WIA to quantify the timing, direction and magnitude of the predominant waves that shape the pressure and velocity waveforms [15,16]. Accurate estimation of PWV is not only important for WIA but is also clinically relevant, because PWV is an important predictor of cardiovascular events [17].
Numerical modelling has been used to assess the following: (i) the ability of WIA to quantify reflection coefficients [18]; (ii) haemodynamic information provided by WIA in a model of aortic coarctation [19] and the fetal circulation [20]; (iii) a modified WIA based on the reservoir-wave separation [21]; and (iv) the performance of several methods for PWV calculation [22][23][24][25]. Numerically generated pressure and flow waveforms are free of measurement errors, and the theoretical values of haemodynamic properties that affect waveforms (e.g. PWV, location of reflection sites, and magnitude of reflected waves produced) are available for comparison with corresponding estimates given by WIA and methods of calculating PWV.
In the present work, we used pressure and flow velocity waveforms measured in vivo in the rabbit or generated numerically in several models of human compliant vessels to (i) show the inability of traditional WIA to identify the important role of peripheral reflections in shaping the pressure waveform; (ii) test the accuracy of the modified P U -loop method of calculating PWV proposed by Mynard et al. [24], which accounts for peripheral reflections originating in previous cardiac cycles; and (iii) propose a new analysis of arterial pulse wave propagation to study the predominant waves that shape pressure and flow waveforms during systole and the contribution to the pressure waveform, over the whole cardiac cycle, of wave reflections originating in previous cardiac cycles, vessel compliances, peripheral resistances, outflow pressures and the flow at the root. We used our new analysis to study the effects of vessel stiffness and peripheral resistance on numerically-generated aortic pressure and flow waveforms.
We generated all numerical data using the nonlinear one-dimensional (1-D) formulation of blood flow in compliant vessels, because WIA is derived from this formulation and, hence, 1-D model pressure and velocity waveforms provide an ideal mathematical framework for our study. Several comparisons against in vivo [26][27][28][29], in vitro [30][31][32][33][34] and 3-D numerical [35] data have shown the ability of the 1-D formulation to capture the main features of pressure and flow waveforms in large human arteries. The nomenclature and abbreviations used in this paper are listed in the supplementary material.

METHODS
We first describe the in vivo (Section 2.1) and numerical (Section 2.2) pressure and flow waveforms used in our work. Next, we summarise the mathematical formulation of traditional WIA (Section 2.3) and show how we quantified the magnitude and timing of wave reflections (Section 2.4). We then introduce the Windkessel model to study arterial haemodynamics during diastole and account for peripheral reflections in WIA and the P U -loop method (Section 2.5).

Rabbit in vivo pressure and flow waveforms
Ten New Zealand white male rabbits (Harlan UK) with the properties shown in Table I were maintained on a standard laboratory diet and housed at 18 ı C on a 12-h light cycle. The rabbits were pre-medicated with Hypnorm (0.1 ml kg 1 ) intramuscularly and anaesthetised with sodium pentobarbitone (35 mg kg 1 ), administered intravenously via the marginal vein of the right ear. They were then placed in the supine position and artificially ventilated with a Harvard Small Animal Ventilator (with a respiratory rate of 50 breaths min 1 ; inflation pressure set for trough 25 cmH 2 O and peak 50 cmH 2 O). Body temperature was maintained near 39 ı C (monitored using a rectal thermistor probe) by placing the rabbits on a heating blanket with a temperature controller (CWE Inc. TC1000). The animal procedures complied with the Animals (Scientific Procedures) Act (1986) and were approved by the Imperial College London local ethical review process.
For each rabbit, a midline thoracotomy was made and the rib cage was retracted to expose the ascending aorta and place a perivascular blood flow probe (6 mm diameter, type MA6PSB, Transonic Systems Inc.) around it, near the aortic valve. The thoracotomy was closed using clamps. The right femoral artery was exposed at the level between the knee and groin for insertion of pressure measurement wires. Two pressure waveforms were simultaneously measured with the aortic flow waveform, from the aortic root to the iliac artery in 1 cm increments (Figure 1(b)) using a dual sensor Millar Mikro-tip catheter transducer (model SPC-721, size 2.5F, sensors 5 cm apart). The aortic pressure waveforms at the root and 5 cm distally (Figure 2(b)) allowed us to calculate the PWV using the 'foot-to-foot' method (c ff ) [25], which we used as the gold standard PWV at the aortic root (Table I). Brief asystole (cessation of heart contraction) were induced by gently tapping the left ventricle (LV) (Figure 2(a,b)). All data were acquired at a sampling rate of 1 kHz using the NOTOCORD-hem acquisition software (NOTOCORD Systems, France).
Prior to euthanasia, heparin (2,000 units) was administered intravenously for anticoagulation of the blood in preparation for production of a resin cast of the arterial system (Figure 2(e)). The thorax was reopened and the LV was cannulated with a polythene tube (filled with saline solution: 9 g/L Table I. Properties of the 10 mature rabbits as follows: age, weight (W), crown-to-rump length (C-R), heart rate (HR), cardiac output (CO), mean pressure (P m ), pulse pressure (PP), outflow pressure (P out ), time constant of the decline in pressure during diastole (R T C T ) and pulse wave velocity at the aortic root calculated using the foot-to-foot method (c ff ), the traditional P U -loop (c PU ) and the modified O P U -loop ( O c PU ), with D 1, 015 kg m 3 in both loop methods.

Rabbit
Age (days) (kg) (cm) (beat/s) (ml/s) (kPa) (kPa) (kPa) (s) (m/s) (m/s) (m/s) Figure 2. In vivo (a) flow rate at the aortic root and (b) pressure at the aortic root (P1) and 5 cm distally in the aorta (P2) with time, measured during an asystole in Rabbit 8. (c) Time constant (R T C T ) and (d) asymptotic pressure (P out ) calculated along the aorta and iliac artery at 5 cm increments (indicated by yellow dots in (e)). R T C T and P out were derived from an exponential function (Equation (19)) fitted to the decline in pressure during an asystole (shaded area in (b)). NaCl) via an apical stab wound. The cannula was flushed through with approximately 5 ml saline. Casting resin (Batson's no. 17, Plastic Replica and Corrosion kits, Polyscience, USA, prepared according to the manufacturer's instructions) was infused into the arterial system via the cannula, at a pressure equivalent to the mean arterial pressure that was recorded in vivo in the thoracic aorta. Infusion was continued until the casting resin ceased to flow because of setting of the resin; the final perfusion volume was approximately 100 ml of resin. The cast was allowed to cure for at least 12 h; then the carcass was completely submerged in an aqueous solution of potassium hydroxide (25% w/v). The alkaline corrosion process was allowed to continue for 14 days. Then the corrosion solution was removed to reveal the arterial cast, which was submerged in a strong warm solution of detergent (Decon 90, Decon Laboratories Ltd., East Sussex, UK) and left for 24 h. The cast was rinsed gently and thoroughly in water and allowed to dry. We used the cast to locate accurately the sites of pressure measurements and determine the luminal cross-sectional area of the aortic root, which is required to convert measured flow rate to flow velocity; the latter, combined with simultaneous pressure at the root, allowed WIA.

Human numerical pressure and flow waveforms
To generate the numerical pressure and flow data for this study, we solved the nonlinear 1-D equations of blood flow in compliant vessels in a single-vessel model of the human thoracic aorta [35] ( Figure 3(a,b,c), Table II) and also in a model of the 55 larger systemic arteries in the human [1] ( Figure 3(d,e,f,g,h), Table III), using a DG scheme with a spectral/hp spatial discretisation [32]. These equations can be derived by applying conservation of mass and momentum to a differential 1-D control volume of the vessel [36,37], 8 < : where t is the time, x is the distance along the vessel, A.x, t/ is the area of the luminal cross section of the vessel, U.x, t/ is the axial blood flow velocity averaged over the cross section, P .x, t/ is the blood pressure averaged over the cross-section, is the density of blood (assumed to be constant), and f .x, t/ D 22 U is the frictional force per unit length, with D 4 mPa s the viscosity of blood. Both models are based on data from young and healthy humans. The single-vessel model is uniform and contains a single peripheral reflection site, which simplifies our initial assessment of WIA. The 55-artery model allows us to assess WIA in the presence of tapering and multiple outflows and reflection sites; all these properties have been shown relevant for analysis of pulse wave propagation phenomena [38][39][40][41].
To account for the fluid-structure interaction of the problem and close the system of equations (1), we obtained the following explicit algebraic relationship between P and A (or tube law). Assuming the arterial wall to be a thin, incompressible, homogeneous, isotropic, Voigt-type visco-elastic membrane, which deforms axisymmetrically, each cross-section independently of the others, we can relate blood pressure (P ) to luminal cross-sectional area (A) through [32] with P e .A,  [1]. They were calculated using the nonlinear purely elastic (a,c) and visco-elastic (d,f,g,h) 1-D equations (1) and (2). At the root of both models, we prescribed a flow rate measured in vivo (labelled 'Root' in (c,g)).
At the outlet of each terminal branch, we coupled a three-element Windkessel model of the perfusion of the microcirculation (b), with R 1 D Z 0 to minimise wave reflections [42]. The properties of the aortic model are shown in Table II. The names and properties of the segments in the 55-artery model are shown in Table III. The uniform Windkessel pressure, p w , given by Equation (17) is shown in red in (a,d,f).
where P e is the elastic component of pressure, h.x/ is the wall thickness, E.x/ is the Young's modulus, '.x/ is the wall viscosity, and A 0 .x/ is the reference area at P D 0 and @A @t D 0. For some simulations, we considered a purely elastic tube law (i.e. P D P e ) as it is assumed in WIA.
We implemented all the boundary conditions of our simulations and solved matching conditions at bifurcations by taking into account the correct propagation of the characteristic information and neglecting energy losses and visco-elastic effects (see [1] for a detailed description).
Both models exhibit the following characteristic features of the pressure and flow waveforms that are observed in vivo under normal conditions. The foot of the pressure and flow waveforms in early systole propagates away from the heart (compare the simulated pressures in Figure 3(a,d) with the in vivo pressures in Figure 1(a,b)). The pulse pressure (the difference between the maximum, systolic, and minimum, diastolic, pressures) increases in the aorta with increasing distance from Table II. Parameters of the single-vessel model of the human thoracic aorta coupled to a three-element Windkessel model of the rest of the systemic circulation ( Figure 3(b)). The resulting wave speed at mean pressure is 5.2 m s 1 .
the heart (Figures 3(a,d) and 1(a,b)), whereas mean pressure gradually decreases § . In the ascending aorta, pressure features a 'shoulder' or point of inflection (Figures 3(a,d) and 1(a,b)). From the ascending aorta to the upper thoracic aorta, a small pressure peak is observed at the start of diastole, which forms the dicrotic notch. This vanishes in the lower thoracic region (Figures 3(d) and 1(a,b)), but it is observed in other proximal arteries such as the common carotid arteries (Figure 3(f)) [43]. Moreover, the initial pressure increase becomes steeper and narrower in time in more peripheral locations (wave steepening) (Figures 3(a,d) and 1(a,b)), and a wide pressure peak appears in diastole from the abdominal aorta to the leg arteries [44,45] and in the brachial artery [46] (Figure 3(d,f)). With increasing distance from the heart, the aortic flow waveform (Figure 3(c,g)) becomes characterised by an increase in width, a reduction in the amount of reverse flow and a decrease in amplitude and mean value ¶ [44,45]. Reversed flow is absent in the suprarenal region of the aorta [29,47], renal arteries [44,48] and carotid arteries [29,49,50]. There is, however, a region of reverse flow in early diastole in the infrarenal region of the aorta and leg arteries [47,50].

Traditional wave intensity analysis
Wave intensity analysis is derived from the system of equations (1). Velocity and pressure waveforms are decomposed into successive wavefronts, with dP and dU as changes in pressure and velocity, respectively, across a wavefront. Fluid viscous losses are assumed to be negligible locally, and A is assumed to depend only upon P through a purely elastic tube law, with uniform and constant properties. Under these conditions, Riemann's method of characteristics applied to the system of equations (1) shows ( Figure 4(a)) that for any point .X , T / in the .x, t/ space, there are two characteristic paths, dt D U˙c, on which W f and W b are constant [15]. The quantities W f and W b are generally known as the characteristic variables or Riemann invariants and satisfy where the pressure-dependent c is the PWV, that is the speed at which pulse wavefronts travel in the absence of convective velocity (U ) (further details are given in Appendix A). For the purely elastic tube law given by Equation (3),  (Figure 3(e)). r in ! r out , mean cross-sectional radii at the inlet and outlet of the arterial segment (radii decrease linearly); c in ! c out , mean wave speed at the inlet and outlet of the segment.  Under physiological flow conditions, c is much greater than the maximum U , so that U C c > 0 and U c < 0 (i.e. the flow is subcritical). Therefore, changes in P and U propagate in the forward and backward directions (we define the forward direction as the direction of mean blood flow, in which x increases) with speeds of U C c and U c, respectively.
In vivo measurements of P and U with time are typically taken at a fixed point x D X , rather than along a characteristic line. Solving the two equations in (4) at x D X for dP and dU yields Wave intensity (dI ) is defined as [15,16] which is the flux of energy per unit area carried by the wavefront as it propagates and has dimensions of power/unit area and SI units W m 2 . Note that dI is calculated with time at a fixed point x D X . According to Equation (7), dI is positive if dW f > dW b and negative if dW f < dW b . Therefore, dI 'measures' the importance with time of changes in P and U in the forward and backward directions at x D X . Whenever dI > 0, forward changes in P and U dominate over backward changes; the flow is accelerated if dP > 0 and decelerated if dP < 0. Whenever dI < 0, backward changes in P and U dominate over forward changes; the flow is accelerated if dP < 0 and decelerated if dP > 0.

Application to measured data.
Given simultaneous measurements of P .t/ and U.t/ at an arterial site of the numerical models or rabbit, we used a Savitzky-Golay filter to smooth where dt is the time between two adjacent sampling points of P or U . This filter is commonly used for WIA because it preserves peaks in dP and dU (and hence dI ) [16,51] [52, p. 650]. We normalised the value of dI given by Equation (7) by dt 2 to make the magnitude of dI independent of the sampling frequency. We used customised MATLAB software (The MathWorks, Inc., MA, USA) for our data analysis.

Forward and backward waveforms.
The measured waveforms P .t/ and U.t/ can be separated into forward-travelling (P f .t /, U f .t /) and backward-travelling Separating dP and dU into changes across the forward (dP f , dU f ) and backward (dP b , dU b ) wavefronts, that is dP D dP f C dP b and dU D dU f C dU b , and using the water hammer equations, yield [15] A derivation of Equations (8) from the system of equations (1) using the method of characteristics is given in [15]. In Appendix A, we provide an alternative derivation of Equations (8) by directly applying conservation of mass and momentum to a control volume moving with the forward or backward pulse wavefronts.
If the PWV (c) is known, Equation (9) allows us to obtain P f,b .t / and U f,b .t / from the measured P .t/ and U.t/ by adding the differences dP We considered the integration constants P 0 and U 0 to be half the pressure and velocity, respectively, at the end of diastole. Thus, wave intensity (dI ) at a fixed point x D X can be separated into forward (dI f > 0) and backward (dI b < 0) components, We can calculate c from simultaneous measurements of P and U using the P U -loop [15,53] method (see Appendix B for more details) or from two measurements of P (or U ) at two different sites using foot-to-foot, least squared difference or cross-correlation techniques [25].

Wave reflections
At sites of impedance mismatch, pulse waves are partly reflected and partly transmitted. Using a linearised version of the 1-D equations (1), we can relate the changes in pressure and flow velocity across the reflected wavefront (dP ref , dU ref ) in the parent and two daughter vessels of an arterial bifurcation to the corresponding changes in pressure and velocity in the incident wavefront (dP inc , dU inc ) through The superscripts p, d1 and d 2 refer to the parent and first and second daughter vessels, respectively. The reflection coefficients for wavefronts propagating in the parent R p f and daughter (R d1 f and R d 2 f ) vessels can be expressed as a function of the characteristic admittance, where c 0 is the PWV at zero pressure (see Appendix C for a detailed derivation of Equations (11) and (12)), At the outlet of each terminal branch coupled to a single resistance (R 1 ) we have [42] dP where P out is the outflow pressure, R f D .R 1 Z 0 /=.R 1 C Z 0 / is the local reflection coefficient and Z 0 D 1=Y 0 is the characteristic impedance of the terminal branch. Note that R 1 D Z 0 yields R f D 0, in which case dP inc and dU inc are completely absorbed by the outflow model. From Equation (11) to Equation (13) with P out D 0, we can relate the wave intensity of the reflected wavefront (dI ref D dP ref dU ref ) to the wave intensity of the incident wavefront (dI inc D dP inc dU inc ) through Given a wave intensity profile separated into dI f and dI b , we defined the initial positive and negative regions as the incident and reflected waves, respectively ( Figure 5). Note that these waves are made of many wavefronts across which wave intensity (dI ) changes. We then calculated the time of arrival (T inc and T ref ) and magnitude (I inc and I ref ) of the incident and reflected waves using either peak values or area-average values; they are indicated by empty squares or stars, respectively, in Figure 5. Given the wavefronts that make up the incident, dI inc , or reflected, dI ref , wave, we defined the area-average values as where the sum is taken from the start to the end of the wave, and t ini and t end are, respectively, the times at the start and end of the wave. These area-average values account for all the wavefronts that make up the incident and reflected waves and not only the wavefronts that define their peaks. Finally, we calculated the transit time (TT) between the incident and reflected waves and the apparent reflection coefficient (R app f ) based on Equation (14) as Thus, we obtained a pair of TT and R app f using peak values and another pair using area-average values. We also calculated TT using two foot values, which we defined as the point when dI inc and dI ref are greater than 1% of the peak I inc and I ref , respectively (they are indicated by empty circles in Figure 5).

Haemodynamics during diastole
We can describe pressure and flow during diastole using a zero-dimensional Windkessel model. At any point in a distributed 1-D model, pressure becomes increasingly well described with increasing time in diastole by a space-independent Windkessel pressure, p w .t / (Figure 3(a,d,f)), given by [39,54] where N is the number of arterial segments, with i D 1 the aortic root and j D 2, : : : , M terminal segments coupled to resistance-compliance-resistance (RCR) Windkessel models (Figure 3(b)), M 1 < N is the number of terminal branches, Q in .t / is the flow waveform at the aortic root, R T is the net peripheral resistance of the arterial network, C T is the total systemic compliance, C c is the total conduit compliance, C i 0D D A i 0 l i =. .c i / 2 / is the compliance of segment i (with luminal area A i 0 , length l i and PWV c i ), C p is the total peripheral compliance, p w .T 0 / is the pressure p w at the reference time t D T 0 (T 0 = 0 in Figure 3(a,d,f)) and q j out .t / is the outflow in the terminal segment j . The parameters of the RCR Windkessel models are R 1 D Z 0 , C , R 2 and P out (Figure 3(b)). During diastole, it is reasonable to assume Q in D 0 and dq j out .t / dt D 0, j D 2, : : : , M , in normal conditions, which reduces p w to Equations (17) and (19) neglect nonlinearities, flow inertia and flow viscous dissipation within the 1-D model arterial segments and assume that wall compliance and fluid peripheral resistance are the dominant mechanisms of blood flow. These are reasonable assumptions towards the end of diastole, when Equation (19) provides accurate predictions of blood pressure, as shown by in vivo studies in dogs [55] and numerical solutions of the three-dimensional Navier-Stokes equations in compliant vessels [35].
The flow rate (q w ) driven by p w P out during diastole is linearly dependent on x in each arterial segment i D 1, : : : where q i in .t / is the flow rate at the inlet of the segment. The wavefronts associated with p w and q i w during diastole at a fixed point x in Segment i, i D 1, : : : , N , are respectively, with dp w .t / uniform in space. The wave intensity dI i w given by dp w and dq i w is where dp w dq i in is the wave intensity at the inlet of Segment i and we have used At the ascending aorta, we have dp w dq in D 0, since we assumed zero flow at the root during diastole. Assuming the aorta to be a uniform vessel without branches, we have dI w D .dp w / 2 x= R T C T .c/ 2 during diastole. Note that .dp w / 2 is exponential with a time constant R T C T =2 (see Equation (19)). During the systolic ejection, we can relate the velocity (dU s ) and pressure (dP s ) wavefronts using dU s D dP s =. c/ [16,53], which yield an early-systole wave intensity dI s D .dP s / 2 =. c/. Finally, dI w and dI s are related through which shows that local (c) and global (R T and C T ) properties of the vasculature are responsible for the change of wave intensity from early systole to late diastole.

RESULTS AND DISCUSSION
We first study the predominant waves identified by traditional WIA in the aorta of our numerical models and rabbits and discuss their role in shaping the pressure and flow waveforms (Section 3.1 to 3.5). In particular, we show that the wave intensity profile is misleading by indicating a reflectionfree period in early systole, which introduces additional error in the estimate of PWV at the aortic root by the traditional P U -loop method (Section 3.2). We then show that WIA does not identify peripheral reflections (Section 3.6) and assess the ability of the Windkessel model to describe haemodynamics during diastole (Section 3.7) and the modified P U -loop method for PWV calculation suggested in [24] (Section 3.8). We describe our new WIA (Section 3.9), which accounts for peripheral reflections using the Windkessel model, and use it to study the effect of vessel compliance and peripheral resistance on blood pressure and flow waveforms (Section 3.10). Lastly, we analyse the sensitivity of our new WIA to sampling frequency and errors in the estimate of PWV (Section 3.11).

Predominant wave intensity waves in the aorta
At any point in the single-vessel aortic model coupled to a matched three-element Windkessel outlet (Figure 3(b)), we identified four predominant waves in the cardiac cycle from the forward (dI f ) and backward (dI b ) wave intensity profiles as follows (Figure 6(c) shows the profiles at the inlet): (i) an initial forward compression wave (i.e. travelling towards the outlet, because dI > 0, and with dP > 0) produced by the increase in blood flow at the inlet (inflow) in early systole; (ii) a backward compression wave indicating reflection of the initial wave at the outlet; (iii) a forward decompression (dP < 0) wave (a 'suction' wave) caused by the decrease in inflow towards the end of systole; and (iv) a forward compression wave due to the short increase in inflow at the end of systole. The flow is accelerated by the first and fourth waves and decelerated by the second and third. These four waves were also observed in the ascending aorta of the 55-artery model (Figure 3(e), results not shown). Furthermore, the four waves were observed in the ascending aorta of the rabbit (Figure 6(d)). Their timing and direction of propagation suggest that the first wave is caused by the contraction of the LV, the second is due to reflections of the initial contraction wave at downstream sites of impedance mismatch, the third is produced by the relaxation of the LV and the fourth results from the closure of the aortic valve. The first three waves have been extensively reported for the human aorta [15,16]. The fourth wave has been observed in human coronary arteries [12,14]; the aorta of the dog [53,57], sheep [58] and swine [59]; and large arteries of the fetal lamb [20].

Left ventricular-contraction wave -an apparent reflection-free period in early systole
Both numerical and in vivo results in Figure 6(c,d) show that there is a period in early systole during which dI b is almost zero. This indicates that changes in pressure and velocity are generated by the contraction of the LV only and not affected by wavefronts reflected in the vasculature. Indeed, the increase in pressure (Figure 6(a,b), orange dashed lines) is made up entirely of forward pressure wavefronts (dP f ) during this apparently reflection-free period (Figure 6(e,f), orange dashed lines). According to Equation (10), if dI b D 0 in early systole, then dP D c dU ; that is pressure (P ) and flow velocity (U ) are proportional and make the linear part of the P U -loop ( Figure 6(g,h), orange dashed lines).
However, Figure 6(e,f) shows clearly that the backward pressure (P b ) and velocity (U b ) waveforms are not constant during the apparently reflection-free period in early systole: P b is made up of backward pressure wavefronts (dP b ) that decrease P , and U b is made up of backward velocity wavefronts (dU b ) that increase U . According to Equation (9), non-zero dP b and dU b indicate that dP D c dU is not satisfied. Indeed, from the slope of the P U -loops (Equation (B.1)) in  (5) with mean area. The mean error of c PU in our ten rabbits relative to c ff is 9.5% (Table I), which is similar to the error reported in [24] using numerical data only. Therefore, the wave intensity profile is misleading by indicating a reflection-free period in early systole, which leads to an error in the estimate of c by the P U -loop method.

Reflection of the left ventricular-contraction wave -pressure augmentation
Reflections of the forward wavefronts that make up the LV-contraction (compression) wave yield a reflected wave that is made up of negative velocity wavefronts (dU ref < 0) and positive pressure wavefronts (dP ref > 0), the wave labelled 'reflection' in Figure 6(c,d). This wave produces a U b that decelerates the net forward flow and a P b that augments pressure (Figure 6(e,f)), thereby generating the 'shoulder' or point of inflexion (Figure 6(a,b)) that defines the pressure augmentation index [60].
In the single-vessel aortic model, reflected wavefronts originate from the outlet, which is coupled to a matched three-element Windkessel model (Figure 3(b)). In the more realistic 1-D distributed model, however, the origin of the reflected wave is more complex. Several studies [38][39][40][41] have identified multiple reflection sites that reflect wavefronts towards the aorta, where they arrive at different times with a magnitude that decreases exponentially with time [39]. The net effect of multiple reflection sites on the TT and R app f given by Equation (16) along the aorta is an apparent reflection site that appears to move away as the measurement location approaches it [40].

Left ventricular-relaxation wave -flow deceleration
The LV-relaxation (forward) wave drops pressure and flow during the deceleration phase towards the end of systole (Figure 6(a,b)). Flow velocity is reduced by both forward (U f ) and backward (U b ) velocities and becomes negative at the end of systole, whereas pressure (P ) is reduced by only the forward pressure (P f ); P b continues increasing P (Figure 6(e,f) Figure 3(e)) using the visco-elastic (Equation (2)) or purely-elastic ( D 0) tube law.
LV-relaxation phase, wavefronts originating at the aortic root decrease both pressure and velocity, whereas wavefronts reflected in peripheral locations increase pressure and decelerate the flow. Two mechanisms have been suggested to describe the origin of the LV-relaxation wave, active myocardial relaxation and flow inertia, which are discussed in [61, p. 95].

Valve wave -dicrotic notch
The closure of the aortic valve produces a forward compression wave that augments pressure, creating the dicrotic notch, and accelerates the flow, from negative to zero at the root (Figure 6(a,b,c,d)). The dicrotic notch is shaped by the forward pressure component (P f ) ( Figure 6(e,f)), that is by wavefronts travelling from the aortic root, and progressively vanishes in the aorta with increasing distance from the heart (Figure 1(a,b)). In the thoracic aorta of the 55-artery model (Figure 3(e)), the dicrotic notch is absent if wall viscosity is modelled and is present if wall viscosity is neglected (Figure 7(a)). Therefore, wall viscosity dissipates the dicrotic notch. Wall viscosity also dissipates all predominant waves in the wave intensity profile as they travel towards distal locations, almost abolishing the valve wave (Figure 7(b)). In vivo measurements in humans show that the valve wave is small or absent in the brachial, radial and femoral arteries [7,11].

Why does wave intensity vanish during diastole? The importance of peripheral wave reflections from previous cardiac cycles
Numerical and in vivo wave intensity profiles in the aorta during diastole do not show any reflection of the four predominant waves discussed previously ( Figure 6(c,d)). Here, we show that WIA fails to identify the important contribution to the pressure waveform of peripheral reflections, using our single-vessel (Section 3.6.1) and distributed (Section 3.6.2) models.

3.6.1.
Single-vessel aortic model. According to this model, wave intensity vanishes during diastole because reflected wavefronts are spread in time by peripheral compliance. To illustrate this, at the inlet of the vessel we prescribed a narrow Gaussian-shaped flow (Figure 8(a)) that approximates a compression wavefront immediately followed by a decompression wavefront. At the outlet, we first coupled a single-resistance with a zero outflow pressure, so that peripheral compliance is absent. This model generates multiple Gaussian-shaped pulse waves (Figure 8(b)). As predicted by Equation (13), both outlet and inlet change the direction of the blood flow; the outlet yields reflected waves travelling towards the inlet with 83% (R f D .R 1 Z 0 /=.R 1 C Z 0 / D 0.83) of the pressure and flow amplitudes of the incident wave, and the inlet behaves as a closed end (R f D 1, because Q in D 0 by the time the first pulse wave is reflected) producing reflected waves with the same amplitude as the incident wave.
The wave intensity profile is able to identify the timing, direction and magnitude of the multiple waves reflected in the model (Figure 8(d)). Equation (16) applied to the wave intensity profile in the midpoint yields errors in the R app f estimates smaller than 2% using peak values and 4% using area-average values (Equation (15)), both relative to the theoretical R f D 0.83. Furthermore, the errors in the estimates of TT in this model are smaller than 1% using foot, peak and area-average values, which shows the ability of Equation (16) to estimate accurately TT and R f in the presence of a single peripheral reflection site without peripheral compliance. Adding peripheral compliance using, for example, a matched RCR Windkessel outflow model (i.e. with R 1 D Z 0 to minimise wave reflections [42], Figure 3(b)) has the effect of smoothing reflected waves (Figure 8(c)), which is necessary for obtaining physiological-like waveforms (e.g. those in Figure 6(a,c)). Although reflected waves alter the pressure and flow waveforms after t D 0.1 s (Figure 8(c)), the wave intensity profile does not reveal them (Figure 8(e)), because the initial pulse wave has a wave energy several orders of magnitude greater than its reflections.

Distributed 1-D model.
We can study the role of peripheral reflections in shaping the pressure waveform using the following methodology. Neglecting nonlinear effects, we can separate the pressure and flow waveforms at an arbitrary point in a given distributed model into a conduit (or arterial) waveform, which is made up of pulse wavefronts propagating from the aortic root and being reflected at the arterial junctions, aortic root and tapered vessels, and a peripheral waveform, which is made up of wavefronts originating from reflections at terminal branches. As detailed in [39], the conduit waveform is obtained by running the simulation with each terminal branch coupled to a single resistance equal to the characteristic impedance of the branch, so that any wavefront leaving the vessel is completely absorbed by the boundary condition; that is R j 1 D Z j 0 so that R j f D 0 in Equation (13), j D 2, : : : , M . The peripheral waveform is the difference between the total and conduit waveforms.
In the aorta and main branches, most of the pressure waveform at the start of systolic ejection is made up of peripheral reflections, as Figure 9(a) shows for the left common carotid artery. These reflections originated in previous cardiac cycles, because they are present before the start of cardiac ejection in the current cardiac cycle. The maximum contribution of peripheral reflections to the pressure waveform occurs in early diastole, when these reflections produce a concave shape in the proximal arteries of our 55-artery model (Figures 3(d) and 9(a)). This feature of the pressure waveform is sometimes present in vivo in the aorta of the human [3,60,62], rabbit (e.g. Figure 6(b)) and dog  Figure 3(e)). [5,55] and in the human carotid, brachial and radial arteries [6][7][8][9]11]. Furthermore, in the carotid artery, peripheral reflections provide a pressure gradient that drives less than 50% of the mean blood flow in systole, but more than 90% of the mean flow in diastole (Figure 9(b)). All wavefronts that did not originate as peripheral reflections and which make up the conduit waveform are responsible for the features of the pressure and flow waveforms in systole (e.g. pulse pressure and dicrotic notch, Figure 9(a); flow amplitude, Figure 9(b)). In diastole, these wavefronts produce a conduit waveform that decreases exponentially and contributes little to the pressure and flow waveforms in the next cardiac cycle.
Although the diastolic pressure and flow waveforms in our distributed 1-D model are mainly made up of peripheral reflections, WIA fails to identify them. The wave intensity profile only reveals the waves that shape the conduit waveform, as Figure 9(c) shows: the conduit pressure and flow waveforms produce a profile similar in shape and magnitude to the 'total' profile, whereas the 'peripheral' wave intensity profile has a considerably smaller magnitude than the 'total' profile, because changes in pressure and velocity across the wavefronts that make the peripheral waveform are smaller than those that make the conduit waveform. As a result, the apparent reflection coefficient given by Equation (16) accounts little for peripheral reflections. Indeed, calculation of reflection coefficients using forward and backward pressures produces smaller errors than using wave intensity profiles [18].
We can use Equation (23) to quantify the reduction in wave intensity in the aorta from systole (dI s ) to diastole (dI w ) because of net peripheral resistance (R T ) and total compliance (C T ) acting on the blood flow through the Windkessel effect. Taking x D 20 cm, c D 5 m s 1 and R T C T = 1.8 s as representative values for the thoracic aorta and dp w D dP s , then dI w is 2% of dI s . However, dp w dP s because the magnitude of pulse wavefronts decreases as they get reflected and attenuated in the arterial network. Therefore, dI w is several orders of magnitude smaller than dI s (Figure 7(b)).
Next, we discuss the validity of the Windkessel model to study haemodynamics during diastole.

Using the Windkessel model to study haemodynamics during diastole
Our rabbit in vivo data show that the space-independent pressure p w .t / given by Equation (19) is able to describe pressure towards the end of diastole. Fitting Equation (19) to the decline in pressure during asystole measured along the aorta and iliac artery at 5 cm increments (at least five measurements were taken for each rabbit) produced a small variability of the time constant (R T C T ) and asymptotic pressure (P out ); for example, R T C T D 318˙10 ms and P out D 4.5˙0.2 kPa in rabbit 8 (Figure 2(c,d) and Table I).
Despite the rabbit having about four times the human heart rate, human and rabbit pressure waveforms are similar in shape when normalised in time (compare Figures 1(a,b)). Therefore, pressure in diastole falls at a faster rate with time in the rabbit, as indicated by the average time constant R T C T in the rabbit (0.29˙0.02 s; Table I) being one order of magnitude smaller than in the human (1.79˙0.13 s [46]). This is due to a smaller total systemic compliance (C T ) in the rabbit (0.21˙0.05 mm 3 Pa 1 ) than in the human (17.0˙1.0 mm 3 Pa 1 [46]); vessels are less compliant (i.e. C 0D D A 0 l=. c 2 / is smaller) for the shorter and thinner rabbit vessels, despite having similar PWV to the human vessels.
In agreement with Equation (20), the flow in the last part of diastole increases linearly from the inlet to the outlet of the single-vessel aortic model (Figure 3(c)) and from the aortic root to the thoracic aorta of the 55-artery model (Figure 3(g)). For a given point x in the aorta, the flow towards the end of diastole decays exponentially (Figure 3(c,g,h)). Furthermore, the P U -loop is approximately linear in the last part of diastole ( Figure 6(g,h), green dashed lines), which is in agreement with p w and q i w being proportional in the last part of diastole, according to Equation (20). These results provide further evidence to support the conclusions in [35,55] that the Windkessel model has application for study of haemodynamics in diastole. Application of this model enables us to account for wave reflections originating from previous cardiac cycles and thereby modify the P U -loop method (Section 3.8) and WIA (Section 3.9). Figure 10(c,d) shows pressure cleared of contributions of waves reflected in previous cardiac cycles O P Á versus flow velocity (U ) at the root of the single-vessel aortic model and in vivo rabbit aorta. In the numerical model, these contributions were obtained by prolonging the decaying p w given by Equation (19) from the previous cardiac cycle into the current cycle, with the time constant (R T C T ) and asymptotic pressure (P out ) given by the theoretical parameters of the model. In the rabbit, R T C T and P out were calculated by fitting an exponential function of the form given by Equation (19) to the decline in pressure during several asystole generated along the aorta and iliac artery at 5 cm increments (Figure 2(b,e)). The values of R T C T and P out were taken to be the corresponding mean values for each rabbit (Figure 2(c,d)). Figure 10(a,b) shows p w prolonged from the previous cardiac cycle and the pressure waveform O P D P p w , which is cleared of contributions of reflected waves originating from previous cardiac cycles, for the numerical and in vivo data. We did not clear U of contributions of reflected waves originating from previous cycles, because U approximates zero at the aortic root. In more distal locations, these contributions could be obtained by prolonging the decaying U from the previous cardiac cycle into the current cycle. Theoretically, the decaying U should be approximated by q w =A d given by Equation (20) Figure 10(e,f), orange dashed lines) and, hence, d O P D c dU is satisfied in early systole (unlike dP D c dU ). As a result, the slope of O P versus U in Figure 10(c,d) (orange dashed lines) yields O c PU D 5.2 m s 1 in the aortic model, in agreement with the theoretical value calculated using Equation (5) with A the mean area, and O c PU D 4.4 m s 1 in the rabbit aorta, in agreement with the value calculated using the foot-to-foot method (Table I). In all ten rabbits, the modified O P U -loop method provided smaller errors than the traditional P U -loop method in the estimates of PWV in the ascending aorta, relative to the foot-to-foot estimates (three last columns of Table I). The mean relative errors for the ten rabbits were 1.3% using the modified O P U -loop method and 9.5% using the traditional P U -loop method, which are comparable to the errors reported in [24] obtained from numerical data.

Modified O P U -loop method
Our results support the use of the modified O P U -loop method when an exponential fit to the decline in pressure during diastole is possible. An exponential fit may be challenging when using in vivo data without asystole and when the pressure decay does not develop into an exponential (e.g. see Figure 6(b)) [63]. Tapering and energy losses are other potential sources of error when using the P U -loop, because they are not accounted for in the derivation of Equation (B.1) (see Appendices A and B).
Lastly, we note that it is not necessary to eliminate the pressure contribution from previous cardiac cycles in order to estimate accurately the local PWV from simultaneous P and luminal area (A) measurements at an arbitrary location in the arterial network. This is because pressure and luminal area wavefronts are reflected with the same reflection coefficient (R f ), whereas velocity wavefronts are reflected with R f and, hence, reflections from previous cardiac cycles have a similar effect on and O U f,b were calculated using the PWV given by (e,g) Equation (5) with A the mean area or (f,h) the foot-to-foot method as described in Section 2.1. For all data contours, the reflection-free period in early systole and the decay in pressure with approximately constant velocity are highlighted in orange and green, respectively.
P and A, unlike on P and U (see Appendix C). Applying Equation (B.3) from Appendix B to P and A at any point of the visco-elastic 55-artery model (Figure 3(e)) yields errors in the estimates of PWV smaller than 2%, relative to the theoretical PWV obtained using Equation (5) with A the mean area over the cardiac cycle. Such a small error, however, remains to be verified using in vivo P and A data.

Novel wave intensity analysis
We propose a modified WIA that consists of (i) the new wave intensity d O using O P and U to study haemodynamics during systole (Section 3.9.1); (ii) the Windkessel model given by Equation (19) to study haemodynamics during diastole (as we showed in Section 3.7) and quantify the contribution to the pressure waveform of wave reflections originating from previous cardiac cycles (Section 3.9.2); and (iii) the Windkessel model given by Equation (17) (with q j out D 0, j D 2, : : : , M , when analysing in vivo data) to analyse the effect on the pressure waveform during both systole and diastole of vessel compliances, peripheral resistances, outflow pressures and the flow at the root (Q in .t /) (Section 3.10).

Haemodynamics during systole. The modified d O
I f,b profiles are not influenced by wave reflections originating from previous cardiac cycles and, hence, allow us to do WIA starting from a 'true' reflection-free period in early systole. In normal conditions, they are very similar to the traditional dI f,b profiles (e.g. compare Figure 6(c,d) with Figure 10(g,h), and 'visco-elastic' in Figure 7(b) with 'M1' in Figure 11(d)) and contain the four predominant waves described in Sections 3.1-3.5. However, using foot values of d O I f and d O I b decreases the error of the TT estimate provided by Equation (16) in the single-vessel aortic model coupled to a matched RCR Windkessel model (Figure 3(b)). From the profiles dI f and dI b in the midpoint of the vessel, calculated with c at diastolic pressure, we obtained TT D 67 ms using peak values, TT D 102 ms using area-average values (Equation (15)) and TT D 65 ms using foot values. Their errors are, respectively, 31%, 100% and 27%, relative to the theoretical 51 ms calculated as l=.2.cCU inc //Cl=.2.c U ref //, according to the characteristics analysis (see Figure 4(b)), where l D 24.1 cm is the vessel length, c D 5.0 m s 1 is the theoretical PWV at diastolic pressure given by Equation (5), U inc D 0.04 m s 1 is the average flow velocity during the propagation towards the outlet of the wavefront that forms the foot of the  Table IV. Contribution from the current cardiac cycle, the first (I), second (II), third (III), fourth (IV) and fifth (V) previous cycles, and the outflow pressure (P out ) to the pressure waveform (P ) at the aortic root of the single-vessel aortic model (Figure 10(a)) and in vivo Rabbit 8 (Figure 10(b)). Each contribution is quantified using p w , as described in Section 3.9.2, and is expressed in % relative to the area under P .
incident wave, and U ref D 0.5 m s 1 is the average flow velocity during the propagation towards the inlet of the wavefront that forms the foot of the reflected wave. Similar relative errors were obtained at any point in the vessel, suggesting that foot values should be used for the calculation of TT. The relative error given by foot values of d O I f and d O I b was reduced to 10%. Conceptually, the new wave intensity profile should be used to study haemodynamics during systole. In practice, however, the new and traditional wave intensity profiles may be very similar due to measurement errors.

Contributions of wave reflections from previous cardiac cycles.
The contribution of wave reflections from several previously occurring cardiac cycles to the pressure waveform (P ) can be studied by prolonging the decay of p w given by Equation (19) at the end of each previous cycle. Figure 10(a,b) shows the contribution to P at the aortic root in the numerical model and rabbit in vivo of the second (II), third (III), fourth (IV) and fifth (V) previous cycles, and Table IV quantifies them. In the rabbit, the pressure at which flow to the microcirculation ceases (P out ) contributes 4.5 kPa to P . In the numerical model, P out was set to zero and does not play a role in shaping P . In both the numerical model and the rabbit, the current cycle contributes to P more than any previous cycle, providing the systolic features of pressure that can be studied using d O I , and, in the rabbit, the diastolic concave shape. Pressure contributions from previous cycles decrease exponentially, with earlier cycles generating a smaller percentage of P than later cycles.
Because conduit pressures vanish almost completely by the end of diastole (Figure 9(a)), pressure contributions from previous cardiac cycles are mostly made up of reflected waves originating from peripheral reflection sites. These reflected waves persist for several cardiac cycles, because they become trapped within the arterial network between the aortic valve and peripheral vessels [39]. In diastole, the valve is closed and so reflects backward wavefronts with a reflection coefficient close to one [39,56].

Effect of vessel compliance and peripheral resistance
We use our modified WIA to study the effects of changes in vessel compliance and peripheral resistance on the pressure and flow waveforms in the thoracic aorta of the visco-elastic 55-artery model (Figure 3(e)), with the parameter values described in Table III (hereinafter referred to as Model M1, normal young) and the following three variations: Model M2, normal old: The elastic modulus E is increased three-fold in all the arterial segments (except for the terminal branches) to decrease the total conduit compliance (C c ) and model the increase in the stiffness of the arterial wall with ageing [64]. Model M3, hypertensive young: The total resistance R 1 C R 2 at the outflow of each 1-D model terminal branch (Figure 3(b)) is increased by 40% to raise blood pressure. Model M4, hypertensive old: The changes introduced in M2 and M3 are combined. Table V. Contribution from the current cardiac cycle, the first (I), second (II), third (III), fourth (IV) and fifth (V) previous cycles, and the outflow pressure (P out ) to the pressure waveform in the thoracic aorta (midpoint of Segment 27 in Figure 3(e)) of the normal young (M1), normal old (M2), hypertensive young (M3) and hypertensive old (M4) models. Using the Windkessel pressure (p w ) as a zero-order approximation to the pressure waveform (P ) provides relevant information on the effects of vessel compliance and resistance. First, we note that the total compliance (C T ) and peripheral resistance (R T ) appear in the term e t 0 =R T C T of the integrand of Equation (17). This term increases the contribution to P of the decaying flow out of the heart (Q in ; Figure 3(g), 'Root') in the last part of systole; the greater C T and R T are, the more they spread the shape of Q in (and hence p w ) in time. Furthermore, the magnitude of P produced by Q in in systole decreases with increasing C T through the term e t=R T C T =C T multiplying the integral. Thus, the pulse pressure decreases with increasing C T . According to Equation (18), C T increases with the increasing compliance of each segment in the arterial network (C 0D ), which is greater in young models (M1 and M3) than in old models (M2 and M4). As a result, the pulse pressure of both p w and P at the thoracic aorta is less in young than old models (Figure 11(a,c)). In vivo measurements in humans show an increase in pulse pressure with age [62] and disease (such as atherosclerosis and diabetes) [65], which can be explained by a decrease in vessel compliance. Lastly, both P and p w decay in diastole at a greater rate in old than in young models (Figure 11(a,c)), because the time constant (R T C T ) is less in old than young models, in agreement with Equation (19).
Using p w to quantify the portion of P at the thoracic aorta originating from previous cardiac cycles reveals similar contributions to P from the first previous cycle in the four models studied (Table V). However, in old models a smaller portion of P is made up of reflected waves originating from earlier cycles because p w falls at a faster rate with decreasing C T , and in hypertensive models a greater portion of P comes from earlier cycles (compare M1 with M3 and M2 with M4 in Table V) because p w falls at a slower rate with increasing R T .
We can study how C T and R T affect systolic haemodynamics using the wave intensity profile d O I . The peak values of d O I are greater with decreasing C T and, hence, the energy carried by pulse wavefronts in old models (M2 and M4) is greater than in young models (M1 and M3) (Figure 11(d)). This result suggests that the LV must produce more energy to propel the same amount of blood flow throughout the vasculature of old models. On the other hand, normal and hypertensive models of the same age (i.e. M1 and M3 or M2 and M4) have similar contours of d O I ; that is, changes in R T have little effect on d O I . Therefore, peripheral reflections have a minor effect on the flux of energy in the thoracic aorta in systole. This result provides further evidence to support our result in Section 3.6 on the inability of the wave intensity profile to identify peripheral reflections. Furthermore, this result suggests that changes in R T have little effect on aortic augmentation index, in agreement with [41].
The foot of the incident wave in the d O I f profile occurs earlier in old models (Figure 11(d)), which feature greater PWV in proximal vessels, but the TT between the incident and reflected waves calculated using Equation (16) with foot values changes by less than 2 ms in all four models. Furthermore, the apparent reflection coefficient R app f is not affected by changes in R T , but it increases with decreasing C T ; R app f is 7% greater in M2 than in M1. Different results are obtained, however, in the single-vessel aortic model; TT decreases with decreasing C T and increasing R T , in agreement with theoretical predictions, and R app f increases with increasing C T and R T (see supplementary  (Figure 12(d)). Lastly, using the local minimum between the fictitious and genuine reflected waves to determine the foot of the reflected wave, the TT estimates given by foot values are more sensitive to PWV errors withiṅ 20% than those given by peak values (Figure 12(e)). Our results suggest a small variability of d O I f,b in the aorta using estimates of PWV with errors up to˙20%, in agreement with results reported for traditional dI f,b using in vivo measurements in human coronary arteries [67]. In Section 3.8, we showed that estimates of local PWV at the ascending aorta with errors smaller than 20% are feasible from simultaneous pressure and flow velocity measurements using the modified O P U -loop method.

CONCLUSIONS
We have presented a novel analysis of arterial pulse wave propagation that combines traditional WIA with identification of Windkessel pressures to account for the effect on the pressure waveform of peripheral wave reflections. Our modified wave intensity profiles allow us to study the predominant waves that shape pressure and flow waveforms during systole, starting from a 'true' reflection-free period in early systole. The Windkessel pressure allows us to (i) quantify the contribution to the pressure waveform over the whole cardiac cycle of reflected waves originating in previous cardiac cycles and (ii) study the buffering effect of the vasculature, which depends on vessel compliances, peripheral resistances, outflow pressures and the flow at the root (Section 3.9). We have shown that traditional WIA identifies the timing, direction and magnitude of the predominant waves that shape aortic pressure and flow waveforms in systole (Section 3.1 to 3.5) but fails to identify the important contribution to the pressure waveform of peripheral reflections. These reflections persist for several cardiac cycles and make up most of the pressure waveform, especially in diastole and at the start of cardiac contraction (Section 3.6). Ignoring the contribution of peripheral reflections to the pressure waveform leads to an erroneous indication of a reflection-free period in early systole and additional error in the estimates of (i) PWV at the ascending aorta given by the P U -loop method (9.5% error using rabbit in vivo data) and (ii) TT to a dominant reflection site calculated from the wave intensity profile (27% error using numerical data). These errors decreased to 1.3% and 10%, respectively, when peripheral reflections were considered in the calculations using the Windkessel pressure (Sections 3.8 and 3.9.1).
We have used our new analysis of wave propagation to study the effects of vessel compliance and peripheral resistance on numerically-generated aortic pressure (P ) and flow (Q) waveforms (Section 3.10). With decreasing compliance, the pulse pressure increases, a smaller portion of P is made up of reflected waves originating from earlier cycles, there is less damping of Q, and pulse wave energy increases, suggesting that the LV must produce more energy to propel the same amount of blood flow throughout a stiffer vasculature. With increasing resistance, the mean pressure raises (but not the pulse pressure), a greater portion of P is made up of reflected waves originating from earlier cycles, and there is little change in Q and wave energy. We have also shown that vessel compliance has a similar effect on reflected waves originating from internal junctions, aortic root and tapered vessels to those originating from the periphery, whereas peripheral resistances only affect reflected waves originating from the periphery.
Lastly, our results suggest a small sensitivity of the forward and backward wave intensity profiles to pressure and velocity data sampled at 0.2 kHz or above or to errors in the estimate of PWV withiṅ 20%, assuming data free of any other error.
It is important to note that our modified WIA differs from the reservoir-wave hypothesis [55,56,66], which has been shown not to be beneficial for WIA [21]. In this hypothesis, all the Windkessel pressure is separated from the measured pressure waveform, whereas our new method separates only the Windkessel pressure from previous cardiac cycles. The latter allows us to do WIA on all components of the pressure waveform generated in the current cardiac cycle, starting from a 'true' reflection-free period.

APPENDIX A: DERIVATION OF THE WATER HAMMER EQUATIONS
Following [15], we consider the area (A), velocity (U ) and pressure (P ) waveforms to be made of successive wavefronts, which travel through the arterial network propagating changes in A, U and P . Across an individual wavefront, we define the changes in area, velocity and pressure as dA, dU and dP , respectively. We take the linear assumption that these changes can be separated into changes across the forward-travelling (dA f , dU f , dP f ) and backward-travelling (dA b , dU b , dP b ) wavefronts ( Figure A.1), that is We define the forward direction as the direction of mean blood flow, in which x increases. Applying conservation of mass in the control volume for the frame moving with the forwardtravelling wavefront ( Figure A.1(c) where is the density of blood assumed to be constant and c is the local PWV. This expression reduces to ignoring the term dU f dA f . Applying conservation of momentum to the same moving frame yields 2) and (A.7) allow us to calculate the local PWV (c) at an arbitrary location in the arteries from simultaneous measurements of P and U , A and U , and P and A, respectively, using the methods described in the following two subsections. In all these methods, c is assumed to be constant; that is c D c.x/, and dA, dU and dP are calculated as the change in the measured A, U and P , respectively, in a sampling period dt (e.g. dA.t / D A.t Cdt/ A.t /) using a Savitzky-Golay filter [16,51] [52, p. 650].

From simultaneous measurements of P (or A) and U using 'loop' methods
The 'loop' methods rely on the existence of a period within the cardiac cycle in which the plots of P versus U (P U -loop) and ln A versus U (lnAU -loop) are approximately linear. In normal physiological conditions, this linearity has been observed in early systole in the P U -loop of the human aorta and main pulmonary artery [16,53] and lnAU -loop of the human carotid artery [68], which suggests that dA, dU and dP in early systole are predominantly made of wavefronts propagating towards distal locations. Thus, dA D dA f , dU D dU f and dP D dP f are reasonable assumptions and Equations (8)  They allow us to calculate c from the slope of the (approximately) linear part in the P U -and lnAU -loops; this slope is c and 1=c, respectively. In the case of the P U -loop, we need to know the density of blood ( ). In Section 3.8, we showed that elimination of the Windkessel pressure generated in previous cardiac cycles from P improves the accuracy of the value of c given by the P U -loop method in the ascending aorta. This equation is similar in form to the so-called 'single-point' method [16,69] for simultaneous P and U measurements, though Equation (B.3) follows directly from conservation of mass and momentum and does not require minimising wave energy over the cardiac cycle as does the 'single-point' method.

APPENDIX C: WAVEFRONT REFLECTIONS AT JUNCTIONS
Consider three arterial segments j , j D p, d1, d 2, joining in a splitting or merging flow junction ( Figure C.1). Three perturbations .a j , p j e , q j / of the initial states .A j , P j e , Q j / D .A j 0 , 0, 0/, j D p, d1, d 2, of luminal cross-sectional area, elastic component of pressure and flow rate, respectively, propagating towards the junction (along each corresponding segment j ) will produce a new wave in each segment, denoted by .A j 0 C ıa j , ıp j e , ıq j /, j D p, d1, d 2, which will propagate away from the junction. This appendix shows how these three new waves can be calculated using the linearised 1-D equations, 8 < : C 1D @p e @t C @q @x D 0, L 1D @q @t C @p e @x D R 1D q, p e D a C 1D , where a, p e and q are the perturbation variables for area, elastic component of pressure and flow rate, respectively, that is .A, P e , Q/ D .A 0 C a, p e , q/, and are the elastic wall compliance, fluid inertia and viscous fluid resistance, respectively, per unit length of vessel. Equations (C.1) follow from linearising (1) and (2) where Y j 0 D 1=Z j 0 , j D p, d1, d 2, is the characteristic admittance of segment j . If we perturb one segment at a time with .a j , p j e , q j /, j D p, d1, d 2, so that p i e D 0 for i ¤ j , then from Equation (C.8) we have , j D p, d1, d 2.
( C . 9 ) Defining the reflection coefficient R j f , j D p, d1, d 2, as the ratio of the change of pressure across the reflected wave to the change of pressure in the incident wave; that is R j f Á ıp j e p j e =p j e , j D p, d1, d 2, and using Equation (C.9) to write ıp j e as a function of p j e yields Equation (12). Combining Equations (C.3), (C.4), (C.10) and (C.11) as described previously for the splitting flow case also yields Equation (12). Note that from p e D a C 1D in (C.1), we have p j e D a j =C j 1D and ıp j e D ıa j =C j 1D , j D p, d1, d 2. Thus, R j f D ıa j a j =a j , j D p, d1, d 2, is satisfied, which shows that the || The expressions given for w f,b follow from applying the method of characteristics to the system of equations (C.1). reflection coefficients also give the ratio of the change of area across the reflected wave to the change of area in the incident wave. However, in terms of changes in the flow rates, we have R j f D ıq j q j =q j , j D p, d1, d 2, which follows from transforming ıq j and q j into pressures using Equations (C.6) and (C.7) for the splitting flow case and Equations (C.10) and (C.11) for the merging flow case. Defining the changes in flow velocity as ıu j D ıq j =A j 0 and u j D q j =A j 0 , we have R j f D ıu j u j =u j , j D p, d1, d 2.