## Abstract

Arterial pulse waves contain clinically useful information about cardiac performance, arterial stiffness and vessel tone. Here we describe a novel method for non-invasively assessing wave properties, based on measuring changes in blood flow velocity and arterial wall diameter during the cardiac cycle. Velocity and diameter were determined by tracking speckles in successive B-mode images acquired with an ultrafast scanner and plane-wave transmission. Blood speckle was separated from tissue by singular value decomposition and processed to correct biases in ultrasound imaging velocimetry. Results obtained in the rabbit aorta were compared with a conventional analysis based on blood velocity and pressure, employing measurements obtained with a clinical intra-arterial catheter system. This system had a poorer frequency response and greater lags but the pattern of net forward-traveling and backward-traveling waves was consistent between the two methods. Errors in wave speed were also similar in magnitude, and comparable reductions in wave intensity and delays in wave arrival were detected during ventricular dysfunction. The non-invasive method was applied to the carotid artery of a healthy human participant and gave a wave speed and patterns of wave intensity consistent with earlier measurements. The new system may have clinical utility in screening for heart failure.

## Key Words

## Introduction

Ventricular ejection causes a wave of increased blood pressure, blood velocity and vessel diameter that propagates along the systemic arteries (

Nichols and O'Rourke, 1998

). When contraction slows and then reverses, a wave of decreased pressure, velocity and diameter propagates in the same direction (Parker et al., 1988

). These waves partially reflect and re-reflect where vessel geometry or structure changes (Segers and Verdonck, 2000

). They carry information concerning cardiac and vascular properties: wave speed (often termed *pulse wave velocity*[PWV]) is related to arterial stiffness (), wave reflections can be altered when vessel tone changes (Weinberg et al., 2001

; Wang et al., 2013

) and wave intensity depends, at least in part, on cardiac performance (Curtis et al., 2007

; Sugawara et al., 2009

). Here we focus on the latter.Curtis et al., 2007

found that the energy of the first wave was markedly reduced in stable compensated systolic heart failure and that the reduction increased with increasing severity on the New York Heart Association scale. The second wave was unaffected. Sugawara et al., 2009

similarly found that the first wave was reduced by approximately 50% in patients with dilated cardiomyopathy (*i.e.,*impaired cardiac contractility), without significant change in the second wave. They also found that the second wave was reduced by approximately 50% in patients with hypertrophic cardiomyopathy (

*i.e.,*impaired diastolic function), without significant change in the first.

Both studies employed the methods of analysis developed by

Parker and Jones, 1990

in which wave intensity, $dI$, is calculated as the product of $dP$ and $dU$—respectively the change in blood pressure ($P$) and velocity ($U$) over a short interval in the artery of interest. If the wave speed is known (and it can be derived from $dP$ and $dU$), the waves can additionally be divided into their forward-traveling and backward-traveling components (Parker, 2009

). Measurements of $P$ and $U$ should be temporally and spatially coincident, and they need to be made with high temporal resolution to capture the rapid changes in wave intensity that occur over the cardiac cycle.These requirements, and particularly the need for frequent pressure measurements, are problematic.

Curtis et al., 2007

used Doppler ultrasound to obtain *U*and tonometry to obtain*P*, both in the carotid artery. Tonometry requires complex calibration to give true pressures, and the measurements of $U$ and $P$ cannot be made simultaneously.Sugawara et al., 2009

used Doppler ultrasound to measure $U$ and M-mode ultrasound to measure arterial diameter $D$, and assumed that $D$ is proportional to $P$. This ignores the well-established non-linearity of arterial stress–strain curves. A technique that avoids these problems is to measure $P$ and $U$ with an intra-arterial catheter that has both a pressure transducer and a Doppler probe at its tip (Davies et al., 2006

). However, the invasiveness of the method limits its clinical utility.Feng and Khir, 2010

developed a formulation in which wave intensity can be derived directly from $U$ and $D$. A subsequent formulation by Biglino et al., 2012a

employs flow rate $Q$ and cross-sectional area $A$. The intensities are not identical to those obtained from $U$ and $P$—even the dimensions are different—but these systems are internally self-consistent. We (Reavette et al., 2020

) and others (Kang et al., 2019

) have determined by numerical modeling that the formulation of Feng and Khir should have clinical utility equivalent to that of the more established method.The advantage of these new systems is that $U$ and $D$, or $Q$ and $A$, can be obtained by non-invasive techniques.

Pomella et al., 2017

, Pomella et al., 2018

) and Kowalski et al., 2019

used Doppler ultrasound to obtain $U$ and B-mode or M-mode ultrasound to obtain $D$, Li et al., 2019

used magnetic resonance imaging (MRI) to obtain $U$ and $D$, and Biglino et al., 2012a

, Biglino et al., 2012b

) used MRI to obtain $Q$ and $A$. Useful wave intensity data were obtained in each case, but limitations are the well-established inaccuracies in Doppler velocity measurements (- Biglino G
- Schievano S
- Steeden JA
- Ntsinjana H
- Baker C
- Khambadkone S
- de Leval MR
- Hsia TY
- Taylor AM
- Giardini A

Modeling of Congenital Hearts Alliance (MOCHA) Collaborative Group. Reduced ascending aorta distensibility relates to adverse ventricular mechanics in patients with hypoplastic left heart syndrome: Noninvasive study using wave intensity analysis.

*J Thorac Cardiovasc Surg.*2012; 144: 1307-1313

Hoskins, 1999

), the difference in optimal beam angles for Doppler and B-mode or M-mode imaging and the high cost and long acquisition time of MRI.In the present work, we describe a non-invasive method in which $D$ and $U$ are both obtained from the same B-mode ultrasound images; singular value decomposition (SVD) is used to separate the weak blood and strong tissue signals, cross-correlation between images to track the moving wall and blood speckles and an ultrafast scanner with plane-wave transmission to adequately resolve the rapid acceleration and deceleration of blood during systole (

Riemer et al., 2022

). The method is validated by comparison with data obtained using the invasive catheter-based method and is found to be capable of detecting ventricular dysfunction in rabbits. Preliminary data indicate that the technique can be translated to human participants.## Methods

### Animal preparation

All experiments complied with the Animals (Scientific Procedures) Act 1986 and were approved by the Animal Welfare and Ethical Review Body of Imperial College London. Twelve New Zealand White rabbits (HSDIF strain, Envigo, Huntingdon, UK; age 81 ± 13 d, weight 2.39 ± 0.38 kg, mean ± SD) were maintained on a 12:12 h light:dark cycle at 18°C and fed a normal laboratory diet and tap water

*ad libitum*. Each animal was pre-medicated with acepromazine (0.5 mg/kg i.m.*,*CVet, London, UK), anaesthetized with fentanyl fluanisone (Hypnorm, 0.3 mL/kg i.m., Janssen, High Wycombe, UK) and midazolam (Hypnovel, 0.1 mL/kg i.v., Roche, Welwyn Garden City, UK), tracheotomized and ventilated at 40 breaths/min with 45 cm H_{2}O peak inflation pressure (Harvard Small Animal Ventilator, Harvard Apparatus, Cambourne, UK). Anaesthesia was maintained with fentanyl fluanisone (0.1 mL/kg) and midazolam (0.1 mL/kg) approximately every 45 min. Body temperature was monitored using a rectal probe and maintained with a heated pad. Blood oxygen saturation was monitored by pulse oximetry.### Ultrasound imaging

A Vantage 64 LE ultrasound research platform (Verasonics, Kirkland, WA, USA) equipped with an L11-4v probe was used to collect high-frame-rate images of the abdominal aorta distal to the renal arteries. The rabbit was tipped from its supine position to elevate its left side, and the probe was angled in the opposite direction and clamped to give a stable, approximately coronal view through the aortic centerline. Sixty-four elements located at the center of the probe were used to transmit and receive, giving a lateral field of view of 19.2 mm, along the axis of the aorta.

B-Mode imaging was performed using a coherent plane wave compounding scheme (

Montaldo et al., 2009

; Riemer et al., 2020

). A pulse sequence consisting of three plane waves (8 MHz, 1 cycle), steered from –5° to 5° in three 5° steps, was transmitted at a pulse repetition frequency of 47.6 kHz for an imaging depth of 20 mm. Successive pulse sequences were separated to give an effective frame rate of 1 kHz after coherent summation of the three angled plane waves. Images were acquired for 2 s. The mechanical index (MI) was 0.1 MPa. Backscattered radiofrequency (RF) signals were sampled at four times the transmit frequency; echoes of each transmission were recorded for post-processing offline.A three-lead electrocardiogram (ECG) signal was recorded concurrently via a PowerLab 26T data acquisition system (AD instruments, Oxford, UK) connected to the Verasonics host computer running LabChart software. An analog trigger signal sent at the start of the image acquisition sequence from the Verasonics to the PowerLab enabled later time alignment of the arterial waveforms and ECG trace.

### Intra-arterial measurements

The proximal right femoral artery was exposed, lidocaine was administered topically (typically three sites, <1 mL total) and the femoral nerve was severed.

A Volcano ComboWire XT wire with pressure and Doppler sensors at the tip, connected to a ComboMap system (Phillips-Volcano, San Diego, CA, USA), was inserted into the artery and advanced into the abdominal aorta, just distal to the imaging site. The wire was rotated until the strongest Doppler signal was obtained. The wall filter was set at 400 Hz and the signal-to-noise threshold was adjusted manually to optimise the real-time Doppler envelope tracking.

The analog pressure and velocity output ports were connected to two analog inputs on the PowerLab data acquisition system (1 V = 100 mm Hg for pressure, 0.5 m/s for velocity, sample rate 1 kHz). The Doppler sampling rate on the ComboMap system is 120 Hz.

### Experimental protocol

At least three data sets were acquired with the image-based (Verasonics) and catheter-based (Volcano) systems. The catheter-mounted flow sensor was disconnected during imaging as it caused interference, and the two techniques were therefore used alternately. Catheter-derived pressure data could, however, be recorded during the ultrasound imaging.

Following these measurements, the effects of ventricular dysfunction were assessed. Esmolol was administered to 10 rabbits; each received three boluses of increasing dose (1.5, 3.0 and 6.0 mg/kg, i.v.) at 10-min intervals. Successive data sets were acquired with both systems every minute, starting just before the first dose. Two rabbits were given an equivalent volume of saline intravenously as a control.

### Post-processing of invasive pressure data

All post-processing was performed in MATLAB R2020b (The MathWorks, Natick, MA, USA).

Signal processing within the Volcano Combowire pressure-measuring system introduces a delay of 20 ms. This was determined in an elastic tube model

*in vitro*by aligning the peak pressure recorded using the Volcano system with that simultaneously recorded using a high-fidelity, zero-delay catheter system (Fig. S1, online only).Even after peak alignment, the foot of the Volcano-derived pressure waveform was less acute and occurred ∼5 ms earlier than the foot in the high-fidelity, catheter-derived waveform when both were employed in the

*in vitro*system (Fig. S1). This difference was also observed when comparing Volcano-derived pressure with diameter measured by the non-invasive ultrasound method. Analysis of the power spectra revealed low-order, frequency-based filtering in the Volcano system. A filter with comparable characteristics was constructed; it was applied to measurements of $D$ and $U$ when processing the data with*PU*and ln(*D*)*P*methods. This adjustment makes data from the different sources directly comparable and avoids anomalous behaviour such as non-linear slopes of the ln(*D*)*P*and*PU*loops in early systole; it does not have a substantial effect on mean wave speeds or intensities (see,*e.g.,*Fig. S2, online only).*In vivo*, there was an additional time delay because of displacement of the Combowire tip from the center of the B-mode field of view. The delay varied between rabbits (from 0.5 to 5 ms) because of differences in both the displacement and the wave speed. This rabbit-specific delay and the system-specific 20-ms lag described previously were corrected together by assuming that reflections and viscoelastic effects dissipate to negligible levels prior to systolic ejection, meaning that $P$ and $D$ should initially rise together. Following

Swalen and Khir, 2009

, the linearity of the systolic portion of the ln(*D*)

*P*loop was optimized by applying arbitrary lags to $P$ and finding the lag that gave maximum correlation. ($D$ had the aforementioned frequency filtering applied first.)

### Post-processing of non-invasive ultrasound data

Parts of the following are based on earlier studies and validations by

Riemer et al., 2020

, Riemer et al., 2022

). The image filtering and particle tracking methods are illustrated in Fig. S3, online only.#### Beamforming and compounding

The RF channel data were reconstructed using delay-and-sum beamforming (assuming a 1540 m/s speed of sound for delay calculations) and then Hilbert transformed. Echoes from the three angled plane waves were coherently summed, resulting in a final ensemble image size of

*Nx*= 128,*Nz*= 416 and*Nt*= 2000 (lateral pixel size = 0.15 mm, axial pixel size = 0.048 mm).#### Singular value decomposition filtering

Singular value decomposition (SVD) filtering was used to separate the blood signal from the tissue signal and noise (

where ${\lambda}_{i}$ are the ordered singular values of $S$, and ${u}_{i}$ and ${v}_{i}$ are the corresponding spatial and temporal singular vectors. Tissue, blood and noise have different spatiotemporal characteristics—the first (and largest) singular values typically correspond to tissue, the next to blood and the smallest to noise. The decomposition can be written as

where $\mathrm{T}{\mathrm{h}}_{1}$ and $\mathrm{T}{\mathrm{h}}_{2}$ are singular value thresholds. These were selected manually, based on the visual appearance of ${u}_{i}$ and the frequency of the temporal vectors ${v}_{i}$ (

Demené et al., 2015

). The image was cropped so that its axial dimension was slightly larger than the vessel diameter. Then the ensemble image $S$ was reshaped into a 2-D spatiotemporal representation (size *Nx*·*Nz*×*Nt*) and decomposed as$S\left(xz,t\right)\phantom{\rule{0.33em}{0ex}}=\phantom{\rule{0.33em}{0ex}}\sum _{i=1}^{{N}_{t}}{\lambda}_{i}\phantom{\rule{0.33em}{0ex}}{u}_{i}\left(xz\right)\phantom{\rule{0.33em}{0ex}}{v}_{i}{\left(t\right)}^{\mathrm{T}}$

(1)

where ${\lambda}_{i}$ are the ordered singular values of $S$, and ${u}_{i}$ and ${v}_{i}$ are the corresponding spatial and temporal singular vectors. Tissue, blood and noise have different spatiotemporal characteristics—the first (and largest) singular values typically correspond to tissue, the next to blood and the smallest to noise. The decomposition can be written as

$S\left(xz,t\right)={S}_{\text{tissue}}+{S}_{\text{blood}}+{S}_{\text{noise}}$

(2)

$\begin{array}{c}\\ S\left(xz,t\right)=\phantom{\rule{0.33em}{0ex}}\sum _{i=1}^{\mathrm{T}{\mathrm{h}}_{1}}{\lambda}_{i}\phantom{\rule{0.33em}{0ex}}{u}_{i}\left(xz\right)\phantom{\rule{0.33em}{0ex}}{v}_{i}{\left(t\right)}^{T}\phantom{\rule{0.33em}{0ex}}+\sum _{i=\mathrm{T}{\mathrm{h}}_{1}+1}^{\mathrm{T}{\mathrm{h}}_{2}}{\lambda}_{i}\phantom{\rule{0.33em}{0ex}}{u}_{i}\left(xz\right)\phantom{\rule{0.33em}{0ex}}{v}_{i}{\left(t\right)}^{T}\\ \end{array}\phantom{\rule{0.33em}{0ex}}\phantom{\rule{0.33em}{0ex}}+\sum _{i=T{h}_{2}+1}^{{N}_{t}}{\lambda}_{i}\phantom{\rule{0.33em}{0ex}}{u}_{i}\left(xz\right)\phantom{\rule{0.33em}{0ex}}{v}_{i}{\left(t\right)}^{T}$

(3)

where $\mathrm{T}{\mathrm{h}}_{1}$ and $\mathrm{T}{\mathrm{h}}_{2}$ are singular value thresholds. These were selected manually, based on the visual appearance of ${u}_{i}$ and the frequency of the temporal vectors ${v}_{i}$ (

Song et al., 2017

). Figure S4 (online only) illustrates how the thresholds were selected and that their precise value is not critical to the derived velocity waveforms. Finally, images were envelope detected and re-sampled using MATLAB's griddedInterpolant() function to obtain isotropic pixel dimensions.#### Velocity measurement

Velocity was measured by tracking the movement of blood speckles between frames in ${S}_{\text{blood}}$. Briefly, an “interrogation window” was defined in the vessel lumen in one image, and the cross-correlation of its speckle pattern with the pattern in a larger window in the next image of the sequence was calculated. The highest correlation identified the most likely location to which the group of scatterers had moved. Velocities were obtained from the displacement and the imaging frame rate.

Correlation was Fourier-based for computational efficiency. Consider a template window ${a}^{\left(k\right)}$ and search window ${b}^{\left(k+1\right)}$ where $k$ is the frame number. Cross-correlation in the spatial domain is equivalent to convolution in the frequency domain eqn (4)

where $A$ and $B$ are the Fourier transforms of the interrogation windows, $\delta x$ and $\delta z$ are pixel shifts in the lateral and axial directions, * is the complex conjugate and ${\mathcal{F}}^{-1}$ is the inverse Fourier transform. Conventionally, the displacement ($\Delta x,\Delta z$) is given by the location of the single Gaussian-shaped peak in the correlation function $r$ (arg max(|

$r\left(\delta x,\delta z\right)=\phantom{\rule{0.33em}{0ex}}{\mathcal{F}}^{-1}\left\{{A}^{\left(k\right)}\left(\delta x,\delta z\right)\xb7{B}^{\left(k+1\right)*}\left(\delta x,\delta z\right)\right\}$

(4)

where $A$ and $B$ are the Fourier transforms of the interrogation windows, $\delta x$ and $\delta z$ are pixel shifts in the lateral and axial directions, * is the complex conjugate and ${\mathcal{F}}^{-1}$ is the inverse Fourier transform. Conventionally, the displacement ($\Delta x,\Delta z$) is given by the location of the single Gaussian-shaped peak in the correlation function $r$ (arg max(|

*r*|)); for this to be evident, velocity gradients within an interrogation window must be negligible.The formulations for wave intensity assume velocity is one dimensional; we were consequently not interested in obtaining the velocity profile across the vessel, and a single interrogation window was therefore used. It was preferred to many smaller windows because it increased computational speed. The window was sized to span the entire lumen unless otherwise stated.

A distribution of velocities within a window causes the correlation peak to become broadened and skewed (

Westerweel, 2008

). The correlation coefficient, $r$, represents the convolution of the “particle image” (autocorrelation of an image) with the probability density function of the velocity field within the window; the location of the maximum thus represents the modal displacement rather than the mean (Adrian, 1998

; Adrian and Westerweel, 2011

). Finding the centroid rather than the peak was used to obtain the mean (Poelma and Westerweel, 2011

).Assuming negligible mean displacement of blood in the $z$ direction, the centroid in the $x$ direction was computed as

where for a window size of $N$, the summation ranges from $-N/2$ to $N/2-1$ and $\delta y=0$. However, the correlation plane contains noise, and thus, integration limits had to be set for eqn (5); for large displacements, the broadened peak is shifted toward the edge of the plane, introducing bias. First the particle image diameter ($\sigma $) was measured from the autocorrelation map (24 pixels, as measured using the 1/

$\Delta {x}_{c}=\phantom{\rule{0.33em}{0ex}}\sum _{\delta x}r\left(\delta x,0\right)\delta x/\sum _{\delta x}r\left(\delta x,0\right)$

(5)

where for a window size of $N$, the summation ranges from $-N/2$ to $N/2-1$ and $\delta y=0$. However, the correlation plane contains noise, and thus, integration limits had to be set for eqn (5); for large displacements, the broadened peak is shifted toward the edge of the plane, introducing bias. First the particle image diameter ($\sigma $) was measured from the autocorrelation map (24 pixels, as measured using the 1/

*e*^{2}rule [Poelma and Westerweel, 2011

]). For non-reversing flows, we assumed that flow at the wall is zero (no-slip condition) and increases toward the center. Thus, the integration limits were defined as $-\sigma /2$ pixels to the maximum displacement ($\Delta {x}_{\text{max}}$) $+\sigma /2$. To determine $\Delta {x}_{\text{max}}$, the modal displacement was found. Finally, velocity ${U}_{\text{raw}}$ was calculated as Δ*xc*/*dt,*where $dt$ is the time between compounded frames (1 ms).Pixel locking is a known problem of centroid calculations (

Forliti et al., 2000

) but is not an issue if $\sigma $ is >3 pixels, as is the case here. Precision of the centroid calculation was increased by interpolating the correlation map; this was done by zero padding in the frequency domain before inverse transformation.Interrogation windows were centered on the part of the vessel where the diameter measurements were made. The window width was set to 80 pixels in the first image, and the search window to 80 plus twice the expected maximum pixel displacement in the second image (we anticipated velocities would not exceed 1 m/s [

Weinberg and Ethier, 2007

]); this ensured the tails of the broadened peak were captured for large displacements.#### Cross-sectional mean velocity estimation

The methods described in the previous section provided a measure of the mean velocity (${U}_{\text{raw}}$) in a 2-D lengthwise slice of the lumen that included the centerline of the vessel. Although the maximum velocity in the 2-D slice is the same as the maximum in the 3-D vessel, that is not necessarily true of the mean velocity. For Poiseuille flow in a cylindrical tube, where the velocity profile is a paraboloid, it is well established that the mean velocity is 0.5 times the maximum. However, in a thin 2-D slice through the centerline of such flow, for which the velocity profile is a parabola, the mean is 0.75 times the maximum. (There is comparatively less of the slow-moving, near-wall fluid.) Real arterial flow is pulsatile and more nearly of the Womersley type than the Poiseuille type.

Wave intensities derived using the 2-D velocity ${U}_{\text{raw}}$ and Womersley-based 3-D velocity ${U}_{\text{Wo}}$ were nearly linearly related (Fig. S5, online only). The former was used in calculations of wave intensity, because only relative values were required to compare

*P-U*and*D-U*methods or ventricular function before and after drug administration. However, absolute numbers were required when calculating wave speeds, to permit comparison with previous values, and ${U}_{\text{Wo}}$ was therefore used. It was calculated as described in Appendix S1 (online only).#### Diameter measurement

The artery was assumed to change dimension only in the radial direction. Diameter waveforms were computed by 1-D cross-correlation of successive frames in ${S}_{\text{tissue}}$. First a region of interest (ROI) with a width equal to 5 A-lines and a height of 30 pixels was selected on the anterior wall of the aorta. The signal values within the ROI were averaged laterally, and the location of the maximum corresponding to the inner layers of the wall was identified from the envelope. This was further refined by subpixel Gaussian fitting.

One-dimensional cross-correlation functions were computed for each lateral location, according to equation 4 with ${A}^{\left(k\right)}$ and ${B}^{\left(k+1\right)}$ vectors of equal size. The resulting correlation functions were averaged, and the maximum was identified to give the axial displacement, again using subpixel Gaussian fitting. This was repeated for the remaining image pairs, with the window offset by the previous displacement each time. The wall motion waveform was given by cumulatively summing the displacements over time. The process was repeated for the posterior artery wall. The diameter change waveform was then given by the difference between the anterior and posterior waveforms; addition of the initial diameter measured in frame 1 gave absolute diameter

*D*.### Wave intensity analysis

Waveforms were smoothed using a second-order Savitsky–Golay filter with window length of 19 ms, equivalent to 19 samples for all data sets. The length was optimized by trial and error. Successive beats in the pressure, diameter and velocity waveforms were ensemble averaged using the R wave of the ECG as a fiducial marker; only complete cardiac cycles were included.

Net wave intensities were calculated as

where the prefix “n” refers to the non-invasive diameter formulation, and $dt$ is the time between two consecutive samples. $dI$ was normalized by $d{t}^{2}$ to make the magnitude of $dI$ independent of the sampling frequency (

$dI=\phantom{\rule{0.33em}{0ex}}\frac{dP}{dt}\phantom{\rule{0.33em}{0ex}}\frac{dU}{dt}$

(6)

${\phantom{\rule{0.33em}{0ex}}}_{\mathrm{n}}dI=\phantom{\rule{0.33em}{0ex}}\frac{dD}{dt}\phantom{\rule{0.33em}{0ex}}\frac{dU}{dt}$

(7)

where the prefix “n” refers to the non-invasive diameter formulation, and $dt$ is the time between two consecutive samples. $dI$ was normalized by $d{t}^{2}$ to make the magnitude of $dI$ independent of the sampling frequency (

Niki et al., 1999

).- Niki K
- Sugawara M
- Uchida K
- Tanaka R
- Tanimoto K
- Imamura H
- Sakomura Y
- Ishizuka N
- Koyanagi H
- Kasanuki H.

A noninvasive method of measuring wave intensity, a new hemodynamic index: Application to the carotid artery in patients with mitral regurgitation before and after surgery.

*Heart Vessels.*1999; 14: 263-271

The peak magnitudes, area under the curve (“wave energy”) and timings with respect to the R wave of the ECG were calculated for the dominant waves. The start and end of the waves were identified by the zero crossing of the wave intensity.

Wave speed $c$ was calculated using single-point “loop methods.” In plots of ensemble-averaged $P$ (or $\text{ln}\left(D\right)$) against $U$, there is a linear portion of the loop in early systole whose gradient can be used to determine $c$ according to

(

(

${c}_{PU}=\phantom{\rule{0.33em}{0ex}}\frac{dP}{\rho dU}$

(8)

(

Khir et al., 2001

; Davies et al., 2006

) or${C}_{ln\left(D\right)U}=\phantom{\rule{0.33em}{0ex}}\frac{1}{2}\phantom{\rule{0.33em}{0ex}}\frac{dU}{dln\left(D\right)}$

(9)

(

Feng and Khir, 2010

), where $\rho $ is the density of the blood (1044 kg/m^{3}). The gradients during the first 10 ms (rabbit) or 50 ms (human participant) of ejection were calculated through linear fitting.The equations are based on the assumption of a reflection-free period in early systole. A third wave speed (

was also calculated when $P$ and $D$ were recorded simultaneously. Not only is this value free of potential velocity errors, but the derivation of the equation does not assume unidirectional wave travel.

Kowalski et al., 2017

)${C}_{ln\left(D\right)P\phantom{\rule{0.33em}{0ex}}}=\phantom{\rule{0.33em}{0ex}}\sqrt{\frac{1}{2\rho}\phantom{\rule{0.33em}{0ex}}\frac{0.5dP}{d\phantom{\rule{0.33em}{0ex}}\text{ln}\left(D\right)}}$

(10)

was also calculated when $P$ and $D$ were recorded simultaneously. Not only is this value free of potential velocity errors, but the derivation of the equation does not assume unidirectional wave travel.

### Human study

Values of $D$ and $U$ were obtained by non-invasive ultrasound of the carotid artery of a healthy young subject, using the parameters described previously except that the MI was increased to 0.4. Data were aligned and ensemble averaged using the R wave of the ECG as a datum.

The experiment was part of a trial approved by the National Health Service (IRAS ID 248724, HRA REC Reference 18/SC/0563) and the Imperial College Healthcare NHS Trust, and registered as ISRCTN 41232. Informed consent was given.

### Statistics

Values of

*p*were obtained using Student's paired*t*-test unless otherwise stated. Limits of agreement in Bland–Altman plots were computed using the SimplyAgree package in R (The R Project for Statistical Computing) to account for repeated measures.## Results

In the following, where waveforms are plotted against time,

*t*= 0 corresponds to the R wave of the ECG. Individual animals are identified by four-letter codes.Figure 1a and 1b illustrate, respectively, invasive pressure measurements and non-invasive diameter measurements for one animal (POEK); data for all animals are provided in Figure S6a and S6b (online only). The two types of waveform clearly have the same shape, but the relation between them is not expected to be perfectly linear: the elasticity of the vessel decreases as diameter increases (“strain stiffening”), and viscous effects cause a lag between changes in pressure and changes in diameter. The non-linearity and hysteresis are visible when $D$ is plotted against $P$ (Fig. 1c and Figure S6c [online only]); data for $D$ were filtered as described previously. Rabbit PRFH had high pressure, anomalous waveforms and a stiff vessel.

Figure 2a provides the corresponding invasive and non-invasive velocity measurements; data for all rabbits are given in Figure S7a (online only). For this comparison, the velocimetry interrogation window was reduced to a small box at the center of the vessel, as the Doppler system identifies the peak velocity, ${U}_{\text{cl}}$. (For the same reason, Doppler-derived velocities can be unsuitable when calculating absolute wave speeds; see

Mynard et al., 2018

). This reduced the spread of velocities within the window, so pixel displacements could be determined using simple Gaussian peak fitting.It is apparent in Figure S7a that the invasive and non-invasive velocity measurements, while in broad agreement, show a systematic difference: the invasive curves resemble smoothed versions of the non-invasive curves. This is particularly apparent at the systolic foot: the change in gradient is visibly sharper in the non-invasive results for many rabbits. The sampling frequency of the Volcano Combowire system and the additional low-pass filtering mean that the system cannot capture the fast early-systolic accelerations in the rabbit. There additionally appears to be a scaling error in rabbits PPEX, PPFD, PRDT and PRFH, which was most likely caused by the difficulty of perfectly aligning the orientation of the catheter with the mean flow direction: the Doppler probe at the catheter tip obtains velocity information only along the beam axis. Because of this issue, B-mode–derived velocity data are used in all subsequent analyses.

Figure 2b and Figure S7b (online only) additionally illustrate velocity waveforms calculated in three ways from an interrogation window spanning the entire vessel: ${U}_{\text{raw}}$ was obtained assuming 2-D Poiseuille flow, ${U}_{\text{Wo}}$ by assuming a 3-D Womersley flow (in both cases the displacement being estimated from the centroid of the cross-correlation peak), and ${U}_{\text{gauss}}$ by using a Gaussian fit to find the cross-correlation peak.

Figure 2c demonstrates a linear relation between velocities derived from catheter-based Doppler and B-mode ultrasound measurements for all rabbits.

Figure 3a–c illustrates calculations of wave speed, $c$, by three different “loop” methods—the

*PU*loop, ln(*D*)*U*loop and ln(*D*)*P*loop—for a single rabbit (POEK); data for all rabbits are given in Figure S8 (online only). In Figure 3a and 3c, $D$ and $U$, obtained from the non-invasive ultrasound methods, were filtered as described previously. The wave speed given for each rabbit is the mean of the wave speed obtained in each repeat, not the wave speed computed from the mean loop. All methods have anomalous results for rabbit PRFH, which had a stiff vessel.Figure 3d–h illustrates the relation between the values obtained by the different methods for all rabbits. Compared with the ln(

*D*)*P*method, which is considered a gold standard, the PU method consistently underestimated wave speed, and the ln(*D*)*U*method consistently overestimated it. However, the mean errors were small: –0.55 and +0.36 m/s, respectively, when compared with a best estimate of 5.02 m/s (the mean value obtained with the ln(*D*)*P*method). The errors do not appear to scale with wave speed but remain relatively constant for all values obtained with the ln(*D*)*P*method.Figure 4a and 4b illustrate wave intensities calculated by both the PU and DU methods for a single rabbit (POEK). In Figure 4a, $U$ is filtered as described previously. Equivalent data for all rabbits are illustrated in Figure S9a and S9b (online only). Note that these intensities are the sum of forward-traveling and backward-traveling wavelets (with positive and negative intensities, respectively) at each time point. Both methods exhibit three waves during the cardiac cycle. There is a large net forward-traveling compression wave (“W1”) followed immediately by a smaller net backward-traveling wave (“R”) and then a further net forward-traveling expansion wave (“W2”) before a return to baseline. The backward-traveling wave comprises reflections from distal sites.

Figure 4c and 4d show comparisons between the two methods for peak intensity and wave energy (the area under the intensity curve). There is a clear positive relation between the two methods for each wave. Rabbit PRFH is again an outlier.

Figure 5 illustrates the effect on wave intensity of increasing doses of esmolol, a short-acting cardioselective β

_{1}-adrenergic blocker that reduces the force and rate of cardiac contraction. Wave intensity was assessed by the invasive and non-invasive methods; illustrative data, plotted in both 2-D and 3-D, are given for rabbit PFLH. Transient dips in wave intensity are evident as a color change in the 2-D plot and as a color change and decreased height in the 3-D plot, shortly after each*dotted white line*representing administration of esmolol. The dips increase in size with increasing dose. A delay in the arrival of the W1 and W2 waves is also evident after the administration of each bolus—in the lefthand panels, there is a rightward shift in the peaks below*each white dotted line*, with the largest delay occurring at the highest dose.In Figure 6a–d are the average peak intensities and energies of the W1 and W2 waves, calculated by both the PU and DU methods. Mean data for nine rabbits are presented ±1 SD to illustrate variability within the sample. (Rabbit PPBT was not included in the analysis as it required an anesthetic top-up during the period of esmolol administration, rendering the results unreliable.)

*DU*-derived intensities and energies are also given for two control rabbits administered vehicle rather than esmolol.Dips in intensity and energy are evident after each esmolol bolus for both waves and both methods. No such trends are seen in the control animals. Statistical tests were performed for the experimental group by comparing the first measurement after each dose of esmolol with the preceding measurement. The drop in W1 peak intensity and energy was significant at all doses with both methods. The drop in W2 wave energy was also significant in all these cases. For W2 peak intensity, the effect was not significant (

*p*> 0.05) for the lowest esmolol dose assessed with the*PU*method and for all doses with the*DU*method. The fact that statistical significance was obtained for wave energies but less consistently with peak intensities is likely due to alignment issues (see the Discussion).Average variation in the time of arrival of the W1 wave, relative to the R wave of the ECG, is illustrated in Figure 7a. (Data for rabbits PPEX and PRDT at

*t*= 25 min were omitted because of a poor ECG signal.) A clear delay is seen after each bolus, and the effect was significant at all doses. Delays of similar magnitude, again all significant, were seen in the interval between the W1 and W2 waves (Fig. 7b). No such changes were seen in the control animals.No significant differences in wave speed were detected after each esmolol bolus with either the

*PU*or ln(*D*)*U*loop method (Fig. S10 [online only],*p*> 0.05). Mean pressure was also not significantly altered (data not shown).Figure 8 illustrates the feasibility of using the non-invasive method in the carotid artery of human participants.

## Discussion

An advantage of calculating wave intensity from arterial diameter and blood velocity is that both can be measured using non-invasive techniques. We developed methods based on ultrasound because its low cost and real-time imaging capability increase clinical utility. Wave intensities are computed from derivatives, which increases the requirement for measurement accuracy and precision. We avoided the use of Doppler ultrasound because of its fundamental limitations in quantifying blood velocity and because accuracy can be further compromised by the difficulty of employing the optimal, orthogonal beam angles for the Doppler and diameter measurements. Instead, changes in diameter and velocity were obtained from the same, successive B-mode images. This precludes the alignment issues apparent when using the current “gold standard” invasive system, which relies on separate pressure and flow transducers.

The new method presents several technical challenges. Blood is a poor ultrasound scatterer at the low frequencies required to penetrate tissue to the depth of conduit arteries. That necessitated the use of SVD to obtain separation of the blood signal from the tissue signal and noise. SVD is computationally intensive, and the manual selection of thresholds was slow. Randomized SVD is faster and has been realized in real-time applications (

Lok et al., 2020

). Threshold selection could be automated by identifying turning points in the singular value magnitude or temporal vector frequency plots (Song et al., 2017

) or by identifying tissue and blood subspaces from the similarity matrix of the spatial vectors (Baranger et al., 2018

). We found that precise selection is not critical: the lower threshold ($\mathrm{T}{\mathrm{h}}_{1}$) could be changed between 115 and 150 and the upper threshold ($\mathrm{T}{\mathrm{h}}_{2}$) could be omitted entirely without visible effect on the resulting velocity waveform (Fig. S4c, online only).A further technical challenge is capturing the fast-moving arterial blood with sufficient temporal resolution. That was achieved by using an ultrafast, plane-wave scanner capable of frame repetition rates of ≈10 kHz. In such systems, the beam is not focused, which accelerates scanning but degrades the image. Image quality is recovered by averaging frames and by compounding scans at different angles. The net frame rate still exceeds that of conventional, focused scanners.

Velocity estimation is also complex. The use of B-mode imaging means that scattering from red blood cells can theoretically be tracked in two dimensions. In conventional ultrasound imaging velocimetry, a small field of scatterers is identified within the vessel in one image, and then a search is made for a field in the succeeding image that has the best-fit pattern of scatterers; the displacement of the field divided by the time interval between the two scans gives the velocity for that region of blood. Repeating the cross-correlation procedure for other fields across the diameter of the artery gives the velocity profile. We modified this procedure: a single field covering the whole diameter of the vessel was used. This reduces noise to a practicable level and increases the speed of computation. It is also consistent with the 1-D formulation of wave intensity. However, the resulting increase in the spread of velocities within each field introduces complexity. Velocity gradients are higher near the wall than near the center of the vessel, and the pattern of scatterers will therefore be more disrupted near the wall as the region of blood moves down the vessel. As a result, the cross-correlation algorithm is biased toward tracking scatterers near the center. Because absolute velocities are higher at the center, this leads to an overestimate of the mean flow velocity. In our method, this effect was ameliorated by using the centroid rather than the peak of the distribution of correlation coefficients

An additional issue is that the ultrasound image approximates a 2-D longitudinal slice through the center of the vessel. (In reality, this slice has finite thickness but that is ignored here.) Considering the simple case of fully developed flow, the velocity profile would be a parabola in the 2-D slice but a paraboloid in the 3-D vessel. The mean velocity for the parabola is three-fourths of the maximum velocity, whereas the mean for the paraboloid is half the maximum velocity. In the real vessel, where Womersley-type flow occurs, the velocity profile is more complex and varies over the cardiac cycle. For calculations of wave speed, where absolute rather than relative magnitudes matter, we used the inverse Womersley solution; elsewhere, to minimise computational cost, the 2-D values were used.

The assessment of arterial diameter over the cardiac cycle also uses the cross-correlation technique but its application in this case is less complex because the wall moves more like a rigid object, with a single velocity at each time step. The method is routinely employed (

*e.g.,*Rabben et al., 2002

). The movement of the wall can be less than a single pixel; again, this problem arises from the need to use relatively low ultrasound frequencies, but subpixel resolution is readily obtained by fitting a curve to the Gaussian cross-correlation function.The calculation of wave speed from the gradient of the

*PU*loop or the ln(*D*)*U*loop assumes that all waves are forward traveling. This condition is most likely to be satisfied during the period when pressure and velocity rise at the start of ventricular ejection. Rabbits and people have broadly similar wave speeds, but the condition is more likely to be violated in the rabbit because of its higher heart rate, meaning there is less time for reflected wave energy to dissipate, and its shorter arteries, which reduce the time before reflections arrive at the measuring site. For these reasons, the present study used the first 10 ms of ejection in rabbits and the first 50 ms in the human participant.Plots of aortic pressure, obtained with an intra-arterial catheter, and aortic diameter, obtained using the B-mode cross-correlation method (Fig. 1; Fig. S6, online only), both revealed in all animals the expected shape of a rapid systolic rise, followed by a slower decay that was clearly biphasic; that is, there was a secondary, diastolic peak or at least an inflection in the gradient of the decay, before the next systole, and not a continuous exponential decay as might be expected in the absence of wave reflections. The relation between pressure and diameter revealed the expected non-linearity and hysteresis caused by strain-stiffening and viscous effects.

Plots illustrating blood flow velocities obtained both by the Doppler ultrasound probe mounted at the tip of the intra-arterial catheter and by non-invasive velocimetry (Fig. 2; Fig. S7, online only) gave the expected triphasic aortic waveform in all animals and with both techniques: there was a period of forward-going flow in early systole, followed by a period of nearly zero flow and, finally a period of low, forward-going flow in diastole. Furthermore, the shapes of the curves obtained by the two methods were similar. A visible discrepancy between the methods was smoothing of the peaks and troughs in the Doppler measurements. The intra-arterial flow measurements were made with a device designed for human use; its sampling rate is probably insufficient to capture fine features of the shorter cardiac cycle in the rabbit. It also has inbuilt, low-pass filtering. A second discrepancy was a difference in the absolute magnitude in four rabbits: the invasive measurements gave lower peak velocities than were obtained with the non-invasive method. We attribute this to the well-known difficulty in aligning the Doppler probe with the predominant flow direction (

Davies et al., 2006

). With this exception, the scatterplot of the two data sets revealed a strong relation, the data lying close to the line of identity.Wave speeds derived from the first 10 ms of the ln(

*D*)*U*loop were higher than those derived from the equivalent period of the*PU*loop (Fig. 3; Fig. S8) but the difference was relatively small: averages were, respectively, 5.36 and 4.48 m/s. Reflections may have been present during this part of ventricular ejection. They introduce errors of opposite direction in the two methods: wave speed is overestimated when using the ln(*D*)*U*loop but underestimated when using the*PU*loop (Alastruey, 2011

; Swillens et al., 2013

; Borlotti et al., 2014

; Segers et al., 2014

; Kowalski et al., 2017

; Campos Arias et al., 2019

). Comparing data obtained using the two methods with data derived from the ln(*D*)*P*loop—thought to be less affected by reflections (Kowalski et al., 2017

)—supports this view as the latter results lie between the other two sets of values (Fig. 3). A similar effect would have been produced by overestimating *U*. Note that wave speed is not required when assessing net wave intensity or net wave energy, but only for separating waves into their forward and backward components and for investigating vessel wall stiffness.The

*PU*and*DU*methods gave comparable patterns of wave intensity (Fig. 4). In all rabbits, both indicated a large net forward-traveling (“W1”) wave in early systole that was followed by a small net backward-traveling (“R”) wave and then a more complex and dispersed net forward-traveling (“W2”) wave. These are thought dominantly to correspond to the initial compression wave caused by ventricular contraction, its reflection and the forward-going expansion wave that derives from an interaction between the inertia of the blood and the reduced and then reversed contraction of the ventricle. The more disperse and complex nature of the W2 wave could be explained at least in part by the longer period between it and the*R*-wave of the ECG, used to time align individual heart beats when creating an ensemble average. This, coupled with beat-to-beat variation in the duration of the cardiac cycle would have the effect of smearing out a wave even it were concentrated at the same relative location in each individual beat. Additionally, the W2 wave may consist of more components than the other waves. Figure 4 reveals a strong correspondence between wave intensities measured by the two methods and also between wave energies (*i.e.,*the integrals of the intensity curves for each wave). The relation in both cases is curvilinear, as expected from the strain-stiffening behaviour of the wall.Ventricular dysfunction was transiently induced with esmolol, a short-acting cardioselective β

_{1}-adrenergic receptor blocker that decreases the force and rate of ventricular contraction. In addition to its clinical use, the drug has been employed in studies investigating the potential performance of left ventricular assist devices in heart failure (*e.g.,*Huang et al., 2003

; Gallagher et al., 2007

). Both the *PU*and*DU*methods were unequivocally able to detect decreased intensity, decreased energy and a delayed arrival of the W1 wave at all the esmolol doses employed (Fig. 5, Fig. 6, Fig. 7). Both could also detect similar effects on the energy and arrival time of the W2 wave (Fig. 5, Fig. 6, Fig. 7). The*DU*method was less sensitive than the*PU*method at detecting the esmolol-induced decrease in peak intensity of the more spread out W_{2}wave (Fig. 6). Neither method achieved a significant result at the lowest dose, suggesting that intensity is less consistent than the energy for this wave.## Conclusions

The new technique, using the diameter–velocity formulation of wave intensity analysis and based on B-mode rather than separate Doppler and M-mode measurements, is sufficiently accurate and non-invasive that it may have clinical utility. The rabbit aorta was chosen as a test bed because its diameter, depth and blood velocity are similar to those of human peripheral arteries, Furthermore, rabbits have wave patterns similar to those of people, unlike small rodents, for example. A proof-of-concept experiment indicated that the non-invasive method can be applied to the human carotid artery. The pattern of wave intensities and the calculated wave speed (Fig. 8) were consistent with data from other non-invasive techniques (

Curtis et al., 2007

). Potential uses include screening for heart failure, tracking its progression and improving the accuracy of prognoses. Although the algorithms required to achieve accurate results have proven complex, they are not difficult to implement now that they have been developed. A recent modeling study (Reavette et al., 2021

) suggests that incorporating artificial intelligence into the analysis improves discrimination to the extent that the method might be able to categorize individual patients.### Acknowledgments

This work was supported by MRC Confidence-in-Concept and Proof-of Concept funds, a BHF Research Excellence Award studentship to K.R. and British Heart Foundation Project Grant PG/18/48/33832.

### Conflict of interest disclosure

P.D.W. is the named inventor on patents filed by Imperial College Innovations, Ltd that describe some of the underlying technology.

### Data availability statement

Data are available from the corresponding author on reasonable request.

## Appendix. Supplementary materials

## References

- Double exposure, multiple-field particle image velocimetry for turbulent probability density.
*Opt Laser Eng.*1998; 9: 211-228 - Particle image velocimetry.Cambridge University Press, New York, NY2011
- Numerical assessment of time-domain methods for the estimation of local arterial pulse wave speed.
*J Biomech.*2011; 44: 885-891 - Adaptive spatiotemporal SVD clutter filtering for ultrafast Doppler imaging using similarity of spatial singular vectors.
*IEEE Trans Med Imaging.*2018; 37: 1574-1586 - A non-invasive clinical application of wave intensity analysis based on ultrahigh temporal resolution phase-contrast cardiovascular magnetic resonance.
*J Cardiovasc Magn Reson.*2012; 14: 57 - Modeling of Congenital Hearts Alliance (MOCHA) Collaborative Group. Reduced ascending aorta distensibility relates to adverse ventricular mechanics in patients with hypoplastic left heart syndrome: Noninvasive study using wave intensity analysis.
*J Thorac Cardiovasc Surg.*2012; 144: 1307-1313 - Experimental evaluation of local wave speed in the presence of reflected waves.
*J Biomech.*2014; 47: 87-95 - Mapping the site-specific accuracy of loop-based local pulse wave velocity estimation and reflection magnitude: A 1D arterial network model analysis.
*Physiol Meas.*2019; 40075002 - Reduced systolic wave generation and increased peripheral wave reflection in chronic heart failure.
*Am J Physiol Heart Circ Physiol.*2007; 293: H557-H562 - Use of simultaneous pressure and velocity measurements to estimate arterial wave speed at a single site in humans.
*Am J Physiol Heart Circ Physiol.*2006; 290: H878-H885 - Spatiotemporal clutter filtering of ultrafast ultrasound data highly increases Doppler and fUltrasound sensitivity.
*IEEE Trans Med Imaging.*2015; 34: 2271-2285 - Determination of wave speed and wave separation in the arteries using diameter and velocity.
*J Biomech.*2010; 43: 455-462 - Bias and precision errors of digital particle image velocimetry.
*Exp Fluids.*2000; 28: 436-447 - Die Elastizität der Blutegefässe.
*Z Biol.*1920; 71: 255-272 - Efficacy and mechanisms of biventricular and left/right direct cardiac compression in acute heart failure sheep.
*Artif Organs.*2007; 31: 39-44 - A review of the measurement of blood velocity and related quantities using Doppler ultrasound.
*Proc Inst Mech Eng H.*1999; 213: 391-400 - HeartPatch implanted direct cardiac compression: Effect on coronary flow and flow patterns in acute heart failure sheep.
*ASAIO J.*2003; 49: 309-313 - On the accuracy of displacement-based wave intensity analysis: Effect of vessel wall viscoelasticity and nonlinearity.
*PLoS One.*2019; 14e0224390 - Determination of wave speed and wave separation in the arteries.
*J Biomech.*2001; 34: 1145-1155 - Robust and practical non-invasive estimation of local arterial wave speed and mean blood velocity waveforms.
*Physiol Meas.*2017; 38: 2081-2099 - Comparison of invasive and non-invasive aortic wave intensity and wave power analyses in sheep.
*Physiol Meas.*2019; 40015005 - Stiffening and ventricular–arterial interaction in the ascending aorta using MRI: Ageing effects in healthy humans.
*J Hypertens.*2019; 37: 347-355 - Real time SVD-based clutter filtering using randomized singular value decomposition and spatial downsampling for micro-vessel imaging on a Verasonics ultrasound system.
*Ultrasonics.*2020; 107106163 - Coherent plane-wave compounding for very high frame rate ultrasonography and transient elastography.
*IEEE Trans Ultrason Ferroelectr Freq Control.*2009; 56: 489-506 - Noninvasive assessment of carotid arterial wave speed and distensibility.
*Am J Physiol Heart Circ Physiol.*2018; 315: H1495 - McDonald's blood flow in arteries: Theoretical, experimental and clinical principles.Edward Arnold, London, UK1998
- A noninvasive method of measuring wave intensity, a new hemodynamic index: Application to the carotid artery in patients with mitral regurgitation before and after surgery.
*Heart Vessels.*1999; 14: 263-271 - An introduction to wave intensity analysis.
*Med Biol Eng Comput.*2009; 47: 175-188 - Forward and backward running waves in the arteries: Analysis using the method of characteristics.
*J Biomech Eng.*1990; 112: 322-326 - What stops the flow of blood from the heart?.
*Heart Vessels.*1988; 4: 241-245 - Generalized displacement estimation for averages of non-stationary flows.
*Exp Fluids.*2011; 50: 1421-1427 - Common carotid artery diameter, blood flow velocity and wave intensity responses at rest and during exercise in young healthy humans: A reproducibility study.
*Ultrasound Med Biol.*2017; 43: 943-957 - Noninvasive assessment of the common carotid artery hemodynamics with increasing exercise work rate using wave intensity analysis.
*Am J Physiol Heart Circ Physiol.*2018; 315: H233-H241 - Ultrasound-based vessel wall tracking: An auto-correlation technique with RF center frequency estimation.
*Ultrasound Med Biol.*2002; 28: 507-517 - Comparison of arterial wave intensity analysis by pressure–velocity and diameter–velocity methods in a virtual population of adult subjects.
*Proc Inst Mech Eng H.*2020; 234: 1260-1276 - Wave intensity analysis combined with machine learning can detect impaired stroke volume in simulations of heart failure.
*Front Bioeng Biotechnol.*2021; 9737055 - Determining haemodynamic wall shear stress in the rabbit aorta in vivo using contrast-enhanced ultrasound image velocimetry.
*Ann Biomed Eng.*2020; 48: 1728-1739 - Contrast agent-free assessment of blood flow and wall shear stress in the rabbit aorta using ultrasound image velocimetry.
*Ultrasound Med Biol.*2022; 48: 437-449 - Role of tapering in aortic wave reflection: hydraulic and mathematical model study.
*J Biomech.*2000; 33: 299-306 - Wave reflection leads to over- and underestimation of local wave speed by the PU- and QA-loop methods: Theoretical basis and solution to the problem.
*Physiol Meas.*2014; 35: 847-861 - Ultrasound small vessel imaging with block-wise adaptive local clutter filtering.
*IEEE Trans Med Imaging.*2017; 36: 251-262 - Clinical usefulness of wave intensity analysis.
*Med Biol Eng Comput.*2009; 47: 197-206 - Resolving the time lag between pressure and flow for the determination of local wave speed in elastic tubes and arteries.
*J Biomech.*2009; 42: 1574-1577 - Comparison of non-invasive methods for measurement of local pulse wave velocity using FSI-simulations and in vivo data.
*Ann Biomed Eng.*2013; 41: 1567-1578 - Alterations in aortic wave reflection with vasodilation and vasoconstriction in anaesthetized dogs.
*Can J Cardiol.*2013; 29: 243-253 - Characteristics of the pulse waveform during altered nitric oxide synthesis in the rabbit.
*Br J Pharmacol.*2001; 133: 361-370 - Twenty-fold difference in hemodynamic wall shear stress between murine and human aortas.
*J Biomech.*2007; 40: 1594-1598 - On velocity gradients in PIV interrogation.
*Exp Fluids.*2008; 44: 831-842

## Article Info

### Publication History

Published online: November 02, 2022

Accepted:
September 24,
2022

Received in revised form:
September 12,
2022

Received:
July 1,
2022

### Publication stage

In Press Corrected Proof### Identification

### Copyright

© 2022 The Author(s). Published by Elsevier Inc. on behalf of World Federation for Ultrasound in Medicine & Biology.

### User License

Creative Commons Attribution (CC BY 4.0) | How you can reuse

Elsevier's open access license policy

Creative Commons Attribution (CC BY 4.0)

## Permitted

- Read, print & download
- Redistribute or republish the final article
- Text & data mine
- Translate the article
- Reuse portions or extracts from the article in other works
- Sell or re-use for commercial purposes

Elsevier's open access license policy