A Computational Method for Analyzing the Biomechanics of Arterial Bruits

The abnormal narrowing of a blood vessel or a “stenosis” can occur at a variety of locations in the vasculature including systemic, pulmonary, and cerebral arteries. The stenoses, especially those occurring in arteries, usually lead to altered local flow patterns, which are known to create distinct bruits (or murmurs). For example, aortic stenosis (narrowing of the aortic valve opening area) can create an ejection murmur during systole [1,2], while coronary artery stenoses generate a distinct murmur during diastole [3]. The bruits are generated as elastic waves and propagate through various organs/tissues, such as lungs, bones, muscles, and fat. As the bruit signals propagate through the tissue, they undergo damping, reflection, and diffraction, and the resulting signals, if strong enough, can be heard on the epidermal surface and used for diagnosis. This detection technique, termed as auscultation, is usually carried out with a stethoscope and has been the primary method for initial screening of heart conditions for centuries [4,5].

Even though this noninvasive and inexpensive technique can provide valuable diagnostic information, the practice of auscultation is virtually unchanged since the invention of the stethoscope and suffers from a lack of specificity and sensitivity. Fortunately, the emergence of low-power, low-cost compact acoustic sensors, advanced signal analysis algorithm, and powerful portable computational devices provides the technological support needed to transform this technique into a more accurate automated diagnostic modality that could potentially play a role in traditional medicine as well as telemedicine [6]. Compared with traditional auscultation, automated auscultation could extract more useful information from these bruits. Moreover, such a modality could be used for longitudinal tracking of patient health. For all these reasons, there is a renewed impetus to transform auscultation into a quantitative and accurate method [3,7–11].

When it comes to the generation of these bruits, studies have converged to the idea that hemodynamic wall pressure fluctuation in the poststenotic region is the main source [12–15]. On the other hand, there is still a lack of understanding of the physics implicated in the propagation of bruits through the intermediate organs and tissues. This lack of understanding at least partially explains the inability to overcome the relative lack of specificity and sensitivity of auscultation-based diagnosis [8,9] and new physical insights are needed to make any significant improvements in this diagnostic modality. Experiments have limited ability to provide such physical understanding due to the fact that most of the data accessible in experiments lies on the epidermal surface and little about sound propagation inside the tissue can be deduced from the surface signal. On the other hand, computational modeling is capable of providing detailed information regarding the bruits at any location inside the domain of interest and this is the primary motivation for developing the current computational framework.

Conventional analysis of wave propagation in elastic media separates the roles of the compression and shear waves. In biological materials, for frequencies in the $O(100 Hz)$ range, the wavelength of the compression wave is usually $O(10 m)$⁠, while that of the shear wave is $O(0.01 m)$ [16–18]. Since the characteristic length of human thorax is $O(0.1 m)$⁠, the long wavelength of the compression wave makes it less affected by structures, such as lungs and bones. Researchers have developed microphones that can separate the compression component from the signal and used it to achieve better directionality [3,19]. On the other hand, the relatively short wavelength of the shear wave makes it suitable for source localization. Owsley and Hull [20] designed experimental studies wherein they found that the shear-wave energy field could accurately locate the source when there were no obstacles in the phantom. Inclusion of “obstacles” (dog ribs in the case of this study), however, complicated the situation, making the source localization more accurate for certain frequency bands, but less accurate for others. Royston et al. [21] used analytical models to study the compression/shear wave generated by a dipole source embedded in the half-space. However, the contribution of the compression/shear component was approximated by doubling the solution obtained in an infinite space, which did not accurately account for the effect of the surface. Ramakrishnan et al. [17] took into consideration the inhomogeneity of the medium and the shear-wave effect when studying the wave propagation in a two-dimensional human thorax model. However, the wave was generated by a prescribed explosive source located at the mitral valve. Seo and Mittal [15] used the flow in a modeled stenosed artery as the bruit source and studied the propagation of the compression wave in the tissue layer. All of these studies point to the need to gain further insights into the propagation mechanism of the compression and shear-wave components of bruits in biological tissues.

In this paper, a computational framework for simulating the generation and propagation of bruits in biological tissues, as well as analyzing the individual roles of the compression wave and the shear wave, is provided. This framework is used to study the generation and propagation of bruits from a modeled stenosed artery. The present study can be viewed as the extension of the study of Seo and Mittal [15], which did not include the effect of shear-wave propagation. This effect is included in the current computational model, and the main objective of the study is the delineation of the effect of the compression and shear-wave components on the propagation and detection of arterial bruits. The problem is tackled with a one-way coupled hemodynamic–acoustic method. First, the flow in a modeled stenosed vessel is simulated. Subsequently, the pressure fluctuations from the flow serve as the input for the acoustic simulation. The medium is treated as a homogeneous linear elastic material, and contributions from the compression wave and the shear wave are separated through a classic decomposition. The accuracy of the computational method is verified by comparison against an analytically derived solution for a point source. Simulation results are used to explore the physics of bruit generation and propagation.

As shown in Fig. 1, the stenosed artery is modeled as a two-dimensional constricted channel similar to the one employed in Ref. [15]. Even though this is not a physiological geometry, it captures the key features of a stenosed artery. The constriction is asymmetric and has the following profile:

where $ρf$ is the density of the fluid, $Umax$ is the maximum centerline velocity at the inlet, and constants $A$ and $B$ are set to 0.75 and 1.5. The maximum and minimum flow rates resulting from this pressure profile are $0.55UmaxD$ and $0.03UmaxD$⁠, respectively. The Strouhal number based on the diameter and the maximum centerline inlet velocity is $St=fD/Umax=0.024$⁠. For typical values in the human aorta where $D=2.0 cm$ and $Umax=0.98 m/s$⁠, this Strouhal number corresponds to a heart rate of $60 bpm$⁠. The Reynolds number in this study is $Re=UmaxD/νf=2000$⁠, where $νf$ is the kinematic viscosity of the blood. It is worth noting that it is the intensity of the stenosis-induced jet that has the dominant effect on the sound generation, rather than the profile of the pressure drop [15].

The bruits generated by the abnormal blood flow propagate through an isotropic, homogeneous tissue layer. Even though the current computational framework can in-principle include the blood vessel wall, this feature is not included here so as to maintain the simplicity of the model. The exclusion of the blood vessel wall is not expected to change the results significantly, since the wavelength of the most energetic components of the bruits is much longer than the wall thickness. The density and compression wave speed of the blood and the tissue are obtained from Ref. [22], and these values are set at 1.05 g/cm³ and 1500 m/s and 1.20 g/cm³ and 1720 m/s, respectively. The shear modulus (⁠$G$⁠) to bulk modulus (⁠$K$⁠) ratio of the tissue layer is set to $2×10−4$⁠, and this is based on the widely used tissue-mimicking silicone gel, Ecoflex-10. The sound signals from various locations on the epidermal surface are monitored and analyzed.

where $Ui$ is the flow velocity vector, and $P$ is the pressure. A sharp-interface immersed boundary method is used to treat the complicated geometry. The fluid region is discretized by a $768×128$ nonuniform grid. The grids are clustered around the stenosis with minimum spacing $Δx=0.01D$ and are stretched toward the inlet and outlet. More details about this immersed boundary formulation and the flow simulation can be found in Refs. [15] and [24].

where $ui$ is the acoustic velocity perturbation vector, $p$ is the acoustic pressure perturbation, $cf$ is the speed of sound in the fluid region, and $D/Dt$ represents the material derivative. The incompressible Navier–Stokes/LPCE method is a two-step, one-way coupled approach suitable for fluid-induced sound at low Mach numbers and more details can be found in Ref. [26].

where $ui$ is the acoustic velocity vector in the tissue, $pij$ is the stress tensor, $ρs$ is the structural density, $δij$ is the Kronecker delta, and $λ and μ$ are the first and second Lame constants of the material, respectively.

As can be seen here, the inhomogeneity of the domain is accounted for by the space-dependent material properties (⁠$ρ(x⇀),λ(x⇀)$⁠, and $μ(x⇀)$⁠), where $λ=ρfcf2$⁠, $μ=0$ in the blood region, as well as the Heaviside function $H(x⇀)$⁠, which equals one inside the blood region and equals zero in the tissue layer. Equation (6) is discretized by a sixth-order compact finite difference scheme [27], and advanced in time by a four-stage Runge–Kutta method. The flow Mach number $M=U/c$⁠, where $c$ is the speed of sound, is set to 0.01 in the current study, instead of the physiological value, which is $O(10−3)$⁠. This tenfold increment is employed to ameliorate the stringent time-step size constraint imposed by the large gap between the velocity scales of the flow and the sound. The same approach is adopted in other studies [15,28] as well and has been shown to have no significant effect on the detected bruits. The stability constraint is more stringent for the acoustic simulations and in order to maintain overall stability of the computational method, the incompressible flow simulation time-step is divided into 16 substeps for the acoustic simulation. The intermediate values of flow variables between two flow time-steps are obtained through a second-order Lagrangian interpolation [26].

A uniform Cartesian grid is used to discretize the acoustic domain with a grid spacing of $0.02D$⁠. The shear-wave length is about $0.14D$ at $St=10.0$⁠, which is resolved by around seven points. This resolution is guided by our previous work [15] and grid refinement studies indicate that this grid provides sufficient resolution to accurately resolve the quantities of interest. The accuracy of the computation is also established in an analytic study, which is described in Sec. 3.2. A bilinear interpolation is used to interpolate the results of the flow simulation onto the acoustic grid. A zero-stress boundary condition is imposed at the epidermal surface, while the nonreflecting energy transfer and annihilating (ETA) boundary condition [29] are applied on the other boundaries.

respectively. Thus, the two components of acceleration can be obtained via a postprocessing of the deformation and stress distributions in the tissue layer. Note that conventional stethoscopes sense the normal acceleration of the skin surface and we therefore focus primarily on this quantity.

To ensure that the results are presented consistently, results from both the flow simulation and the acoustic simulations are nondimensionalized by the same characteristic parameters: velocity scale, $Umax$⁠; length scale, $D$⁠; time scale, $D/Umax$⁠; and pressure/stress scale, $ρUmax2$⁠.

Figure 2 shows the snapshots of the vorticity and pressure field in the modeled artery. As the flow starts to accelerate through the stenosis, a vortex street consisting of fairly evenly spaced counter-rotating vortices is formed and this is convected into the poststenotic region. As the flow decelerates, the shear layers in the poststenotic region become unstable and complex patterns of vortex mergers can be observed. The pressure distribution in the poststenotic region is driven by the vortices with the centers of these vortices coinciding with the regions of low pressure. It has been shown in Ref. [15] that arterial bruits are well correlated with the time-derivative of the integrated pressure force on the arterial wall, instead of the pressure fluctuation in the blood flow region. Thus, the time histories and the corresponding spectra of the time-derivative of the pressure on the upper lumen at different downstream locations are plotted in Fig. 3. The lower limit of the $St$ in the spectra corresponds to the fundamental frequency of the pulsatile flow, i.e., $St=0.024$⁠. The general trend of the spectra shows initial decline after the fundamental frequency, but more energy is generated in the high-frequency region due to the stenosis-induced pressure fluctuation. The strongest pressure fluctuation is observed around $St=0.9$⁠. Based on the spectra, the intensity of the time-derivative of the pressure is calculated and plotted along the upper lumen, as shown in Fig. 4. This plot clearly shows that the source of the bruits is located around $4D$ downstream of the stenosis.

The computational model for the acoustics as well as the method for computing the compression and shear waves associated with arterial bruits presented here are new, and the fidelity of this method needs to be verified. This is particularly important since the relatively short wavelength of the shear wave increases the resolution requirement. It is therefore desirable to establish that the current method accurately resolves the propagation of all the relevant waves in the tissue layer. In order to accomplish this verification in the most comprehensive manner possible, we employ the exact solution associated with a point source in the infinite domain and then incorporate the time history of the actual sound source obtained from the above flow simulation. This allows us to examine the accuracy of the computational method for the entire range of relevant frequencies.

where $H(1)$ is the Hankel function of the first kind, and $kp$⁠, $ks$⁠, and $r$ are defined as $kp=ω/cp,ks=ω/cs, and r=|x⇀|$⁠. In the above equations, $cp$ and $cs$ are the compression and shear-wave speeds, respectively. It is also noted that the first term of the Green's function represents the compression wave, while the second term is solely related to the shear wave. Finally, the acceleration at any location is obtained as $ai(x⇀,ω)=−ω2Gij(x⇀,ω)〈sj〉(ω)$⁠.

As shown in Fig. 5, the acceleration point source given by $si=(0, s2(t))$ is placed at the origin, $(0D,0D)$⁠, and we monitor the response at two points located at $(2D,0D)$ and $(0D,2D)$⁠. The signal at one location in the flow simulation (⁠$4D$ downstream of the stenosis on the upper lumen wall) is used as the point source. Its spectrum is shown in Fig. 3(b). The shear to bulk modulus ratio here is also set to $2×10−4$⁠, and a uniform grid with the same spacing, $0.02D$⁠, is employed for the spatial discretization of the domain. At this resolution, the shear wave at the highest frequency, $St=10$⁠, is resolved by seven points per wave length. Consistent with the exact solution in free-space represented by Eq. (10), the ETA boundary condition is applied at all the boundaries to allow for the transmission of the waves.

Figure 6 plots the spectra of the vertical acceleration signal at locations (a) and (b). Also, the exact solutions obtained from the Green's function are plotted in this figure along with the numerical solutions. The numerical solutions match well with the exact solutions even at $St=10$⁠, which suggests that seven points per wave length provides adequate resolution of the acoustic wave. However, there are some noticeable differences between the numerical solutions and the exact solutions for low frequencies around $St∼10−1$⁠. These discrepancies can be attributed to the inherent limitations of the ETA boundary treatment for low-frequency band [29] and not to any resolution issues. Fortunately, the frequencies of bruits are usually orders of magnitude higher than the fundamental frequency, i.e., the heart rate [32]. Thus, the simulations still possess high fidelity in the frequency range of interest here.

With the bruit propagation modeling and simulation verified, we turn to examining the characteristics of the propagation of the bruits. Figure 7 shows the computed instantaneous vertical acceleration field when the intensity of the hydrodynamic source (⁠$DP/Dt$⁠) reaches its maximum. The compression component and the shear component are obtained through the aforementioned decomposition. As shown in the figure, the wave fronts of the compression component align parallel to the arterial lumen and the high intensity region is located immediately downstream of the stenosis. On the other hand, the shear component is transmitted radially away from a region localized in the vicinity of the stenosis. The contours clearly indicate that the energy of the shear component propagates through the tissue in an oblique angle to the lumen, and its effect on the surface signal will be felt more at locations that are significantly downstream of the stenosis. This phenomenon is related to the characteristics of the shear wave. Since the fluid cannot sustain shear movements, shear wave can only be created at the interface of the blood region and tissue layer. The pressure force exerted by the blood on the interface acts as a point source aligned in the vertical direction. Thus, the generated shear wave will mainly propagate in the downstream direction. Its region of influence will grow as the wave travels, forming the oblique wave pattern shown in Fig. 7(b).

Despite the different wave patterns, both plots in Fig. 7 suggest that the main source of the bruits is located in the poststenotic region, where the wall pressure fluctuation is the strongest. This observation is consistent with the results of the flow simulation. Since one issue of interest here is to localize the source of the bruits using auscultation, the following analyses focus on the acoustic signal measured on the epidermal surface after the stenosis.

To study the effect of shear modulus (⁠$G$⁠) on wave propagation, a separate simulation in which the shear modulus is set to zero is also conducted. Figure 8(a) shows the comparison of the surface acceleration at $4.5D$ downstream of the stenosis from simulations with and without the shear modulus. This specific location is chosen because it has the highest signal intensity, which is shown later. Assuming the influence of the fundamental frequency is negligible above its tenth harmonic, the analyses focus on the high-frequency range (⁠$St∈[0.24 ,10]$⁠), which carries most of the energy of the bruits. Interestingly, the presence of shear-wave propagation mechanisms has little effect on the vertical acceleration spectrum at this location. To further investigate the effect of the shear modulus, the decomposition described before is applied to separate the shear component from the compression component and the result is shown in Fig. 8(b). When $St$ is below 0.9, the compression wave and the shear wave share similar amplitudes, but the phase difference between them results in a total signal with a much lower amplitude. When $St$ is above 0.9, the shear-wave energy decreases sharply, leaving the compression wave as the dominant component. It is noted that the compression wave spectrum in Fig. 8(b) is different from the signal spectrum for $G=0$ in Fig. 8(a) in the low-frequency range. This signifies the fact that the shear and compression wave mechanisms are not linearly additive. Instead, the inclusion of the shear modulus induces complex interaction between the compression wave and the shear wave at the surface, which can only be captured by including both effects simultaneously. The importance of including the shear modulus into these bruit propagation models is further demonstrated by the signal at $x/D=15D$ in Fig. 9. As shown in Fig. 9(a), the influence of the shear wave is more apparent at this location. Significant differences still exist when $St>1.0$ in Fig. 9(a), while the corresponding decomposition in Fig. 9(b) shows clear compression dominance in this region.

Motivated by the notion of using quantitative auscultation to localize the sound source, we examine the spatial distribution of the surface signal, and the band-limited spectral energy of the total vertical acceleration, which is calculated from the spectrum using $f2=∑St=0.2410〈f〉$⁠, is plotted along the epidermal surface in Fig. 10. As shown here, the intensities from both cases overlap very well up to $8D$ after the stenosis, indicating a trivial contribution from the shear component in this region. However, the case with shear modulus produces a stronger surface signal further downstream, due to the increasing contribution from the shear wave. In Fig. 10, both simulations with and without the shear modulus put the source of the bruits around $4.5D$ after the stenosis, which is in good agreement with the actual source location from the flow simulation. It is worth noting that the prediction here does not contradict with the observation in Ref. [15], where the peak is located between 5D and 6D, since the latter is predicted by the vertical velocity fluctuation on the epidermal surface.

The simulation results (Figs. 7 and 10) show that the effect of the shear wave on the auscultation signal is felt further downstream from the stenosis. This implies the thickness of the tissue between the lumen, and the epidermal surface should have a differential effect on the contribution of the compression and shear wave. The thickness of the intervening tissue may vary significantly based on individual anatomy (body size and body-mass index) as well as the artery of interest. For instance, peripheral arteries such as the carotid and femoral lie closer to surface, whereas other vessels commonly implicated in arterial occlusive diseases such as the coronary and iliac arteries, as well as the aorta, lie deeper in the body.

To explore the effect of tissue thickness on the bruit propagation and detection, another set of simulations with and without the shear modulus are carried out, but with the thickness of the tissue layer reduced from $9D$ to $6D$⁠. Figure 11 shows the spectral energy distribution along the epidermal surface for both cases. A single peak at $4.5D$ after the stenosis is observed when the shear modulus is not included, while two peaks are observed when the shear modulus is included. For this latter case, the first peak is again at $4.5D$ after the stenosis, which is in line with the peak generated in the previous case. The other peak, which is at $10.5D$ downstream of the stenosis, is related to the contribution from the shear wave. Based on the wave patterns in Fig. 7, we can see that the decrease of thickness has little effect on the propagation of the compression wave, and the maximum signal intensity is observed at the same location. On the other hand, reducing the thickness of tissue layer also reduces the distance between the stenosis and the location where the shear waves with high intensity interact with the surface. Hence, the contribution from the shear wave comes into effect earlier, creating a second peak slightly downstream of the first peak. This second peak also has an intensity level that is comparable to the first peak and would therefore appear as an equally strong auscultatory signal. The appearance of this second peak could therefore complicate the task of source localization.

A coupled hemodynamic–acoustic simulation method and a wave-decomposition based analysis approach are developed to study the generation and propagation of arterial bruits. A key feature of the model is that it includes compression as well as shear-wave contributions to the bruits. The model and wave-decomposition method are successfully verified against a canonical problem solution and are then used to study the effect of shear waves on the propagation and detection of bruits generated by a modeled stenosis. The decomposition reveals that the compression wave propagates perpendicular to the arterial lumen, while the shear wave is mainly transmitted in an oblique direction through the tissue. As a result, the shear wave contributes more significantly to the surface signal further downstream of the stenosis. For a relatively thick tissue layer between the lumen and tissue surface, the shear wave has a limited effect on the auscultatory signal. For this case, the peak location of the surface signal corresponds reasonably well to the source location. However, a reduction in tissue thickness amplifies the contribution from the shear wave, creating a second peak in the auscultatory signal slightly downstream of the first peak. This “false” peak could confound the task of source localization and this issue will be explored in future simulations. Besides the thickness of the tissue, there are other factors that can affect the generation and propagation of the bruits, such as curvature of the artery and heterogeneity of the intervening tissue. The current method can be extended to study these issues in the future.

Even though it is difficult to separate compression and shear components from the surface signal in the clinic with available stethoscopes, the present computational framework can potentially facilitate the development of more sophisticated auscultation tools and signal processing algorithms. For example, the simulation results can be used to fine tune the source localization algorithm to take account of potential false peaks and generate more accurate information for diagnosis. Such analysis could also be “inverted” to obtain information about the tissue properties, and in doing so, provide a noninvasive means for tissue biopsy. Finally, it is noted that while the focus of this paper is on arterial bruits, the current methodology is equally applicable to heart murmurs where the source location lies within the cardiac chambers.

The authors would like to acknowledge the financial support from the National Science Foundation (NSF) Grant Nos. IIS-1344772 and CBET 1511200, computational resources by XSEDE NSF Grant No. TG-CTS100002, the Maryland Advanced Research Computing Center (MARCC), and the Institute for Computational Medicine (ICM) cluster.

a_y =

vertical acceleration

b =

severity of the stenosis

c_f =

speed of sound in the fluid

c_p/c_s =

compression/shear-wave speed in the tissue

D =

diameter of the modeled artery

Re =

Reynolds number

St =

Strouhal number

U_max =

maximum center line velocity at the inlet

λ, μ =

Lame constants

ν_f =

fluid kinematic viscosity

ρ_f/ρ_s =

fluid/structure density