Shear wave elastography plaque characterization with mechanical testing validation: a phantom study

Determining plaque vulnerability is critical when selecting the most suitable treatment for patients with atherosclerotic plaque. Currently, clinical non-invasive ultrasound-based methods for plaque characterization are limited to visual assessment of plaque morphology and new quantitative methods are needed. In this study, shear wave elastography (SWE) was used to characterize hard and soft plaque mimicking inclusions in six common carotid artery phantoms by using phase velocity analysis in static and dynamic environments. The results were validated with mechanical tensile testing. In the static environment, SWE measured a mean shear modulus of 5.8  ±  0.3 kPa and 106.2  ±  17.2 kPa versus 3.3  ±  0.5 kPa and 98.3  ±  3.4 kPa measured by mechanical testing in the soft and hard plaques respectively. Furthermore, it was possible to measure the plaques’ shear moduli throughout a simulated cardiac cycle. The results show good agreement between SWE and mechanical testing and indicate the possibility for in vivo arterial plaque characterization using SWE.

Keywords: shear wave elastography, plaque characterization, atherosclerosis, carotid artery, ultrasound, phantom, mechanical testing (Some figures may appear in colour only in the online journal)

Introduction
It is estimated that cardiovascular diseases (CVDs) are responsible for over 17.3 million deaths per year and are the leading causes of death in the world (World Health Organization 2011). CVDs due to atherosclerosis include ischemic heart disease or coronary artery disease, cerebrovascular disease (stroke), and arterial diseases such as hypertension and peripheral vascular disease. Atherosclerotic disease of the carotid arteries has been associated with the development of cerebrovascular events due to rupture of carotid plaques (Salem et al 2013). The disease is influenced by several risk factors such as diet, age, genetics, smoking habits, and alcohol consumption where the slow progression of thrombotic atheroma plaque formation can take up to 20 years (Libby 2001). Carotid artery plaque formation is a complex process that does not follow a single pathway, making it possible to develop several different types of plaques. Moreover, plaque formation evolves through the progression of atherosclerosis from fatty necrotic, to mixed, to fibrous core (Libby 2001). Physicians are interested in identifying vulnerable plaques, which often have a necrotic core and a thin fibrous cap (Naghavi et al 2003b), as they are most likely to result in thrombosis. Furthermore, neovascularization, found in advanced atherosclerotic lesions, is considered an element of vulnerability because microvessels are prone to rupture causing intraplaque haemorrhage, which has been shown to accelerate plaque formation (Kolodgie et al 2003).
In current clinical practice, ultrasound duplex scanning is typically the first level of screening for atherosclerotic plaques (Saba et al 2012). The degree of stenosis typically has been considered the parameter of choice to determine the therapeutic approach (Barnett et al 1998), but several investigations (Naghavi et al 2003a, 2003b, Schwarz et al 2013 have demonstrated that the degree of luminal stenosis is only an indirect indicator of the atherosclerotic process and that direct assessment of the plaque structure and composition may be key to predict the development of future cerebrovascular ischemic events. Moreover, it is critical to characterize plaques to determine the most suitable treatment (endarterectomy, angioplasty, or medication) for the patient.
Currently, clinical non-invasive ultrasound-based methods for plaque characterization are limited to visual assessment of plaque morphology, hypoechoic area, and echo reflection in the plaque (Gray-Weale et al 1988, Kakkos et al 2013. There have been quantitative plaque characterization attempts by grey-scale median (GSM) measurements of the plaque (Biasi et al 1999, Kanber et al 2013, Salem et al 2013 as well as using contrast-enhanced ultrasound imaging in combination with Tissue Doppler imaging (Deyama et al 2013), however both techniques suffer from various limitations. Thermal Strain Imaging (TSI) has also been used in vivo to detect lipids in atherosclerotic plaques in the femoral arteries of rabbits (Mahmoud et al 2013) but the technique is not ready for clinical use. A promising technique to characterize plaque is by measuring the plaque strain throughout the cardiac cycle via block matching speckle tracking techniques on B-mode images (Shi et al 2008, Widman et al 2015. Small successful in vivo pilot studies have been conducted, but more data is needed from a larger patient population. Several successful plaque characterization studies have been completed using intravascular ultrasound (IVUS) (de Korte et al 2000, Schaar et al 2003, Zhang et al 2011. It should be noted that non-invasive vascular elastography has also been performed with ultrasound (Hansen et al 2013, Naim et al 2013.
Another promising technique is acoustic radiation force impulse (ARFI) imaging, where ultrasound radiation force excites the tissue followed by measuring the tissue's axial displacement (<10 µm). Pilot studies to characterize carotid plaques have been performed (Dumont et al 2006, Behler et al 2009, Dahl et al 2009, Allen et al 2011, Doherty et al 2013, but ARFI is a relative displacement measurement, which limits the evaluation of plaque stiffness to surrounding proximal tissue. This implies that plaques on opposing anterior and posterior arterial walls cannot be compared to each other. Despite previous efforts, quantitative methods are still needed to characterize plaque without user-subjectivity influence. Shear wave elastography (SWE) is an emerging imaging modality that uses the shear modulus of tissue as the contrast mechanism in images. SWE uses ultrasound radiation force to non-invasively generate shear waves in the tissue and ultrasonic tracking methods measure the shear wave propagation speed, which is directly proportional to the material properties. To date, SWE has been applied in tissues such as breast (Tanter et al 2008), liver (Muller et al 2009), arterial wall (Couade et al 2010), kidney (Gennisson et al 2012), and muscle (Gennisson et al 2010) where it successfully detects tissue changes related to disease or physiological processes that may not be detected with B-mode ultrasound alone. For large uniform tissues such as breast or liver, the tissue's stiffness is estimated assuming that the generated shear wave propagates in an incompressible, homogeneous, isotropic, linear elastic material (Bercoff et al 2004). Plaques and arteries do not have these properties and therefore phase velocity analysis of the shear wave is required to accurately assess the stiffness of the material, which will be discussed in greater detail in the section 2: wave propagation in arteries.
In this paper, we attempted to use SWE to characterize hard and soft plaque phantoms in a static setup using phase velocity analysis of the propagating shear wave and validate our results against mechanical testing. The shear wave propagation through the plaque mimicking inclusion was modelled as a lamb wave traveling through a thin plate to which experimental phase velocity data was fitted to estimate the plaque's shear modulus. Additionally, we measured stiffness changes in the plaque phantom throughout a simulated cardiac cycle to demonstrate feasibility of an in vivo application. Finally, we compared the average phase velocity, group velocity, as well as the mean shear wave speed calculated by a commercial ultrasound machine to determine if phase velocity analysis is required to characterize plaque stiffness.

Wave propagation in arteries
Ultrasonic radiation force gives rise to two types of waves in the excited medium: compressional waves associated with the medium's bulk modulus; and shear waves correlated to the shear modulus of the medium (Sarvazyan et al 1998, Nightingale et al 2003. Compressional waves are in the MHz frequency range, are millimetres in wavelength, travel very fast through the medium, and there is little variation in the bulk modulus between different biological media. Shear waves on the other hand have a typical bandwidth <2 kHz, have wavelengths in centimetres, travel slowly through biological tissue, and the shear modulus varies greatly between different biological media. For elastography of breast and liver, an incompressible, homogeneous, isotropic, infinitely large medium is assumed where the shear wave speed is equal to the group velocity (c g ). With these assumptions, the shear modulus (µ) and Young's modulus (E) is where ρ is the mass density of the material. However, these assumptions do not hold for arterial and plaque applications. Arteries can be considered thin hollow cylinders, whereas arterial plaques can be of various geometrical configurations often a few millimetres in size. Hereafter, the plaque will be considered part of the arterial wall to simplify the analysis, which is a realistic geometric assumption. An analytical solution for a hollow cylinder with appropriate boundary conditions is needed to characterize the motion of the propagating shear wave in the artery. However, the equations and boundary conditions describing the shear wave motion through the hollow cylinder are complex and difficult to solve analytically, but many wave propagating motions in complex geometries have successfully been approximated by a fundamental antisymmetric Lamb wave mode traveling through an infinite plate (Couade et al 2010, Bernal et al 2011, Nenadic et al 2013. The approximation is appropriate when the phase velocity curve of the fundamental mode of the shear wave is similar in shape to that of a fundamental antisymmetric Lamb wave. In recent years, attempts have been made for arterial applications to incorporate theories of guided waves into SWE post-processing algorithms. For vascular tissue, Couade, et al (2010) used values from literature for vessel geometry, mass density, and shear wave speed to suggest that the cut-off frequencies of all non-zero-order modes was beyond the frequency spectra of an induced shear wave. This, in addition with the principal detection of antisymmetric wave modes, made Couade, et al propose a solution for the wave motion based on the zero-order antisymmetric (A 0 ) wave mode alone. The phase velocity (c p ) was expressed in terms of the frequency ( f ), plate thickness (h), and apparent shear wave velocity (c s ), as given by with the derived apparent shear wave velocity representing an average velocity of the travelling shear wave and a correction factor K derived from curve-fitting the shear wave phase velocity curve to an antisymmetric Lamb wave propagating in a plate. It was also suggested that for measurements on hollow cylindrical structures, the curve fitting procedure was to be performed solely on frequencies above 1 kHz, but this frequency limit was somewhat arbitrarily chosen based merely on empirical data indicating that the differences in the phase velocity curves between a wave propagating through a plate and a hollow cylinder will mainly be visible at low frequencies.
Other groups studying shear wave propagation in hollow cylinders have shown that it is possible to measure zero-order and first-order symmetric as well as antisymmetric shear waves in hollow cylinders by taking the first and second partial derivatives of the k-space with respect to frequency and wave number (Bernal et al 2011). Furthermore, Bernal et al provides a comprehensive derivation of a Lamb wave dispersion equation for a homogeneous elastic plate submerged in incompressible water like fluid which simplifies to where ρ, ω, c L , k L , k s , h, and µ are the mass density, angular frequency, frequency dependent Lamb wave velocity, Lamb wave number, shear wave number, half-plate thickness, and shear modulus respectively. Equation (4) can be fit to experimental data such that the shear modulus µ can be estimated.
It is important to note that the outlined work has been based on Lamb wave plate theory. So far, no analytical expression for a strict hollow cylindrical structure submerged in a non-viscous liquid has been attempted in SWE setups due to the complexity of the solutions (Zhang et al 2003). For hollow cylindrical structures, some initial trials have been made using the modified Moens-Korteweg (M-K) equation to directly relate mechanical properties of the hollow cylinder to the propagating wave velocity (Bernal et al 2009, Widman et al 2012. The modified M-K equation however represents a purely empirical expression and the results indicate a poor performance when compared to mechanical testing (Widman et al 2012).

Method
Carotid artery phantoms with soft and hard plaque mimicking inclusions were created, placed in an enclosure, and the shear modulus was estimated in the plaque and vessel wall using SWE in a static setup. The shear wave propagation through the plaque mimicking inclusion was modelled as a lamb wave traveling through a thin plate to which experimental phase velocity data was fitted to estimate the plaque's shear modulus. The shear moduli were compared to values calculated from tensile mechanical measurements. Thereafter, the phantom was connected to a programmable pump where a dynamic test was performed and the shear modulus was measured throughout the simulated cardiac cycle using SWE.

Phantom construction
Six in vitro carotid artery phantoms with a plaque mimicking inclusion (three soft and three hard) were constructed from a mixture (w/w) of 87% deionized water, 10% polyvinyl alcohol (PVA) with a molecular weight of 56.140 g mol −1 (Sigma-Aldrich, St Louis, MO, USA) and 3% graphite powder with particle size < 50 μm (Merck KGaA, Darmstadt, Germany). The solution was heated to 50 °C and stirred until the mixture thickened and was fully dissolved. It was subsequently poured in the vessel phantom mould (figure 1(b)), which consisted of a hollow acrylic block (cylindrical shape with a 12 mm diameter, 100 mm length, and moulded fixing collars) and a bronze rod with a 6 mm diameter placed in the centre of the hole. The metallic rod was 50 mm long and had a 1.5 mm thick extension attached to create a cavity in the phantom wall (figure 1(a)). Once the solution was poured in the mould, it stayed confined between the acrylic block and the rod, resulting in a hollow cylindrical-shaped phantom. The phantom was frozen for 24 h at approximately −23 °C and then thawed for 24 h at room temperature (≈22 °C). This freeze-thaw cycle was repeated twice for three phantoms to construct phantoms with soft plaque mimicking inclusions. Subsequently, the phantom wall cavity was filled with a PVA-graphite solution, prepared the same way as previously described, and the entire phantom was exposed to one extra freeze-thaw cycle, resulting in three freeze-thaw cycles for the vessel wall and one for the cavity mimicking the plaque. Since PVA stiffens with increasing freeze-thaw cycles, the phantom wall became stiffer than the inclusion, mimicking the vessel wall and the soft plaque, respectively. A photo of the acrylic mould and rod, a schematic representation of the vessel phantom, and B-mode long-axis image of the soft plaque phantom is shown in figure 1. As the plaque cavity was filled with PVA, the thin vessel wall (0.75 mm) extruded into the lumen creating a partial lumen occlusion shown in figure 1(c).
Three phantoms with hard plaque mimicking inclusions, as shown in figure 2(a), were created by pouring the same PVA mixture into an acrylic mould to create small plaque inclusions (partial cylindrical shape with 120° circumference, 50 mm length, and 1.5 mm thickness), as shown in figure 2(b). The mould underwent seven freeze-thaw cycles creating hard plaque inclusions. The inclusions were placed in the acrylic vessel phantom mould (figure 1(b)) and filled with a new batch of PVA mixture, which encompassed the plaque inclusions. The mould was exposed to three freeze-thaw cycles, resulting in three freeze-thaw cycles for the vessel wall and ten for the hard plaque inclusion. A long-axis B-mode image of the hard plaque phantom can be seen in figure 2(c).
Two control vessel phantoms without plaque inclusions were created for both the soft and hard plaque phantoms by using the acrylic vessel phantom mould with a rod without an extrusion (creating a symmetrical vessel phantom without a plaque inclusion). The phantoms were created for mechanical testing which assumed geometric symmetry and material incompressibility in order to calculate the shear modulus. For validation of the soft plaque and vessel wall, the control phantoms had one and three freeze-thaw cycles respectively and were created from the same PVA batch as the soft plaque phantoms. Likewise, the control vessel phantoms for the hard plaque and vessel wall had ten and three freeze-thaw cycles respectively and were also made from the same batch of PVA as the hard plaque phantoms.

Static experimental setup
A static test was performed so that the conditions from the SWE test could be reproduced in the mechanical testing setup, making it possible to compare shear moduli values by SWE and mechanical testing. The PVA phantoms were attached to a cylindrical acrylic plastic enclosure (100 mm length, 140 mm diameter) by placing two plastic discs around the fixing collars of the phantoms and tightening them to the enclosure with screws (figure 3). The enclosure had an adjustable short end that allowed the length of the cavity to be adjusted, which was used to pre-stretch the phantom to avoid curvature of the vessel wall when attached to the enclosure. The phantoms were surrounded by water and orientated such that the plaque was positioned in the anterior wall with its cross-sectional centre facing the transducer. A thin sheet (5 mm thickness) of rubber was placed at the bottom of the enclosure to reduce ultrasonic reflections from the acrylic plastic.
The phantoms were connected to two hoses with a 6 mm inner diameter which was subsequently filled with water. No air bubbles were trapped in the phantom lumen which was verified by using B-mode ultrasound scanning. Thereafter, the movable wall of the enclosure was adjusted 8 mm to remove bulging from the water filled phantoms.

Data acquisition
Ultrasound in-phase quadrature (IQ) data were collected on an Aixplorer ultrasound machine (Supersonic Imagine, Aix-en-Provence, France) with a custom shear wave research package using a fixed SL15-4 linear array transducer. The transducer was placed 1 cm above the plaque inclusion along the long-axis of the phantom (figure 3(a)) with a region of interest (ROI) placed on the anterior wall such that the second ultrasound pushing beam in the supersonic push sequence consisting of three pushes (Bercoff et al 2004) was located in the centre of the plaque. The pushing beam centre frequency was 6 MHz and the burst duration for each pushing beam was 150 µs, resulting in a total burst time of 450 µs. After exciting the phantom  vessel wall, the machine switched to ultrafast plane wave imaging to track the shear wave motion at a frame rate of 8 kHz at a centre frequency of 7.5 MHz with no angle compounding at an imaging depth of 3 cm. Similarly, when collecting data from the vessel wall, the ROI was moved such that the push was located in the anterior wall 2.5 cm from the end of the plaque with the ROI covering the vessel wall. Three consecutive data acquisitions, each consisting of a supersonic push sequence, spaced 100 ms apart were performed at each location. Finally, the diameter of the plaque and wall was measured from the B-mode image (centre frequency 6.63 MHz) to later be used as input data for the post-processing algorithm.

Data post processing
The unprocessed IQ data was exported from the Aixplorer ultrasound machine and postprocessed in MATLAB 2012b (MathWorks, Natick, MA, USA). For each pushing location, the three data acquisitions were post-processed with the following algorithm and the estimated shear moduli were averaged together. A flow-chart of the post processing algorithm can be seen in figure 4. The IQ data was divided such that the left and right shear wave could be studied independently. The left shear wave was analysed first followed by the right wave and the results were averaged together. First, a 2-D autocorrelator (Loupas et al 1995) measured the plaque axial displacement by finding the phase-shifts between consecutive pulse-echo measurements. The plaque motion was measured along five and three lateral pixel lines in the centre of the plaque and averaged together for the soft (figure 5) and hard plaque respectively. The hard plaques where thinner (approximately 1.9 mm in cross-sectional diameter) than the soft plaques (approximately 5 mm in cross-sectional diameter) than the soft plaques (approximately 5 mm in cross-sectional diameter), hence fewer lines were selected to analyse the shear wave motion in order to avoid wave reflection effects from the plaque boundaries. Thereafter, an axial displacement colour map was created displaying the shear wave propagation as time (ms) versus distance (mm). Once the shear wave had fully dissipated, data along the time axis in the spatiotemporal matrix were cropped out to remove late wall reflections and noise from the signal. Similarly, the analysis was performed on the vessel wall where displacement was measured along five lateral pixel lines in the centre of the wall and averaged together. Second, the axial displacement signal was converted from the time domain to the frequency domain by calculating the 2D fast Fourier transform of the signal where d x t ( , ) z is the axial displacement with respect to the transducer (z direction) as a function of distance (x) and time (t). Plotting H k f ( , ) is often referred to as the k-space representation of the propagating shear wave, where frequency ( f ) is represented on the x-axis and the wave number ) on the y-axis. The k-space plot was cropped to the upper right (right propagating wave) or lower right (left propagating wave) quadrant depending on the direction of propagation for the analysed wave. Thereafter, noise was removed by a threshold filter set to retain 90% of the maximum signal intensity.
Third, a phase velocity (c p ) map was created by multiplying each paired frequency ( f ) and wave number (λ = k 1/ L or 1/y-axis) for every position in the 2D k-space velocity data, cut-off frequencies were introduced at each end of the frequency spectrum. The high cut-off frequency was defined as the frequency for which the intensity has decreased to less than 30% of the maximum k-space intensity, whereas the low cut-off frequency was  defined as the frequency for which the relative change in phase velocity between neighbouring data points were greater than 50%. The low cut-off frequency removed high-intensity background noise signals near 0 Hz.
Fourth, a mathematical model for an antisymmetric Lamb wave propagating in a plate (4) was fitted to the phase velocity curve using a least squares fit (Couade et al 2010, Bernal et al 2011 with the shear modulus (µ) as the only varying parameter with other dependent variables such as wall thickness, mass density, and speed of sound given as input to the curve fit. The mass density and speed of sound of PVA were measured in a pre-study while the plaque wall thickness was measured from B-mode images using the Aixplorer ultrasound machine. Following the reasoning of Couade et al (2010), a second shear modulus was calculated using the same least squares fit procedure above the half-bandwidth of the phase velocity spectrum as cylindrical-and plate-based dispersion curves converge at high frequencies and deviate at low frequencies. The half-bandwidth criterion was selected in order to retain as much valuable dispersion data as possible.
To evaluate the quality of the curve fits, the root mean square error (RMSE) between the shear wave phase velocity curve and Lamb wave curve fit was calculated as where c (n) p is the shear wave phase velocity, c n ( ) L is the Lamb wave propagation speed from the curve fit, and N is the number of data samples. The RMSE was calculated for both the full spectrum and upper half-bandwidth curve fits.

Comparison of shear wave speed and shear moduli
A comparison of shear wave speeds and shear moduli using the phase velocity (c ) p , group velocity (c g ), and mean velocity from the Aixplorer ultrasound machine colour maps was made. The phase velocity was approximated once the shear wave phase velocity curve had plateaued by calculating the mean phase velocity on the upper half-bandwidth of the shear wave spectrum. Thereafter, the frequency the mean phase velocity occurred at was graphically interpreted and defined as the phase velocity frequency.
The group velocity (the speed of the envelope of the propagating wave) was calculated using the Radon transformation on the spatiotemporal shear wave data. The Radon transformation sums the axial displacement in the spatiotemporal matrix along line trajectories rotating around the centre of the matrix and identifies the trajectory which best characterizes the wave motion, which is identified by the peak Radon sum. The group velocity can be calculated from the slope of the peak Radon sum wave trajectory (Rouze et al 2010, Hah et al 2012, Urban and Greenleaf 2012. The shear moduli calculated from the group velocity was found by using (1).
The mean velocity from the Aixplorer ultrasound machine SWE colour maps was calculated using the machine's Q-box function, where the mean shear wave speed was calculated in a user-defined ROI. Two ROIs were placed equidistance (approximately 5 mm from the push location) on either side of the supersonic push line, ensuring that the ROI only covered the plaque or vessel wall depending on the measurement location. The means of the left and right shear waves were averaged together. The shear moduli calculated from the mean velocity from the Aixplorer system was found by using (1).

Mechanical testing
The shear modulus estimated by phase velocity analysis was validated by mechanical testing. The control vessel phantoms without plaque inclusions with the same number of freeze-thaw cycles as the plaques and vessel walls were used for the validation test. The American Society for Testing and Materials suggest the use of angular twist for a hollow cylindrical sample when evaluating elastic shear properties (ASTM 2013). However, this test method is developed for structural materials and is not suitable for soft tissue or rubber samples as twisting may induce torsional buckling complicating the assessment of mechanical properties. Hence, a tensile test was developed that would mimic the loading conditions imposed on the phantom from the static shear wave experiment.
A fixture was created (figure 6(a)) which allowed the phantoms to be simultaneously filled with water while a tensile test was performed by an Instron 5567 (Instron Worldwide, Norwood, MA, USA) using a 100 N (calibrated 0.05% maximum error) electro-mechanical single range load cell. The phantoms were attached to the fixture by placing two plastic discs around the fixing collars of the phantoms and tightening them to the end fixtures with screws followed by filling the phantom with water until an equilibrium pressure state was reached as in the static experimental setup. The phantom was pre-stretched the same distance as in the static test (8 mm) to remove slack in the system and to mimic the test conditions from the static SWE experiment. The phantoms were axially loaded up to 10% strain at a rate of 10 mm min −1 . The strain rate was empirically selected from initial testing where the results indicated a strain rate independence which was in agreement with previous PVA mechanical characterization studies (Stammen et al 2001). The 10% strain kept the Young's modulus response to a finite linear region of the hyper-elastic material. The test was repeated 3 times and performed blinded by a scientist who did not partake in the static experiment. The applied axial force (F z ) and deformation (δ z ) were exported to MATLAB 2012b. Using the dimensions of each sample, the stretch (λ z ) and applied Cauchy stress (σ z ) were calculated as where L i and A i are the initial length and cross-sectional area of the sample. Several models were considered for calculating the Young's modulus. A pre-study found that non-linear models were not required and a local average least-squares fit linearization was selected. According to Hooke's law, the axial stress of a homogeneous, linear material may be expressed as where E is the material's linear elastic modulus and ε the applied Cauchy strain (ε λ = −1 z ). A mean linear elastic modulus (E ) was calculated using thirty equidistant sub-segments along the strain axis of the stress-strain curve. For each sub-segment, an elastic modulus (E i ) was derived by a local least-squares fit linearization (Fromageau et al 2003). The shear modulus was calculated as assuming a perfectly incompressible material (Poisson's ratio = v 0.5). This assumption was believed to be valid as the Poisson's ratio for the PVA material was measured in a pre-study to be = ± v 0.48 0.2 for a PVA sample with 3 freeze-thaw cycles.

Dynamic experimental setup
For the dynamic experiment two phantoms, each with a hard or soft plaque inclusion, were connected to a programmable pump (CompuFlow 1000MR, Shelley Medical Imaging Technologies, Ontario, Canada) by the previously described hoses and plastic discs. The pump was programmed to simulate a carotid flow profile at a rate of 72 cycles/min. A bloodmimicking solution of 40% glycerine (Merck KGaA, Darmstadt, Germany) and 60% deionized water mixture was pumped through the phantom at a peak flow rate of 30 mL s −1 . The pump was purged before the experiments to remove air bubbles from the solution. The interlumen pressure was measured with a pressure sensor element (Samba Preclin 420 LP/Samba 202 Pressure Measurement System, Samba Sensors, Västra Frölunda, Sweden) connected to an optic fibre which was placed approximately 2 cm from the plaque in the phantom lumen. The pressure measurement system measured pressure relative to the atmospheric pressure, and was calibrated prior to insertion in the phantom.
The same pushing and image acquisition parameters were used as in the static experiment with the exception of the number of data acquisitions. Eight consecutive data acquisitions spaced 100 ms apart spanning the simulated cardiac cycle (period = 828 ms) were performed in the phantom plaques. The acquisition sequence was triggered by a simulated electrocardiogram (ECG) signal sent from the pump to the ultrasound machine triggering on the R-wave of the signal. Pressure and ultrasound data were acquired simultaneously, although it was not possible to synchronize the pressure measurement system with the simulated ECG signal. The location and placement of the ROI on the plaque was the same as in the static experiment.

Data analysis
For the dynamic data analysis, a similar post-processing procedure was employed as in the static experiment with a few minor variations. For the soft plaque phantom, the plaque thickness varied throughout the simulated cardiac cycle due to radial strain forces (along the vessel radius) from the pressure wave induced by the programmable pump. To account for this variation, the plaque and wall diameter were manually measured using in-house software from the B-mode image for every pushing sequence in the simulated cardiac cycle. For the hard plaque phantom the plaque thickness was measured from the diastolic B-mode image and held constant for all shear wave estimations as the thickness did not change throughout the simulated cardiac cycle.
The radial motion induced by the pressure wave from the simulated cardiac cycle consisted of low frequencies. Thus, prior to estimating the shear wave axial displacement, the data was filtered with a zero-phase bandpass filter with cut-off frequencies at 100 and 2000 Hz to The slope of the shear wave propagation corresponds to the shear wave speed, where a slow travelling wave (soft plaque) has a steeper slope compared to a fast travelling wave (hard plaque). The tissue displacement's dynamic range is approximately −1 (blue) to 1 (red) µm. The corresponding k-space maps for the soft plaque, hard plaque, and vessel wall can be seen in panels (d)-(f) respectively. The viscoelastic effects of the soft plaque limit the shear wave frequency spectrum to lower frequencies compared to the hard plaque and phantom wall. Lamb wave model curve fits (blue lines) and the phase velocity of the shear wave (red dots) for the soft (g, j) and hard (h, k) plaques as well as the phantom vessel wall (i, l) for full spectra (third row from top) and upper halfbandwidth (bottom row) respectively. remove motion artefacts. Moreover, since the phase velocity spectrum changed throughout the simulated cardiac cycle, the half-bandwidth was automatically selected for each pushing sequence for the upper half-bandwidth Lamb wave curve fit. After post-processing, the left and right propagating shear waves were averaged together. Finally, the retrieved signal from the pressure sensor was filtered with a median filter of 10 samples in length to remove noise.

Static experiment
Sample shear wave axial displacement colour maps can be seen in figure 7 showing the shear wave propagation through soft (figure 7(a)) and hard (figure 7(b)) phantom plaques as well as the phantom vessel wall (figure 7(c)) of the soft plaque phantom. The slope of the shear wave propagation in the axial displacement map indicates the speed of the shear wave traveling through different media. The k-space plots of the left propagating shear waves can be seen in figures 7(d)-(f) where the soft plaque phantom (figure 7(d)) has a narrower bandwidth than both the hard plaque phantom (figure 7(e)) and vessel wall (figure 7(f)). Figures 7(g)-(i) and figures 7(j)-(l) show the corresponding Lamb wave model curve fits (blue lines) and the phase velocity of the shear wave (red dots) for the full frequency spectrum and upper half-bandwidth respectively. Table 1 shows the thickness of the plaque and vessel wall for all hard and soft phantoms measured from B-mode images. These values were used as input for the plate thickness in the model of a Lamb wave propagating through a plate. The results indicated good agreement between SWE and mechanical testing shear modulus estimation for both the soft and   The mean shear modulus in the soft plaque measured by SWE was 5.8 ± 0.3 kPa and 3.3 ± 0.5 kPa measured by mechanical testing. In the vessel wall, the mean shear modulus measured by SWE and mechanical testing was 25.0 ± 1.2 kPa and 29.7 ± 0.5 kPa respectively. Similarly, the shear modulus of the left and right shear waves and the shear modulus calculated from mechanical testing for the hard plaque and vessel wall are shown in figure 8(b). The mean plaque shear modulus measured by SWE and mechanical testing was 106.2 ± 17.2 kPa and 98.3 ± 3.4 kPa respectively and the mean shear modulus in the vessel wall measured by SWE was 23.3 ± 2.1 kPa and 29.7 ± 0.5 kPa by mechanical testing. The results for the mean shear moduli and standard deviations calculated by SWE over the full shear wave bandwidth, the upper half-bandwidth, and by mechanical testing are shown in table 2. Table 3 shows that the RMSE for the upper half-bandwidth was less compared to the full spectrum RMSE for both the soft (p < 0.001) and hard (p < 0.001) phantom plaques. A comparison of the means and standard deviations of the shear wave phase velocity (c p ), group velocity (c g ), and mean velocity computed from the Aixplorer ultrasound machine SWE colour map (c Aix ) as well as the respective corresponding shear moduli for the plaques and vessel walls of the hard and soft plaque phantoms can be seen in table 4. The shear wave velocities were similar for the three techniques. However, the estimated shear moduli were significantly different due to invalid assumptions when calculating the shear modulus using the group velocity.   Means and standard deviations of shear wave phase velocity (c p ), phase velocity measurement frequency, phase velocity shear moduli (μ p ), group velocity (c g ), group velocity shear moduli (μ g ), Aixplorer mean colour map velocity (c Aix ), and aixplorer shear moduli (μ Aix ) for the soft and hard plaque phantoms (n = 9).

Dynamic experiment
The normalized carotid artery flow profile curve (blue line) from the programmable pump is shown in figure 9(a) with the ECG triggered shear wave data acquisitions (red dots) highlighted along the curve. The corresponding shear wave axial displacement colour maps (marked 1-8) for the soft plaque phantom can be seen in figure 9(c). The hard and soft plaque phantom lumen pressures during the simulated cardiac cycle can be seen in figure 9(b) with 62 mmHg and 82 mmHg peak transient pressures and a 14 mmHg and 15 mmHg static baseline pressures for the soft plaque and hard plaque phantom respectively. Table 5 shows the means and standard deviations of the measured plaque thickness, phase velocity, frequency at which the phase velocity was evaluated, and SWE estimated shear modulus calculated on the upper half-bandwidth of the shear wave spectrum for the soft plaque phantom. Similarly, table 6 shows the means and standard deviations of the phase velocity c p , frequency at which the phase velocity was measured, and SWE estimated shear modulus calculated on the upper half-bandwidth of the shear wave spectrum for the hard plaque phantom.

Discussion
Non-invasive quantitative characterization of carotid artery plaque is essential in determining the optimal treatment for patients with atherosclerosis. Current clinical standards limit physicians to visual assessment of ultrasound B-mode images and carotid angiograms allowing observer error to influence the results (Kan et al 2012). The goal of this study was to use SWE to characterize hard and soft plaque mimicking inclusions in a phantom setup using phase velocity analysis of the propagating shear wave and compare the results to mechanical testing, as well as measure stiffness changes in the plaque phantom throughout a simulated cardiac cycle. To the best of our knowledge, no previous attempts have been made to assess plaque shear modulus via SWE using phase velocity analysis or to validate the results with mechanical testing. The results indicated good agreement between SWE and mechanical testing for estimating the shear moduli of both hard and soft plaques as well as the vessel walls. In both soft (figure 8(a)) and hard (figure 8(b)) plaques, the SWE measurement slightly overestimated the shear moduli, whereas in the vessel wall the SWE measurement slightly underestimated the shear moduli in comparison to the mechanical testing (table 2). The RMSE in table 3 as well as the sample curve fits (figure 7(g) through figure 7(l)) indicate that a better fit between the Lamb wave model and the SWE phase velocity curve in the soft (p < 0.001) and hard (p < 0.001) plaques as well as the vessel walls (p < 0.001) was obtained by limiting the frequency range to the upper half-bandwidth. In the hard plaques, the shear waves propagated faster compared to the soft plaques and arterial walls (table 4). Thus for hard plaques, less ultrasound B-mode frames recorded the shear wave propagation making the shear wave more difficult to track with the 2D autocorrelator algorithm resulting in an increased error in the shear moduli estimates. This can be seen in figure 8(b), where there is greater variation in the shear modulus estimation for hard plaques compared to the phantom walls and soft plaques ( figure 8(a)). The relative error between SWE using phase velocity analysis and mechanical testing was large, particularly for the soft plaque. However, the goal of this work was to use SWE for plaque characterization and despite large relative errors, this study demonstrated the feasibility of SWE to distinguish between hard and soft plaques. This distinction is possible even with a shear modulus margin of error of a few kPa. Table 4 shows that the phase and group velocities for the soft plaques as well as the vessel wall are quite similar, indicating that the phase velocity curves plateau near the shear wave group velocities. This is not true for the hard plaque as the phase velocity never plateaus (figure 7(h)), making the hard plaques' phase velocity shear wave speed in table 4 less meaningful. The group velocity was in good agreement with the mean velocity measured from the Aixplorer SWE colour map. However, it is not known how the shear wave velocities in the Aixplorer are calculated. Despite similar shear wave speeds for phase and group velocities, the resulting estimated shear moduli calculated from phase and group velocities are different for both hard plaques and phantom walls. The assumptions of an incompressible, homogeneous, isotropic, infinitely large medium when using (1) to calculate the shear modulus using group velocity are not valid resulting in significant underestimates indicating that phase velocity analysis is necessary. The shear moduli calculated by phase velocity and group velocity for the soft plaque were similar, indicating that the soft phantom plaques were large enough for the homogeneous, isotropic, infinitely large medium assumptions to hold, proving that the shear modulus calculated by group velocity is valid in the given experimental setup. Nevertheless, these conclusions do not translate to in vivo soft plaques where a similar study comparing different types of shear wave speed calculations in plaques of varying consistency must be performed.
The axial displacement maps show a reflected wave from the phantom walls at approximately three and four milliseconds for the soft plaque (figure 7(a)) and phantom wall (figure 7(c)) repetitively, while two reflected waves are present in the hard plaque displacement map (figure 7(b)) at approximately 1.5 and 3.5 ms. The shear wave travels faster in the hard plaque allowing for more time to measure reflections during the 6.5 ms acquisition time window. These reflections are not present in the k-space plots as they are cropped out prior to the 2D Fourier transformation to reduce noise in the spectrum.
The post-processing algorithm was able to estimate shear moduli in both soft and hard plaques throughout a simulated cardiac cycle despite arterial wall and plaque motion (approximately 2 mm radial motion) due to the pulse-pressure wave induced by the pump. Although there was some variation in the shear modulus throughout the simulated cardiac cycle, it was less than expected. Couade et al (2010) measured the shear modulus in vivo in a human common carotid arterial wall throughout the cardiac cycle and found a dynamic range of approximately 50 kPa. The small variation in the dynamic experiment can be explained by the fact that the pressure provided by the pump was considerably lower than the pressure in carotid arteries in vivo, as well as differences in PVA's elastic properties to that of in vivo vascular tissue. Despite this, the dynamic test indicated the possibility to conduct future in vivo measurements in a dynamic arterial vessel.
Compared to ARFI (Dahl et al 2009) and quasi-static elastography where strain is estimated by applying a static deformation and measuring tissue displacement (Varghese 2009), SWE is not a relative measurement and not dependent on the applied force. Moreover, SWE offers the advantage of reproducibility (Ramnarine et al 2014) and can potentially be used to track the progression of atherosclerotic plaque development.
Some limitations of the study deserve to be pointed out. The accuracy and appropriateness of modelling the shear wave propagation through the plaque as a Lamb wave propagating through a plate can be questioned, but the results indicate that the model is in agreement with mechanical testing values. Given the complexity of an analytic expression for a wave propagating in a strict hollow cylindrical structure submerged in a non-viscous liquid, the plate model is the current best choice when curve fitting the phase velocity curve to find the shear modulus.
A difficulty when comparing SWE and mechanical testing shear moduli estimations is that the same phantoms are not used for both measurements. The mechanical test required symmetry for the geometry and elastic properties in order to estimate the shear modulus. Since these conditions were not met by the plaque phantoms, it was not possible to use them in the mechanical setup. A biaxial mechanical test was considered, but attaching the PVA material to a biaxial setup proved to be difficult, thus a second set of geometrically symmetrical phantoms were constructed with the same number of freeze-thaw cycles from the same batch of PVA. However, even though different phantoms were tested in the SWE and mechanical test setup, special care was taken in the design of the mechanical test fixture in order to create an identical mechanical stress state between SWE and the mechanical setup. Additionally, it should be noted that although the mechanical test was used as a reference value, it was also subject to error. Measuring the shear modulus for the plaque control phantom with one freeze-thaw cycle in the mechanical setup was difficult due to the soft consistency of the phantom. The phantom sagged once attached to the mechanical fixture and filled with water; thus, the geometrical symmetry assumption was no longer valid. Therefore, the soft plaque shear modulus found by mechanical testing in table 2 is believed to be underestimated as it was not perfectly straight in the mechanical testing fixture, thus imposing less force on the load cell as the phantom was displaced. Adding a frequency in the range of 100-1200 Hz on top of the axial stretch and intraluminal pressure was considered, however the consistency of the results between mechanical testing and SWE indicates that excluding the dynamic SWE push has a negligible effect. The linear behavior of the PVA in the tested range indicates that the effects induced by an SWE frequency would be minor compared to the effects introduced by the intraluminal pressure and the axial stretch, which our mechanical test compensated for.
In this study, only the fundamental (A 0 ) antisymmetric Lamb wave in a plate was studied and a significant amount of data was cropped out of the axial displacement maps in order to reduce noise to achieve a better curve fit between the experimental data and model. It is possible that wave reflections could provide information regarding plaque composition. Moreover, it is feasible to study higher order modes for both symmetric and antisymmetric waves which might reveal clinically relevant information (Bernal et al 2011).
The static test setup was not pressurized due to the sensitivity of the plaque phantoms. Pressurized experiments had been attempted in a pilot study, but only low static pressures (approximately 20 mmHg) were possible before the PVA phantom would mechanically fail. Similarly, the baseline pressure in the dynamic test was lower than observed in vivo diastolic blood pressure values (Hao et al 2013, Zhu et al 2013. The baseline pressure was limited by the PVA phantoms' sensitivity as well as the pump's capabilities.
In this study, the carotid artery phantoms were orientated such that the plaques were only located in the anterior phantom wall. A study with plaques in the anterior and posterior vessel walls must be performed to evaluate the effects of depth when characterizing plaques using phase velocity analysis.
A shortcoming of the plaque phantoms is that they were uniformly hard or soft, while real atherosclerotic plaques are often of inhomogeneous composition consisting of calcified, fibrous, and necrotic regions (Libby 2001). However, the vessel walls were of comparable stiffness to porcine arteries. Bernal et al (2011) measured the shear modulus using a similar SWE technique in excised porcine carotid arteries and found shear moduli ranging from 24-45 kPa for 10-100 mmHg static pressure. Characterization of the biomechanical properties of human plaques have been attempted using tensile and compressional testing with shear modulus values ranging from 7-100 kPa (Heiland et al 2013). This compares favourably with the stiffness of our plaque mimicking phantoms.
The tests were performed in a controlled static setup and the authors acknowledge that it would be more difficult to measure shear moduli in arterial plaques in a clinical setting. Reproducibility requires the transducer and patient to be motionless during image acquisition. In practice this is difficult for elderly patients, especially those with additional ailments and illnesses. Stabilization of the ultrasound probe in a fixture would be beneficial in clinical applications and could potentially improve results and reproducibility. Furthermore, it is essential to focus the ultrasonic push in the arterial wall in proximity to the plaque to ensure shear wave propagation through the plaque. An interactive interface would be required to allow the clinician to click on the B-mode image where the ultrasonic push should be focused. Despite these difficulties, SWE has the potential to be of clinical value when characterizing carotid plaques. To assess the likelihood of rupture, plaque vulnerability must be correlated to the plaque shear modulus in future studies. It is still unknown how the model of a lamb wave propagating through a thin plate will perform on in vivo plaque data, although Couade et al (2010) used the same plate model to measure arterial stiffness in vivo. An in vitro reproducibility study using SWE to evaluate plaque stiffness with a commercial ultrasound machine has been performed, but phase velocity analysis was not considered and the measured stiffness values can be questioned (Ramnarine et al 2014). The present approach using phase velocity analysis is a step towards a more accurate plaque shear modulus estimation, as shown in tables 2 and 4.
Patient safety must be taken into consideration when exciting in vivo plaques with ultrasound radiation force to avoid plaque rupture, which could cause a transient embolic stroke. A simulation of exciting various types of plaques with acoustic radiation force showed that the stress caused by the acoustic radiation is approximately three orders of magnitude less than the force imposed by blood pressure (Doherty et al 2013). ARFI imaging has been used in vivo in a pilot study to characterize carotid plaques (Dahl et al 2009, Allen et al 2011. Furthermore, an in vivo reproducibility study measuring carotid plaque stiffness with SWE using a commercial ultrasound machine has been performed (Ramnarine et al 2014).
Future studies should consider using a model of a wave propagating through a hollow cylinder to further improve the curve fitting procedure. Obtaining B-mode long-axis images of straight in vivo carotid arteries can be difficult and an algorithm that tracks shear wave propagation along a curvature should also be developed. Ex vivo tests on excised carotid arteries with plaques and validation with biaxial mechanical testing should be performed as well as in vivo studies on patients undergoing endarterectomy procedures followed by histological validation.

Conclusion
This study shows that it is possible to measure the shear modulus using SWE in plaque phantoms with phase velocity analysis. The results were validated with mechanical testing and good agreement was found for both the soft and hard plaque phantom as well as the vessel wall. Moreover, it was shown that phase velocity analysis was necessary in the hard plaque phantom and vessel wall, unlike the soft phantom plaque where the shear modulus calculated by group velocity gave similar results as the shear modulus calculated by phase velocity analysis. Additionally, the study shows that it is possible to measure the phantom plaque shear modulus in a dynamic setup throughout a simulated cardiac cycle. The results indicate the possibility to apply the technique successfully in vivo with confidence in the obtained results.