Influence of wall thickness and diameter on arterial shear wave elastography: a phantom and finite element study

Quantitative, non-invasive and local measurements of arterial mechanical properties could be highly beneficial for early diagnosis of cardiovascular disease and follow up of treatment. Arterial shear wave elastography (SWE) and wave velocity dispersion analysis have previously been applied to measure arterial stiffness. Arterial wall thickness (h) and inner diameter (D) vary with age and pathology and may influence the shear wave propagation. Nevertheless, the effect of arterial geometry in SWE has not yet been systematically investigated. In this study the influence of geometry on the estimated mechanical properties of plates (h  =  0.5–3 mm) and hollow cylinders (h  =  1, 2 and 3 mm, D  =  6 mm) was assessed by experiments in phantoms and by finite element method simulations. In addition, simulations in hollow cylinders with wall thickness difficult to achieve in phantoms were performed (h  =  0.5–1.3 mm, D  =  5–8 mm). The phase velocity curves obtained from experiments and simulations were compared in the frequency range 200–1000 Hz and showed good agreement (R2  =  0.80  ±  0.07 for plates and R2  =  0.82  ±  0.04 for hollow cylinders). Wall thickness had a larger effect than diameter on the dispersion curves, which did not have major effects above 400 Hz. An underestimation of 0.1–0.2 mm in wall thickness introduces an error 4–9 kPa in hollow cylinders with shear modulus of 21–26 kPa. Therefore, wall thickness should correctly be measured in arterial SWE applications for accurate mechanical properties estimation.

obtained from experiments and simulations were compared in the frequency range 200-1000 Hz and showed good agreement (R 2 = 0.80 ± 0.07 for plates and R 2 = 0.82 ± 0.04 for hollow cylinders). Wall thickness had a larger effect than diameter on the dispersion curves, which did not have major effects above 400 Hz. An underestimation of 0.1-0.2 mm in wall thickness introduces an error 4-9 kPa in hollow cylinders with shear modulus of 21-26 kPa. Therefore, wall thickness should correctly be measured in arterial SWE applications for accurate mechanical properties estimation.
Keywords: arterial geometry, arterial stiffness, finite element method, phantom, phase velocity, wall thickness, shear wave elastography (Some figures may appear in colour only in the online journal)

Introduction
Changes in arterial mechanical properties strongly affect cardiovascular function and blood pressure levels (Hamilton et al 2007, Chirinos et al 2012, Palatini et al 2011. As arteries become stiffer, the arterial compliance throughout the cardiac cycle decreases, increasing the work on the heart to pump blood throughout the vascular tree (O'Rourke 2007, Maksuti et al 2016a. Notably, increase in arterial stiffness has proven to be an independent predictor for many cardiovascular diseases (Laurent 2006, Hamilton et al 2007, Palatini et al 2011, Scuteri et al 2014, which are the leading cause of death in the world (WHO 2011). Quantitative, noninvasive assessment of arterial mechanical properties could therefore be highly beneficial in the clinics for both early diagnosis of arteriosclerosis and follow up of treatment.
Current commercially-available methods aim at measuring global arterial stiffness by detecting the pulse wave velocity (PWV) between two arterial sites. This technique is at the moment the most simple and reproducible method (Laurent et al 2016), but suffers from several limitations. These are mostly due to the fact that PWV measurements are performed over a very long segment of the arterial tree, which can introduce large bias errors, up to 30%, because of inexact knowledge of the true length of the arterial segment (Nichols et al 2011, Davies et al 2012. Several research methods based on arterial motion tracking in ultrasound imaging have the potential to overcome these problems and measure arterial stiffness accurately and locally (Teixeira et al 2016).
Shear wave elastography (SWE) is an ultrasound-based technique for quantitative and local estimation of tissue stiffness introduced in the late 1990s (Sarvazyan et al 1998) and since then has been applied to many clinical areas, such as detection of breast cancer (Chang et al 2011), characterization of liver fibrosis (Zheng et al 2015) and stratification of follicular-patterned thyroid nodules (Samir et al 2015). The SWE technique is based on generating shear waves by means of highly-focused ultrasound beams and then measuring their propagation speed in the tissue. The shear wave speed is physically related to the shear modulus, which can therefore be derived and used as a measure of the tissue stiffness. Large isotropic organs, whose smallest dimension is much larger than the shear wavelength (order of magnitude of millimetres), can be considered of infinite extent and the so-called group velocity of the shear wave can be used to estimate the shear modulus. The infinite-extent assumption can be considered valid in large organs, however, in smaller structures, such as arteries, this assumption is no longer valid since geometrical wave velocity dispersion occurs. In such cases, a frequency-dependent velocity, the phase velocity, must be analysed to retrieve an accurate shear modulus of the tissue (Maksuti et al 2016b) (see section 2 for details). Dispersion is also seen in large organs such as the liver due to viscosity (Chen et al 2013. Similarly to the viscosity-dependent dispersion, also geometrical dispersion can be studied by analysing phase velocity-frequency relations, i.e. dispersion curves. Phase velocity analysis based on guided wave propagation in a plate has been previously applied to quantify arterial stiffness in vitro (Couade et al 2010, Bernal et al 2011, Widman et al 2015, Maksuti et al 2016b, ex vivo (Bernal et al 2011, Widman et al 2016 and in vivo (Couade et al 2010) showing the capability of the technique to measure local shear modulus multiple times throughout the heart cycle. Phase velocity analysis using a plate model to derive the shear modulus has also been validated against mechanical testing in pressurized arterial phantoms (Maksuti et al 2016b).
In order to reach clinical applicability of arterial SWE, additional aspects such as the effect of anisotropy, viscoelasticity and arterial geometry need to be addressed. A feasibility study in an ex vivo horse aorta sample showed the capability of SWE to detect anisotropy (Shcherbakova et al 2014), but phase velocity analysis was not performed. Geometry (Bernal et al 2011) and viscoelasticity (Chen et al 2004) can both generate dispersion of the wave velocity with frequency. However, simulation studies using plate models showed that viscosity did not substantially affect the dispersion curve (Nguyen et al 2011, Caenen et al 2015. Arterial geometry varies with age (O'Rourke 2007, Hickson et al 2010 and pathology (Krejza et al 2006, Lim et al 2008, Mitchell et al 2008 and may influence wave velocity dispersion and, thus, the estimated shear modulus. Nevertheless, the influence of arterial geometry, in terms of diameter and arterial wall thickness, has not yet been thoroughly investigated. Therefore, the primary aim of this study was to assess how changes in wall thickness and inner diameter influence dispersion curves of guided waves in arterial phantoms and by finite element analysis. The secondary aim was to evaluate the magnitude of the error that arises by using a plate model to derive the shear modulus in arterial phantoms if thickness is not correctly measured.

Shear waves propagation in free space and in waveguides
Shear waves propagate in a homogenous, isotropic, linear and elastic bulk medium at a constant and unique velocity, which is a function of shear modulus μ and density ρ (Achenbach 1973): Equation (1) is currently used in commercial SWE systems to estimate the shear modulus in bulk organs such as liver and breast. However, the bulk medium assumption is not correct for shear wave propagation within arterial walls. In fact, it has been shown that deriving mechanical properties based on group velocity, corresponding to c T in (1), significantly underestimates the shear modulus (Maksuti et al 2016b). Arteries can be considered as tubes containing blood and being surrounded by soft tissue. When waves propagate in the arterial wall there are multiple reflections at the interfaces wallblood and wall-soft tissue. Mechanical energy tends to remain confined within the arterial vessel, which therefore behaves as a waveguide. Shear waves propagate in the vessel walls only at low frequencies (<2 kHz) possibly because of viscous and geometric-based attenuation, and due to shear to longitudinal wave conversion. These frequencies correspond to wavelengths of few millimetres that are of the same order of magnitude as the arterial wall's thickness. Lamb (1917) and Gazis (1959) extensively studied the theory behind guided wave propagation in plates and hollow cylinders, respectively. Both studies highlighted that elastic wave propagation in a waveguide is affected by dispersion. This means that more than one wave mode can propagate at the same time and that each wave mode velocity is a function of frequency, i.e. waves propagating at a given frequency travel at a specific phase velocity, which results in a specific wavelength. The function between phase velocity and frequency is usually referred to as dispersion relation. The effect of dispersion in a waveguide can be seen in figure 1.
Because of dispersion, equation (1) is not suitable for SWE in confined geometries such as arteries, and another approach must be used. Moreover, even though the dispersive behaviour in plates and hollow cylinders appears to be different, the fundamental antisymmetric A(0) mode for the plate is comparable to the flexural F(1, 1) mode for the hollow cylinder, at higher frequencies. This suggests that an analytic function derived from a plate could also be used for hollow cylinders.

Numerical solutions of guided wave propagation
Currently, there are several software routines that are able to evaluate plates and cylinders dispersion relations by numerically approximating the Lamb's and Gazis' equations. Among these, DISPERSE (Imperial College, London) initially described in Pavlakovic et al (1997) is a commercial software based on the global matrix method developed by Lowe (1995). DISPERSE is able to retrieve dispersion curves for plates and cylinders with an arbitrary number of layers and under the assumption of a waveguide in vacuum or surrounded by liquids or solids. Figure 1 shows the dispersion curves (fundamental modes) for an elastic plate and hollow cylinder in vacuum and in water with equal thickness and elastic properties. Dispersion curves for an elastic plate and a hollow cylinder in vacuum and in water, obtained using the software DISPERSE. Material properties for both plate and tubes were Young's modulus E = 66 kPa, density ρ = 1060 kg m −3 , Poisson's ratio ν = 0.4999, resulting in a longitudinal wave speed of c L = 1500 m s −1 . Wall thickness was 1 mm for both geometries. The hollow cylinder's inner diameter was 6 mm. The mode definition is the one adopted by Silk and Bainton (1979) and only fundamental modes are shown. The Lamb's A(0) mode in the plate is comparable to the Gazis' F(1, 1) mode in the hollow cylinder. It can be observed that A(0) and the F(1, 1) modes tend to the same phase velocities at high frequencies.

Wave propagation in plates immersed in an incompressible fluid
In order to take the geometrical wave velocity dispersion into account while reducing the complexity, Bernal et al (2011) modified the Lamb wave propagation model in a plate in vacuum by adding the surrounding fluid and approximated the complex propagation in the artery with this modified zero-order antisymmetric Lamb wave mode. This simplification is particularly valid at high frequencies, where the dispersion curves of the plate and the hollow cylinder converge (figure 1). Being able to use a plate model for arterial SWE applications could allow for fast inversion to obtain the shear modulus in nearly real time during an ultrasound examination. The antisymmetric Lamb wave mode for a solid elastic plate surrounded by a nonviscous incompressible fluid of similar density is described by the following equation (Bernal et al 2011: is the shear wavenumber, k k L 2 s 2 β = − , c p being the frequency-dependent Lamb wave velocity, ω is the angular frequency, h is the plate thickness. Because k s and β are functions of the shear modulus μ, it is possible to estimate μ by fitting experimental SWE data to (2).

Methods
In order to study the effect of geometrical dispersion in arterial SWE, multiple steps were performed. Specifically, geometrical dispersion was first studied in a simple confined geometry, i.e. a plate, and then in a more realistic geometry, i.e. a hollow cylinder. Both geometries were surrounded by water in both experiments and simulations. Water was also present inside the cavity of the hollow cylinder. This section describes the methods and materials used to achieve this goal and is divided in four main parts. The first part describes the experimental measurements on plate and hollow cylinder phantoms, the second part describes FEM simulations of equivalent phantoms as the ones used in the experiments, the third part describes the quantitative comparison between experiments and simulations and the fourth part calculated the thickness-dependent error that can be generated in arterial SWE applications. Plate and hollow cylinder phantom geometries were chosen in order to reproduce dimensions of the human carotid arteries with intima-media thickness in the range 0.5-1.5 mm (Lim et al 2008) and inner diameter in the range 5-8 mm (Krejza et al 2006). Thicker phantoms of 2-3 mm were also considered in order to evaluate the robustness of the technique and the effect of geometry on wave dispersion. For the thinnest wall thickness only FEM simulations were performed due to the complex manufacturing of thin arterial phantoms.

Experimental setup
3.1.1. Phantoms construction. Two sets of phantoms were constructed: plates and hollow cylinders. Multiple samples were manufactured with varying thicknesses of the plate (0.5, 0.8, 1, 1.5, 2, 2.5, 3 mm) and of the hollow cylinder wall (1, 2, 3 mm). Complete information about the geometry of the phantoms is reported in table 1. All phantoms were made out of a solution of 10% in mass poly(vinyl alcohol) (PVA) cryogel (molecular weight: 56.140 g mol −1 , density: 1.269 g cm −3 , Sigma-Aldrich, St. Louis, MO, USA), 87% deionized water and 3% graphite powder (molecular weight 12.01 g mol −1 , density 5 g cm −3 , particle size <50 μm, Merck KGaA, Darmstadt, Germany). The solution was continuously stirred and heated up to a temperature of approximately 55 °C on an electrical stove. Thereafter, the solution was cooled down to approximately 45 °C and poured into customized moulds.
In order to obtain homogeneous mechanical properties, multiple plate phantoms were manufactured simultaneously using the customized mould in figure 2. The mould was constructed as an open box divided in sections by panels, which could slide along slots engraved in the box. The panels were needed to separate the plate phantoms from each other as well as to facilitate the phantom removal from the mould. The distance between the panels was the desired plate thickness. The panels were covered with a thin layer of car wax before usage in order to prevent the plate phantoms from attaching and breaking during removal from the mould. The PVA solution was poured into the box and successively the panels were inserted. Next, the mould with the PVA solution was put in a freezer at approximately −21 °C for 12 h and then let thaw at room temperature for an additional 12 h. This freeze-thaw (F/T) cycle was repeated three times to increase the phantom stiffness.
For the hollow cylinder phantoms, three cylindrical acrylic moulds with a length of 10 mm and different outer diameters of 7, 8 and 9 mm were constructed. The inner diameter of the  5, 5.5, 6, 6.5, 7, 7.5, 8] 100 phantoms remained constant as the wall thickness increased. This was achieved by inserting a metallic rod with a fixed diameter equal to 6 mm in the middle of the cylindrical mould (figure 2). The rod was used to create a cavity in the middle of the cylinder, representing the arterial lumen, and was successively removed. Using this set of moulds, three types of arterial phantoms were constructed, all of which had a constant inner diameter of 6 mm and varying wall thickness of 1, 2 and 3 mm. The moulds were covered with a thin layer of petroleum jelly to enable an easy removal of the phantoms and then the PVA solution was poured into the moulds and underwent three F/T cycles of 24 h each. For the thinnest phantoms (wall thickness of 1 and 2 mm), the mould was placed on the vibrating plate of an homogenizer (UltraTurrax, IKA Werke GmbH, Staufen, Germany) after the PVA solution was poured into the mould in order to aid the solution to sink and allow air bubbles to escape before freezing the material. Due to the viscosity of the solution, air bubbles are very likely to be trapped in the thin cavities, creating holes in the arterial phantoms. An additional sample phantom (solid cylinder, height = 2 cm and diameter = 3 cm) was manufactured in the same way as the plate phantoms and served as a sample for measuring the attenuation coefficient of PVA.

Attenuation coefficient measurement.
The PVA attenuation coefficient was measured using single-crystal transducers with centre frequencies equal to 5, 7.5 and 10 MHz (models V309/320/311-SU, Olympus NDT, Waltham, MA, USA), covering the SL15-4 frequency bandwidth in SWE mode (SuperSonic Imagine, Aix-en-Provence, France). The attenuation coefficient was calculated by recording both the spectrum of the ultrasound signal between the transducer and a reflector at a known distance and the spectrum of the signal when inserting the PVA phantom between the transducer and the reflector. The second spectrum was subtracted to the first one in order to obtain the attenuation in decibels at the central frequency of each transducer, given the traveling distance as twice the distance between the transducer and the reflector (Grishenkov et al 2009). Each measurement was repeated three times, after removing the phantom and placing it again in water. The attenuation coefficient of PVA was needed as an input to the simulations.
3.1.3. Shear wave elastography measurements. The phantoms were placed in a water bath, one by one, and ultrasound data were acquired using the Aixplorer ® system (SuperSonic Imagine, Aix-en-Provence, France) with a SuperLinear ™ SL15-4 transducer and a customized research package allowing for export of the in-phase and quadrature (IQ) data.
A schematic representation of the plate phantom measurements is shown in figure 3. The plates were gently placed on a V-shaped block such that the plate surface would lay parallel to the transducer. The plate phantom was kept still by fixing it with an elastic band placed around the V-shaped block. Attention was paid to not stretch the plates, which was easily achievable since the PVA plates were partially floating once immersed in water.
The hollow cylinder phantoms were attached to a customized enclosure, filled and surrounded by water but not pressurized, as shown in figure 3.
The transducer was placed approximately 10 mm from the plate or hollow cylinder upper wall. The Aixplorer system generates the acoustic radiation force (ARF), or push, at three focal points at different depths (Bercoff et al 2004). The region of interest for SWE data acquisition was chosen such that the second focal point of the supersonic shear imaging push would hit the plate or cylinder wall. After ARF generation, the shear wave propagation was imaged with a pulse repetition frequency of 8000 Hz and IQ data were exported to MATLAB (R2014a, Mathworks, Inc., Natick, MA, USA). The axial incremental displacement generated by the shear wave propagation was tracked using a 2D auto-correlation approach (Loupas et al 1995) according to the incremental displacement definition used in Pinton et al (2006).The obtained axial displacement map wave was cropped in the time direction after 2.5 ms in order to exclude circumferential waves from the field of view. The effect of cropping the data in time can be seen in figure 4. A secondary wave appears in the axial displacement map few milliseconds after the first wave front (figure 4(a)). By cropping the axial displacement map before this wave appears in the field of view, it is possible to better isolate the main wave front and therefore the F(1, 1) mode. If the axial displacement map is not cropped early in time, multiple modes merge together in the phase velocity map (figure 4(c)), which can then result in a local maximum in the dispersion curve as shown in a previous study (Maksuti et al 2016b). When the axial displacement map is cropped in time, this does not occur and the local maximum does not appear in the dispersion curve. After being cropped in time, the axial displacement map was analysed as described in our previously study (Maksuti et al 2016b). In summary, data were converted by means of the Fourier transform to the frequency domain (k-space) in order to obtain a phase velocity intensity map (i.e. phase velocity versus frequency map). The dispersion curve of the A(0) mode for the Lamb model, respectively F(1, 1) mode for the Gazis model (see section 2), was obtained by identifying the maximum intensity of the phase velocity map at each frequency. By excluding the secondary wave from the axial displacement map, it was possible to isolate the F(1, 1) mode in the frequency analysis in a more robust way since the secondary wave did not appear and merge with the F(1, 1) mode in the phase velocity map (figures 4(c)-(d)).
Finally, the shear modulus was retrieved by fitting the dispersion curve for frequencies above 200 Hz to equation (2) by maximizing the R 2 coefficient obtained for different values of μ. The experimentally-estimated μ was then used as input to the FEM models.

Finite element method simulations
Simulations were conducted using COMSOL Multiphysics Acoustics Module (version 5.0, COMSOL, Sweden).
3.2.1. Acoustic radiation force simulation. The ARF generation was modelled by adopting a previously proposed method (Caenen et al 2015, Palmeri et al 2005. The body load distribution F N m 3 ( ) − generated by the ARF application was calculated from the temporally averaged acoustic intensity I W m 2 ( ) − (Nyborg 1965, Torr 1984, Palmeri et al 2005 as where c the longitudinal speed of sound and α is the attenuation coefficient of PVA. This latter was measured as described in section 3.1.2 to be equal to 0.9 dB cm −1 • MHz. The body load distribution in three different planes can be seen in figure 5. The acoustic intensity I can be derived from the acoustic pressure p by means of where Z is the acoustic impendence and T the time duration of the ARF application. The body force was assumed parallel to the Poynting vector, which defines the direction of acoustic energy flux density (Nyborg 1965). It was assumed that this vector is directed purely in the The red line represents the maximum intensity at each frequency, i.e. the final dispersion curve. If secondary waves are present in the axial velocity map, different modes appear in the phase velocity map, which might not be correctly distinguished in experimental measurements. By cropping the axial velocity map before the secondary wave appears in the field of view, it is possible to better isolate the F(1, 1) mode. The fusion of multiple modes can appear as a local maximum in the dispersion curve as shown in a previous study (Maksuti et al 2016b). When the axial displacement map is cropped in time, this does not occur and the local maximum does not appear in the dispersion curve.
axial direction (i.e. parallel to z axis of the models) for locations within 10% of the focal distance.
The acoustic intensity field was estimated using the simulation tool FOCUS (version 0.905, Michigan State University) (McGough 2004, Chen andMcGough 2008), which takes into account both the technical specifications of the ultrasound probe and the acoustic properties of water (c L = 1540 m s −1 , α = 0.04 dB cm −1 • MHz, ρ = 1000 kg m −3 ). Supersonic Imagine's SL15-4 probe specifications were used as input (kerf = 0.02 mm, width = 0.2 mm, height = 4 mm, 256 elements, central frequency f 0 = 6 MHz). Focal distance was set to F A = 10 mm in the azimuth plane, and the equivalent curvature radius of the elevation focusing lens to R c = 30 mm. The number of activated elements was 30, resulting in an F-number of approximately 1.6.
FOCUS was used to generate a three-dimensional (3D) acoustic intensity field I x y z , , ( ) under continuous wave emission. The main lobe was fit to a 3D Gaussian distribution. This was achieved by separately fitting the acoustic intensity along each dimension with a onedimensional (1D) Gaussian using a least square approach. Specifically, where I max , b i and c i are the parameters to be estimated. I max was then removed from each function. The 3D Gaussian distribution I x y z , , ( ) was then evaluated as the product of the 1D Gaussian profiles and imported into COMSOL as an analytical function, Given I(x, y, z), equation (3) was used to define the body force distribution, which was imported into COMSOL.
The shear waves source was imposed in the red volume of the COMSOL model in figure 6 as a body load distribution.

Geometry and mesh.
In order to reduce computational costs, only a quarter of the complete geometry reported in table 1 was taken into consideration for both plate and tube models, when implemented in COMSOL. This was achieved by imposing symmetry boundary conditions, as will be explained in section 3.2.5.
For plates, the lateral extent of the model was set to L 2 min λ + * . The parameter min λ is defined as the smallest wavelength to be resolved in the simulation, and occurs when the problem is solved at the maximum frequency of interest (1 kHz in this study). Plate thickness was set to h. The height of the model was set to h 2 1 mm min ( ) λ * + + . Therefore two layers of water were prescribed, one above and one under the plate, both 1 mm min λ + thick. A plate's corner (highlighted in red in figure 6) was designated as the ARF application region.
For tubes, axial length of the model was set to L 2 min λ + * . The tube's inner diameter was set to D and its wall thickness to h. The surrounding water domain was a rectangular box. The height and width of the model were set to D h 2 1 mm min ( ) λ + * + + . The ARF application domain was defined as an ellipsoid (highlighted in red in figure 6) located inside the tube's upper wall.
Unstructured (i.e. tetrahedral) elements were used to mesh the entire geometry but the outer frame (light blue in figure 6), where semi-structured (i.e. swept) elements were used instead. Mesh size was set according to the 10-12 degrees of freedom per wavelength criterion suggested in the literature (Harker 1984, Alleyne 1991), corresponding to 5-6 second order Lagrange elements per wavelength. The minimum wavelength λ min parameter was used to define mesh dimension. Attention was paid to assure that there were at least 5 mesh elements throughout the perfectly matched layer (PML) thickness because this condition is necessary for the PML to work correctly (see section 3.2.6 for more details). The ARF application domain was meshed using unstructured elements of size 0.25 mm, consistent with (Palmeri et al 2005), in order to more accurately resolve the ARF distribution.

Model's material properties.
A linear elastic material was chosen to describe the plates' and tubes' elastic behaviour. The material was characterized by the Young's Modulus E, Poisson's ratio ν and density ρ. For each phantom to be simulated and assuming incompressibility of the medium, the shear modulus μ obtained by SWE was used to calculate the corresponding E according to For both tube and plate phantoms, density was set to ρ = 1060 kg m −3 . Poisson's ratio was defined in such a way that the longitudinal wave speed c E 1 1 1 2 L ( )/ ( )( ) ν ρ ν ν = − + − would be equal to 1500 m s −1 . This condition resulted in ν values greater than 0.4999 (i.e. a nearly incompressible material) and therefore a hybrid formulation for equation (8) was needed in order to cope with the amplification of numerical errors in the evaluation of volumetric strain caused by the high bulk modulus (Bower 2009).
In order to study the effect of small wall thickness difficult to achieve in phantoms, additional simulation in hollow cylinders with different thickness and diameter were performed according to table 1. For these additional simulations, a constant μ = 22 kPa (E = 66 kPa) was chosen for all geometries.

Waves propagation equations.
The acoustic problem was evaluated in the frequency domain by considering a harmonic source with the same amplitude at all frequencies. Different sets of equations were implemented for wave propagation simulations in the solid and fluid domain. In order to evaluate shear wave propagation in solids, the following displacement equation of motion was used (Achenbach 1973 The linearized acoustic equation was implemented to evaluate wave propagation in fluids (Temkin 1981) where Q m and q d are contributions related to monopolar and dipolar sources, if present, and ρ f is the density of the fluid.
3.2.5. Boundary conditions. In addition to equation (8) and (9), a third set of equations was implemented to couple the two independent variables at the boundaries between fluid and solid domains n p , σ = (11) where n is the vector normal to the boundary between fluid and solid domains. Figure 6 shows an overview of the additional boundaries conditions. The regions highlighted in light and dark blue indicate that symmetry conditions were applied. For the regions belonging to the solid domain, this condition was equivalent to assuming the displacement component normal to the symmetry plane equal to zero, i.e. For the regions belonging to the fluid domain, it was equivalent to prescribing infinite acoustic impendence, i.e.
At the remaining boundaries, zero stress condition was imposed in the solid domain and infinite acoustic impendence condition was imposed in the fluid domain.
3.2.6. Perfectly matched layer. Perfectly matched layers (PMLs) were initially developed for simulations of electromagnetism phenomena (Berenger 1994) in order to prevent reflections at the computational domain borders and can also be applied to mechanical problems. A PML assures strong attenuation of waves propagating through the region where it is applied, by locally scaling the coordinate system according to where f ( ) ξ is a complex function of real variable. The vector n ξ defines the direction of coordinate stretching and is always normal to the outer boundaries. The normalized variable ξ is defined as where w ∆ is the PML region thickness along the direction n ξ and x 0 is a scalar value that represents the PML coordinate along the direction n ξ .
By activating PML equations in the domains highlighted in light blue in figure 6 it was possible to cancel back propagation of reflected waves from the borders of the computational region. (8) and (9) define the waves propagation problem in the frequency domain. A solution can be interpreted as the steady state response of the system to a sinusoidal excitation of given angular frequency ω. Each case was solved in the range 10-1000 Hz with a frequency step of 10 Hz, yielding a total of 100 solutions per case. Each solution was retrieved using the parallel sparse direct solver (PARDISO), which directly inverts the stiffness matrix using LU decomposition. Therefore, there was no need to set a tentative initial condition or to perform preconditioning like in iterative methods such as gradient descent-based solvers.

Solving the waves equations. As previously stated, equations
3.2.8. Simulations post processing. The post-processing used to analyse the FEM simulations results was conceptually the same as the one used for the processing of SWE experimental data (section 3.1.3) and first adopted in Bernal et al (2011) and Maksuti et al (2016b). Some changes were made in order to cope with the differences between FEM and experimental datasets and are summarized in this section.
The displacement component parallel to the z-axis was evaluated at equally spaced points along a line of sight lying in the middle of the wall.
The corresponding frequency-space displacement matrix ( ) u f x , was then imported in MATLAB. A 1D discrete Fourier transform was applied to u f x , ( ) along the columns (i.e. along the spatial dimension). U f k , ( ) is the obtained Fourier transform. The distribution of ( ) U f k , was then converted to a phase velocity-frequency map U c f , p ( ) by relating the given wavenumber k to the corresponding phase velocity c p for each pixel in the matrix, by means of equation ( ) for each column (i.e. for each frequency the problem was solved for). Not to accidentally extract other modes than the fundamental one, only the points below the maximum velocity for the fundamental mode were considered.
Dispersion curves were, first, compared with the ones retrieved by means of DISPERSE software and, then, with experimental measurements in phantoms as described in the following section.

Experiments and simulations comparison
The dispersion curves obtained by SWE for all plate and hollow cylinder phantoms were compared with the corresponding FEM simulations. The quality of fit between the experimental and simulated dispersion curve was measured by the R 2 coefficient and the absolute error was measured by the root mean square error (RMSE).

Thickness-dependent error
In the clinical setting, arterial thickness might not always be measured accurately. Therefore, the error induced by inserting an incorrect thickness in the plate model (equation (2)) was quantified. To evaluate this thickness-dependent error, the dispersion curve obtained by FEM simulations for a specific hollow cylinder phantom was fit to the plate model equation and the estimated shear modulus was derived, similar to the phantom experiments (see section 3.1.3). Initially, the input thickness was correct. Thereafter, the input thickness was altered from the true values with ±0.1 mm steps, corresponding to the smallest detectable distance in the Aixplorer system. For each input thickness, a new shear modulus estimation was obtained and compared with the shear modulus estimate obtained with the correct input thickness. The thickness-dependent error e h was finally calculated as where e h is the thickness-dependent error in kilopascals, μ(h ± Δμ) is the shear modulus retrieved with the thickness value altered by Δh and μ(h) is the shear modulus retrieved with the correct thickness value.

Experiments and simulations comparison
4.1.1. Plates. An example of the phase velocity map and the corresponding A(0) mode obtained with FEM and by SWE on a plate phantom can be seen in figure 7. In the FEM phase velocity map, the energy of the signal is more confined around the A(0) mode compared with experimental data, where the intensity is higher in a broader region. There was good agreement between A(0) from FEM simulation and experiments for all thicknesses, as shown in table 2.

Hollow cylinders.
The comparison between the Gazis mode with n = 1 obtained with DISPERSE and the phase velocity map obtained by FEM simulations for a hollow cylinder of thickness 1 mm and diameter 6 mm is shown in figure 8(a). The low-order modes from DIPSERSE overlap well with the high intensity regions of the phase velocity map obtained by FEM simulations. Figure 8(b) shows a comparison between the dispersion curve for the plate and the cylinder with the same mechanical properties and thickness. The two dispersion curves merge after about 400 Hz for this specific geometry and mechanical properties. Figure 9 shows the phase velocity maps obtained by FEM simulations and exper imental SWE measurements for different hollow cylinder phantom wall thickness. The maps obtained by FEM simulations can differentiate the propagating modes. On the contrary, only F(1, 1) can clearly been distinguished in the experimental phase velocity maps. A weak higher order mode is also visible in the range of frequencies 200-400 Hz, which moves to higher    frequencies with increasing thickness. The intensity of the higher order mode is low in the experimental map also as a result of the time crop of the ultrasound data described in figure 4. The F(1, 1) mode obtained by FEM simulations agreed well with the F(1, 1) mode obtained by experiments, as can be seen in table 3. The comparison was done for frequencies above 200 Hz and 400 Hz. When increasing the frequency of the lower cut off from 200 to 400 Hz, the absolute error in meters per seconds calculated in terms of RMSE was reduced, but the quality of fit quantified by the R 2 coefficient was also reduced, since R 2 is a normalized error and at higher frequencies the variance of the data decreased more than the residuals did.

Diameter and thickness changes
Thickness changes, while keeping mechanical properties and inner diameter of the hollow cylinder unchanged, influence the dispersion curve ( figure 10(a)), whereas changes in diameter, while keeping thickness constant, do not considerably affect the dispersion curves above 400 Hz.

Thickness-dependent error
When fitting the FEM dispersion curve from a hollow cylinder to the plate model in (7) and inserting the true thickness (1, 2 or 3 mm), the two curves showed a good agreement (figure 11). The true shear modulus, i.e. the value used as input in the FEM simulations, was underestimated by 2-3 kPa.   The additional thickness-dependent error is shown in figure 12. The input thickness-error plots report the additional error in kilopascals that is generated by inserting the wrong wall thickness as input in the plate model used for the fit.

Discussion
Simulation and phantom experiments showed that arterial geometry influences dispersion by affecting the phase velocity curve (figures 9 and 10) and consequently the shear modulus Figure 12. Error in kilopascals resulting from using an incorrect input wall thickness in the plate model used to estimate the shear modulus in arterial shear wave elastography. The error is reported for hollow cylinders of different wall thickness. estimation by SWE ( figure 12). Specifically, wall thickness had a larger effect than inner diameter on phase velocity at all frequencies and diameter did not substantially affect the dispersion curve above 400 Hz ( figure 10(b)). An underestimation of 0.1 mm or 0.2 mm in a 1 mm thick arterial wall can introduce an error of 4 kPa or 9 kPa, respectively. Due to the nonlinearity of the relation between geometrical parameters and dispersion, the wall thickness measurement error in millimetres does not linearly correspond to the stiffness estimation error in kilopascals. For instance, it can be noticed that the magnitude of the stiffness error is smaller if the wall thickness is overestimated rather than underestimated (figure 12).
The plate geometry was first analysed because the plate equation is often used in arterial SWE applications to retrieve the shear modulus (Bernal et al 2011, Widman et al 2015, Maksuti et al 2016b. This was a preliminary validation of the FEM model with physical experiments. A plate model could also be useful in other SWE applications, such as e.g. assessment of skin elasticity (Lee et al 2015). Figure 8(a) shows that the A(0) mode in the plate and the F(1, 1) mode in the hollow cylinder merge at higher frequencies, confirming that a plate model is a valid approximation for the hollow cylinder if considered in a suitable frequency range.
There was good agreement between FEM and experimental results for both plates (table 2) and hollow cylinders (table 3), as also shown in figures 7 and 9. Additionally, there was good agreement between FEM and numerical results of Gazis equations solved in DISPERSE (figure 8), even though the comparison was limited to the lowest propagating modes. This comparison aimed at validating the FEM model in a simple geometry, for which analytical solutions are available and can be simulate by means of DISPERSE. Once the model has been validated, it can be used for more complex cases in order to study e.g. the effect of surrounding media, anisotropy, viscosity or the presence of a plaque. These more complex cases cannot be studied with the use of DISPERSE.
From both simulations and experiments it can be concluded that A(0) and F(1, 1) modes show similar dispersion relations at frequencies above 400 Hz, or a frequency-thickness product of 0.4, for geometry of interest in arterial applications (figure 9). Using a plate model to fit dispersion curves in a hollow cylinder generated an error of 2-3 kPa in models with shear modulus of 21-26 kPa, which is within the overall accuracy of the technique of 10% when compared with mechanical testing (Maksuti et al 2016b). The fact that diameter did not substantially affect the dispersion curves is of particular relevance since diameter is clearly not a parameter that can be modified in the plate model (a plate can be seen as a hollow cylinder of infinite inner diameter). Therefore, it is feasible to fit the experimentally measured F(1, 1) mode of a hollow cylinder to the A(0) mode of a plate in order to estimate the shear modulus. The use of a simpler model for wave propagations in the arterial wall is advantageous compared with solving the wave equations in a hollow cylinder, which would require longer time and higher level of complexity. On the contrary the plate model equation can be easily implemented in a commercial scanner and used to estimate the shear modulus in real time.
Previous studies have analysed wave dispersion either experimentally in tubes with a constant geometry (Bernal et al 2011, Maksuti et al 2016b or numerically and experimentally in plates in vacuum (Caenen et al 2015). To our knowledge, no studies have investigated wave dispersion for SWE applications in plates and hollow cylinders in water both numerically and experimentally, assessing the sensitivity of the technique to geometry. Therefore, these results represent a further important step toward the clinical application of arterial SWE, having characterized the sensitivity of the technique to wall thickness and diameter, and assessed the validity of applying a plate model to a hollow cylindrical geometry.
Tissue-mimicking phantoms can be a useful tool in the development of ultrasound-based techniques such as SWE and have extensively been used for this purpose (Fromageau et al 2003, 2007, Couade et al 2010, Caenen et al 2015, Widman et al 2015, Maksuti et al 2016b, since they allow to have a controlled experiment and at the same time take the most relevant characteristic of soft tissues into account. However, repeatable and uniform arterial phantoms with small wall thickness (1 mm) are difficult to construct. In such cases, FEM simulation can be a valuable substitute given that geometry is not a limitation. Moreover, FEM simulations can help in the understanding of complex wave propagation, which is not always possible to fully characterize experimentally. In return, phantom experiments can be used as a reference for validating the FEM simulation. Therefore we considered the two methods as complementary to each other in this study. An example of the advantage of combining the two methods is that simulated phase velocity maps in figure 9 clearly shows that the local maximum often present in dispersion curves of hollow cylinder (figure 4(c)) and reported in previous studies (Maksuti et al 2016b) derives from secondary waves merging into the F(1, 1) mode and is not a feature of the F(1, 1) mode itself. The fusion of the two modes was not obvious when observing experimental phase velocity maps but it clearly appeared in the FEM phase velocity maps (figure 9). Knowing that the local maximum derives from the fusion of the two modes, one can remove the effect of these additional waves by cropping the axial displacement map earlier in time in the SWE post processing (figures 4(b)-(d)). The secondary wave could be explained by the generation of shear waves in the posterior wall of the phantom by the supersonic push (Bercoff et al 2004) of the Aixplorer system. These waves travel along the circumference and can reach the anterior wall during the ultrafast imaging acquisitions time. Circumferential waves have been previously observed and reported (Couade et al 2010, Urban et al 2013, Hansen et al 2015. Phase velocity at low frequencies (below 200 Hz) could not be reliably acquired in most phantoms measurements and therefore the fit between the plate model and the experimental dispersion curve was performed above 200 Hz.
FEM models were solved in the frequency domain, in the 10-1000 Hz range. This approach offers several advantages over solving the problem in the time domain as previously done (Palmeri et al 2005, Caenen et al 2015. In the frequency domain each frequency component is evaluated as an independent simulation step, whereas in the time domain the previous solution is used as the starting point in the Newton-Raphson algorithm, which needs a convenient time step to achieve convergence (i.e. retrieving a solution within a pre-specified error tolerance). When solved in the frequency domain, the wave propagation problem is linear, enabling the possibility of solving complex simulations using the superposition principle. Solving the models in the frequency domain also means retrieving steady-state solutions for the acoustic problem. The modes are a result of reflections of waves between wall and water. These reflections take time to come back (similarly to what shown in figure 4) and to create a steady-state guided wave. Therefore, truncating the data in the time direction can help differentiate the fundamental mode.

Limitations
Some limitations of the study deserve to be addressed. Figure 1 shows dispersion curves obtained with DISPERSE for a plate and a hollow cylinder with similar mechanical properties as for the PVA phantoms. However, extracting modes at these low frequencies can be difficult for soft materials, which results to a not uniform sampling of the curves as visible in figure 1. The ARF excitation imposed in the FEM model was a simplification of the real experimental set-up. As previously mentioned, the source generating the ARF was assumed to vibrate with constant amplitude at all frequencies, which does not occur in reality since the push is not an ideal impulse, i.e. a Dirac delta function. This implies that the spectral ampl itude in the k-space looks brighter where the modes resonate with highest intensity, i.e. when the modes appear in the k-space at their specific cut-off frequency (vertical white bands in figure 8). Another implication of this approach is that if the frequency data were converted in the time domain, the displacement magnitude could not be directly compared with experimental measurements, acquired as a function of time. This comparison could be done if one applies a window to the source in the frequency domain, but this was not possible since the exact frequency content of the shear wave source was not known. However, the constant amplitude of the source at all frequencies does not influence our results and conclusions since the analysis was based on the phase velocity dispersion curve, which depends on material properties and geometry and not on the frequency content of the source.
The ARF from the Aixplorer system is likely to hit both walls in real experiments and not only the anterior wall as simulated in the FEM models. If both walls are excited simultaneously, multiple waves might be generated, which can interfere with each other. The early time crop performed in the SWE post processing aimed at limiting the impact of these interferences and the agreement between experiment and simulations indicates that it had minor effect on the F(1, 1) mode detection.
FEM simulations were conducted up to 1000 Hz for computational reasons. It is probable that for stiffer hollow cylinders, such as arteries, a larger bandwidth would be needed to obtain reliable shear modulus estimations, since previous studies have shown that a large bandwidth is important to obtain adequate data for evaluating the modulus (Widman et al 2016).
A linear elastic material model was used in all simulations, while PVA has shown viscoelastic and hyperelastic behaviour (Karimi et al 2014). However, the experimental and numerical phantoms used in the current study were not pressurized and therefore there was no need to use a hyperelastic model, since PVA behaves linearly for small deformations such as the ones induced by ARF (Fromageau et al 2003). In addition, previous studies have shown that the effect of viscoelasticity on dispersion curve in gelatine plate phantoms was minor (Nguyen et al 2011, Caenen et al 2015. Therefore, the costs and benefits of increasing complexity of the model's material properties should be carefully evaluated.
The shear wave modulus values obtained in the phantoms and used as an input in the FEM simulations are likely to be lower than the one present in the human carotid artery (Deng et al 1994, Kamenskiy et al 2014. The purpose of this study was to assess whether geometry had a significant impact in arterial SWE and therefore the stiffness of the PVA phantoms was the reference stiffness. The comparison with the experiments demonstrated that the simulation framework was validated and could be used for other studies with varying stiffness. Preliminary simulations showed similar trends as those in figure 10 when simulating thickness and diameter changes on stiffer models. Shear modulus values similar to the in vivo values will be further investigated in future studies.
The presented framework is a simplification of the in vivo common carotid artery case, since arteries are constituted of anisotropy tissue and surrounded by soft tissue and another large vessel such as the jugular vein. The main aim of this study was to assess the influence of geometry alone on dispersion curves and on shear modulus estimation. Future studies will evaluate the combined effect of e.g. anisotropy, viscosity and surrounding tissue.

Clinical significance and future directions
The geometry of the hollow cylinder phantoms was based on real dimensions of the human carotid artery for both wall thickness, measured in clinics as intima-media thickness (Lim et al 2008), and inner diameter (Krejza et al 2006).
Arterial thickness is already clinically measured and could therefore be easily added as an input to the ultrasound system when performing arterial SWE by phase velocity analysis. Because of the high sensitivity, the thickness needs to be measured accurately. A single measurement throughout the heart cycle would be sufficient since the variation amplitude of the mean intimamedia complex during the heart cycle is less than 0.1 mm in healthy volunteers and patients (Zahnd et al 2014). Additionally, many medical ultrasound manufacturers provide tools for automatic intima-media thickness measurements, which would reduce user-dependent errors.
Reference values for changes in arterial shear modulus that are of clinical relevance have not been established yet, therefore the thickness-dependent error could be more or less significant depending on how small changes should be detected in order to differentiate normal from pathological subjects.
Our results show that diameter does not influence dispersion curves substantially, especially at higher frequencies. This is a positive result for the implementation of the SWE by phase velocity analysis in the clinics for two main reasons. First, the plate model is applicable because only thickness is needed in the phase velocity analysis. Second, the common carotid diameter is not only changing among individuals but also throughout the cardiac cycle, with changes up to 0.8 mm in young subjects (Giannattasio et al 2008) and therefore taking diameter changes into account could be challenging.
Future studies should investigate the effect of anisotropy while preserving the tubular geometry of the artery and perform measurements in both the long-axis and short-axis view of the arteries, since multiple modes can be detected (Urban et al 2013) and might provide useful clinical information. Furthermore, a critical review of the major error sources in arterial SWE should be performed in order to identify the factors that need to be addressed and the factors that can be neglected given that they generate a minor error source, before the technique can be reliably applied in the clinics.

Conclusions
Experiments and FEM simulations in plate and hollow cylinder phantoms showed good agreement. Additionally, they showed that a plate model can be a good approximation of the hollow cylinder at higher frequencies. Moreover, simulations in hollow cylinders indicated that wall thickness influences wave dispersion within the phantom's wall and consequently the shear modulus estimation by phase velocity analysis. An underestimation of 0.1-0.2 mm introduces an error of 4-9 kPa in the estimated shear modulus of 21-26 kPa. Therefore, wall thickness should be accurately measured in arterial SWE applications. On the contrary, diameter did not substantially influence the dispersion curve after 400 Hz. These results suggest that is feasible and justifiable to use a plate model in an arterial geometry in order to derive the shear modulus by SWE, given that arterial wall thickness is accurately measured.