Physics of Within-Tissue Wave Propagation Generated by Pulse Propagation in the Carotid Artery

: (1) Background: We aimed to assess the validity of laser Doppler vibrometry (LDV) as an emerging method to measure the local pulse wave velocity (PWV) from skin displacement generated by the pressure pulse inside an underlying artery. (2) Methods: A finite element model representing a simplified common carotid artery embedded within a soft tissue mimicking material was used to reproduce how tissue motions due to a wave propagation along the artery radiates to the skin surface. A parametric study was set up, varying: (i) the pressure conditions inside the artery (shock and traveling pressure impulse), (ii) the arterial depth and (iii) the geometry in a patient-specific artery model. (3) Results: under all conditions, the arterial pulse induced primary and secondary waves at the skin surface; of which the propagation speed deviated from the imposed PWV (deviations between − 5.0% to 47.0% for the primary wave front). (4) Conclusions: the propagation of a short pressure impulse induced complex skin displacement patterns revealing a complicated link between PWV and measured propagation speeds at the skin surface. Wave propagation at the skin level may convey information about arterial PWV, however, advanced signal analysis techniques will be necessary to extract local PWV values.


Introduction
Cardiovascular disease (CVD) remains the leading cause of death worldwide and contributes substantially to morbidity, disability and death in patients [1]. Besides continued efforts to find more effective treatment, identifying subjects at risk of CVD-especially among those categorized as having a low risk according to classical risk factors (age, gender, blood pressure, lipid levels, smoking status . . . )-has received a lot of attention in pre-clinical research. A particular aim of such studies is to prevent the development of CVD at a young age [2,3]. There is compelling evidence that arterial stiffness is a biomarker that carries prognostic power above and beyond these classical risk factors [4].
The clinical "gold standard" method for evaluating arterial stiffness is carotid-femoral pulse wave velocity (PWV c-f ), calculated as the ratio of the distance between the carotid and femoral measuring sites, and the time taken for the arterial pulse to propagate from the carotid to the femoral artery [5][6][7]. Different devices are available to measure these pulse waveforms, but current methods have methodological issues (there is no unequivocal path between the carotid and femoral artery, and the measurement basically excludes the most distensible proximal aorta). Furthermore, they are relatively demanding to use (e.g., arterial applanation tonometry, ultrasound), requiring trained operators and, in the case of magnetic resonance imaging, are associated with high costs [8,9]. These limitations justify the continued search for an accurate, low-cost and non-contact measurement device for PWV assessment that can easily be used in daily clinical practice. In the early stages of atherosclerotic disease, alteration of the arterial elastic properties starts in localized regions to then become more widespread. Thus, a local PWV value represents a more precise indicator of the ongoing arterial stiffening compared with the regional PWV value calculated by the most common commercial devices [8]. Therefore, an easy to use device which allows arterial PWV assessment locally is needed.
Within this context, we are exploring laser Doppler vibrometry (LDV) as a method to locally measure PWV in the common carotid artery (CCA), a superficial large elastic artery known to be prone to remodel and stiffen with age and disease [10][11][12]. LDV allows optical measurement of out-of-plane skin displacements as a function of time, in which the mechanical perturbation caused by the arterial pressure pulse can be identified. This measurement is then performed in at least two different spatial locations along the carotid artery (CA) and PWV can consequently be calculated as the ratio of the distance between the laser beams and the time shift of the recorded displacement signals. Due to the limited distance between the laser beams and the total length of the CCA (a few centimeters), LDV provides a local estimation of PWV compared to the carotid-femoral PWV measurement.
The technical feasibility of LDV as a technique for PWV assessment has been explored in some in vitro and in vivo studies [13][14][15][16][17]. Recent studies have also shown good correspondence between PWV evaluated with LDV and ultrasound-based pulse wave imaging (PWI) in a phantom, and demonstrated the potential of LDV to detect arterial wall stiffening [18]. Even though previous studies have shown promising results for LDV as screening tool for CVD, the understanding of the generated wave propagation mechanics between arterial wall and skin surface is still limited. Indeed, the waves are subjected to reflections, attenuation and interaction with the surrounding structures and organs, potentially leading to wave interference or wave mode conversions, complicating the relationship between the arterial pressure pulse and the corresponding skin displacement.
A computer modelling environment able to reproduce the complex wave phenomena associated with the propagation of the mechanical perturbation caused by the arterial pulse provides a flexible way to perform a complete analysis of the wave propagation. The goal of this study is to provide deeper insights into the naturally occurring wave physics in non-stenotic carotid vessels using numerical simulations. These insights can ultimately improve current LDV analyses for PWV assessment and pinpoint limitations of the LDV technique.
In this work, a model representing the CCA embedded within soft tissues was created using the finite element method (FEM). Its realism in accurately representing wave propagation was assessed through modelling specific modes of bulk and surface waves, where a known analytical solution is at hand. Finally, we used this model to perform a parametric study investigating the correlation between arterial pressure and displacement pattern at the skin surface. The effect of the following parameters was examined: (i) pressure conditions within the artery, (ii) artery depth below the skin surface and (iii) arterial geometry (patient-specific carotid bifurcation).

Theoretical Background: Bulk and Surface Wave Physics
The propagating pressure pulse inside the CA generates a combination of wave modes inside the tissue (longitudinal and transverse waves) and on the skin surface (e.g., Rayleigh waves). All these types of wave obey Navier's governing equation of motion as follows, when approximating tissue as an isotropic linear elastic material: where λ is the Lamé material constant, µ the shear modulus, u the displacement vector, ∇∇· u = grad(div u) and ∇ 2 is Laplace operator.
For an infinite medium, Equation (1) can be further simplified using Helmholtz decomposition, resulting in two equations describing longitudinal or compressive and shear or transverse waves. These equations result in the following formula for longitudinal wave propagation speed c L and transverse wave propagation speed c T : where E is the Young's modulus, ν the Poisson's ratio and ρ the tissue density [19]. For a semi-infinite medium (solid with a single plane surface that extends to infinity in all directions), the free surface will couple the longitudinal and transverse waves to form a surface wave. For the well-known Raleigh wave which travels on the surface, Equation (1) transforms into the following normalized particle displacements describing an elliptical trajectory: where u and w are the displacements in the longitudinal (z) and transverse (y) direction respectively, A the amplitude, k the wavenumber, c R the Rayleigh wave speed and parameters q = The elliptical orbit of the particles decreases in amplitude with increasing depth below the surface and the direction of polarization changes from clockwise to counterclockwise. The surface wave speed can be estimated as [20]:

Materials and Methods
In Section 3.1, we discuss the basic settings of the computational model, including: (1) geometry and boundary conditions, (2) wave propagation and sensitivity mesh analysis, and (3) material model. In Section 3.2, an overview is provided of the different simulations that were carried out to investigate: (1) the effect of different pressure conditions (shock impulse and traveling impulse) in the artery on wave travel and vibration patterns, (2) relevance of the artery depth below the skin surface and (3) a patient specific geometry. In Section 3.3, the data analysis is presented.

Geometry and Boundary Conditions
The anatomical setting of the common carotid artery surrounded by soft tissues was abstracted as a cuboidal domain (width 10 mm, length 50 mm, height 18 mm) with a cylindrical cavity of a diameter representative of the CCA (internal diameter 6 mm). The top surface of the model (skin surface) was located 10 mm above the modelled artery, indicated with parameter d in Figure 1a, but the skin layer was not explicitly accounted for. The model was created using the FEM software Abaqus 14.1 (Abaqus Inc., Providence, RI, USA). The free (skin) surface is in dark grey, while infinite mesh elements placed on the outer boundaries are in light grey. In red are highlighted the node sets at the skin and lumen surface, where displacement data were extracted for the wave analysis. (b) Geometry of the patient-specific carotid artery CA embedded within a cuboidal domain. Sub regions common carotid artery, CCA; bifurcation, BIF; internal carotid artery, ICA and external carotid artery, ECA are indicated. Again, the top external surface (dark grey) is a free surface, while all the other boundaries in light grey are infinite elements. The lines in red show the paths along which the data for the wave analysis were exported. For both panel (a,b), the coordinate system in red indicates the actual position of the model, but for visualization purposes a view triad (located in the right corner of each panel) will be used further in the manuscript to display the view orientation.
The model was meshed with 3D solid hexahedral elements with reduced integration (C3D8R). The skin surface was modelled as a free surface, whereas the other boundaries of the cuboidal domain were covered with infinite elements (CIN3D8) to prevent wave reflection, as shown in Figure 1a.

Wave Propagation and Mesh Sensitivity Analysis
Numerical settings such as time increment size, numerical damping and mesh size were carefully chosen to ensure computationally efficient and reliable wave propagation simulations. The equations of motion for wave propagation analyses were solved in Abaqus using an explicit direct time integration scheme, which is computationally the most efficient for our considered problem. This method is conditionally stable, in which the size of the stable time increment is automatically calculated as the minimum time that a longitudinal wave takes to propagate across any mesh element in the model. This means that the magnitude of the time increment is dependent on the chosen mesh size and material parameters. When approximating soft tissue as a linear elastic material, the longitudinal wave speed expression (Equation (2)) is highly sensitive to the choice of Poisson's ratio ν: if ν approaches 0.5 (incompressibility-as is often assumed for soft tissue), the numerical time increment size gets smaller leading to longer simulation times. Therefore, Poisson's ratio should be cautiously selected based on a trade-off between realistically modeled compressional speeds for soft tissue (typically 1540 m/s) and reasonable computational time. The finally chosen material model together with its parameters is further discussed in Section 3.1. 3. In order to improve the modeling of high-speed dynamic events, the explicit integration scheme in Abaqus automatically introduces bulk viscosity as a numerical damping term. The bulk viscosity model consists of two parameters, linear and quadratic bulk viscosity, which generates a bulk viscosity pressure that is linear and quadratic with the volumetric strain rate. Linear bulk viscosity damped "ringing" at the highest mesh element frequency, while quadratic bulk viscosity introduced a resisting pressure that prevented mesh elements from collapsing. For more information about the implementation and physical meaning of the bulk viscosity, we refer to the Abaqus manual [21], but it is important to note that the bulk viscosity is intended as a numerical effect only and it is thus not part of the material's constitutive response. Both parameters retained their default value in the simulations (0.06 and 1.2 for linear and quadratic bulk viscosity respectively).
A mesh sensitivity analysis was performed to determine the optimal mesh size, taking into account the literature-suggested criterion of 10 to 15 elements per wavelength λ [22]. As longitudinal wavelengths (order of meters) are much larger than shear or Rayleigh wavelengths (both in the order of centimeters) in biological materials, the mesh criterion is more restrictive for shear/Rayleigh waves than for longitudinal waves. Therefore, we considered only shear and Rayleigh wave modeling in the mesh sensitivity analysis. The model with the largest mesh size to which the simulated speeds of both wave types still corresponded well with the theoretical speeds of Equations (3) and (6), gave the optimal mesh size for the parametric study in Section 3.2. Furthermore, the accuracy of the computational model was further verified by comparing the simulated displacement patterns of the Rayleigh wave in the far field (x > λ) to the analytical solutions in Equations (4) and (5).
For the mesh sensitivity analysis, we implemented an elastic material model-similar to that of the theoretical derivation (Equations (3)-(6))-with a density ρ of 1045.5 kg/m 3 , a Young's modulus E of 73.0 kPa and a Poisson's ratio of 0.49999 [23]. Five mesh element sizes were considered: 0.15, 0.20, 0.25, 0.30 and 0.5 mm. Shear waves were modelled by applying a localized, impulsive displacement in the y-direction with a magnitude of 0.01 mm for a duration of 0.1 ms on a single mesh node located in the middle of the spatial domain between the arterial lumen and the skin surface, as shown in Figure 2a. For the Rayleigh wave simulations, the geometry of Figure 1a was considered without the arterial lumen and reduced to a 2D geometry, with a length of 50 mm and a height of 20 mm. Rayleigh waves were induced by applying the theoretical displacement profile described in Equations (4) and (5) (C L = 1078.77 m/s, C T = 4.82 m/s and C R = 4.59 m/s at a frequency of 500 Hz) to the cross sectional plane at z = 0 mm. The simulation duration was set at 5.1 ms and 10 ms for the shear and Rayleigh wave simulations respectively.

Material Model
For the parametric study in Section 3.2, the elastic material definition of Section 3.1.2 was expanded with time-dependent behavior using a 2-term Prony series of the shear modulus G as follows [24]: where G 0 is the initial value of the shear modulus, g i are the ratio of the considered and initial shear modulus, and τ i are the material time constants. The material coefficients were derived from previous uniaxial mechanical testing of tissue-mimicking polyvinyl alcohol (PVA) phantoms in [23] and are summarized in Table 1. Table 1. Summary of the viscoelastic material characteristics of polyvinyl alcohol (PVA) [23].

Material Characteristics
It should be noted that Abaqus allowed only an elastic material definition for the infinite structural elements modeled at the edges of the domain of interest.

Parametric Study
Using the global model definitions as described in Section 3.1, we assessed the effect of different pressure conditions, arterial depth, and the use of a patient-specific geometry on simulated wave physics.

Pressure Conditions in the Artery
We considered two kinds of pressure loads, one causing free wave propagation in the model (called 'shock impulse') and the other mimicking the propagation of the pressure pulse traveling inside the artery (called 'traveling impulse'). Both loads were uniformly distributed over a ring-shaped region 1 mm wide in the arterial lumen and had a magnitude of 40 mmHg, as shown in Figure 2a. Their spatial and temporal profile is the following: • Shock impulse. At z = 1 mm in Figure 1a, the load was gradually applied within 0.1 ms, constantly applied for 0.2 ms and gradually released in 0.1 ms (see left panel of Figure 2b).

•
Traveling impulse. The pressure load traveled along the lumen's long axis from z = 1 mm to z = 50 mm (z-direction in Figure 1a) with a propagation velocity within the range reported for in vivo PWV [25]. We considered three propagation velocities: 4, 7 and 10 m/s (see right panel of Figure 2b) which corresponds to simulating an increase of PWV due to arterial stiffening.
Wave propagation was modeled for a physical time of 20 ms in total, which is more than sufficient to ensure that a pulse propagating at speeds of 4 to 10 m/s has reached the end of the simulation domain (length of 50 mm).

Artery Depth Below the Skin Surface
As carotid depth varies substantially between subjects due to variability in body weight and anatomy, we altered the arterial depth from 10 mm to 5 and 15 mm. Both models were used to simulate the traveling impulse (velocity of 4, 7 and 10 m/s), as described in Section 3.2.1, within a total simulation time of 20 ms.

Patient-Specific Carotid Artery
To study the geometrical effect of a realistic artery on the wave propagation analysis, a patient-specific carotid artery (CA) model was investigated (see Figure 1b). The geometry consisted of the CCA (length 34 mm) and the bifurcation region (BIF; length~16 mm), where the CCA splits into the internal (ICA) and external (ECA) carotid artery. The patient-specific CA was embedded in a cuboidal domain (width 18 mm, length 50 mm, height 26 mm and averaged CCA depth d~15 mm). Due to the complexity of the geometry, we used software developed in-house Pyformex (Pyformex 1.0, Ghent University, Belgium; http://pyformex.org) to model and mesh the geometry using both hexahedral and tetrahedral elements. The boundary conditions were the same as for the CCA model: infinite elements were placed on the external boundaries of the domain to avoid wave reflections except for the free top surface, as shown in Figure 1b. The element size definition of the mesh varied depending on the vicinity of the BIF, but the largest overall element size was chosen based on the outcome of our mesh sensitivity analysis of the CCA model. The accuracy of the selected mesh was verified by modeling shear wave propagation in a similar manner to that described for the CCA model in Section 3.1.2. After this verification step, a pressure impulse traveling along the centerline of the artery with a velocity of 4, 7 and 10 m/s (magnitude and spatial profile similar to that defined in Section 3.2.1) was imposed on the model for a time of 13 ms.

Data Analysis to Estimate Wave Speeds
Wave speeds can be determined using a time-of-flight (TOF) method that analyzes wave position as a function of time. TOF methods assume a priori a known direction of wave propagation and characterize wave position in time by identifying the main excited displacement component of the medium. For the simulated shear and Rayleigh wave propagation in the mesh sensitivity analysis (Section 3.1.2), the main excited displacement component corresponded to the direction of the applied loading (y-direction), while wave propagation occurred perpendicular to this direction i.e., the z-direction. This is visualized in Figure 2a. For the shock and traveling impulse pressure conditions considered in Section 3.2.1, we analyzed the y-component of the displacements propagating in the z-direction at the level of (i) the free surface at x = 0 mm, representing the skin surface directly above the artery (y = 13 mm for CCA model with d = 10 mm)-mimicking actual LDV measurement-and (ii) the lumen surface at x = 0 mm (y = 3 mm for CCA model with d = 10 mm), as visualized in Figure 1a,b.
Concerning TOF-method specifications, we used the Radon transform to estimate the wave speed [26]. Therefore, we first extracted the main component of the particle displacement along the selected wave propagation path at a sample rate of 10 kHz at the FEM grid. In the case of the patient-specific artery, it was necessary to interpolate the extracted displacements from the FEM grid to an equally spaced rectangular grid (0.21 mm) due to the non-uniformity of the mesh grid. To eliminate any high-frequency numerical noise, a 5th order Butterworth low pass filter with cut-off frequency of 1000 Hz was applied to the raw displacement data in Abaqus. Next, the temporal spatial map of displacements was Radon transformed using projection angles from 0 • to 179.75 • in steps of 0.25 • . The angle with the maximum of the Radon transform represented the linear wave trajectory with the largest displacements, and thus its tangent was taken as the wave speed of the main wave front (indicated with V 1 ). Using this method, the wave speed of potential other (secondary) wave fronts can also be identified by searching for the local maxima of the Radon transformed displacements. These are indicated with V 2 (and V 2 * ) . All post-processing procedures were performed in Matlab R2015a Research (The Math Works Inc., Natick, MA, USA).

Mesh Sensitivity Analysis: Shear and Rayleigh Wave Simulations
For the shear wave simulations, shear wave speeds were calculated for both the forward (positive z-direction) and backward (negative z-direction) directions with respect to the loading location as visualized in Figure 2a. Their averaged value is shown in Table 2 for all considered mesh grids together with the computational time. All meshes reproduced the shear wave propagation with a good degree of approximation with respect to the theoretical shear wave velocity value of 4.82 m/s (deviations varied between −0.6% for the finest mesh to −10.0% for the coarsest mesh). The mesh with 0.20 mm element size was chosen as the optimal mesh for our study purpose due to the combination of result accuracy, reasonable computational time and minimized residual computational noise after load release. The accuracy of the selected mesh size of 0.20 mm for wave propagation analysis was further verified by simulating Rayleigh waves. Figure 3a compares the spatial characteristics for simulation and theory, whereas the temporal characteristics of the simulated Rayleigh wave at z = 5 mm are shown in Figure 3a for three depth values (y-direction) together with the theoretical solution. This figure demonstrates that both the spatial and temporal characteristics of the Rayleigh waves were well reproduced in the simulation. Furthermore, the estimated surface wave velocity is 4.63 m/s in the simulation, deviating only 0.9% from the theoretical value (4.59 m/s).

Load Pressure Conditions: Shock Impulse and Traveling Impulse
Wave propagation is illustrated in Figure 4a Figure 4a, where the index r refers to the reflected wave front). Throughout the manuscript, we will focus on the generated positive displacements, as these have the same direction as the applied load.  Figure 4b demonstrates the simulated wave propagation when a propagating 40 mmHg pressure pulse is applied to the lumen (velocity V z of 7 m/s). At the first time instance, the imposed travelling load is visible near the initial position of the lumen (z = 0 mm); however, note that the displacement magnitude (maximum of 0.03 mm) is not yet fully developed at this time point. At t = 2 ms, the induced displacement magnitude of the main wave front in the lumen is about 0.12 mm. At 4 ms, the load has travelled about half the length of the lumen, while the wave front at the skin surface is visible at about one third of the length of the model. The main wave front clearly has an oblique trajectory throughout the thickness of the model, causing a spatial lag along the z-location between the z-location of the main wave front at the lumen and skin level (about 10.48 mm). Maximal displacement magnitude of the main wave front has increased to 0.13 mm at the level of the lumen and 0.10 mm at the level of the skin. At 6 ms, we again observe a secondary wave front at the lumen surface (indicated with f r ) due to reflection at the skin surface.
Out-of-plane displacement (y-direction) as a function of time and propagation path (z-direction) is shown in Figure 5 for the top of the lumen and skin surface under shock and traveling impulse conditions. The propagation of a single shock impulse is shown in Figure 5a, demonstrating one main wave front and one secondary wave front propagating in the lumen displacement plot with speeds V 1 and V 2 of 4.66 m/s and 11.6 m/s, respectively. At the skin level, the main wave front propagated with a speed V 1 of 9.74 m/s and two secondary waves travelled with speeds V 2 and V 2 * of 4.79 and 15.7 m/s, respectively. The speed of secondary wave V 2 deviated only 2.8% from the main wave speed V 1 at the lumen surface. For the traveling impulse condition, Figure 5b demonstrates one single broad wave front propagating at the lumen surface traveling with a velocity V 1 of 6.99 m/s. This is in good agreement with the imposed loading conditions (V z of 7 m/s). However, at the skin level, two wave fronts are visible. The main wave front travels with a velocity V 1 of 6.99 m/s, in very close agreement with the imposed V z . The secondary wave traveled with a velocity V 2 of 11.70 m/s. Note that in both panels of Figure 5 wave reflections are not completely absent due to the non-ideal behavior of the infinite elements. For example, wave reflection is visible at a time (7 ms (lumen) and 9 ms (skin) in the TOF plots of Figure 5b. The right panels of Figure 5 show the temporal profile of the y-displacement at two locations along the length of the model (z = 20 mm and z = 30 mm) for the top of the lumen and skin surface, to replicate the signal measured by LDV. Under shock impulse loading conditions, the lumen displacement profile is characterized by multiple peaks and it has a different shape for z = 20 mm compared to z = 30 mm. For the traveling load condition, the displacement patterns have clearly one traveling peak, allowing an easy identification for LDV analysis. For both loading conditions represented in Figure 5, the tracked trajectory using the Radon transform corresponded well with the main displacement peak. Figure 6 shows the spatio-temporal displacement (y-direction) plots at the skin surface when considering traveling pressure pulse velocities V z of 4 and 10 m/s. The main wave front propagated with a speed V 1 of 5.05 m/s (26.3%) and 9.75 m/s (−2.5%) for an imposed V z of 4 and 10 m/s, respectively. For all loading conditions, the secondary wave fronts propagated with a higher speed V 2 than the main wave front, as indicated on the top of each plot in Figure 6.
Comparing the results of all imposed traveling pulses in Figures 5 and 6, the maximum amplitudes of the main wave front (indicated with speed V 1 ) were 0.06, 0.10 and 0.08 mm for V z of 4, 7 and 10 m/s respectively, whereas the maximal amplitude of the secondary wave fronts (indicated with speed V 2 ) remained constant (0.030 mm for V z of 4, 7 and 10 m/s respectively). The estimated spatial lags between the wave front propagating at the lumen and skin surface increased when V z increased (1.20, 10.48 and 19.05 mm for V z of 4, 7 and 10 m/s at 5 ms, respectively).

Depth of the Carotid Artery
Another parameter that might affect the wave propagation from the arterial lumen to the skin surface is the depth of the artery. To investigate this effect, two additional depths of 5 and 15 mm were considered while simulating a traveling pressure with velocities of 4, 7 and 10 m/s. Figure 7 shows the spatio-temporal displacement (y-direction) plots for the skin, revealing a general decrease in the maximum displacement amplitude of the main wave front (indicated with speed V 1 ) with increasing depth of the artery. Indeed, the maximum amplitude of the main wave front varied from 0.12, 0.11 and 0.08 mm in the model with d = 5 mm for a PWV of 4, 7 and 10 m/s, respectively, to 0.03, 0.06 mm and 0.05 mm in the models with d = 15 mm for a PWV of 4, 7 and 10 m/s, respectively. Nevertheless, the variation of arterial depth did not affect the wave speed estimation of the main wave front at the skin surface for V z of 7 and 10 m/s, where the agreement was optimal with the imposed conditions (maximal deviation of −5% between estimated wave speed at skin surface and imposed PWV). However, for an imposed V z of 4 m/s, estimated V 1 at the skin surface deviated 4% and 33.5% from imposed velocity for 5 and 15 mm arterial depth, respectively. The secondary wave fronts propagated with a variable velocity (V 2 ) which is indicated on the top of the plots in Figure 7. This figure also clearly demonstrates a different arrival time of the main wave front at the skin surface for both arterial depths: ±4.5 ms for d = 15 mm vs. ±2 ms for d = 5 mm due to the difference in traveling distance for the main wave front to reach the skin. Furthermore, the spatial lag between the main wave front location at lumen and skin surface increased with increasing arterial depth and increasing V z : the spatial lag is 0.15, 5.50 and 9.36 mm in the models with d = 5 mm for imposed V z of 4, 7 and 10 m/s respectively and increased to 4.28, 15 and 22.50 mm in the models with d = 15 mm for imposed V z of 4, 7 and 10 m/s respectively.

Patient-Specific Carotid Artery
The chosen mesh for the patient-specific carotid artery model consisted of 2,095,200 elements (average maximum edge length 0.40 mm and average minimum edge length 0.15). Simulated shear waves propagated with a velocity of 4.72 m/s (−2.07% deviation compared to the theoretical value of 4.82 m/s). Figure 8 shows the radial displacement inside the patient-specific CA model for four time instances, when imposing a travelling pressure impulse (V z of 7 m/s) inside the arterial lumen. At the time 2 ms, the traveling pressure is localized near the proximal end of the artery (z = 0 mm) and the induced displacement has a maximum amplitude of 0.03 mm. At 4 ms, the induced lumen displacement due to the applied load is visible at the position z = 20 mm and the wave has almost reached the skin level. When the pulse passes through the bifurcation (time 6 and 8 ms), the wave propagation evolves into a more complex displacement pattern: secondary perturbations are generated and the wave propagation direction from lumen to skin surface has a quasi-horizontal trajectory due to the bending of the external carotid artery (ECA). Maximum displacement magnitude is then 0.10 mm at the lumen level and 0.06 mm at the skin level.  Figure 9 shows the TOF plots of the out-of-plane displacement (y direction) along the centerline distance at the level of the top of the lumen and along the artery long axis for the skin of the patient-specific CA model. The travelling wave front at the lumen surface is composed of a wave front which changes speed along the length of the model due the geometrical changes of the arterial upper wall compared to the centerline of the artery. Three different speeds were calculated for the main wave front, i.e., 7.09, 10.1 and 14.5 m/s, as reported in Figure 9. At the skin surface, the main wave front has a speed of 6.78 m/s, which is in good agreement with the imposed V z (−3.1%). A secondary wave front (maximum amplitude of 0.03 mm compared to 0.06 mm for the main wave front) with propagation velocity of 9.79 m/s was also observed. As with the CCA model, the temporal profile of the y-displacement at locations z = 20 and 30 mm is shown in the right panel of Figure 9, where the red markers indicate the wave trajectory determined using the Radon transform. Both displacement signals show a pattern in which the main displacement peak, representing the main wave front, is easy identifiable.   Comparing the results of all imposed traveling pulses in Figure 10, it is clear that the maximum amplitudes of the main wave front at the skin surface (indicated with speed V 1 ) were 0.05 mm for V z of 4 and 10 m/s and 0.06 mm for V z of 7 m/s, whereas the maximal amplitude of the secondary wave fronts (indicated with speed V 2 ) was 0.03 mm for V z of 4, 7 and 0.04 mm for V z of 10 m/s.

Discussion
In the last decade, few studies have investigated the potential of LDV-based methods to assess local PWV estimation from skin vibrations generated by the pressure pulse traveling inside the underlying artery. LDV measurements have been performed with a two-beam LDV both in vitro (phantom set up mimicking arterial pulse propagation), and in vivo (CCA in healthy volunteers) [16,17]. Preliminary results have assessed the technical feasibility of LDV to accurately measure the displacement waveform at two skin locations, but also emphasized that the complexity of the skin motion prevented accurate estimation of the pulse transit time between the two measurement sites, resulting in an uncertainty in the PWV results. The transit-time error depends on many factors: shape differences between the two recorded pulses, perturbations of the LDV outputs due to a non-optimal or too weak a degree of light reflection, the existence of spurious reflections, artifacts due to secondary sources of motion (e.g., swallowing, any involuntary body movement, movement of the LDV and vibrations coming from the environment) as well as limitations of the algorithm used to estimate pulse transit time [16]. An in vivo study on 48 patients [27] has shown that signal analysis is still the most difficult challenge to validate local PWV through an LDV-based method. As expected pulse transit times are small (ranging from 2.5 to 1 ms for a beam separation of 10 mm when assuming a PWV between 4-10 m/s). Thus, a small inaccuracy in the pulse transit time estimation can lead to a considerable error in PWV estimation. Several algorithms are currently available to determine the pulse transit time. Some of them use the first or second derivative of the displacement signal (e.g., second derivative method, tangent-intersect method) but these methods are sensitive to noise [28]. The cross-correlation method measures the similarity of two waveforms as a function of a time-lag applied to one of them but fails when comparing signals with shapes that are too different [17]. Common algorithms to estimate the pulse transit time from in vivo measurements are based on the identification of a reference point, typically the dicrotic notch or the systolic foot [29]. However, the identification of the reference point can often not be assessed unambiguously in a substantial portion of the data [27]. Further studies have tried to improve local PWV estimation by increasing the number of measuring points from two to four with 15 mm spacing between them. PWV obtained by such LDV measurements in an in vitro set-up consisting of an elastic tube embedded in tissue-mimicking materials, showed a good agreement with measurements performed by the more established ultrasound-based pulse wave imaging [18]. The LDV technique was also able to distinguish phantoms with greater stiffness (higher PWV) but not to determine the precise location of a stiffer inclusion due to the low spatial resolution of the measurements.
In this work, we have presented a finite element model to gain a more fundamental understanding of wave propagation from mechanical perturbations initiated by arterial pressure pulses in a healthy carotid artery model, to the skin surface (measured by the LDV-technique). As explained in Section 3.1.2, the choice of numerical settings such as time increment size, numerical damping and mesh size, together with the material parameters (e.g., Poisson's ratio) highly influenced the computational efficiency and the reliability of the wave propagation simulations. Therefore, the appropriateness of the chosen numerical settings was initially verified by comparing the simulated displacement patterns of waves where a theoretical solution was at hand (shear and Rayleigh waves, as described in Section 2). It should be noted that these analytical solutions assume an infinite (shear waves) and semi-infinite (Rayleigh waves) isotropic linear elastic material, whereas our final model did not meet these geometrical requirements. Despite this limitation, the simulation results showed that the chosen numerical settings made it possible to reproduce a shear wave propagation speed with a deviation of −2.0% from the analytical value, whereas Rayleigh wave speed propagation deviated only 0.9%. Additionally, a good match between the spatial and temporal displacement characteristics of simulation and theory was observed (see Figure 3). The geometry of the patient-specific model introduced a higher degree of complexity which considerably complicated the test of the computational settings. Only shear waves were simulated in this set-up, resulting in a shear wave speed propagation that differed by −2.0% from the analytical value.
Consecutively, a parametric study was set up to study the impact of different arterial pressure loads and geometrical conditions on the induced skin displacement in a controlled environment. Tracked displacement patterns showed the coexistence of a primary (higher amplitude) and secondary wave front travelling with different velocities both at the boundary of the lumen and the skin surface. Figure 11 summarizes the estimated wave velocities for the primary (V 1 ) and secondary (V 2 ) wave fronts at the skin level for all simulated cases. Initially, a localized, non-propagating short pressure impulse was used to induce free wave propagation in the soft-tissue mimicking material surrounding the CCA model. Tracked displacement patterns in the CCA model showed a main wave front propagating with a velocity V 1 of 4.66 m/s at the lumen surface and a similar wave speed was found for the first secondary wave front at the skin level (V 2 = 4.79 m/s). These values approached the analytical value of shear (4.82 m/s) and Rayleigh (4.59 m/s) wave speed, however, Figures 4 and 5 also clearly showed a complex wave pattern with multiple wave fronts due to phenomena such as wave recombination, refraction and reflection, and therefore the relation between the observed wave patterns and the ideal free wave propagation modes is unclear. Additionally, it should be noted that we are only analyzing one displacement component in the TOF-plots (component relevant for LDV-analyses), whereas the traveling waves will have displacement components in multiple directions.
When considering a traveling impulse (wave speeds varying within the physiological range 4-10 m/s), it can be observed that for imposed propagation speeds of 7 and 10 m/s, all tested conditions in the CCA model yielded an estimated propagation speed V 1 close to the imposed value (differences −5.0 to −0.4; Figure 11). However, the lower imposed wave speed of 4 m/s is overestimated at the skin surface by 4%-33.5% in the model, in which the deviation increased with increasing arterial depth. It should be noted that the selected model's geometry and boundary conditions play a role in reaching a steady state regime for wave pattern development at the skin surface. However, since the tracked skin speeds for the higher imposed PWV corresponded fairly well with the imposed value for all considered arterial depths, we assume that the present transition phenomena in current modeling environment are minimal. As such, these results suggest that observed changes in wave velocity at the skin level may not fully reflect the effective changes in the arterial pulse wave velocity, as shown in Figure 11.
For the patient-specific CA model, the estimated lumen wave speed was composed of a split main wave front characterized by different velocities along the length of the model, as shown in Figure 9. This split wave front is caused by the differences in length and orientation of the upper wall compared to that of the centerline at x = 0 mm. The estimated wave speed corresponded best with the imposed PWV (1.3%) in the proximal part of the model (z < 20 mm), since the geometry of the upper wall and centerline is then very similar (see Figure 1b). At the skin level, almost no effect of this split wave front was observed and traveling wave fronts V 1 and V 2 were tracked by a single line. Similarly to the CCA model, estimated wave speed at the skin corresponded best with the imposed V z when considering a PWV of 7 m/s (−3.1%). The deviation was higher for a V z of 4 m/s (47.0%) and a V z of 10 m/s (36.0%). Figure 11 shows a large range in estimated speeds for the secondary waves in all simulation cases. It should be noted that in some cases, i.e., the CCA model with shock impulse loading condition and the patient-specific CA model with imposed PWV of 10 m/s, the speed of the secondary wave corresponded better with the expected or imposed velocity than the velocity of the main wave at the skin surface. These results suggest that the main wave front at the lumen level not necessarily translates into a main wave front at the skin level. However, these results again need to be interpreted with care as the modeled geometry and boundary conditions also have an effect. In other simulation cases, the nature of these secondary waves might be related to the propagation of the free vibrations of the surrounding tissue that is excited during the pulse transmission; However future work should include a more in-depth study of all displacement components and the application of advanced wave analysis techniques to better understand the nature of these secondary waves.
Results supported the conclusion that the skin displacement is not a simple overlapping of free and induced waves, but the result of a combination of multiple reflections and re-combinations which iterate along the model's length in time, and are thus hard to predict and interpret. Therefore, we can conclude that the transmission of a single pulse creates a complex skin displacement pattern that manifests with changes in displacement amplitude, spatial lag and propagation velocities (V 1 and V 2 ).
Our results can be considered as an initial step towards a better understanding of the LDV-measured skin vibration signals by studying the wave propagation physics between the lumen and the skin. The model showed that the correspondence between the actual PWV in the lumen and the PWV measured at the skin surface depended on the magnitude of the PWV in the lumen and the studied subject (arterial depth and geometry). To the best of our knowledge, this work represents the first computational study that has investigated in detail the induced skin displacement in a structural model of a healthy artery for the estimation of pulse wave propagation speed. Indeed, until now in this context, a fluid-structure finite element approach has been adopted only to study the sound generated in a stenosed artery and transmitted through biological tissues to the skin, to investigate the possibility of using the LDV technique for arterial stenosis detection [30][31][32].
Although the developed numerical model has allowed us to study the propagation of a single traveling impulse under prescribed conditions in a simplified environment, it rests on several assumptions. For instance, the artery was approximated as a straight cylinder in the CCA model without any wall thickness. The absence of modeling the arterial wall did not allow the control of arterial stiffness in a direct way (e.g., increasing the Young's modulus of the material), thus arterial stiffening was indirectly simulated by varying the speed of the propagating pulse. Moreover, the effect of the surrounding biological structures (e.g., muscle, fat, bones) was not reproduced and the soft tissue was modeled as an isotropic homogenous viscoelastic medium without any biological discontinuities. It is generally known that soft tissue is more viscoelastic than PVA-material and this will have an extra damping effect on the wave patterns reducing all high frequency patterns. Also, the boundary conditions to minimize wave reflections did not work optimally due to the nature of the infinite elements (see Figure 1). The loading conditions did not replicate the shape of a physiological pulse nor its periodic course, therefore the actual mechanical interaction of the pressure pulse and the lumen wall was not realistically simulated. In addition, we did not account for the fluid-structure interaction or any vibrations that may arise from instabilities in the flow field, propagating to the wall and skin. It should be kept in mind that the accuracy of the presented results depends on the chosen finite element mesh size, interpolation grids in Matlab and the projection angle of the Radon transform. This is also related to the available computer power and performance requirements, which were very demanding for the FEM simulations of the patient-specific model (72 h on a Linux cluster using a single computing node containing two Intel E5-2670 processors each with 8 cores, 2.60 GHz processor base frequency, 64 GB of memory and Fourteen Data Rate InfiniBand). The main limitation of the current modeling environment is the lack of validation. However, we believe that the good correspondence between simulation and theoretical results for shear and surface wave propagation (see Table 1 and Figure 3) lends confidence to the current modeling results.
The developed finite element model offers the advantage of visualizing and investigating the within-tissue wave propagation mechanics between lumen and skin, generated by the traveling pulse wave, which is difficult to obtain from experiments. Furthermore, it allows one to analyze multiple components of the displacements in 3D and gives access to other variables such as stresses. Compared to the LDV technique, the simulations provided a higher spatial resolution and adjustable temporal resolution (10 kHz and 0.2 mm over a 50 mm vs. 500 kHz and 15 mm over 45 mm [18]) and are not contaminated with any measurement noise. This high spatial resolution allowed us to employ an alternative method to the traditional LDV algorithms to estimate PWV, i.e., the Radon transform, an approach which is commonly used in the field of elastography to estimate wave speeds [33]. In contrast to the available methods for the LDV technique, this method provides a robust way of estimating wave speeds by yielding a linear wave trajectory based on the maximal Radon-projected displacements generated by the propagating wave. As this method requires a higher spatial and temporal resolution in combination with more measurement points, it might be more productive in the context of the recently developed LDV system consisting of 6-beams spaced 5 mm apart (total length 25 mm) by Li et al. [34]. Nevertheless, the developed simulation framework in this study can act as a validation tool for newly proposed LDV-processing techniques, since the simulated displacement data are not contaminated by any of the complexities present in an actual experiment (noise, limited spatial and temporal resolution, etc.) and it offers direct access to the actual PWV in the lumen.
In addition to patient-specific intrinsic factors affecting the wave propagation from lumen to skin (e.g., carotid depth-as considered in this study), there are also some technical and operator-dependent factors influencing the LDV measured PWV value such as laser signal strength and the localization of the measurement points. The latter issue might be improved by providing a real-time data visualization of the spatio-temporal plots to the operator that demonstrates wave propagation at the skin, as shown in the results section of our model. Indeed, this real-time tool may help the user to identify LDV positions where the LDV device is correctly aligned with the artery's axis and without any signal contamination caused by other vibrating structures (e.g., jugular vein). Combining this real-time tool together with palpation can thus improve positioning and alignment of the LDV device with the artery, potentially resulting in a more robust and accurate PWV estimation.

Conclusions
We have developed a versatile computational model which has made it possible to study induced wave patterns under prescribed arterial pressure conditions in a controllable environment. The model gives access to wave propagation along 3D wave paths and permits both a visual and quantitative analysis which is impossible to reproduce in any experimental set up. Our study has revealed the coexistence of multiple wave fronts at the skin surface that may affect the accuracy of local PWV estimation. Increasing the number of measuring points may provide a more detailed understanding of the skin displacement. In spite of its limitations, the computational model developed here provides a starting point for further analysis of the complex interaction between the displacement of the arterial wall due to the pressure pulse and the way that the energy due to this displacement is propagated through the surrounding soft tissue to the skin surface.