A Fast Method to Estimate the Wall Shear Stress Waveform in Arteries

Abstract

Introduction

Ultrasound has been applied to measure vessel diameter and blood flow velocity to compute the wall shear rate (WSR) in arteries. This paper describes a fast technique to assess the WSR waveform using an image of a pulsed Doppler waveform downloaded from a modern clinical ultrasound scanner.

Methods

A walled vascular phantom has been developed to mimic the physiological condition of brachial arteries, from where measurements were made. A MATLAB program has been developed and used to compute the WSR waveform in a flow phantom from a pulsed Doppler image. The mean WSR obtained from the WSR waveform was compared with the mean WSR derived from the flow rate obtained from a timed collection method. Measurement errors in Doppler velocity estimates from ultrasound scanners were also investigated and used to determine correction factors in WSR calculations.

Results

For three different flow phantom depths, 9.5,14.5 and 19.5 mm, the mean percentage errors between the true and measured WSR were found to be 4.5% (SD = 4.0), 7.4% (SD = 5.1) and 14.2% (SD = 4.1) respectively.

Conclusions

The results demonstrated the feasibility of calculating WSR based solely on an image of the Doppler spectrum and arterial diameter measurement, which opens up the possibility of obtaining WSR estimates from generic scanners.

Keywords

Doppler ultrasound wall shear stress flow phantom

Wall shear stress (WSS) in arteries is the viscous drag of blood on the vessel wall. The WSS is very small (typical value 1–20 Pa) compared with values of the order of 10,000 Pa for blood pressure. However, WSS is important because the endothelium, the layer of cells lining arteries, is responsive to alterations in WSS. Mechanosensors, located on the endothelial cell surface and within the cell, stimulate a variety of physiological effects with different timescales, including release of nitric oxide (seconds), re-alignment of endothelial cells with the WSS vector (hours) and remodelling of the artery (weeks).¹

WSS also has a role in cardiovascular disease. It has long been known that atherosclerosis is initiated at regions of low WSS,² and it has been hypothesized that the high WSS in atherosclerotic plaque leads, via a sequence of events involving local inflammation, to thinning of the fibrous cap, rendering the plaque at risk of rupture.^{3, 4}

Measurement of WSS in advanced disease requires the use of patient-specific modelling (integration of three-dimensional [3D] imaging with computational fluid dynamics) in order to account for the complex haemodynamics.⁵ In early disease, normal vessel geometry is mostly preserved, enabling the use of simplified methods for measurement of WSS, which have mostly involved ultrasound imaging.⁶ The procedure involves measurement of the vessel diameter from B-mode images and measurement of the time-averaged maximum velocity using a spectral Doppler, from which mean wall shear rate (WSR) γ can be calculated using Equation (1):

γ = \frac{2 (2 + n) v}{d}

(1)

where d is the internal arterial diameter, v is the time-averaged mean blood velocity and n represents the shape of the velocity profile. For a fully developed parabolic profile, n = 2; this is the normal assumption when estimating the shear rate. The mean WSS can then be calculated as the product of the viscosity of blood multiplied by the mean WSR.

Estimation of the change in WSR or WSS through the cardiac cycle is more challenging. Most work in this area involves patient-specific modelling, which is resource intensive, requires specialist knowledge of computational modelling and hence is not available for clinical practice. The WSR waveform has been measured using dedicated multigate Doppler ultrasound systems.⁷ However, these are specialist systems not available to the wider user-community. In order to overcome these limitations, Blake et al.⁸ described a methodology for estimation of the WSR and WSS waveform, which used commercial ultrasound technology. This method assumed fully developed flow and used the Womersley equations.⁹ The software required entry of the measured vessel diameter from the B-mode image and of the maximum velocity waveform from the spectral Doppler. In the original implementation, ultrasound data were transferred digitally from an HDI 5000 (Philips Medical Systems, Bothell, WA, USA) to an offline programme, which measured the vessel diameter and extracted the maximum velocity waveform, from which the WSR waveform was calculated. This procedure was optimized for the HDI 5000 scanner, for which the software was developed over many months, requiring detailed knowledge of the HDI file format and the writing of programmes to decode the Doppler data. While Blake et al.⁸ described a simple procedure, the requirement for detailed knowledge of file format data does not make this procedure generalizable. The problem, which this paper attempts to resolve, is to make the procedure for estimation of the WSR/WSS waveform generally applicable.

The methods developed in this paper are relevant to different arteries. However, the work arose from a study on flow-mediated dilation, for which a brief review of the background is given here. Endothelial dysfunction is an early event of atherosclerosis preceding formation of plaques,¹⁰ which impairs vasodilation in response to increases in blood flow associated shear stress. Over the last two decades, many studies have attempted to non-invasively test the endothelial function by investigating the shear stress on the vascular endothelium that causes endothelium-dependent flow-mediated vasodilation (FMD).^{10, 11} After a five-minute occlusion of the brachial artery with a cuff, removal of the cuff leads to an increase of the WSS, stimulating the release of nitric oxide, a vasodilator, from the endothelial cells into the smooth muscle.

Methods

Calculation of the WSR waveform

The software developed by Blake et al.⁸ was modified to allow extraction of the ultrasound data from image frames taken from any ultrasound scanner, provided that the B-mode image of the vessel and spectral Doppler waveform are shown.

The WSS is the force exerted on the artery wall by blood moving past it, which is given by

τ_{ω} = - μ {\frac{\partial v}{\partial r} |}_{r = R}

(2)

where τ_w is the WSS, μ the dynamic viscosity and $\frac{\partial v}{\partial r} |_{r = R}$ the velocity gradient, or WSR at the artery wall. The estimation of WSS from Equation (2) relies on the measurement of the velocity gradient and knowledge of the blood viscosity.

The velocity profile can be obtained from the Womersley equation as below.⁹ For the k = 0 component, which corresponds to steady flow, the velocity profile is parabolic; for all other components, for which k > 1, Equation (3) gives the velocity profile:

v (y, ​ t) = \sum_{k = 0}^{\infty} Re {V_{mean, k} e^{j (k ω t - Ψ k)} [\frac{J_{0} (τ_{k}) - J_{0} (τ_{k} y)}{J_{0} (τ_{k}) - 2 J_{1} (τ_{k}) / τ_{k}}]}

(3)

where Re{.} represents the real part of a complex function, J₀ and J₁ are the zeroth and first-order Bessel functions of the first kind, j is the imaginary number, y is the normalized radial coordinate, ω is the angular frequency, t is the time and Ψ represents the phase of each harmonic. For the condition k = 0, which corresponds to steady flow, the velocity profile is parabolic. The component τ_k is given by

τ_{k} = R \sqrt{\frac{k ω}{v_{k}}} j^{3 / 2}

(4)

where υ_k is the kinematic viscosity of the fluid and R is the vessel radius.

In order to make use of the maximum velocity waveform, as measured using Doppler ultrasound, the mean velocity in Equation (3) is expressed in terms of maximum (centre-line) velocity, as shown in Equation (5) (see Reference¹²):

V_{mean, k} e^{- j Ψ k} = V_{centre, k} e^{- j ϕ k} [\frac{J_{0} (τ_{k}) - 2 J_{1} (τ_{k}) / τ_{k}}{J_{0} (τ_{k}) - 1}]

(5)

Substituting Equation (5) into Equation (3) then gives an expression relating the centre line velocity to the velocity profile:

v (y, ​ t) = \sum_{k = 0}^{\infty} Re {V_{centre, k} e^{j (k ω t - ϕ k)} [\frac{J_{0} (τ_{k}) - J_{0} (τ_{k} y)}{J_{0} (τ_{k}) - 1}]}

(6)

The WSR γ_w′, is then calculated using the expression below

γ_{w} = {\frac{\partial v}{\partial r} |}_{r = R}

(7)

Ultimately, the WSS can be derived from the product of blood viscosity and γ_ω. This method may be used to estimate the variation in WSR as a function of time during the cardiac cycle, from which maximum and mean WSR can be computed. The centreline velocity in Equation (6) can be measured by using any commercial clinical ultrasound scanner.

Common file format types used as output for ultrasound systems include DICOM and tiff. For analysis of the duplex image, the Blake et al.⁸ software was modified to allow isolation of the region containing the B-mode image and conversion of this into a matrix array. The vessel centre was identified by placement of a cursor somewhere in the vessel and the diameter measured automatically as described in Blake et al.⁸ The image required entry of scaling markers for the purposes of distance calibration. The spectral Doppler waveforms were similarly isolated and converted to a matrix. The modified software identified the baseline corresponding to zero Doppler frequency shift. The user interactively defined the time and velocity scale, as well as the region of interest on the Doppler trace. To calculate an ensemble average from the corrected centreline velocity waveform, each individual cycle must first be isolated, aligned and averaged. The beginning of each cycle was determined using a gradient threshold. The resulting cycles were aligned by cross-correlating each cycle with the first and shifted in time by the resulting correlation peak. The mean of the aligned signals was calculated; the length of the average cycle was taken as the median of the lengths of the individual cycles. Following correction for geometric spectral broadening (see below) the ensemble averaged peak velocity waveform and the arterial diameter were used to generate the Womersley velocity profiles from which WSR was estimated according to the procedure described by Blake et al.⁸ The WSS waveform was then calculated from the WSR waveform. Data were saved to a disc as a delimited file for uploading to spreadsheet software such as Excel (Microsoft, Seattle, WA, USA) for further analysis and display.

Preparation of the flow phantom

The original methodology as described by Blake et al.⁸ was extensively validated using flow phantoms. However, as modifications to the original software were made, some validation using flow phantoms was performed as part of this study. This was done with respect to the brachial artery, due to the current clinical interest in FMD, a recognized marker of endothelial function and cardiovascular risk.¹³

A walled flow phantom was constructed to mimic the brachial artery. C-flex tubing (Cole Palmer, London, UK) was used with inner and outer diameters of 0.094 and 0.156 inches respectively. The C-flex tubing was mounted within a Perspex box and a line was marked on the wall of the box above the centreline of the tubing at distances of 9.5,14.5 and 19.5 mm respectively, which was representative of the depth of the brachial artery in the patients. Tissue mimicking material (TMM) was prepared, based on an agar-based formulation described by Teirlinck et al.,¹⁴ and poured into the Perspex box until it reached the line marked on the wall, and then allowed to set. Once set, a small layer of 9% glycerol in water solution was poured onto the surface of the TMM and sealed using cling-film to prevent drying out. The phantom was connected into a flow loop driven by a gear pump (Micropump Model 132; Michael Smith Engineers Ltd, Surrey, UK) powered by an electric motor and an external controller. A plain text programme was written to be loaded into LabVIEW (National Instruments Corporation Ltd, Newbury, UK) to drive the gear pump to output a flow with similar patterns to those of the brachial artery.⁸ Blood-mimicking fluid (BMF) was prepared based on the recipe described by Ramnarine et al.¹⁵ and the flow phantom was primed with the BMF.

Acquisition and processing of data from the flow phantom

The flow phantom was set to give pulsatile flow using a brachial waveform (Figure 1). Ultrasound data were acquired using an HDI 5000 (ATL Ultrasound, Bothell, WA, USA). The machine settings were optimized for vascular applications using a linear array transducer (L12-5 38, central B-mode frequency 12 MHz, central Doppler frequency 5 MHz). The gain, brightness and contrast were set as recommended in the clinical application protocol. The sample volume was positioned at the centre of the lumen to ensure that the centreline velocity was measured. The angle correction cursor was aligned with the vessel wall. The ultrasound data were saved as a DICOM image and transferred to the PC for analysis as described above.

Figure 1

Spectral Doppler trace (a) and its WSR waveform (b) of the vascular phantom

Evaluation of the estimated mean WSS requires knowledge of the true mean WSR. This was calculated from flow rate (Equation (8)) measured using timed collection with a measuring cylinder and stopwatch. For this purpose the flow rate Q was set to 70 mL s⁻¹, which is within the normal physiological range in the brachial artery.

γ_{mean} = \frac{32 Q}{π d^{3}}

(8)

where d is the diameter of the C-flex tubing and Q is the flow rate.

The true mean WSR was compared with the time-averaged WSR obtained from the WSR waveform. It is noted that both the Womersley method and Equation (8) assume fully developed flow.

Correction factor

Under standard conditions, the maximum velocity is overestimated in many clinical scanners by 0–29%,¹⁶ due to geometric spectral broadening. A BBS string phantom with an O-ring rubber filament (BBS Medical Electronics, Hagersten, Sweden) was used to investigate this source of error and to establish a correction factor for the machine used in the present study. The string phantom was contained within a tank filled with 9.0% glycerol in water solution by volume, which has an acoustic velocity of 1540 m s⁻¹ at 20°C. The indication speed on the phantom controller was set to 57 cm s⁻¹. Two short sections of thin nylon filament were wrapped around the O-ring rubber separated by 6 cm. The nylon filaments produced brief spikes on the Doppler spectrum. The true filament speed was calculated by dividing 6 cm by the time separation of the two spikes, as measured using the ultrasound system on-board measurement package. The Doppler-estimated maximum velocity was compared with the true velocity. A look-up table of correction factors was created by varying the target depth and angle between the ultrasound beam and the moving filament.

Ultrasound data acquisition from patients

Ultrasound B-mode and spectral Doppler waveforms were acquired from a patient undergoing FMD in order to demonstrate the feasibility of the technique in patients. Ethical permission was obtained and the patient gave written informed consent. A screen shot showing the ultrasound images was captured and stored in DICOM format, which was then converted to BMP format for analysis.

Results

The correction factors for Doppler velocity measurements were found to be 1.2, 1.1 and 1.0 at depths of 9.5, 14.5 and 19.5 mm respectively, with a beam-target angle of 70°.

A velocity waveform, mimicking blood flow in the brachial artery, was generated by the gear pump under LabVIEW control. The spectral Doppler traces, saved as an image file, were transferred to the MATLAB script for the WSR computation. Examples of the spectral Doppler trace and the WSR waveform are given in Figure 1.

Figure 2 shows comparisons between the MATLAB calculations and mean WSRs obtained from timed collection measurements for depths of 9.5, 14.5 and 19.5 mm. The mean percentage errors are 4.5% (SD = 4.0), 7.4% (SD = 5.1) and 14.2% (SD = 4.1), respectively. Figure 3 gives the residual values of the original data. The phantom with the depth of 14.5 mm gave the smallest residual values, the other two phantoms showed similar residual values.

Figure 2

Comparison of the measured and computed mean WSR in a vascular phantom, for (a) depth = 9.5 mm, (b) depth = 14.5 mm and (c) depth = 19.5 mm. The mean percentage errors between the true and measured WSR were found to be (a) 4.5% (SD = 4.0), (b) 7.4% (SD = 5.1) and (c) 14.2% (SD = 4.1) respectively

Figure 3

The residuals of the MatLAB calculation of WSR at: (a) depth = 9.5 mm, (b) depth = 14.5 mm and (c) depth = 19.5 mm

Discussion

Errors in mean WSR at 4.5–14% are similar to those found by Blake et al.⁸ which were 9% for the brachial artery phantom. This gives confidence that the new technique operates with a similar level of accuracy to the original method described by Blake et al.⁸ The new method, which can accept B-mode and spectral Doppler data in file formats which are universally available, is now much more generalizable than the original Blake et al.⁸ method, which used proprietary file formats from Philips Medical Systems. This may enable future clinical work to be performed with data from any ultrasound scanner.

The measured errors in flow rate are considered further. Maximum velocity measurements using Doppler ultrasound may be subject to very large errors (up to 60%), which the correction method above attempts to reduce. We have not undertaken an evaluation of the residual error following correction using this method. However, this approach is similar to angle correction from the edge of the Doppler aperture, for which the error is 2%.¹⁷ For diameter measurement, provided that the interface between fluid and vessel wall is clearly seen, the diameter is estimated in an unbiased manner from the B-mode image, with an error of 2–3% from the cursor placement and pixel size. However, if the interface cannot be clearly visualized, diameter is underestimated by typically half the pulse length,¹⁸ equivalent to 10–15% error. WSR in its simplest form is related to velocity divided by diameter. Expected errors in WSR are therefore in the range of 4% for best-case conditions to 10–15% for worst-case conditions. This corresponds with the errors found in this study of 4–14%, and as noted above, agrees with errors found in the study by Blake et al.⁸

The methods described in this paper assume fully developed flow. For steady flow this means that the velocity profile is parabolic. For pulsatile flow, this means that the velocity profile, at a certain time during the cardiac cycle, does not change with distance along the vessel. During pulsatile flow, the velocity profile is not parabolic and changes through the cardiac cycle. However, the time-averaged velocity profile (averaged over the cardiac cycle) is parabolic. This justified the use of simplified equations such as Equation (2), which can be derived by assuming a parabolic velocity profile. These issues are explored in more detail in the article by Hoskins.¹⁹

The current study was based on the use of a single ultrasound scanner, a Philips HDI 5000. In future studies, it would be worth performing the WSR calculations using different ultrasound systems. In the current study, only three phantom vessels with different depths, which are close to the physiological range of depths relevant to the brachial artery, were used. In future studies, a larger range of depths and diameters relevant to other arteries could be explored.

The current study included just one example from a patient for the purposes of demonstration of feasibility in vivo. It is noted that a larger data-set involving FMD is the subject of a separate study.

Conclusion

A method for the calculation of the WSR from an image of a Doppler spectral waveform has been developed. The method has been tested using a vascular phantom. It should be possible to apply this technique to other clinical ultrasound scanners to investigate cardiovascular conditions, for instance, WSR alternation in FMDs.

Declarations

Competing interests: None.

Funding: None.

Ethical approval: Not applicable.

Guarantor: XY.

Contributorship: XY and PH researched literature and conceived the study. FA and FK were involved in protocol development, gaining ethical approval. LH was involved in experiment and data analysis. All authors reviewed and edited the manuscript and approved the final version of the manuscript.

Footnotes

Acknowledgements

None.

References

Davies

. Flow-mediated endothelial mechanotransduction. Physiol Rev 1995; 75: 519–60.

Caro

, Fitz-Gerald

, Schroter

. Atheroma and arterial wall shear. Observation, correlation and proposal of a shear dependent mass transfer mechanism for atherogenesis. Proc R Soc Land B 1971; 177: 109–59.

Slager

, Wentzel

, Gijsen

, Schuurbiers

, van der Wal

, van der Steen

. The role of shear stress in the generation of rupture-prone vulnerable plaques. Nat Clin Pract Cardiovasc Med 2005; 2: 401–7.

Slager

, Wentzel

, Gijsen

, Thury

, van der Wal

, Schaar

. The role of shear stress in the destabilization of vulnerable plaques and related therapeutic implications. Nat Clin Pract Cardiovasc Med 2005; 2: 456–64.

Hoskins

, Hardman

. Three-dimensional imaging and computational modelling for estimation of wall stresses in arteries. Br J Radiol 2009; 82(Spec No 1):S3–17

Stoner

, Sabatier

. Use of ultrasound for non-invasive assessment of flow-mediated dilation. J Atheroscler Thromb 2012; 19: 407–21.

Brands

, Hoeks

, Hofstra

, Reneman

. A noninvasive method to estimate wall shear rate using ultrasound. Ultrasound Med Biol 1995; 21: 171–85.

Blake

, Meagher

, Fraser

, Easson

, Hoskins

. A method to estimate wall shear rate with a clinical ultrasound scanner. Ultrasound Med Biol 2008; 34: 760–74.

Womersley

. Method for the calculation of velocity, rate of flow and viscous drag in arteries when the pressure gradient is known. J Physiol 1955; 127: 553–63.

10.

Celermajer

, Sorensen

, Gooch

. Non-invasive detection of endothelial dysfunction in children and adults at risk of atherosclerosis. Lancet 1992; 340: 1111–5.

11.

Pyke

, Tschakovsky

. Peak vs. total reactive hyperemia: which determines the magnitude of flow-mediated dilation? J Appl Physiol 2007; 102: 1510–19.

12.

Holdsworth

, Norley

, Frayne

, Steinman

, Rutt

. Characterization of common carotid artery blood-flow waveforms in normal human subjects. Physiol Meas 1999; 20: 219–40.

13.

Khan

, Galarraga

, Belch

. The role of endothelial function and its assessment in rheumatoid arthritis. Nat Rev Rheumatol 2010; 6: 253–61.

14.

Teirlinck

, Bezemer

, Kollmann

. Development of an example flow test object and comparison of five of these test objects, constructed in various laboratories. Ultrasonics 1998; 36: 653–60.

15.

Ramnarine

, Nassiri

, Hoskins

, Lubbers

. Validation of a new blood-mimicking fluid for use in doppler flow test objects. Ultrasound Med Biol 1998; 24: 451–9.

16.

Hoskins

. Accuracy of maximum velocity estimates made using Doppler ultrasound systems. Br J Radiol 1996; 69: 172–7.

17.

Hoskins

. A comparison of single and dual beam methods for maximum velocity estimation. Ultrasound Med Biol 1999; 25: 583–92.

18.

, McDicken

, Hoskins

. Blood vessel diameter measurement by ultrasound. Physiol Meas 1993; 14: 291–7.

19.

Hoskins

. A review of the measurement of blood velocity and related quantities using Doppler ultrasound. Proc Inst Mech Eng H 1999; 213: 391–400.