Advertisement

Non-invasive Assessment by B-Mode Ultrasound of Arterial Pulse Wave Intensity and Its Reduction During Ventricular Dysfunction

Open AccessPublished:November 02, 2022DOI:https://doi.org/10.1016/j.ultrasmedbio.2022.09.016

      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 WW
      • O'Rourke MF.
      McDonald's blood flow in arteries: Theoretical, experimental and clinical principles.
      ). When contraction slows and then reverses, a wave of decreased pressure, velocity and diameter propagates in the same direction (
      • Parker KH
      • Jones CJ
      • Dawson JR
      • Gibson DG.
      What stops the flow of blood from the heart?.
      ). These waves partially reflect and re-reflect where vessel geometry or structure changes (
      • Segers P
      • Verdonck P.
      Role of tapering in aortic wave reflection: hydraulic and mathematical model study.
      ). They carry information concerning cardiac and vascular properties: wave speed (often termed pulse wave velocity [PWV]) is related to arterial stiffness (
      • Frank O.
      Die Elastizität der Blutegefässe.
      ), wave reflections can be altered when vessel tone changes (
      • Weinberg PD
      • Habens F
      • Kengatharan M
      • Barnes SE
      • Matz J
      • Anggård EE
      • Carrier MJ.
      Characteristics of the pulse waveform during altered nitric oxide synthesis in the rabbit.
      ;
      • Wang JJ
      • Bouwmeester JC
      • Belenkie I
      • Shrive NG
      • Tyberg JV.
      Alterations in aortic wave reflection with vasodilation and vasoconstriction in anaesthetized dogs.
      ) and wave intensity depends, at least in part, on cardiac performance (
      • Curtis SL
      • Zambanini A
      • Mayet J
      • McG Thom SA
      • Foale R
      • Parker KH
      • Hughes AD
      Reduced systolic wave generation and increased peripheral wave reflection in chronic heart failure.
      ;
      • Sugawara M
      • Niki K
      • Ohte N
      • Okada T
      • Harada A.
      Clinical usefulness of wave intensity analysis.
      ). Here we focus on the latter.
      • Curtis SL
      • Zambanini A
      • Mayet J
      • McG Thom SA
      • Foale R
      • Parker KH
      • Hughes AD
      Reduced systolic wave generation and increased peripheral wave reflection in chronic heart failure.
      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 M
      • Niki K
      • Ohte N
      • Okada T
      • Harada A.
      Clinical usefulness of wave intensity analysis.
      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 KH
      • Jones CJ.
      Forward and backward running waves in the arteries: Analysis using the method of characteristics.
      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 KH.
      An introduction to wave intensity analysis.
      ). 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 SL
      • Zambanini A
      • Mayet J
      • McG Thom SA
      • Foale R
      • Parker KH
      • Hughes AD
      Reduced systolic wave generation and increased peripheral wave reflection in chronic heart failure.
      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 M
      • Niki K
      • Ohte N
      • Okada T
      • Harada A.
      Clinical usefulness of wave intensity analysis.
      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 JE
      • Whinnett ZI
      • Francis DP
      • Willson K
      • Foale RA
      • Malik IS
      • Hughes AD
      • Parker KH
      • Mayet J.
      Use of simultaneous pressure and velocity measurements to estimate arterial wave speed at a single site in humans.
      ). However, the invasiveness of the method limits its clinical utility.
      • Feng J
      • Khir AW.
      Determination of wave speed and wave separation in the arteries using diameter and velocity.
      developed a formulation in which wave intensity can be derived directly from U and D. A subsequent formulation by
      • Biglino G
      • Steeden JA
      • Baker C
      • Schievano S
      • Taylor AM
      • Parker KH
      • Muthurangu V.
      A non-invasive clinical application of wave intensity analysis based on ultrahigh temporal resolution phase-contrast cardiovascular magnetic resonance.
      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 RM
      • Sherwin SJ
      • Tang M
      • Weinberg PD.
      Comparison of arterial wave intensity analysis by pressure–velocity and diameter–velocity methods in a virtual population of adult subjects.
      ) and others (
      • Kang J
      • Aghilinejad A
      • Pahlevan NM.
      On the accuracy of displacement-based wave intensity analysis: Effect of vessel wall viscoelasticity and nonlinearity.
      ) 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 N
      • Wilhelm EN
      • Kolyva C
      • González-Alonso J
      • Rakobowchuk M
      • Khir AW.
      Common carotid artery diameter, blood flow velocity and wave intensity responses at rest and during exercise in young healthy humans: A reproducibility study.
      ,
      • Pomella N
      • Wilhelm EN
      • Kolyva C
      • González-Alonso J
      • Rakobowchuk M
      • Khir AW.
      Noninvasive assessment of the common carotid artery hemodynamics with increasing exercise work rate using wave intensity analysis.
      ) and
      • Kowalski R
      • Mynard JP
      • Smolich JJ
      • Cheung MMH.
      Comparison of invasive and non-invasive aortic wave intensity and wave power analyses in sheep.
      used Doppler ultrasound to obtain U and B-mode or M-mode ultrasound to obtain D,
      • Li Y
      • Hickson SS
      • McEniery CM
      • Wilkinson IB
      • Khir AW.
      Stiffening and ventricular–arterial interaction in the ascending aorta using MRI: Ageing effects in healthy humans.
      used magnetic resonance imaging (MRI) to obtain U and D, and
      • Biglino G
      • Steeden JA
      • Baker C
      • Schievano S
      • Taylor AM
      • Parker KH
      • Muthurangu V.
      A non-invasive clinical application of wave intensity analysis based on ultrahigh temporal resolution phase-contrast cardiovascular magnetic resonance.
      ,
      • 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.
      ) 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 PR.
      A review of the measurement of blood velocity and related quantities using Doppler ultrasound.
      ), 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 K
      • Rowland EM
      • Broughton-Venner J
      • Leow CH
      • Tang M
      • Weinberg PD.
      Contrast agent-free assessment of blood flow and wall shear stress in the rabbit aorta using ultrasound image velocimetry.
      ). 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 H2O 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 G
      • Tanter M
      • Bercoff J
      • Benech N
      • Fink M.
      Coherent plane-wave compounding for very high frame rate ultrasonography and transient elastography.
      ;
      • Riemer K
      • Rowland EM
      • Leow CH
      • Tang MX
      • Weinberg PD.
      Determining haemodynamic wall shear stress in the rabbit aorta in vivo using contrast-enhanced ultrasound image velocimetry.
      ). 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 MJP
      • Khir AW.
      Resolving the time lag between pressure and flow for the determination of local wave speed in elastic tubes and arteries.
      , 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 K
      • Rowland EM
      • Leow CH
      • Tang MX
      • Weinberg PD.
      Determining haemodynamic wall shear stress in the rabbit aorta in vivo using contrast-enhanced ultrasound image velocimetry.
      ,
      • Riemer K
      • Rowland EM
      • Broughton-Venner J
      • Leow CH
      • Tang M
      • Weinberg PD.
      Contrast agent-free assessment of blood flow and wall shear stress in the rabbit aorta using ultrasound image velocimetry.
      ). 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é C
      • Deffieux T
      • Pernot M
      • Osmanski BF
      • Biran V
      • Gennisson JL
      • Sieu LA
      • Bergel A
      • Franqui S
      • Correas JM
      • Cohen I
      • Baud O
      • Tanter M.
      Spatiotemporal clutter filtering of ultrafast ultrasound data highly increases Doppler and fUltrasound sensitivity.
      ). 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(xz,t)=i=1Ntλiui(xz)vi(t)T
      (1)


      where λi are the ordered singular values of S, and ui and vi 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(xz,t)=Stissue+Sblood+Snoise
      (2)


      S(xz,t)=i=1Th1λiui(xz)vi(t)T+i=Th1+1Th2λiui(xz)vi(t)T+i=Th2+1Ntλiui(xz)vi(t)T
      (3)


      where Th1 and Th2 are singular value thresholds. These were selected manually, based on the visual appearance of ui and the frequency of the temporal vectors vi (
      • Song P
      • Manduca A
      • Trzasko JD
      • Chen S.
      Ultrasound small vessel imaging with block-wise adaptive local clutter filtering.
      ). 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 Sblood. 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(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)
      r(δx,δz)=F1{A(k)(δx,δz)·B(k+1)*(δx,δz)}
      (4)


      where A and B are the Fourier transforms of the interrogation windows, δx and δz are pixel shifts in the lateral and axial directions, * is the complex conjugate and F1 is the inverse Fourier transform. Conventionally, the displacement (Δx,Δ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 J.
      On velocity gradients in PIV interrogation.
      ). 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 R.
      Double exposure, multiple-field particle image velocimetry for turbulent probability density.
      ;
      • Adrian RJ
      • Westerweel J.
      Particle image velocimetry.
      ). Finding the centroid rather than the peak was used to obtain the mean (
      • Poelma C
      • Westerweel J.
      Generalized displacement estimation for averages of non-stationary flows.
      ).
      Assuming negligible mean displacement of blood in the z direction, the centroid in the x direction was computed as
      Δxc=δxr(δx,0)δx/δxr(δx,0)
      (5)


      where for a window size of N, the summation ranges from N/2 to N/21 and δ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 (σ) was measured from the autocorrelation map (24 pixels, as measured using the 1/e2 rule [
      • Poelma C
      • Westerweel J.
      Generalized displacement estimation for averages of non-stationary flows.
      ]). 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 σ/2 pixels to the maximum displacement (Δxmax) +σ/2. To determine Δxmax, the modal displacement was found. Finally, velocity Uraw 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 D
      • Strykowski P
      • Debatin K.
      Bias and precision errors of digital particle image velocimetry.
      ) but is not an issue if σ 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 PD
      • Ethier CR.
      Twenty-fold difference in hemodynamic wall shear stress between murine and human aortas.
      ]); 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 (Uraw) 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 Uraw and Womersley-based 3-D velocity UWo 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 UWo 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 Stissue. 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
      dI=dPdtdUdt
      (6)


      ndI=dDdtdUdt
      (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 dt2 to make the magnitude of dI independent of the sampling frequency (
      • 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.
      ).
      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
      cPU=dPρdU
      (8)


      (
      • Khir A
      • O'Brien A
      • Gibbs J
      • Parker K.
      Determination of wave speed and wave separation in the arteries.
      ;
      • Davies JE
      • Whinnett ZI
      • Francis DP
      • Willson K
      • Foale RA
      • Malik IS
      • Hughes AD
      • Parker KH
      • Mayet J.
      Use of simultaneous pressure and velocity measurements to estimate arterial wave speed at a single site in humans.
      ) or
      Cln(D)U=12dUdln(D)
      (9)


      (
      • Feng J
      • Khir AW.
      Determination of wave speed and wave separation in the arteries using diameter and velocity.
      ), where ρ is the density of the blood (1044 kg/m3). 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 R
      • Beare R
      • Willemet M
      • Alastruey J
      • Smolich JJ
      • Cheung MMH
      • Mynard JP.
      Robust and practical non-invasive estimation of local arterial wave speed and mean blood velocity waveforms.
      )
      Cln(D)P=12ρ0.5dPdln(D)
      (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.
      Fig 1
      Fig. 1(a) Pressure and (b) diameter waveforms, and (c) the relation between them for a single rabbit (POEK). Dashed lines represent the ensemble average of several heartbeats for each repeat, and the solid line is the average of all (n = 3) repeats. In (c), mean pressure and mean diameter have been subtracted from instantaneous values. Time during the cardiac cycle runs clockwise around the loop. Data for all rabbits are given in Figure S6 (online only).
      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, Ucl. (For the same reason, Doppler-derived velocities can be unsuitable when calculating absolute wave speeds; see
      • Mynard JP
      • Kondiboyina A
      • Clarke MM
      • Kowalski R
      • Cheung MMH
      • Smolich JJ.
      Noninvasive assessment of carotid arterial wave speed and distensibility.
      ). This reduced the spread of velocities within the window, so pixel displacements could be determined using simple Gaussian peak fitting.
      Fig 2
      Fig. 2(a) Waveform of centerline blood velocity (Ucl) 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 (nUcl, 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: Uraw (black line) was obtained assuming 2-D Poiseuille flow, UWo (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 Ugauss 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 B-mode 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.
      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: Uraw was obtained assuming 2-D Poiseuille flow, UWo by assuming a 3-D Womersley flow (in both cases the displacement being estimated from the centroid of the cross-correlation peak), and Ugauss 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.
      Fig 3
      Fig. 3(a) PU, (b) ln(D)U and (c) ln(D)P loops for a single rabbit (POEK). Data for all rabbits are given in 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.
      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.
      Fig 4
      Fig. 4Wave 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.
      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.
      Fig 5
      Fig. 5Wave intensity in a single rabbit (PFLH) during administration of esmolol, determined by the (a) invasive PU method and (b) non-invasive DU method. The x-axis indicates a single cardiac cycle, while the y-axis represents the 30-min duration of the experiment. Esmolol was administered at increasing doses immediately after recordings were made at 0, 10 and 20 min (dashed white lines). Wave intensity is represented by color (left) or color and height (right).
      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.
      Fig 6
      Fig. 6Effect 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.
      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.
      Fig 7
      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.
      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.
      Fig 8
      Fig. 8Carotid artery wave intensity measured in a healthy human participant. In all panels, dashed lines represent the ensemble average for each repeat, and the solid line represents the average of all repeats. (a) Diameter. (b) Velocity calculated assuming 2-D Poiseuille flow (Uraw, black line) and 3-D Womersley flow (UWo, crosses). (c) ln(D)U loop, with wave speed computed for the first 50 ms of ejection (blue). (d) Net wave intensity.

      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 UW
      • Song P
      • Trzasko JD
      • Daigle R
      • Borisch EA
      • Huang C
      • Gong P
      • Tang S
      • Ling W
      • Chen S.
      Real time SVD-based clutter filtering using randomized singular value decomposition and spatial downsampling for micro-vessel imaging on a Verasonics ultrasound system.
      ). Threshold selection could be automated by identifying turning points in the singular value magnitude or temporal vector frequency plots (
      • Song P
      • Manduca A
      • Trzasko JD
      • Chen S.
      Ultrasound small vessel imaging with block-wise adaptive local clutter filtering.
      ) or by identifying tissue and blood subspaces from the similarity matrix of the spatial vectors (
      • Baranger J
      • Arnal B
      • Perren F
      • Baud O
      • Tanter M
      • Demene C.
      Adaptive spatiotemporal SVD clutter filtering for ultrafast Doppler imaging using similarity of spatial singular vectors.
      ). We found that precise selection is not critical: the lower threshold (Th1) could be changed between 115 and 150 and the upper threshold (Th2) 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 SI
      • Baerum S
      • Sørhus V
      • Torp H.
      Ultrasound-based vessel wall tracking: An auto-correlation technique with RF center frequency estimation.
      ). 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 JE
      • Whinnett ZI
      • Francis DP
      • Willson K
      • Foale RA
      • Malik IS
      • Hughes AD
      • Parker KH
      • Mayet J.
      Use of simultaneous pressure and velocity measurements to estimate arterial wave speed at a single site in humans.
      ). 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 J.
      Numerical assessment of time-domain methods for the estimation of local arterial pulse wave speed.
      ;
      • Swillens A
      • Taelman L
      • Degroote J
      • Vierendeels J
      • Segers P.
      Comparison of non-invasive methods for measurement of local pulse wave velocity using FSI-simulations and in vivo data.
      ;
      • Borlotti A
      • Li Y
      • Parker KH
      • Khir AW.
      Experimental evaluation of local wave speed in the presence of reflected waves.
      ;
      • Segers P
      • Swillens A
      • Taelman L
      • Vierendeels J.
      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.
      ;
      • Kowalski R
      • Beare R
      • Willemet M
      • Alastruey J
      • Smolich JJ
      • Cheung MMH
      • Mynard JP.
      Robust and practical non-invasive estimation of local arterial wave speed and mean blood velocity waveforms.
      ;
      • Campos Arias D
      • Stergiopulos N
      • Rodríguez Moliner T
      • Segers P
      Mapping the site-specific accuracy of loop-based local pulse wave velocity estimation and reflection magnitude: A 1D arterial network model analysis.
      ). Comparing data obtained using the two methods with data derived from the ln(D)P loop—thought to be less affected by reflections (
      • Kowalski R
      • Beare R
      • Willemet M
      • Alastruey J
      • Smolich JJ
      • Cheung MMH
      • Mynard JP.
      Robust and practical non-invasive estimation of local arterial wave speed and mean blood velocity waveforms.
      )—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 Y
      • Gallagher G
      • Plekhanov S
      • Morita S
      • Brady PW
      • Hunyor SN.
      HeartPatch implanted direct cardiac compression: Effect on coronary flow and flow patterns in acute heart failure sheep.
      ;
      • Gallagher GL
      • Huang Y
      • Morita S
      • Zielinski RR
      • Hunyor SN.
      Efficacy and mechanisms of biventricular and left/right direct cardiac compression in acute heart failure sheep.
      ). 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 W2 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 SL
      • Zambanini A
      • Mayet J
      • McG Thom SA
      • Foale R
      • Parker KH
      • Hughes AD
      Reduced systolic wave generation and increased peripheral wave reflection in chronic heart failure.
      ). 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 RM
      • Sherwin SJ
      • Tang MX
      • Weinberg PD.
      Wave intensity analysis combined with machine learning can detect impaired stroke volume in simulations of heart failure.
      ) 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

        • Adrian R.
        Double exposure, multiple-field particle image velocimetry for turbulent probability density.
        Opt Laser Eng. 1998; 9: 211-228
        • Adrian RJ
        • Westerweel J.
        Particle image velocimetry.
        Cambridge University Press, New York, NY2011
        • Alastruey J.
        Numerical assessment of time-domain methods for the estimation of local arterial pulse wave speed.
        J Biomech. 2011; 44: 885-891
        • Baranger J
        • Arnal B
        • Perren F
        • Baud O
        • Tanter M
        • Demene C.
        Adaptive spatiotemporal SVD clutter filtering for ultrafast Doppler imaging using similarity of spatial singular vectors.
        IEEE Trans Med Imaging. 2018; 37: 1574-1586
        • Biglino G
        • Steeden JA
        • Baker C
        • Schievano S
        • Taylor AM
        • Parker KH
        • Muthurangu V.
        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
        • 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
        • Borlotti A
        • Li Y
        • Parker KH
        • Khir AW.
        Experimental evaluation of local wave speed in the presence of reflected waves.
        J Biomech. 2014; 47: 87-95
        • Campos Arias D
        • Stergiopulos N
        • Rodríguez Moliner T
        • Segers P
        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
        • Curtis SL
        • Zambanini A
        • Mayet J
        • McG Thom SA
        • Foale R
        • Parker KH
        • Hughes AD
        Reduced systolic wave generation and increased peripheral wave reflection in chronic heart failure.
        Am J Physiol Heart Circ Physiol. 2007; 293: H557-H562
        • Davies JE
        • Whinnett ZI
        • Francis DP
        • Willson K
        • Foale RA
        • Malik IS
        • Hughes AD
        • Parker KH
        • Mayet J.
        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
        • Demené C
        • Deffieux T
        • Pernot M
        • Osmanski BF
        • Biran V
        • Gennisson JL
        • Sieu LA
        • Bergel A
        • Franqui S
        • Correas JM
        • Cohen I
        • Baud O
        • Tanter M.
        Spatiotemporal clutter filtering of ultrafast ultrasound data highly increases Doppler and fUltrasound sensitivity.
        IEEE Trans Med Imaging. 2015; 34: 2271-2285
        • Feng J
        • Khir AW.
        Determination of wave speed and wave separation in the arteries using diameter and velocity.
        J Biomech. 2010; 43: 455-462
        • Forliti D
        • Strykowski P
        • Debatin K.
        Bias and precision errors of digital particle image velocimetry.
        Exp Fluids. 2000; 28: 436-447
        • Frank O.
        Die Elastizität der Blutegefässe.
        Z Biol. 1920; 71: 255-272
        • Gallagher GL
        • Huang Y
        • Morita S
        • Zielinski RR
        • Hunyor SN.
        Efficacy and mechanisms of biventricular and left/right direct cardiac compression in acute heart failure sheep.
        Artif Organs. 2007; 31: 39-44
        • Hoskins PR.
        A review of the measurement of blood velocity and related quantities using Doppler ultrasound.
        Proc Inst Mech Eng H. 1999; 213: 391-400
        • Huang Y
        • Gallagher G
        • Plekhanov S
        • Morita S
        • Brady PW
        • Hunyor SN.
        HeartPatch implanted direct cardiac compression: Effect on coronary flow and flow patterns in acute heart failure sheep.
        ASAIO J. 2003; 49: 309-313
        • Kang J
        • Aghilinejad A
        • Pahlevan NM.
        On the accuracy of displacement-based wave intensity analysis: Effect of vessel wall viscoelasticity and nonlinearity.
        PLoS One. 2019; 14e0224390
        • Khir A
        • O'Brien A
        • Gibbs J
        • Parker K.
        Determination of wave speed and wave separation in the arteries.
        J Biomech. 2001; 34: 1145-1155
        • Kowalski R
        • Beare R
        • Willemet M
        • Alastruey J
        • Smolich JJ
        • Cheung MMH
        • Mynard JP.
        Robust and practical non-invasive estimation of local arterial wave speed and mean blood velocity waveforms.
        Physiol Meas. 2017; 38: 2081-2099
        • Kowalski R
        • Mynard JP
        • Smolich JJ
        • Cheung MMH.
        Comparison of invasive and non-invasive aortic wave intensity and wave power analyses in sheep.
        Physiol Meas. 2019; 40015005
        • Li Y
        • Hickson SS
        • McEniery CM
        • Wilkinson IB
        • Khir AW.
        Stiffening and ventricular–arterial interaction in the ascending aorta using MRI: Ageing effects in healthy humans.
        J Hypertens. 2019; 37: 347-355
        • Lok UW
        • Song P
        • Trzasko JD
        • Daigle R
        • Borisch EA
        • Huang C
        • Gong P
        • Tang S
        • Ling W
        • Chen S.
        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
        • Montaldo G
        • Tanter M
        • Bercoff J
        • Benech N
        • Fink M.
        Coherent plane-wave compounding for very high frame rate ultrasonography and transient elastography.
        IEEE Trans Ultrason Ferroelectr Freq Control. 2009; 56: 489-506
        • Mynard JP
        • Kondiboyina A
        • Clarke MM
        • Kowalski R
        • Cheung MMH
        • Smolich JJ.
        Noninvasive assessment of carotid arterial wave speed and distensibility.
        Am J Physiol Heart Circ Physiol. 2018; 315: H1495
        • Nichols WW
        • O'Rourke MF.
        McDonald's blood flow in arteries: Theoretical, experimental and clinical principles.
        Edward Arnold, London, UK1998
        • 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
        • Parker KH.
        An introduction to wave intensity analysis.
        Med Biol Eng Comput. 2009; 47: 175-188
        • Parker KH
        • Jones CJ.
        Forward and backward running waves in the arteries: Analysis using the method of characteristics.
        J Biomech Eng. 1990; 112: 322-326
        • Parker KH
        • Jones CJ
        • Dawson JR
        • Gibson DG.
        What stops the flow of blood from the heart?.
        Heart Vessels. 1988; 4: 241-245
        • Poelma C
        • Westerweel J.
        Generalized displacement estimation for averages of non-stationary flows.
        Exp Fluids. 2011; 50: 1421-1427
        • Pomella N
        • Wilhelm EN
        • Kolyva C
        • González-Alonso J
        • Rakobowchuk M
        • Khir AW.
        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
        • Pomella N
        • Wilhelm EN
        • Kolyva C
        • González-Alonso J
        • Rakobowchuk M
        • Khir AW.
        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
        • Rabben SI
        • Baerum S
        • Sørhus V
        • Torp H.
        Ultrasound-based vessel wall tracking: An auto-correlation technique with RF center frequency estimation.
        Ultrasound Med Biol. 2002; 28: 507-517
        • Reavette RM
        • Sherwin SJ
        • Tang M
        • Weinberg PD.
        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
        • Reavette RM
        • Sherwin SJ
        • Tang MX
        • Weinberg PD.
        Wave intensity analysis combined with machine learning can detect impaired stroke volume in simulations of heart failure.
        Front Bioeng Biotechnol. 2021; 9737055
        • Riemer K
        • Rowland EM
        • Leow CH
        • Tang MX
        • Weinberg PD.
        Determining haemodynamic wall shear stress in the rabbit aorta in vivo using contrast-enhanced ultrasound image velocimetry.
        Ann Biomed Eng. 2020; 48: 1728-1739
        • Riemer K
        • Rowland EM
        • Broughton-Venner J
        • Leow CH
        • Tang M
        • Weinberg PD.
        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
        • Segers P
        • Verdonck P.
        Role of tapering in aortic wave reflection: hydraulic and mathematical model study.
        J Biomech. 2000; 33: 299-306
        • Segers P
        • Swillens A
        • Taelman L
        • Vierendeels J.
        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
        • Song P
        • Manduca A
        • Trzasko JD
        • Chen S.
        Ultrasound small vessel imaging with block-wise adaptive local clutter filtering.
        IEEE Trans Med Imaging. 2017; 36: 251-262
        • Sugawara M
        • Niki K
        • Ohte N
        • Okada T
        • Harada A.
        Clinical usefulness of wave intensity analysis.
        Med Biol Eng Comput. 2009; 47: 197-206
        • Swalen MJP
        • Khir AW.
        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
        • Swillens A
        • Taelman L
        • Degroote J
        • Vierendeels J
        • Segers P.
        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
        • Wang JJ
        • Bouwmeester JC
        • Belenkie I
        • Shrive NG
        • Tyberg JV.
        Alterations in aortic wave reflection with vasodilation and vasoconstriction in anaesthetized dogs.
        Can J Cardiol. 2013; 29: 243-253
        • Weinberg PD
        • Habens F
        • Kengatharan M
        • Barnes SE
        • Matz J
        • Anggård EE
        • Carrier MJ.
        Characteristics of the pulse waveform during altered nitric oxide synthesis in the rabbit.
        Br J Pharmacol. 2001; 133: 361-370
        • Weinberg PD
        • Ethier CR.
        Twenty-fold difference in hemodynamic wall shear stress between murine and human aortas.
        J Biomech. 2007; 40: 1594-1598
        • Westerweel J.
        On velocity gradients in PIV interrogation.
        Exp Fluids. 2008; 44: 831-842