NON-INVASIVE ASSESSMENT BY B-MODE ULTRASOUND OF ARTERIAL PULSE WAVE INTENSITY AND ITS REDUCTION DURING VENTRICULAR DYSFUNCTION

—Arterial pulse waves contain clinically useful information about cardiac performance, arterial stiff- ness and vessel tone. Here we describe a novel method for non-invasively assessing wave properties, based on measuring changes in blood ﬂow velocity and arterial wall diameter during the cardiac cycle. Velocity and diam- eter 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 proc- essed 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 clin- ical 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. (E-mail


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 (Frank 1920), 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. (2017Pomella 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. (2012aBiglino 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 (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.

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 highfidelity, 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 Volcanoderived pressure with diameter measured by the noninvasive 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. (2020Riemer 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 (Demen e 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 where λ 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 Th 1 and Th 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 blood . Briefly, an "interrogation window" was defined in the vessel lumen in one image, and the crosscorrelation 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 ðkÞ and search window b ðkþ1Þ 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, dx and dz are pixel shifts in the lateral and axial directions, * is the complex conjugate and F À1 is the inverse Fourier transform. Conventionally, the displacement (Dx; Dz) 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 dy ¼ 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 (s) 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 À s=2 pixels to the maximum displacement (Dx max ) þs=2.
To determine Dx max , the modal displacement was found. Finally, velocity U raw was calculated as Dx c /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 s is >3 pixels, as is the case here. Precision of the centroid 476 Ultrasound in Medicine & Biology Volume 49, Number 2, 2023 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 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 slowmoving, 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 raw and Womersley-based 3-D velocity U 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 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 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 ðkÞ and B ðkþ1Þ 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 dt 2 to make the magnitude of dI independent of the sampling frequency (Niki et al. 1999). 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 lnðDÞ) against U , there is a linear portion of the loop in early systole whose gradient can be used to determine c according to (Khir et al. 2001;Davies et al. 2006) or (Feng and Khir 2010), where r 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 (Kowalski et al. 2017) 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 ttest unless otherwise stated. Limits of agreement in BlandÀAltman plots were computed using the Simply-Agree 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 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 noninvasive 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 raw was obtained assuming 2-D Poiseuille flow, U 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 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 Bmode 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 noninvasive 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 backwardtraveling 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 Fig. 2. (a) Waveform of centerline blood velocity (U cl ) for a single rabbit (POEK), derived from catheter-based Doppler (crosses) and non-invasive B-mode (black line) ultrasound measurements, with root mean square error (RMSE) indicating the difference between them. (b) Blood velocity waveforms calculated in four different ways from non-invasive B-mode ultrasound measurements: the centerline velocity (nU cl , grey line) is repeated from (a), and the remaining three waveforms are derived from an interrogation window spanning the entire vessel rather than just the central region but using different methods for estimating the significant displacement of speckles between frames: U raw (black line) was obtained assuming 2-D Poiseuille flow, U Wo (circles) was obtained by assuming 3-D Womersley flow (in both cases the displacement being estimated from the centroid of the cross-correlation peak (see Fig.  S3) and U gauss was obtained using a Gaussian fit to find the cross-correlation peak. In (a) and (b), data are the averages of all (n = 3) repeats. Data for all rabbits are given in Figure S7 (online only). (c) Relation between the peak and mean centerline velocities derived from catheter-based Doppler and Bmode ultrasound measurements (the dashed line being the line of identity) for all rabbits. Colors and four-letter codes indicate the individual animals. The non-invasive method is indicated by the prefix "n" in the axis label.
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 b 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.  Figure S8 (online only). Dashed lines represent the ensemble average for each repeat, and the solid line is the average of all (n = 3) repeats. Wave speed, c, was calculated over the first 10 ms of ejection, shown in blue. This duration was determined by trial and error, with the aim of using the same duration in all rabbits and with both P-U and D-U methods. (dÀf) Comparisons between the loop methods for all rabbits. The solid black line is the orthogonal ("Deming") regression, and the dotted line is the line of identity. (g, h) BlandÀAltman plots with repeated measures revealing discrepancies with the PU and ln(D)U methods, assuming the ln(D)P method is correct. Solid lines represent the means, and dashed lines the 95% limits of agreement. 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 Fig. 4. Wave intensities derived from (a) PU and (b) DU data in a single rabbit (POEK). Data for all rabbits are given in Figure S9 (online only). Dashed lines represent the ensemble average for each repeat, and the solid line is the average of all (n = 3) repeats. Three waves are visible: a net forward-traveling compression wave, "W1," followed immediately by a small net backward-traveling reflected "R" wave and then a further net forward-traveling expansion wave, "W2." (c, d) Comparisons between the two methods for (c) the peak intensity and (d) the energy of each wave in all rabbits. The non-invasive method is indicated by the prefix "n" in the y-axis label.
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 (Th 1 ) could be changed between 115 and 150 and the upper threshold (Th 2 ) could be omitted entirely without visible effect on the resulting velocity waveform (Fig. S4c, online only).
A further technical challenge is capturing the fastmoving 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 Bmode 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 Fig. 6. Effect of esmolol on wave intensity and energy. All plots illustrate the mean (solid line) § standard deviation (dashed lines) for n = 9 rabbits. Intensity of (a) the W1 wave and (b) the W2 wave, and energy of (c) the W1 wave and (d) the W2 wave, were determined by the invasive PU method (left) and the non-invasive DU method (center). Esmolol was administered at increasing doses immediately after the readings at 0, 10 and 20 min. Data were normalized by the mean value of the intensity or energy throughout the 30-min period in the same rabbit. In the right-hand plots are individual data for two rabbits (RNLX, RNNE) administered vehicle only, determined by the non-invasive diameterÀvelocity method. Asterisks indicate significance: *p < 0.05; **p < 0.01; ***p < 0.005. 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 intraarterial 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 Fig. 7. (a) The interval between the R-wave of the electrocardiogram and the arrival of the W1 wave, and (b) the interval between the arrival of the W1 and W2 waves, determined by the non-invasive diameter-velocity method. Left hand plots show the mean (solid line) § SD (dashed lines) for n = 9 rabbits; esmolol was administered at increasing doses immediately after the readings at 0, 10 and 20 min. The average interval throughout the 30-min period in the same rabbit was subtracted from each reading. The right-hand plots provide individual data for two rabbits (RNLX, RNNE) administered vehicle only. Asterisks indicate significance: **p < 0.01; ***p < 0.005. 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 loopthought 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 strainstiffening behaviour of the wall.
Ventricular dysfunction was transiently induced with esmolol, a short-acting cardioselective b 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 (Figs. 5À7). Both could also detect similar effects on the energy and arrival time of the W2 wave (Figs. 5À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 Bmode 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-ofconcept 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.