Finite element modeling of acoustic wave propagation and energy deposition in bone during extracorporeal shock wave treatment

Starting in the 1980 s, shock wave lithotripsy (SWL) has initially exhibited great success in clinical treatments by pulverizing hardened calcific deposits in the human body (e.g., kidney stones).^1–3 Recently, extracorporeal shock wave treatment (ESWT) has been investigated as a novel application of shock wave (SW) therapies for treating a variety of musculoskeletal disorders, such as plantar fasciitis, bone fractures and non-unions, calcified tendonitis of the shoulder, and later epicondylitis of the elbow.^4–9 An increasing body of literature suggests that ESWT is capable of providing a non-surgical and relatively pain free alternative treatment modality for patients who do not respond well to conservative treatments.

During SW treatments, cavitation bubbles are generated from the tensile phase of each SW pulse. When the driving SW pressure is sufficiently high, pre-existing vapor or gas voids undergoes rapid growth during the negative pressure phase of the SW field and then collapses violently and implosively during the compression part of the field. This phenomenon is called inertial cavitation (IC). Many researchers have indicated that IC plays an important role in both SWL and ESWT by fragmenting kidney stones,^10–12 removing dead cells from the surface of the injury site and induce acute injury, triggering the body's natural healing process, or stimulating revascularization in the targeted region.^13–16 

Most ESWT devices are based on SWL technology but scaled down in size. Treatment sites are relatively close to the skin and the penetration depths are normally less than 40 mm. Although ESWT has already been applied clinically since the 1980 s, there are few specific studies clearly addressing questions of physical mechanisms of SW actions (e.g., inertial cavitation); practical ESWT parameters and protocols are usually simply determined by clinicians according to personal clinical experience. Sometimes, these subjective ESWT strategies may generate unforeseen damage in parts of the patient body, as the SW pulses might not ideally strike the area of interest due to complicated wave propagation and interactions in biological tissues/structures.¹⁷ Therefore, better physical understandings of SW propagation characteristics and energy deposition during ESWT processes may be a prerequisite for standardizing treatment dosage, optimizing future protocols, and even providing guidance to manufacturers regarding machine designs.

Earlier physical studies regarding the SW treatments have been made by simply measuring the acoustic field characteristics (e.g., peak pressure, intensity, etc.) in a water bath, then extrapolating the experimental results to soft tissues.^18–21 Later on, some simple numerical models were developed to predict SW propagation in water or even to some degree in homogeneous soft tissue.^22,23 However, these numerical models can neither predict wave reflections due to the assumption of axis symmetry nor contain the mode conversion to account for the elastic properties of tissue structures. More importantly, for the purpose of simplification, the presence of bone is often neglected in both experimental measurements and numerical simulation, although it has been explicitly reported that the presence of bony structures will significantly change the shape of SW focal zone, alter the SW propagation path, shield areas from SW exposure and even excite complex vibrations in the bone and along the surface at tendon bone junctions.²⁴ 

In the present paper, a Finite Element Method (FEM) model was developed based on linear elasticity equations and acoustic propagation equations to examine SW propagation and deflection in soft tissue with the presence of a mimic musculoskeletal ankle bone. The maximum negative pressure produced in the wave propagation path, which would result in the generation of cavitation bubbles, was numerically predicated as a function of variables that can be controlled by the clinician (e.g., the orientation of the SW source). High-speed photography experiments were also performed to record the cavitation bubbles generated by SW pulses with the presence of mimic bone under conditions similar to numerical settings. After verifying the proposed model by comparing the simulated data with experimental results, the strain energy density distribution in the bone was also predicted according to numerical simulations, which should be helpful for better understanding where and how the energy was deposited in the musculoskeletal structure during ESWT. The long term goal is to develop an effective tool based on integration of actual patient magnetic resonance imaging/computed tomography (MRI/CT) scans and numerical modeling so that patient-specific treatment strategies can be provided for ESWT in clinic.

In the present work, the propagation of a SW pulse in the presence of a bony structure (e.g., a talus bone) is simulated using an FEM model based on linear elasticity and acoustic wave propagation equations.

It is well known that a SW pulse is characterized by a sudden rise from the ambient pressure to the maximum pressure in the system. Meanwhile, as the velocity of wave components increases with increasing pressure, the wave will deform and steepen into a shock. Electrohydraulic lithotripsy, one of the most common technologies to generate SW pulses in water, can be modeled as an underwater explosion system. As illustrated in Fig. 1, the electrohydraulic lithotripter has a spark plug source sitting at the F1 focus of a half ellipsoidal reflector to generate spherical SW pulses, which will reflect away from the reflector and then propagate through the water bath into body tissues. In practical applications, coupling gel is used to minimize the attenuation between the water bath membrane and the skin. Here, the difference between the acoustic impedances of water and soft tissues is regarded as small enough to be neglected. In other words, the bony structure is assumed to sit in a liquid environment. Ideally, SW pulses should eventually collide with the intended target (marked as point T on the bone surface) centered at the second focus of the ellipse (marked as F2). The geometric origin point is set at the F2 focus.

The typical form of a therapeutic SW pulse that would appear at the F2 focus and the associated physical parameters are illustrated in Fig. 2. There is a sudden positive rise in pressure (about 50–100 MPa) with a short life cycle of approximately 5–10 ns, which is followed by a negative pressure wave (about −10 MPa) that lasts for a few microseconds and may cause cavitation activity. In the numerical model, the shock wave source sent out at F1 will be simply simulated as a homogeneous plane wave

with $α=3.5×105 , f=50kHz$ and P_a = 50 MPa. The tensile region of the SW would be of particular interest, as it may result in cavitation bubbles that are thought to contribute to the angiogenic responses involved in ESWT.

The dynamic behavior of SW pulses propagating through water/tissues and bony structure can be modeled with the acoustic wave equation

where $ρ$ is the density of propagating medium, $v$ is velocity vector of the acoustic wave, and $λ,μ$ are Lame constants in the medium. The mimic talus bone is assumed to be made from isotropic solid materials, e.g., acrylonitrile butadiene styrene (ABS), which can be described using the linear elastic equations. Hence, the following linear elastic equations are used to describe the properties of the mimic bone:

where $σij(i,j=x,y,z)$ is the stress tensor and $εij(i,j=x,y,z)=12(uij+uji)$ is the strain tensor. The displacement $uij(i,j=x,y,z)$ is related to the velocity vector with the formula $v=∂u∂t$⁠. The bone surface profile used in the FEM model is extracted from CT images of the bony structure. In simulations, densities and Lame parameters of water and mimic bone are set to be $ρwater=1000 kg/m3$⁠, $λwater=2190.4 MPa$⁠, $μwater=0 MPa$⁠, $ρABS=1210 kg/m3$⁠, $λABS=3264 MPa$⁠, and $μABS=1032MPa$⁠. Each FEM grid cell will be assigned values of ρ, λ, and μ depending on the material. Adaptive mesh refinement is used so that the finest grids can be adopted near the SW propagating regions to show detailed dynamic behavior.

In the FEM model, a rigid boundary condition is set for the half ellipsoid reflector, and the rest parts in water are assumed to have absorbing boundaries. At the interface between the water and bone, continuity boundary conditions should be satisfied for the normal velocity and stress tensors, which means that $vtissue,n=vbone,n$ and $σtissue=σbone$⁠.

In Sec. IV, numerically calculated regions of maximal tension indicate where cavitation bubbles would occur in laboratory experiments. Thus, this is the basis for the comparison between simulated and measured results. In order to investigate the influence of the bone placement on SW propagation and acoustic energy deposition in the bone, the vertical distance between the SW geometric focus (F2) and the intended target point (T) one the bone surface is varied within a range between −12 mm (the T point below F2 focus) and 9 mm (T point above F2 focus) with a step size of 3 mm.

As presented in the above section, an FEM model could be used to predict SW propagation properties in the presence of a mimic talus bone. In order to verify the effectiveness of the proposed model, some experimental measurements were also performed to photograph cavitation bubble clouds generated in SW fields in the presence of a mimic talus bone using a high-speed imaging system.

Figure 3 illustrates the picture of the experimental system. According to CT images of a human talus bone, a piece of mimic bony structure was constructed with the material of ABS by using Rapid Prototyping technology. The mimic bone was then mounted on an adjustable stage and immersed in the water tank. The orientation of the mimic bone with respect to shock waves was similar to that used in clinical procedures. An electrohydraulic spark discharge device, a replica of the Dornier HM-3 Lithotripter, was used as SW source.²⁰ The SW device was operated at 22 kV with a 2-Hz pulse repetition frequency. Similar to the numerical studies, a series of observations was performed by varying vertical distance between the intended target point on the bone surface and the SW F2 focus with a step size of 3 mm and five replicated measurements were perform at each position.

To build a stable cavitation field, the spark was fired continuously for at least 20 shots before every measurement. The cavitation fields generated at the negative phase of SW pulses were photographed using a high-speed camera (Imacon 200, DRS technologies, Inc., NJ). This digital camera is capable of frame rates up to 2 × 10⁷ frames/s and can collect 14 high-resolution frames with a streak option, which allows real-time bubble motion measurements. The time between frames can be varied from 50 ns to 20 ms. Eventually, the photographed cavitation fields were analyzed using matlab image processing programs (Mathworks, Natick, MA, USA).

To investigate the propagation properties of a SW as it interacts with bony structure, numerical simulations were performed first based on an FEM model. Since the simulation results have to be verified by comparing them with experimental observations and the high-speed photograph technology used for current experiments takes pictures in 2D space, only a series of 2D FEM simulation results in the x-z plane (viz., y = 0) are illustrated in Fig. 4. Each image shows the deflection of the SW's tensile component with the bone location shifting up along the z-direction. The pentagram marker in the image indicates the location of the F2 focus and the solid dot represents the intended treatment point (T point) on the bone surface. The dark bands correspond to the calculated regions of maximum tension accumulated for 150 images that indicate where cavitation bubbles would occur during SW propagation for 1.5 μs. When the bone surface is relatively far away from the F2 focus (viz., z = −12 mm), the SW pulse propagates more like in a free field, and the maximum tension region exhibits an hour glass contour that corresponds to the usual focal zone of a focused SW device. As the mimic talus bone moves closer and closer relative to the focal zone of the lithotripter along the z-direction, it is clearly shown that the deflection of the SW field becomes more significant.

In order to compare with numerical simulation results, the cavitation fields generated by SWs sending from the Dornier HM-3 lithotripter were recorded using a high-speed imaging system, which can provide a visual impression about how the wave fields changed as a bone is placed in the SW path. In the series of photos shown in Fig. 5, the bubble clusters generated by the negative pressure of SW pulses (coming from the left side of the image) clearly indicate the SW's propagation paths. As the bone height relative to the F2 focus of the lithotripter is raised from z = −12 mm to 9 mm with a step size of 3 mm (images labeled a to h), the shape of the bubble cloud region actually bends upwards, away from the bone surface. This bending phenomenon agrees with the numerical prediction, which suggests that the acoustic energy is deflected to a different region when the incident SW pulse gets closer to the bone.

In the above sections, the propagation properties of SW pulses in the presence of bony structure are investigated both numerically and experimentally. Equally instructive information is illustrated in Figs. 4 and 5, which suggests that the SW field will be deflected in the presence of the bony structure, and varying deflection angles can be observed as the bone is shifted up in the z-direction relative to the F2 focus. However, because of the linear nature of the wave propagation and elastic material property equations and the computational limitations of mesh refinement level, the calculated pressure wave field shown in Fig. 4 is more smeared out than those photographed in high-speed imaging experiments. As a result, the simulated profiles of deflected SW beam do not fully agree with the laboratory results, although the general behavior of SW propagation is similar in that there is an upward deflection of the cavitation path.

Thus, in order to further confirm the validity of the proposed FEM model for predicating the SW propagation while considering its interaction with the bone, a quantitative comparison was performed between the simulated and measured deflection angles of SW-induced cavitation fields. Figure 6(a) is a sample image illustrating how the SW deflection angle was quantified after image processing procedures done by matlab programs. First of all, histogram equalization and image contrast enhancement were applied to each image by adjusting its gray-scale intensity and expanding its histogram, which is helpful for identifying and emphasizing the region of interest (ROI; viz., the cavitation bubble region in Fig. 5 or the maximum tension region in Fig. 4). Then, by setting an appropriate intensity threshold, the image was converted to a binary image (viz., image in black and white) and edge detection was applied to ROI. As a result, the incident SW pulse and the outgoing deflecting beam are identified as the traces ending at the F2 focus and starting from the F2 focus, respectively. Finally, the deflection angle was calculated as the intersection angle (α) between the central line of the outgoing beam (solid line) and the extended line (dashed dotted line) of the central incident SW path (dashed line).

After applying the above processing scheme to all the images, the deflection angle of the cavitation field can be quantified as a function of the vertical distance between the intended target point on the bone surface (T point) and the geometric F2 focus of the lithotripter. The experimental data are plotted as Mean ± STD for 5 replicated measurements. As shown in Fig. 6(b), it is obvious that the numerical data agree well with the experimental results, which suggests that the proposed FEM model should be capable of effectively predicting the propagation properties of SW even taking into account the interaction between acoustic waves and elastic bony structures.

One can easily notice that an essential property of the electrohydraulic ESWT system described above is that, due to the SW focusing, the target treatment area (point T in Fig. 1) should be overlapped with the second focus of the ellipsoid (F2 in Fig. 1) where the maximal energy should be deposited. However, both the numerical simulations and experimental observations have shown that as SW pulses enter the body they might be deflected and dissipated depending on the material parameters of the tissues/bones in the wave path, and the deflection angle will be altered as the location of the bone changes with respect to the geometry focus (F2 focus) of the lithotripter. According to the quantified data shown in Fig. 6(b), when the distance between the F2 focus and the expected target point (T point) on the bone surface shifts by about 21 mm, the maximum deflection angle of SW field can reach ∼35°. In some cases, as SW-induced cavitation fields deviate far from the intended treatment region, undesired damage might be caused in bone and/or surrounding tissues because of unforeseen acoustic energy deposition.^17,24 Therefore, the strain energy distribution should be taken seriously in SW treatments.

According to the theoretical model, because the incident SW meets the bone/tissue interface with different acoustic impedances, the SW energy will be released and shear stress will be generated within the bone. Assuming the bony structure is composed of linear elastic material, the strain energy density $u$ in the bone can be calculated by $u=σ22E$⁠, where $E$ is the Young's modulus and Lame constants can be related to the Young's modulus (E) and Poisson's ratio (⁠$δ$⁠) through the equations of $λ=Eδ(1+δ)(1−2δ)$ and $μ=E2(1+δ)$⁠. Here, since the validity of the proposed FEM model has been verified in the above section, the distribution of strain energy density in the bone can be numerically simulated. Figure 7(a) gives a sample of calculated strain energy density distribution in bone, as the vertical distance between SW's F2 focus and the T point on the bone surface is −6 mm. It is noticed that higher energy is deposited on the bone surface facing to SW source and the energy will gradually dissipate with the increasing propagating depth in bone.

In order to investigate how the acoustic energy deposited in the bone is affected by the placement of the bony structure, the total strain energy $U$ deposited in a certain area (the circle area in Fig. 7(a)) around the intended target point (T point) is integrated through $U=∫∫su dS$⁠, for different situations. As illustrated in Fig. 7(b), the temporal evolution curves of the strain energy U integrated in the circle region around T point with a radius of 5 mm are plotted for varied vertical distance between F2 focus and T point. It clearly shows that, for the current electrohydraulic spark discharge device, it will take about 0.18 ms for incident SW pulses to reach the bone region after they are send out from the lithotripter, which can make the acoustic energy deposited in the bone achieve peak values. As the bone is placed far below F2 focus (e.g., z = −9 mm), the acoustic energy propagating into the bone is pretty low as the T point misses the real focus. As the bone is shifted up along z-direction, the T point gets closer and closer to the SW geometric focus, so that the more and more acoustic energy is delivered into the bone. When the T point almost overlaps with F2 focus (e.g., z is around 0 mm), a nearly maximum energy is achieved. However, if the bone height is set to be much higher than F2 focus, the energy deposition in the intended treatment region might also be impaired, because the mimic talus bone is an irregular structure so that most SW energy could be blocked by other parts of the bone. All of the above results suggest that, with the combination of the FEM model and CT images, more detail information about SW propagation and interaction with bone/tissues during ESWT applications can be obtained, which is valuable for clinicians to determine a patient-specific treatment strategy.

Although the comparison between the simulated and measured data suggests that the proposed FEM model should be capable to predict the propagation and deflection of SW pulses with the presence of bony structures, it is also mentioned in Sec. IV C that the simulated pressure wave field is more smeared out than the measured one due to the computational limitations of FEM mesh refinement level. In FEM calculations, it is usually recommended that the grid should be smaller than 1/6 of the wavelength. Thus, smaller grids will be applied to simulate wave motions with higher frequency (viz., shorter wavelength), which requires larger amount of computational efforts. Meanwhile, the sharper waveform distortion induced by higher driving pressure and poorer convergence caused by wave nonlinearity will also dramatically increase the computation time.

Moreover, in the present work, the simulated data were only calculated for one set of parameters and the experiments were performed by setting a mimic bony structure in the water tank. Although the water is degassed and the bony structure was constructed based on a CT scan image of a real patient, there are still discrepancy between laboratory study conditions and realistic human tissues/structures. Therefore, some efforts would be made in the future to improve the FEM model. For instance, the absorption, attenuation, the inhomogeneity of surrounding tissues or even the scattering of bony structures might be considered in future models.

Without the knowledge of wave propagation properties and energy deposition during ESWT processes, it is impossible to better understand the physical mechanisms of SW actions (i.e., cavitation and shear stress). Moreover, the efficacy evaluation of ESWT treatments cannot be well established, nor can the comparison between various ESWT devices and applications. The major objective of the current work is to investigate how a SW field would change if a bony structure exists in the path of an acoustic wave. An FEM model is developed to predict the propagation properties of SW pulses in the presence of a mimic bone. The results suggest that a SW-induced cavitation field will be significantly deflected and the deflection angles will vary with the vertical distance between the intended target point on the bone surface and the geometric focus of a lithotripter. The numerical simulations are successfully confirmed by experimental observations based on high-speed photographic technology, which indicates that the FEM simulation is a fast and easy way to model SW propagation in biological structures. Furthermore, while considering different placement schemes of the bony structures, strain energy distributions in bone are also numerically analyzed using the proposed FEM model, which will give us a better physical understanding of ESWT treatments and should be helpful for optimizing ESWT applications in clinic.

This work was partially supported by the National Basic Research Program 973 (2011CB707900), National Natural Science Foundation of China (81127901, 81227004, 11074123, 11174141, 11274170, 11161120324, and 11104140), the Natural Science Foundation of Jiangsu Province of China (BK2011543 and BK2011812), the Six Top Talents of Jiangsu Province of China (2010-WS-082), NCET-11-0236, and the Project Funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions. The authors are grateful to Brian MacConnaghy in the Mechanical Engineering Department at the University of Washington for creating the mimic ABS bone.