Pulsatility Index as a Diagnostic Parameter of Reciprocating Wall Shear Stress Parameters in Physiological Pulsating Waveforms

Arterial wall shear stress (WSS) parameters are widely used for prediction of the initiation and development of atherosclerosis and arterial pathologies. Traditional clinical evaluation of arterial condition relies on correlations of WSS parameters with average flow rate (Q) and heart rate (HR) measurements. We show that for pulsating flow waveforms in a straight tube with flow reversals that lead to significant reciprocating WSS, the measurements of HR and Q are not sufficient for prediction of WSS parameters. Therefore, we suggest adding a third quantity—known as the pulsatility index (PI)—which is defined as the peak-to-peak flow rate amplitude normalized by Q. We examine several pulsating flow waveforms with and without flow reversals using a simulation of a Womersley model in a straight rigid tube and validate the simulations through experimental study using particle image velocimetry (PIV). The results indicate that clinically relevant WSS parameters such as the percentage of negative WSS (P[%]), oscillating shear index (OSI) and the ratio of minimum to maximum shear stress rates (min/max), are better predicted when the PI is used in conjunction with HR and Q. Therefore, we propose to use PI as an additional and essential diagnostic quantity for improved predictability of the reciprocating WSS.


Introduction
Arterial wall shear stresses (WSS) are known as a main regulator of endothelial cell function and vascular structure and condition. Numerous in-vitro and in-vivo studies have shown that vascular regions with disturbed flow accompanied by turbulent flow, low, oscillatory or instantaneous negative WSS and high WSS gradients are strongly correlated with vascular pathologies, cardiovascular diseases and disorders [1][2][3].
WSS parameters were also found as indicators for advanced congestive heart failure, where the combination of reduced cardiac output compensated by increased HR promotes retrograde flow and negative WSS in the aorta over a major part of the cardiac cycle [19].
Clinical evaluation of time-dependent WSS in-vivo is not yet the part of routine cardiovascular exams. This is mainly because it requires accurate measurements of either the direct flow velocity gradient near the wall-using complicated MRI or Doppler ultrasound technique [20][21][22], or the indirect estimates using intrusive pressure gradient measurements. Pressure gradient measurements in conjunction with the classical Womersley solution for pulsatile flow in rigid tubes [23] may provide time-dependent WSS according to: where j is the imaginary unit, D is the tube inner diameter; J 0 and J 1 are the zero-and firstorder Bessel functions of the first kind. The solution is a sum of N+1 harmonics of the basic frequency ω 0 , where k sn are the amplitudes of oscillatory pressure gradient and α n = 0.5D(ωn/v) 1/2 is the Womersley number of the n th harmonic (ω n = ωn), respectively. v is the kinematic viscosity. A third option is to estimate in-vivo time-dependent WSS using the Doppler ultrasound flow waveform measurements and "Inverse Womersley method" [24]. This method allows to extracts the time-dependent WSS from the measured waveform Q(t), by estimating Womersley velocity profiles (u(t,r)) based on the time-dependent flow in a straight tube, and extracting the shear stress from the estimated velocity gradients: where μ is dynamic viscosity. Thus, it should be possible, at least in theory, to estimate the time dependent WSS by applying either straight Womersley method on the time-dependent measured pressure or the inverse Womersley method on the time-dependent measured flow waveform. However, these methods are still too complex to be applied in clinical practice since they require accurate time-resolved measurements. In addition, the accuracy of these methods relies on the validity of Womersley solution, approximating the blood vessel to a straight uniform and rigid tube. The reported discrepancies of WSS measured using the inverse Womersley method as compared to more complex non-linear approximations, are in the range of 23% [25].
Given the above, it is desirable to establish clinical evaluation based on simpler and more accessible measurements that would provide improved estimation of the reciprocating WSS parameters, such as OSI, P[%] or min/max.
Some studies tried to correlate reciprocating WSS parameters with the basic integral parameters of averaged or maximal flow rate (Q) and heart rate (HR), which are the simplest to measure. For example, Finol et al. [14] found that the maximal values of non-dimensional mean WSS and non-dimensional WSS gradient increase with the Reynolds number (non-dimensional representation of the flow rate according to Re = 4Q/πDv). Other researchers [26,27] have shown correlations between increased WSS magnitude and increased Womersley number (non-dimensional representation of the HR according to α = 0.5D(2π Ã HR/60v) 1/2 ) during exercise. Increased HR is considered a major risk factor for cardiovascular disease [27,28], and when medically treated, was found to reduce heart failure [29].
However, as we show in this study, the correlations based only on HR and Q cannot be sufficient due to the nature of the pulsating flow waveforms. In some flow cases, for the same HR (Womersley number) and averaged flow rate (Reynolds number) the flow waveform and associated local WSS parameters may differ drastically due to reverse flow rate during a portion of the cardiac cycle. For example, the aortic wave steepening downstream of the aorta, caused by wave reflection and diameter reduction, results in different flow waveforms at different locations along the artery, and thus lead to different WSS parameters for the same HR and Q [30,31]. Under such conditions, the peak of the maximum flow rate increases but it is compensated over the cardiac cycle by the stronger reverse (negative) flow rate, causing a lower net positive flow rate.
It is important to distinguish between the reverse flow which is local (with negative velocity values only in the near wall region but with the overall positive flow rate), and global reverse flow, when the flow rate at the given cross-section of the artery is negative (Q < 0). In both local and global reverse flow cases there is negative WSS, but on very different time scales [32].
We suggest that the pulsatility index (PI) could be used as a third dimensionless quantity to predict WSS parameters. PI was proposed by Goslin and King [33] to measure the variability of blood velocity in a vessel, and it is defined as the difference between the peak systolic flow and minimum diastolic flow rates divided by the mean flow rate (PI = (Q max − Q min )/Q mean ). PI is readily accessible in-vivo [34]. It has been shown to correlate with presence of arterial stenosis [35], and been recently recognized as an important factor in aneurysmal flows [7,15,18,36] or acute intracerebral hemorrhage [37]. In a different context, the PI is called reverse/forward flow index (measured in percent, see Hashimoto and Ito [38], for instance) and it shows correlations with the pulse pressure amplification and arterial stiffness.
In this study we present in a simulation and verify experimentally that in a straight tube, addition of PI allows better prediction of WSS parameters than Q and HR alone and would potentially improve prediction of negative WSS during the cardiac cycle. To the best of our knowledge, prediction of WSS parameters utilizing PI, as developed in this study was not published elsewhere.
The manuscript is organized as follows. Section 2 presents the methods used in the study, namely the simulation based on the Womersley model and an experimental study using particle image velocimetry (PIV). Section 3 presents the results of the observed correlations and followed by discussion of the results along with the limitations of the present studies and the future perspectives of this research.

Materials and Methods
In order to demonstrate the capabilities of PI in improving prediction of WSS parameters for waveforms with global reverse flow and negative WSS and distinguish those from cases of negative WSS with positive flow rate, we demonstrate in this manuscript the correlation of WSS parameters with Reynolds, Womersley and PI for 14 representative pulsating waveforms (cases). The waveforms were chosen such that they could be realized in our experimental setup and implemented in both simulations and the experiment. Chronologically, the experiments were performed before the simulation. The experimentally accessible waveforms were analyzed to obtain the amplitudes and phases of the flow rate waveforms and simulated using inverse Womersley solution [31]. WSS was estimated using Eq 2.
Based on time-dependent WSS, τ(t), we quantify the WSS-related parameters (using Eq 3), shown to correlate with clinical evidence of vascular dysfunction, namely a) the ratio of min/ max WSS, b) the time interval of the negative WSS during the cycle (P[%]), and c) the measure of overall oscillation (the so-called oscillation shear index, OSI). The first two parameters, min/ max and P[%] were defined by Gharib and Beizaie [19] and were shown to be correlated with the clinical findings. OSI is an index defined by Ku et al. [4] to describe the degree of deviation of the WSS from its average direction. These parameters are defined as follows: where τ min and τ max are the minimal and maximal values of time-dependent WSS in the cycle, t| τ<0 is the total time during which WSS is negative, and T is the time period of the pulsation (inverse of the HR).
A custom-design flow test rig was developed to create controllable pulsatile flows, as shown in Fig 1a. Pulsatile flow waveforms were developed by three DC voltage-driven computer-controlled gear pumps: two pumps connected in parallel drove the flow in the streamwise (forward) direction, and the third pump was installed in the inverse direction to control the flow reversals (Fig 1a). The flow was created in a slightly distensible tube (made of Tygon B-44-4X, Saint Gobain, 80 cm long, inner diameter of 1.9 cm and wall thickness of 0.32 cm, L/d ! 40, elasticity modulus of 12 MPa, estimated distensibility 1.36 × 10 −6 Pa -1 ). The wave speed is approximately 30 m/s and wavelength λ is (using the pulsation period of 1-4 seconds) between 30 and 120 meters. For such a short tube, L/λ ( 1, the effect of distensibility on the flow is negligible [39]. The experimental cases examined had typical physiological characteristics composed of 30-50% "systolic" (forward) part followed by a "diastolic" component. The cases differ by the mean flow rate, presence/absence of the negative WSS and presence/absence of the total flow reversal (i.e. negative flow rate) during a certain part of the cycle. Examined flow waveforms were with relevant parameters at the relevant range for large vessels with Reynolds number ranged between Re max = 70 Ä 750 (associated with maximal peak systolic flow rate), Womersley number α = 6 Ä 13 and PI = 1.5 Ä 9 [40]. The dimensionless parameters of each experimental run are listed in Table 1.
The measurements included pressure transducers at the inlet and outlet of the tube (EW-68075-02, Cole Parmer), flow rate (magnetic flow meter MAG 1100, Danfoss) and velocity distribution using particle image velocimetry (PIV) as shown in Fig 1b. In order to reduce optical reflections and distortion, the tube was placed in a 800 × 300 × 200 mm glass tank filled with a 60%w glycerin-water solution, which refractive index matches to the Tygon tube. Due to power limitations of our pumping system we could not pump the viscous 60%w glycerinwater solution through the tube at the velocity providing physiologically relevant Reynolds numbers. Therefore, the trade-off was set to use for the working fluid inside the tube, a 40%w glycerin-water solution. The mismatch in the index of refraction on only one side of the tube causes a negligible effect on the PIV results (S1 Fig presents a photo of the PIV setup and a zoom-in view into a pipe with 40%w versus 60%w glycerin-water solutions). This Newtonian and homogenous fluid is a commonly used solution to match the blood viscosity [41]. Its dynamic viscosity is 3.72 cP and density is 1.09 g/cm 3 as compared to blood with 3-4 cP and 1.06 g/cm 3 . The small refractive mismatch between the tube and the working fluid resulted with some light reflection from the pipe walls that affect PIV measurement resolution at the near wall region.
We implemented PIV technique using a dual Nd:YAG laser (532 nm, 120 mJ/pulse, Solo 120XT, New Wave Research), a high-resolution CCD camera (12 bit, 4008×2672 pixels, TSI Inc.). The field of view was located 0.6 m from the inlet, downstream of the inlet region of the tube, in order to assure the fully developed flow at the measurement site. The field of view was set to 60 x 40 mm with a magnification ratio of 15 μm/pixel. The PIV analysis was performed on rectangular windows of 64 x 16 pixels. Silver-coated hollow glass spheres (14 μm, 1.05 g/ cm3, TSI Inc.) were used as seeding particles. An analog pressure measurement located at the inlet was used to trigger the PIV synchronizer at fixed phase instants. Seven double-laser pulses and PIV images were taken at fixed time instants during each cycle (at t/T = 0, 0.125, 0.25, 0.375, 0.5, 0.625 and 0.75 equally spaced sections of the normalized period). For each of the seven phases during the cycle, 100 images were acquired, resulted with 50 instantaneous flow realizations. Velocity fields were estimated using a commercial software (Insight3G, TSI Inc.) and verified with an open source software (www.openpiv.net) [42]. First, the 50 flow realizations at each phase were averaged to obtain the phase averaged PIV flow realizations. Second, the velocity was averaged along the streamwise (x) direction, arriving at the phase averaged velocity profiles, u(r,t). The phase averaged velocity profiles are shown in Fig 2 as symbols and colors are used to distinguish between different phases. In some cases, the velocity profiles overlap at different phases, and only a part of those are shown for the sake of clarity. The flow rate waveforms for each of these cases are shown in Fig 3 (dashed lines). The experimental results were analyzed for quantification of the WSS and its relation to the presence of local (near the wall) or global reverse flow (negative flow rate) during a part of each waveform cycle. The most accurate method was the extended version of inverse Womersley solution. The PIV system provided a set of seven time phases along the cycle but with high spatial resolution in the tube cross-section. We applied constrained non-linear least-squares method to find the best fit of Womersley profiles to the phase-averaged measured velocity profiles, according to: j 3=2 a n J 0 ðj 3=2 a n Þ À j 3=2 a n J 0 j 3=2 a n 2r D À Á j 3=2 a n J 0 ðj 3=2 a n Þ À 2J 1 ðj 3=2 a n Þ e jont The non-linear least squares fit provides the amplitudes and phases of flow rate harmonics. Using the amplitudes and phases we evaluated velocity gradients according to Eq 4 and then calculated WSS according to Eq 2 and the respective parameters defined in Eq 3. We then studied their correlations with PI, Re and α. Fig 4a shows an example of the analysis for two example cases: one case with only positive flow throughout the entire cycle (Run 10, left) and the second case with a total flow reversal (Run 12, right). It shows the comparison between the flow rate (Q) waveform obtained directly from integration of PIV velocity profiles in radial direction (symbols with error bars at seven phases) and the time-dependent flow rate obtained from the inverse Womersley method (solid line). Similarly, the WSS estimated derivative of the velocity profiles near the wall as measured by PIV (symbols in Fig 4b) are compared with the WSS obtained from the inverse Womersley method according to Eqs 4 and 2 (solid lines).
The differences emphasize intrinsic limitations of both the direct WSS measurements from PIV data near the wall and the fact that we fit a linear Womersley solution in a rigid tube to the real flow case in a slightly distensible tube. Although it is theoretically possible to include nonlinearities arising from tube's distensibility, curvature, tube-end conditions and nonuniformity using the modified Womersley solution for thin wall elastic tubes [39], this method would significantly complicate the simulations.

Results
Using the simulation and experimental results, we estimate WSS parameters defined in Eq 3 and present them as functions of the Reynolds, Womersley, and pulsating index numbers for each waveform, as listed in Table 1.
The main result of this study is shown in Figs 5-7 and summarized in Table 2. In each figure, panels (a-c) present the WSS parameters of interest, namely OSI, P[%] and min/max, versus each of the three hemodynamic parameters of Reynolds, Womersley and PI, separately. We add to each figure a best-fit curve, the R-square value of which is listed in Table 2.
We limit our fit to the 2 nd order polymonial, for the sake of simplicity and robustness. Higher order polynomials, as well as other types of functions have not shown better fit, though increased complexity and noise amplification during the following, predictive step. As can be seen in Next step after analysing the separate dimensionless quantities, was to find a non-linear fit of two dimensionless quantities together and construct the WSS-related parameters as functions of two arguments, for instance of Reynolds and Womersley numbers (blue squares in panels d). These attempts are compared in  to the models of Re and PI as the arguments (magenta triangles) and also the three-variable model that we focus on in this work (red circles). The models are constructed by non-linear least squares fit of the 2 nd order polynomials in the form ofŶ ¼ c 0 þ c 1 Re þ c 2 Re 2 þ c 3 a þ c 4 a 2 þ c 5 PI þ c 6 PI 2 , where c i are the  Table 2 strengthen the central result of this work. Inclusion of the pulsatility index in the list of significant hemodynamic parameters (in addition to the Reynolds and Womersley numbers) that characterize the Womersley-like flow in straight and rigid pipes provides the best prediction for all three WSS-parameters.

Summary and Discussion
In this study we show that the use of the flow rate (Q) and heart rate (HR) is not sufficient for prediction of WSS parameters, and we suggest to improve the diagnostic evaluation by adding a third dimensionless parameter, the pulsatility index (PI).
Figs 5-7 and Table 2 demonstrate the limitations of non-dimensional Reynolds and Womersley numbers alone and the benefit of adding the PI for better prediction of WSS parameters.
The results show that Re, α and PI separately have poor correlation with the WSS parameters, with the Womersley number lacking any correlation as a single hemodynamic characteristic to the WSS parameters. This is expected, as the definitions of the WSS parameters are mostly normalized to the total cycle period and, therefore, the basic frequency has very little When coming in pairs, the correlation of the dimensionless parameters does not improve dramatically in predicting WSS parameters. Only when the three parameters are used together, the correlations are significantly improved.
Among the three parameters of interest, min/max seems to be the parameter, which is poorly predicted by the hemodynamic quantities. Apparently this is due to the nature of its definition of the ratio of two peak values unrelated to the underlying fluid dynamics. High (or low) peaks can appear in weak (or strong) flows with weak (or strong) pulsations or with high (low) flow rates. Thus, prediction of min/max is problematic and if this parameter is of interest, direct measurement of the peaks is required.
The P[%] and most importantly, the OSI parameters are very well related to the 2 nd order polynomial multi-variable model that involves the three hemodynamic parameters. Large OSI  and P[%] values for large PI values imply that the global flow reversal is a major contributing mechanism to WSS variability during the cycle as well as to the duration of the negative WSS. The results undoubtedly reflect the need for an additional, third, parameter that influences WSS parameters in pulsating flows independently of the HR or Q.
The application seems to be straightforward: by means of sampling flow waveform (e.g. using ultrasound Doppler) at several phases during the cardiac cycle, one can readily obtain PI in addition to the traditional HR and flow rate measurements. Joint analysis of these three values will then provide the necessary prediction of negative WSS and severity of associated WSSrelated parameters.
Addition of PI to the measurement list may improve significantly prognosis and prediction of vascular disease development due to vascular remodeling and variations in flow pulsatility. Future research could implement clinically the above suggestion and investigate correlations of local PI, HR and Q, extracted from Doppler measurements with vascular dysfunctions.  Table 1

Study limitations and further research
The inverse Womersley method in the present study considers a simplified geometry of a straight long circular tube. We do not take into account complicated anatomical geometries, arterial stenoses, arterial curvatures, bifurcations and nonuniform cross sections-which may contribute vortical or helical 3D flow, skewed flow profiles, separation regions, etc. In addition, wall motion may affect WSS parameters in up to 15% [39]. Therefore, the described analytical  method approximates flow profiles and corresponding WSS only for the specific case described. Such simplifications were used in other studies (e.g. [25,43]). Another limitation relates to the working fluid. The working fluid used as a model for blood is not perfectly matched with the Tygon tube. The small refractive mismatch resulted with some light reflection from the pipe walls. Static and dynamic calibrations were implemented to ensure accurate measurement. However, measurement resolution at the near wall region was limited.
Due to limited time resolution of the PIV system, the velocity profiles were measured only at seven time instants during the cycle. It could be beneficial to achieve a better accuracy with the time-resolved PIV systems [42]. Nevertheless, using the high spatial resolution of the PIV method and the inverse Womersley method, we overcame this limitation and succeeded identifying an important improvement in correlation between the WSS and WSS-related parameters, and nondimensional numbers, using the addition of the pulsatility index.
The inverse Womersley method is based on Womersley's simplified model of a straight, rigid and infinitely long tube. It was shown in previous studies that the goodness of fit of the models may be affected by the geometry and boundary conditions [44], Sahtout and Ben Salah [45] have shown that Womersley linear theory is accurate enough for the long rigid vessel with no taper or distensibility. Yet, the present experimental case is neither rigid, nor a straight infinite pipe with a fixed inner diameter. The Tygon tube, clamped at its end points was slightly curved (due to its distensibility), its internal diameter is known only approximately along the length. Moreover, the pipe is elastic and reacts to the pressure pulses by small, yet non-negligible changes in diameter. These local changes are different for different pressure pulses and flow rates. Nonlinearities arising from distensibility and/or tubes nonuniformity may affect the accuracy of the results by up to 18% [25]. Although the model we used for rigid straight tubes has limited accuracy, it is the simplest model that allows for extending our analysis to the full cardiac cycle, as well as to explain in future work cases that we did not examine in this study.
Recent studies show the feasibility of using the inverse Womersley method in conjunction with the Doppler ultrasound measurements to assess the flow rates and the velocity profiles in rigid tubes [43]. Our contribution is the comprehension of the inconsistency in the analysis based on the mean flow rate and the frequency only, and the necessity to include the third parameter, the pulsatility index.