Cavitation induced shear and circumferential stresses on blood vessel walls during photo-mediated ultrasound therapy Open Access

Cavitation may significantly impact the function of a blood vessel. In the 1980s, several studies found that the vapor cavity inside blood vessels was responsible for vascular injuries during shock wave lithotripsy treatment.^1,2 Since then, many studies have shown possible use of microbubbles in the presence of ultrasound for therapeutic applications. Many of these investigations are based on the use of ultrasound activated microbubbles for drug and gene delivery in vivo^3–8 and in vitro.^4,9–11 Other potential applications include temporary blood–brain barrier opening^12,13 and lysis of blood clots and cell membrane.^14–17 Inertial and non-inertial cavitation inside microvessels can generate significant shear and circumferential stresses through liquid jets and microstreaming, which will lead to various biological responses in blood vessels.^18–25 For example, the shear stress induced by oscillating bubbles on the blood vessel wall can lead to the activation of the ion channel, reversible perforation of the membrane, and cell detachment and lysis.²⁶ Bubble-induced circumferential stress, on the other side, is often related to the rupture of vessels. One study showed that ultrasound contrast agents inside the vessel can cause axial direction rupture in vessel phantoms.²⁷ This suggests that circumferential stress that is tangential to the vessel results in axial opening in the vessel. Other studies have shown that the changes in circumferential stress were related to parameters such as the vessel diameter, initial bubble size, ultrasound frequency, and amplitude.^12,28,29

In recent years, studies using various numerical methods have been conducted to study the interactions between the bubble and the blood vessel wall when a bubble was placed inside a blood vessel. Many of these studies have been focused on the natural frequency of microbubbles and stresses on the vessel wall that resulted from the oscillations of microbubbles inside a blood vessel. Some studies have used the boundary element method for confined bubble modeling. Ye and Bull³⁰ developed a two-dimensional (2D) model using the boundary element method to study stresses on the vessel wall. The stresses were produced due to gas bubbles generated through perfluorocarbon droplet vaporization inside a rigid and flexible tube. They found that vessel wall flexibility and the initial bubble diameter to vessel diameter ratio affected wall stresses significantly. Later, Miao et al.³¹ and Wang et al.³² used the boundary element method coupled with the finite element method (FEM) and potential flow theory to study the 2D asymmetric bubble oscillation inside an elastic vessel. Miao et al. concluded that the circumferential stress would increase with decreased vessel thickness, vessel radius, and acoustic frequency. On the other hand, Wang et al. found that the resonance frequency of a bubble would increase and its oscillation amplitude would decrease when the bubble was inside an elastic vessel. Numerical studies with the lumped parameter model have also been used to study the natural frequency of a confined bubble. Qin et al.³³ and Martynov et al.³⁴ used the lumped parameter model and lumped parameter model coupled with the Navier–Stokes equation to calculate the natural frequency of a 2D asymmetric bubble. Qin et al. found that the natural frequency of a microbubble in a compliable vessel would increase with decreased vessel size but decrease with increased vessel rigidity. Martynov et al. also showed that a bubble confined inside a long elastic vessel had a higher natural frequency than an unconfined bubble. The finite volume model and FEM are also among the numerical methods used for calculating bubble-produced stresses in a blood vessel. Gao et al.³⁵ developed a 2D finite volume model to study the asymmetric bubble motion and the resulting circumferential stress inside a deformable pseudo-elastic vessel. They found that vessel confinement forced the bubble to take ellipsoidal shape and exert larger stresses than spherical shape. All these studies assumed free microbubbles without shells. Hosseinkhah et al.^29,36 developed a 3D FEM bubble model to study the dynamics of ultrasound contrast agents in a confined vessel. They used a modified Rayleigh–Plesset equation that included the effect of the elastic shell to define the bubble motion and its coupling with the surrounding fluid. They found that the resonance frequency would increase and the oscillation amplitude would decrease with increased vessel rigidity. They also calculated circumferential and shear stresses as functions of the initial vessel size, initial bubble size, ultrasound frequency, and amplitude. All above numerical studies concluded changes in the natural frequency and decrease in the oscillation amplitude of microbubbles due to vessel confinement. Some experimental studies had also shown substantial decrease in the microbubble amplitude under vessel confinement.^37,38 In all of these studies, the radius of the microbubble was on the same order of magnitude as that of the blood vessel, and hence, the confinement effect of the blood vessel was significant.

We have developed a novel technique, termed photo-mediated ultrasound therapy (PUT), which utilizes combined ultrasound and laser irradiation to enhance cavitation inside blood vessels. Through the enhanced cavitation effect, PUT can remove microvessels and reduce local blood perfusion, a process that can benefit many medical conditions involving neovascularization such as cancer, macular degeneration, and diabetic retinopathy.^39–42 PUT uses concurrently applied laser and ultrasound pulses for enhanced cavitation.^43–46 This technique does not require ultrasound contrast agents and depends on photoacoustic cavitation, which is induced by the photospallation effect.^47,48 The photospallation effect can result in high rarefaction pressure at the center of a cylindrically shaped optical absorber such as a blood vessel and induce cavitation activity.⁴⁹ A recent numerical study on PUT showed considerable decrease in ultrasound pressure needed for bubble growth if nanosecond laser pulses were applied concurrently with ultrasound waves.^43,50 Several experimental studies showed PUT-induced bioeffects on blood vessels such as the reduced perfusion rate and rupture of blood vessels as a function of ultrasound and laser parameters.^45,46 Circumferential and shear stresses are thought to be responsible for the bioeffects on vessels during PUT. However, how exactly cavitation bubbles impact blood vessels during PUT has not yet been studied.

The current study focuses on the cavitation produced shear and circumferential stresses under different ultrasound and laser parameters during PUT. We developed a 2D axisymmetric numerical model of a bubble confined in an elastic vessel to study the induced shear and circumferential stresses on the blood vessel wall during PUT. The bubble was assumed to be filled with air, and no elastic shell was considered as ultrasound contrast agents were not used in PUT. The model was then solved by using the FEM in COMSOL Multiphysics to calculate the produced stresses on the blood vessel wall. Throughout the simulation, unless otherwise indicated, the sizes of bubbles and blood vessels are given in radius.

In our model, an air-filled microbubble was represented by a hollow cavity, similar to the Hosseinkhah model,³⁶ as velocity and pressure calculations inside a bubble were not required for calculating stresses on the vessel wall. The microbubble was assumed to be at the center of a vessel and surrounded by blood; the vessel was further surrounded by a thick elastic tissue, as shown in Fig. 1. The semi-circular boundary of the hollow cavity represented the bubble–blood interface and was moved according to the velocity calculated from the Keller–Miksis equation for bubble dynamics, which was similar to the Hosseinkhah model.³⁶ To include the effect of laser irradiation, the laser-produced PA wave was calculated using the PA wave model presented by Li et al.⁵⁰ and was applied in the Keller–Miksis equation along with ultrasound waves. Blood was modeled as a Newtonian and incompressible fluid, while the vessel and tissue were modeled as linear elastic solids. The vessels with radii of 50 μm, 100 μm, and 150 μm with a vessel wall thickness of 4 μm were used in this study. Bubbles with the equilibrium radii of 50 nm and 100 nm were used for most of simulations. These parameters were selected because they can best demonstrate the impact of combining laser and ultrasound, and it is also consistent with that reported in our earlier studies.^43,50 This model calculated stresses on the vessel wall that resulted from bubble oscillations, both linear and nonlinear, with bubbles expanding up to 100 times their initial sizes in highly nonlinear oscillations during PUT, i.e., combined laser and ultrasound application. In this study, micro-jets due to bubble collapse were not considered as bubbles were placed at the center of a vessel and far from the boundary. Our model is, therefore, only valid for bubbles at the center of a vessel as it is based on 2D axisymmetry, which is consistent with that reported in previous studies, where bubbles have been assumed to be at the center of a blood vessel for PA cavitation during PUT.^43,50 In addition, the confinement effect of a vessel on bubble dynamics was not considered as the initial bubble size was less than 1/1000th of the vessel size (to be discussed in Sec. II B).

To study the effect of combined laser and ultrasound on cavitation inside a blood vessel, the PA wave due to the thermo-elastic relaxation process was superposed on an ultrasound wave near its negative peak. When a nanosecond laser pulse is applied to an optically absorbing medium, under stress and thermal confinement, a thermo-elastic pressure wave is generated due to instantaneous heating. This pressure wave will converge inside cylindrical absorbers (such as blood vessels), resulting in a significantly high rarefaction pressure near the center.⁴⁹ The initial pressure inside a vessel due to instantaneous delta heating is given by⁵¹ 

where Γ is the Gruneisen constant, μ_a is the optical absorption coefficient, and F(r) is the laser fluence. For blood, μ_a is 250 cm⁻¹ (assuming 532 nm laser wavelength) and Γ is 0.2. For a laser fluence [F(r)] of 20 mJ/cm², the initial pressure in the vessel will be 1.0 MPa.

The pressure wave inside a vessel at different positions resulting from a nanosecond laser pulse with the Gaussian temporal profile can be calculated with the following equation:⁴⁹ 

where p_δ is the pressure wave for delta heating and T is the laser pulse width defined at the full width of half maximum. This pressure wave p(r, t) depends on the size of the absorber (the vessel) and distance from the center of the absorber. In a cylindrical vessel, the strength of the converging PA wave increases with the size of the vessel and is maximum at the center of the vessel.⁵⁰ 

The Keller–Miksis equation was used to calculate the bubble radius as a function of time t and has the following form:^52,53

The dots represent the time derivative, R is the bubble radius, c is the speed of sound in the surrounding medium, ρ is the density of the surrounding medium, P_b is the pressure in the surrounding medium side at the bubble and medium interface, P_∞ is the pressure in the surrounding medium at infinity, σ is the surface tension coefficient, μ is the viscosity of the medium, and k is the polytropic index. The process was assumed as adiabatic (k = 1.4), and the values of ρ, P_∞, σ, μ, and c were assumed as 1055 kg/m³, 104.6 kPa, 0.072 N/m, 0.005 Pa s, and 1500 m/s, respectively. $pt+Rc$ represents the pressure due to ultrasound and PA wave. The Keller–Miksis equation calculates the bubble radius by assuming a bubble is surrounded by infinite liquid, without considering vessel confinement. Vessel confinement will result in a smaller amplitude of bubble oscillation and asymmetric oscillation such as the ellipsoidal bubble shape. In this study, however, the effect of vessel confinement was negligible as the vessel was 10³ times larger than the initial bubble size. Even for strong nonlinear bubble oscillations, the maximum bubble radius was less than 1/10th of the vessel radius that was considered in this study. Therefore, the use of the Keller–Miksis equation for bubble dynamics was justified.

Blood was modeled as a fluid domain and assumed as incompressible and Newtonian. The transient Navier–Stokes equation (5) and continuity equation (6) for an incompressible fluid were used to solve the pressure and velocity fields in the blood domain,

where ρ_o is the density, v is the velocity vector, p is the pressure, and F is the volume force vector. The blood viscosity (μ) was assumed to be 0.005 Pa s,⁵⁴ and the blood density (ρ_o) was constant and assumed to be 1055 kg/m³.⁵⁴ 

The vessel and tissue were modeled as solid domains and were considered isotropic linear elastic solids. The deformation in vessel and tissue was given by

where ρ_s is the solid density, u is the displacement vector, σ is the stress tensor, and F is the volume force vector. The vessel was given Young’s modulus of 1.5 MPa⁵⁴ and a density of 1070 kg/m³.⁵⁴ The blood vessel was assumed to be surrounded by an elastic muscle tissue with Young’s modulus of 0.5 MPa⁵⁴ and a density of 1055 kg/m³.⁵⁴ Both solids were given Poisson’s ratio of 0.49.

The finite element model was solved in COMSOL Multiphysics 5.5 (Burlington, USA) using Laminar Flow and Solid Mechanics modules. The ambient pressure was given as 104.6 kPa. In the 2D axisymmetric model, the exterior part of tissue surrounding the vessel was a perfectly matched layer (PML). The PML acts as an absorbing boundary and prevents reflections by virtually stretching the tissue (PML) toward infinity, giving rise to infinite tissue. In the horizontal direction, the ends of the fluid domain were an open boundary and the ends of the solid domain were fixed, as shown in Fig. 1(c). The open boundary condition was applied using Eqs. (8) and (9); it assumes zero gauge pressure, and the fluid can freely enter and exit through the boundary. The fixed boundary condition in the solid domain was applied using Eq. (10). The semi-circular boundary of the hollow cavity representing the bubble–blood interface was a no slip boundary condition using Eq. (11), where the blood velocity (u) was equal to the wall velocity (u_tr). The wall velocity was calculated using the Keller–Miksis model,

The adjacent boundary to the bubble in the fluid domain representing the axis of symmetry in Fig. 1(c) was the axial symmetry boundary condition. A fully coupled fluid–solid interaction (FSI) boundary condition along with no slip boundary condition was applied at the blood and vessel interface using the following equations:

where v_fluid is the blood velocity at the boundary, u_solid is the vessel displacement at the boundary, n is the normal vector to the FSI boundary, σ is the stress tensor in the vessel, and Γ is the stress tensor in blood. A two-way coupling was applied at the vessel and blood interface, such that the pressure and viscous forces in blood shown by Eq. (14) were given as loads to the vessel using Eq. (13), which caused deformation (u_solid) of the vessel. In response, the deformation of the vessel caused changes in the blood domain using Eq. (11).

In the fluid domain (blood), a moving mesh based on the Arbitrary Lagrangian–Eulerian (ALE) formulation was used. Generally, the Lagrangian formulation is used for simulating solids, and the Eulerian formulation is used for simulating fluids. The Lagrangian formulation based on the material coordinate system allows moving boundaries, while the Eulerian formulation based on the spatial coordinate system does not allow moving boundaries. Therefore, a combination of both Lagrangian and Eulerian formulations known as ALE was used to account for the movement as a function of time at the vessel–blood interface and bubble–blood interface. The hyperelastic smoothing method, which uses minimum mesh deformation energy criteria, was used for mesh deformation. It is highly non-linear and gives robust results for highly stretched meshes also. The automatic remeshing method using distortion criteria of 2 was used to maintain the quality of mesh throughout the simulation. It works by remeshing the geometry each time for the increase in distortion beyond a limit of 2.

Triangular elements, which used piecewise linear interpolation for velocity and pressure fields in the fluid domain and quadratic Lagrange interpolation for the displacement field in the solid domain, were used for generating meshes. The mesh was fine near the bubble–blood interface with 32 elements on the interface and coarse away from the interface. The total mesh elements depended on the model size with ∼17 000 elements for a 50 μm vessel and 30 000 elements for a 150 μm vessel. A free time step with a maximum time step of 1 ns was used for solving the model, and a maximum time step of 0.1 ns was used when bubble acceleration was high during the compression phase after maximum expansion. The length of the vessel was ten times the radius of the vessel, i.e., for a 150 μm vessel, the vessel length was 1500 μm.

After solving the model, the circumferential and shear stresses on the vessel wall due to bubble oscillation were calculated. The circumferential stress was calculated on the vessel wall by assuming the vessel as a thick cylinder using

where P_i is the pressure on the inner radius of the vessel and r_i and r_o are the inner and outer radii. For a thin vessel with an effective vessel radius at least 20 times more than the vessel thickness, the equation reduces to

The shear stress due to blood velocity gradients on the vessel wall was calculated using

where v is the fluid velocity in the y direction and w is the fluid velocity in the z direction.

The model was first validated by comparing it with the Hosseinkhah model. We have reproduced the Hosseinkhah model in COMSOL Multiphysics 5.5 by following the original paper.³⁶ We were running this model using parameters from the Hosseinkhah paper, and a bubble of 2 μm inside a 5-μm vessel, which had a vessel thickness of 2 μm and Young’s modulus of 5 MPa, was used. The oscillation of the 2 μm bubble at an ultrasound amplitude of 261.5 kPa at 1 MHz is shown in Fig. 2(a), which is identical to the bubble oscillation in the original paper.³⁶ In our model with the same parameters as in the Hosseinkhah paper, a bubble of 2 μm was forced to follow the bubble oscillation, as shown in Fig. 2(a). The calculated vessel displacement, circumferential stress, and shear stress from the Hosseinkhah model and our model were compared. The vessel displacement, circumferential stress, and shear stress obtained in our model matched very well with those in the Hosseinkhah model, as shown in Figs. 2(b) and 3. The RMSD (root mean square deviation) between our model and the Hosseinkhah model for vessel displacement, circumferential stress, and shear stress was 0.001 μm, 1.29 kPa, and 60.6 Pa, which have about 2% deviation. Overall, this validation verified that our model was correct and accurate.

Figure 4 shows the waveform of a PA wave in a blood vessel and when it is combined with an ultrasound wave. The PA wave was calculated using Eqs. (1) and (2) ⁵⁰ at a location that is 25 nm from the axis of a 50 μm vessel [Fig. 4(a)]. A significantly high rarefaction pressure of 2.55 MPa was generated in the PA wave. When the PA wave was combined with the ultrasound wave during PUT, the negative peak of the PA wave was superposed on the ultrasound wave, which had a peak negative pressure of 1550 kPa, at a phase angle of 261° (725 ns) [Fig. 4(b)]. The default laser fluence of 20 mJ/cm² and ultrasound frequency of 1 MHz were used in all simulations throughout this study unless otherwise indicated. The laser fluence of 20 mJ/cm² is within the safety standard recommended by the American National Standards Institute⁵⁵ for skin exposure and has been used in our earlier experiments.^43,45,50

Figure 5(a) shows the bubble oscillation when a 50 nm bubble was driven by combined PA wave and ultrasound wave application that is shown in Fig. 4(b), while Fig. 5(b) shows the bubble oscillation when the same bubble was driven by the ultrasound wave only. Both results were calculated with the Keller–Miksis equation. The bubble expanded to a maximum radius of 108 times its initial radius with the combined PA wave and ultrasound application [Fig. 5(a)]. In contrast, for ultrasound-alone, the bubble only expanded to a maximum radius of 1.8 times its initial radius [Fig. 5(b)]. The ultrasound wave phase angle at which the negative peak of the PA wave was imposed plays an important role in bubble oscillation.⁵⁶ In this case, the maximum bubble radius was found when the negative peak of the PA wave was superposed on the ultrasound wave at a phase angle of 261° (at 725 ns). To best demonstrate the impact of combining laser and ultrasound, this phase angle will be used in most of the following studies, unless otherwise indicated.

In our FEM model, the 50 nm bubble inside the 50 μm vessel was forced to move following the oscillation pattern, as shown in Fig. 5(a) for combined laser and ultrasound application, i.e., PUT, and Fig. 5(b) for ultrasound-alone application. The bubble movement will generate pressure and a velocity gradient inside the blood, resulting in the circumferential and shear stresses on the vessel wall. Figure 6(a) shows the produced circumferential stress as a function of time. For ultrasound-alone, the maximum and minimum of circumferential stresses of 0.002 kPa and −0.005 kPa were obtained. During PUT, where the PA wave and ultrasound wave were combined, the maximum and minimum amplitudes of 192 kPa and −392 kPa were obtained. The positive and negative circumferential stresses in Fig. 6(a) correspond to vessel compression and expansion. The maximum magnitude of the negative circumferential stress is used for maximum circumferential stress (MCS) comparison in the current study because it is the stress responsible for causing vessel rupture. Figure 6(b) shows the produced shear stress as a function of time. For ultrasound-alone, the maximum shear stresses (MSSs) of −0.000 35 Pa and 0.000 39 Pa were obtained. During PUT, the maximum amplitudes of −61 Pa and 100 Pa were obtained. The positive and negative shear stresses in Fig. 6(b) are due to the different direction of the fluid velocity during bubble compression and expansion phases. The maximum magnitude among the positive and negative shear stresses is used for maximum shear stress (MSS) comparison in the following results. The maximum circumferential stress on the vessel wall was obtained at a point directly above the center of the bubble as in the other studies.^30,36,57 The value of shear stress was, however, 0 at a point directly above the bubble center.^30,36,57 The maximum shear stress in Fig. 6(b) was found at a point 30 µm directly above the bubble center on the vessel wall.

The stresses shown in Fig. 6 were calculated until 1242 ns. The Keller–Miksis equation used for solving the bubble dynamics did not produce accurate results during the late compression phase for strong inertial oscillations. The solutions from the Keller–Miksis equation for combined ultrasound and laser application [Fig. 5(b)] resulted in a very high velocity at the bubble–blood interface during the late compression stage (around 1244 ns). The velocity of the bubble–blood interface exceeded the speed of sound in blood for a few nanoseconds, leading to inaccurate results from the Keller–Miksis equation. Therefore, we limited our study to solve the model only up to a time step where the bubble–blood interface velocity was less than the speed of sound in blood. The other consideration was that during this collapsing phase, bubbles were compressed to a size that was much smaller than the size of the blood vessel. Therefore, we assumed the effect on the vessel wall would be small when the bubble was located at the center of a blood vessel. We agree that if a bubble is close to the vessel wall, the strong collapsing phase must be considered as it is when the strongest impact occurs. However, it is out of the scope of the current study as we do not consider the effect of micro-jets.

Figure 7 shows the Maximum Normalized Radius (MNR) of a 50 nm bubble as a function of the ultrasound amplitude for ultrasound-alone application and PUT in 50 μm, 100 μm, and 150 μm vessels. For ultrasound-alone application, the MNR was close to 1 until 1550 kPa and jumped to 80 at 1600 kPa ultrasound amplitude. For PUT application in a 50 μm vessel, the MNR was close to 1 up to 1350 kPa and increased to 65 at 1400 kPa ultrasound amplitude. The threshold for the jump in the MNR was reduced from 1600 kPa to 1400 kPa during PUT for the bubble inside the 50 μm vessel. The threshold was further reduced to 1150 kPa for the 100 μm vessel and 1000 kPa for the 150 μm vessel during PUT. The threshold for ultrasound-alone was 1600 kPa. The threshold decreased during PUT when the vessel radius was increased because of the increased PA wave amplitude.

Figure 8 shows the MCS and MSS due to the oscillation of a 50 nm bubble for ultrasound-alone application on the wall of a 50 μm vessel and PUT application on the wall of 50 μm, 100 μm, and 150 μm vessels. For ultrasound-alone in a 50 μm vessel, the MCS and MSS were almost 0 up to 1550 kPa ultrasound amplitude and jumped to 215 kPa and 43 Pa, respectively, at 1600 kPa ultrasound amplitude. The MCS and MSS were close to 0 below 1550 kPa ultrasound amplitude because the bubble oscillation amplitude was small and only generated weak pressure and velocity gradients in blood and did not produce significant circumferential stress and shear stress on the vessel wall. At 1600 kPa ultrasound amplitude, the bubble oscillation amplitude was large and resulted in a significant increase in the MCS and MSS on the vessel wall. For PUT application in the 50 μm vessel, the MCS and MSS started rapidly increasing at 1400 kPa ultrasound amplitude due to a large bubble oscillation amplitude. Similarly, for the 100 μm and 150 μm vessels, the MCS and MSS started rapidly increasing at 1150 kPa and 1000 kPa ultrasound amplitudes, respectively.

To better understand the results, we defined a threshold peak ultrasound pressure (TPUP) for the MCS and MSS. The TPUP for the MCS and MSS was when their gradient exceeds 100 and 0.001, respectively, for the first time while the ultrasound amplitude was increased from 950 kPa to 1600 kPa. The TPUP for the MCS and MSS for ultrasound-alone was 1600 kPa, while the TPUP for PUT was 1400 kPa, 1150 kPa, and 1000 kPa for vessels of 50 μm, 100 μm, and 150 μm, respectively. The TPUP was decreased with increased vessel radius because the PA wave amplitude inside a vessel was increased with increased vessel radius. However, when the vessel size was increased, the rate of increase in both the MCS and the MSS as a function of the ultrasound amplitude slowed down. The slope of the increase in both the MCS and the MSS as a function of the ultrasound amplitude was largest for the 50 μm vessel and smallest for the 150 μm vessel beyond the TPUP. While the slope of the increase in the MCS reduced slightly with increased vessel size, the slope of the increase in the MSS reduced more drastically with increased vessel size.

Figure 9 shows the MCS and MSS on the wall of a 50 μm vessel as functions of the ultrasound amplitude due to the oscillation of a 50 nm bubble for ultrasound-alone and PUT at 1 MHz and 0.7 MHz ultrasound frequencies. The negative peak of the PA wave was superposed on the ultrasound wave at a phase angle of 261°, which corresponded to 725 ns for 1 MHz and 1036 ns for 0.7 MHz ultrasound frequency, measured from the beginning of an ultrasound cycle. The MCS and MSS were larger for 0.7 MHz ultrasound frequency as compared with 1 MHz ultrasound frequency for both ultrasound-alone and PUT. The TPUPs of both the MCS and the MSS were 1400 kPa for PUT at both frequencies. For ultrasound-alone, the TPUPs of the MCS and MSS were 1600 kPa at 1 MHz frequency and 1550 kPa at 0.7 MHz frequency.

When the laser fluence changes, the amplitude of the resulting PA wave will change. As a result, the MCS and MSS resulting from an oscillating bubble will be affected. Figure 10 shows the MCS and MSS on the wall of a 50 μm vessel as functions of the laser fluence due to the oscillation of a 50 nm bubble for PUT. Both the MCS and the MSS decreased nearly linearly as the laser fluence decreased and reached a near zero value at 4 mJ/cm² for the 50 μm vessel at an ultrasound amplitude of 1550 kPa.

Figure 11 shows the MCS and MSS on the wall of a 100 μm vessel as functions of the bubble radius for ultrasound-alone and PUT. The ultrasound amplitude of 1200 kPa and laser fluence of 20 mJ/cm² were used for all simulations. The negative peak of the PA wave was superposed at different ultrasound phase angles for each size of bubbles in order to obtain the maximum amplitude for PUT. The MCS and MSS produced during PUT were significantly higher than those produced during ultrasound-alone when the size of a bubble was small. The difference can be observed from the gap between the two curves. However, with the increased bubble size, the gap between stresses produced during PUT and ultrasound-alone was reduced and the curves almost overlaid each other for large size bubbles. For example, the MCS was 1 kPa for ultrasound-alone and 100 kPa during PUT for a 50 nm bubble. For a 200 nm bubble, it became 500 kPa for ultrasound-alone and 590 kPa during PUT. Similarly, the MSS was 0.1 Pa during ultrasound-alone and 5 Pa for PUT for a 50 nm bubble, while it was 35 Pa during ultrasound-alone and 40 Pa for PUT for a 200 nm bubble. The effect of the laser diminished with increased bubble size, which could be due to the combination of the PA wave amplitude and bubble dynamics. Particularly, smaller bubbles tend to exhibit inertial cavitation behavior, which is very sensitive to changes in driving pressures. A small change in driving pressure, such as the superposition of the PA wave, could induce a large change in bubble oscillation and, thus, large changes in the MCS and MSS. Our results showed that the effect of PUT was only pronounced for smaller size bubbles with the radius less than 100 nm.

The magnitude of shear stress has been measured in normal blood vessels. Peak and mean shear stress values in blood vessels have been measured to be less than 5.5 Pa and 1.6 Pa for large arteries with diameters from 2.7 mm to 8.5 mm.⁵⁸ The mean shear stress measured in arterioles that are larger than 15 μm diameter in animals such as rabbits, cats, and rats was no more than 5 Pa.⁵⁹ The value of shear stress derived from theoretical calculation was 1.5 ± 50% Pa^58,59 and was assumed constant along the arterial tree. In another study, by combining experimental measurement and theoretical simulation, the shear stress was estimated to be slightly more than 10 Pa in capillaries that were smaller than 10 μm in diameter.⁶⁰ All these results indicated that under normal conditions, the shear stress in a blood vessel does not exceed 10 Pa in large vessels and may slightly exceed 10 Pa only in small capillaries. In our simulation, when the laser and ultrasound were combined, i.e., during PUT, the calculated shear stresses produced by an oscillating bubble can easily exceed 10 Pa, while the calculated shear stresses were nearly zero with ultrasound-alone. As increased shear stresses will promote the activation of endothelium cells in the blood vessel wall and result in various vessel biological responses, this result could partially explain the enhanced bio-effects of the reduced blood perfusion rate^45,46,61 and thrombolysis⁶² observed during PUT in our previous studies.

The circumferential stress is responsible for the rupture of the vessel wall and can result in oval opening in the vessel wall with its larger axis along the vessel axis.⁶³ A study with frog mesenteries estimated a circumferential stress of 800 kPa for vessel rupture,⁶³ while, in another study involving rabbits, capillary rupture was observed for a circumferential stress above 80 kPa.⁶⁴ In one study using abdominal aortic aneurysm wall tissue, the circumferential stress for vessel rupture ranged from 200 kPa to 2000 kPa depending upon the diameter, thickness, and other vessel properties.⁶⁵ In the current study, we found that the calculated MCSs are mostly less than the values needed for vessel rupture for the vessel sizes we simulated, while vessel rupture may still be possible depending on the selected laser and ultrasound parameters during PUT.

In all our simulations, when ultrasound was applied alone, the MCS and MSS on the vessel wall were nearly zero. Only during PUT, significant stresses can occur and may activate bio-effects that are triggered by circumferential and shear stresses. On the other side, since shear stress and circumferential stress play different roles in the vessel response, one needs to be very careful to select the proper laser and ultrasound parameters during PUT in order to achieve the desired bio-effects. Based on our simulation, it is possible to find a parameter range where the MSS is above the normal value, but the MCS is less than the value that may cause vessel rupture. Hence, PUT has the potential to induce bio-effects through the elevated shear stress without causing vessel ruptures.

We have examined various parameters that affect the bubble oscillations and resulting stresses on the vessel wall. The parameters such as the ultrasound amplitude, ultrasound frequency, laser fluence, superposition between PA and ultrasound waves, size of the blood vessel, and bubble size played an important role in PUT. Particularly, as the ultrasound amplitude increases, the MCS and MSS values increase, as shown in Figs. 8 and 9. At the TPUP, circumferential stress and shear stress increase significantly from the near zero value. The TPUP depends on the peak rarefactional pressure and the time period of the rarefactional phase. The TPUP decreases for PUT as compared to ultrasound-alone, as shown in Figs. 8 and 9, due to high peak rarefactional pressure resulting from the PA wave, which reduces the peak rarefactional pressure requirement from the ultrasound wave. The TPUP for ultrasound-alone reduces in Fig. 9 with the decrease in the ultrasound frequency due to the increase in the time period of the rarefactional phase. For a very high value of the ultrasound amplitude beyond the TPUP for ultrasound-alone, the stresses from ultrasound-alone and PUT will be similar. The stresses are near zero below the TPUP for PUT. For laser parameters, when the laser fluence is decreased, the produced stresses will decrease, as shown in Fig. 10, and fall to zero below a certain laser fluence. Previous experimental studies have shown vascular changes only occur when the applied ultrasound amplitude and laser fluence were above a certain level, and further increased ultrasound pressure and laser fluence could induce vessel rupture.^45,66 These results can be further validated with the simulation results in Figs. 8–10. Furthermore, bubble oscillations are also strongly affected by the ultrasound phase angle at which the negative peak of the PA wave is superposed. Figure 12 shows the MNR of an oscillating bubble as a function of the ultrasound phase angle at which the negative peak of the PA wave is superposed under different ultrasound pressures. The optimal phase angle, which is defined as the phase angle at which the maximum of the MNR occurs, shifts to the left as the ultrasound pressure increases. The effect of the phase angle has been demonstrated experimentally in our previous studies.⁵⁶ 

Small size cavitation nuclei are present in human blood in vivo. When an ultrasound wave is applied, due to rectified diffusion, these cavitation nuclei can grow. A previous study has shown that a 100 nm bubble can grow to 10 μm due to rectified diffusion during PUT.⁵⁰ As a bubble grows, the resulting stresses on the vessel wall can increase many folds for the same ultrasound amplitude, as shown in Fig. 11. For example, the MCS is 99 kPa for a 50 nm bubble and 1395 kPa for a 5 μm bubble during PUT with the same parameters. Similarly, the MSS is 5 Pa for a 50 nm bubble and is 90 Pa for a 5 μm bubble during PUT with the exact same parameters. However, for large size bubbles, the effect of the laser is negligible.

Our study is not without limitations. The main limitation is that we assumed that a bubble always remained at the center of a blood vessel. The assumption has been adopted in previous studies^43,50 and can be true when a bubble is at the center during the beginning phase of PUT. PUT utilizes concurrently applied laser pulse and ultrasound burst. In our previous study, a nanosecond laser pulse and a 10 ms ultrasound burst were combined during each PUT cycle. The mechanism contained two steps: the first step was the substantial increase in the equilibrium bubble size through the combined laser pulse and ultrasound wave at the center of a blood vessel during the beginning phase, and the second step was the oscillation of the formed bubble that was driven by the subsequent ultrasound cycles. The current study is limited to the simulation of the first step. During the second step, as a bubble may very likely drift off the vessel center and move closer to the vessel wall, the magnitude of the induced MCS and MSS can be significantly increased for ultrasound-alone during the subsequent ultrasound cycles in PUT. When a bubble is off the center in a blood vessel, the model is not symmetric and the axis of symmetry cannot be defined, which restricts the use of the 2D axisymmetric approach to solve the model. In this case, a 3D numerical model must be employed, which is particularly a computation challenge for inertial cavitation oscillations. The second limitation is that our simulation was only performed when the bubble wall velocity was less than the velocity of sound in the fluid due to limitations from the Keller–Miksis equation and COMSOL model. During the collapsing phase, a bubble becomes very small and far from the vessel wall, leading to negligible force on the vessel when the bubble is placed at the center of a large blood vessel. Therefore, we expect that the impact of this limitation on our conclusion should be insignificant. We agree that when a bubble is closer to the vessel wall or a vessel size is small, this effect must be considered. The use of an elastic model for the vessel and tissue in our study is also a limitation. Although most of the past studies have used an elastic model, a viscoelastic model could be better if the most appropriate model can be identified.

In this study, a 2D axisymmetric FEM model was developed to simulate the interaction between blood vessels and oscillating bubbles during the beginning phase of PUT, and the resulting circumferential and shear stresses on the vessel wall were calculated. Both the circumferential and shear stresses increased when the ultrasound pressure and laser fluence were increased but decreased when the ultrasound frequency and vessel size were increased. In addition, the produced stresses were found to be affected by the ultrasound phase angle at which the negative peak of the PA wave was superposed on the ultrasound wave. The effect of the laser during PUT was found to be significant only for bubbles smaller than 100 nm in radius. The findings here can be used to select appropriate parameters for maximizing efficiency in experiments using PUT. In the future, a 3D model can be developed to investigate the stresses on the vessel wall for off-center bubbles during PUT.

This work was supported in part by the National Institutes of Health (NIH) under Grant No. R01EY029489.

The data that support the findings of this study are available from the corresponding author upon reasonable request.