摘要:声空化机械效应是聚焦超声治疗的重要物理机制. 以脂类包膜微泡/纳米相变液滴为空化核可显著地增强空化效应, 本文耦合空化动力学、组织和脂类包膜黏弹性模型, 构建了组织内脂类包膜微泡声空化动力学模型, 数值分析了微泡声空化动力学行为以及周围组织内机械应力的时空分布规律, 并探究了包膜材料、组织黏弹性和驱动声压等关键参数的影响. 包膜和组织黏弹性都将抑制微泡振动, 但组织黏弹性的抑制作用更大. 组织内机械应力在膨胀阶段为挤压应力, 而在收缩阶段和反弹初始阶段为拉伸应力, 且应力局部分布于微泡壁附近, 随着距离增大而显著减小, 其中拉伸应力衰减率明显更大. 包膜黏弹性可减小应力, 但声压较大时, 应力减小可忽略不计. 应力随着组织弹性增大而减小, 随着组织黏度增大而先增大后减小, 随着声压增大而增大. 本研究可为进一步阐释聚焦超声治疗中组织机械损伤的内在机制奠定重要理论基础.
关键词: 气泡动力学/
包膜微泡/
组织内声空化/
空化机械效应

English Abstract

全文HTML

--> --> -->

超声空化是指在超过一定声压阈值的超声波作用下, 液体中压力变化引发的气泡产生、生长、振荡和溃灭等一系列动力学过程, 它是强超声在液体中引发的一种特殊的瞬态物理现象^[1-3]. 超声空化可导致血管内皮细胞损伤、毛细血管破裂、血脑屏障开放、细胞声致穿孔和组织毁损等力生物学效应^[4-10], 促进药物/基因细胞内递送, 甚至直接机械破坏肿瘤细胞, 对肿瘤乏氧区内细胞杀伤具有重要意义^[6]. 声空化在生物组织或细胞周围产生的微流、应力、冲击波和微射流等机械效应在上述过程中发挥重要作用, 其中声空化产生的组织应力与应变已被认为是重要机制之一^[11-14]. 由于空化气泡的微米尺度、高速振动和组织的不透明性, 实验研究组织内的声空化动力学行为以及量化组织内产生的应力(如时间演化与空间分布等)变得十分困难^[11,14]. 基于建模仿真的理论研究可以规避上述实验限制, 阐释组织内声空化动力学瞬态物理过程和量化组织内应力的时空分布规律, 进一步深入理解声空化造成细胞/组织机械损伤的机制, 有助于规划和监控超声治疗过程以实现精准高效的超声肿瘤治疗.
生物软组织是一种复杂的黏弹性介质, 组织内声空化动力学行为因需要考虑组织黏弹特性而变得更复杂. 目前一般可通过耦合牛顿流体水中的空化动力学方程和组织黏弹性模型来实现^[15]. 例如, Yang和Church^[16]将空化动力学Keller-Miksis (KM)方程与组织黏弹特性的Kelvin-Voigt模型结合, 得到了软组织内声空化动力学模型, 目前该模型已被广泛使用. Hua和Johnsen^[17]通过耦合KM方程与组织黏弹性Zener模型, 研究了组织弹性和应力松弛时间对空化泡大幅度振动的影响, 深化了黏弹性组织内声空化动力学的研究. 近年来, 为了阐释高强度聚焦超声进行组织毁损治疗的内在机制, Mancia等^[11,13]通过耦合KM方程和非线性弹性的Kelvin-Voigt模型构建了组织内声空化动力学模型, 分析了组织内声空化产生的机械应力分布. 然而在上述理论模型中, KM方程只适用于微泡振动速度$ \dot{R} $较小的情况(即马赫数$ \dot{R}/c < 1 $)^[18]. 在实际的超声治疗中, 较高的超声声压通常会使得空化泡在坍塌阶段时收缩速度的马赫数$ \dot{R}/c>1 $; 另外, 常用的Kelvin-Voigt模型也未考虑组织的应力松弛过程, 而在极短的超声辐照时间内(微秒级), 组织应力松弛过程会变得十分重要^[18]. 对此, Zilonova等^[18]和Pahk等^[19]提出选择适用于更高马赫数($ \dot{R}/c < 2.2 $)的空化动力学Gilmore方程与考虑组织弹性、黏性和应力松弛的Zener模型耦合, 得到适用范围更广且更加符合实际情况的组织内声空化动力学模型.
组织内声空化阈值通常较高, 为了降低空化阈值以增加空化治疗时的安全性和可控性, 需要引入人工空化核, 如造影微泡、回声脂质体、纳米相变液滴和介孔固体颗粒等微纳粒子^[20,21]. 其中造影微泡和纳米相变液滴通常由脂类、白蛋白或聚合物等包膜封装以增加其稳定性、声散射特性和分子靶向性等^[21]. 包膜是影响包膜微泡声空化动力学的重要因素之一, 包膜微泡比无包膜微泡更硬、更黏, 将限制微泡振动幅度, 降低空化非线性效应^[22-27]. 引入微泡的包膜参数对声空化动力学方程进行修正可开发出多种包膜微泡的声空化动力学模型^[25,26]. 脂类包膜微泡因包膜厚度更薄、弹性更大, 近年来得到了更广泛的关注和应用. 实验研究发现, 脂类包膜微泡存在“仅压缩”和“剪切变稀”现象^[28-30]. 对此, Marmottant等^[29]将微泡包膜表面张力分为皱缩坍塌、弹性和破裂三个阶段, 有效模拟了脂类包膜微泡的“仅压缩”行为. Doinikov等^[31]考虑包膜黏性与包膜剪切速率之间的函数相关性解释了包膜“剪切变稀”的现象. 国内章东研究团队^[32-34]认为包膜弹性和黏性都是变量, 建立了一个既可描述“仅压缩”现象, 又可表示“剪切变稀”行为的动力学模型. 通过拟合实验数据, 发现构建的新模型不仅可以很好地模拟脂类包膜微泡的“仅压缩”行为, 同时还可减小包膜弹性和黏度对微泡初始半径的依赖性, 使结果更符合包膜材料的物理本质. 然而, 上述研究均是在牛顿流体水中进行的, 在脂类包膜微泡/纳米相变液滴增效超声治疗的实际应用中, 组织内包膜微泡空化动力学行为以及空化造成周围组织机械损伤的机制尚需进一步研究.
本文耦合Gilmore方程、组织黏弹性Zener模型和脂类包膜的非线性黏弹性模型, 构建了组织内脂类包膜微泡的声空化动力学模型, 定量研究了微泡的声空化动力学行为以及微泡周围组织内机械应力的时空分布规律, 最后探讨了微泡包膜、组织弹性与黏性以及超声声压的影响.

2

-->

2.1.空化动力学方程

考虑单个微泡在黏弹性生物组织内的声空化动力学行为以及组织内应力的时空分布$ {\tau }_{rr}(r, t) $, 如图1所示. 声空化过程中, 假设微泡仍然保持球形、微泡内部压力和温度均匀分布、微泡内气体满足理想气体定理以及微泡与周围介质间无物质交换. 本文选用适用性更广的Gilmore方程描述声空化动力学^[18,19]:
图 1 黏弹性组织内微泡动力学与应力示意图
Figure1. Schematic of the bubble dynamics and stress in the viscoelastic tissue.

$\begin{split} & R\ddot R\left(1 - \frac{{\dot R}}{c}\right) + \frac{3}{2}{{\dot R}^2}\left(1 - \frac{{\dot R}}{{3c}}\right) \\=\;& H\left(1 + \frac{{\dot R}}{c}\right) + \frac{{R\dot H}}{c}\left(1 - \frac{{\dot R}}{c}\right) + \left(1 + \frac{{\dot R}}{c}\right) \\ &\times \left[\frac{1}{\rho }\int_R^\infty \left(\frac{{\partial {\tau _{rr}}}}{{\partial r}} \right.\right. +\left.\left.\frac{{3{\tau _{rr}}}}{r}\right) {\rm{d}}r\right] \\ & + \frac{R}{c}\frac{{\rm{d}}}{{{\rm{d}}t}}\left[\frac{1}{\rho }\int_R^\infty {\left(\frac{{\partial {\tau _{rr}}}}{{\partial r}} + \frac{{3{\tau _{rr}}}}{r}\right)} {\rm{d}}r\right] ,\end{split} $

$H = \frac{1}{\rho }\frac{n}{{n - 1}}\left({p_0} + {p_{\rm{a}}} \!+\! B\right)\left[\!{\left(\frac{{{p_{\rm{b}}} + B}}{{{p_0} + {p_{\rm{a}}} + B}}\right)^{\tfrac{{n - 1}}{n}}} \!\!- \!1\right], $

$c = c{}_0{\left(\frac{{{p_{\rm{b}}} + B}}{{{p_0} + {p_{\rm{a}}} + B}}\right)^{\tfrac{{n - 1}}{{2n}}}},$

式中, R为空化泡半径; $ \dot{R} $, $ \ddot{R} $分别为R对时间的一阶导数和二阶导数; ρ为介质密度; H为介质的焓; $ \dot{H} $为H对时间的一阶导数; n和B为Gilmore方程常数; p₀为静态压力; p_a为超声声压; p_b为介质内微泡壁处压力; c为微泡壁处局部声速; c₀为介质中声速; τ_rr为介质内r方向的应力.
2 -->

2.2.空化动力学方程与组织黏弹性模型耦合

黏弹性是生物组织一个非常重要的物理特征参数, 可显著影响组织内空化动力学行为和空化力生物学效应. 目前已构建多种黏弹性本构方程来阐释生物组织应力-应变-时间关系^[35]. 选用可同时描述组织弹性、黏性和应力松弛特征的Zener模型, 其本构方程为^[35]

${\tau _{rr}} + {\lambda _1}{\dot \tau _{rr}} = 2G{\gamma _{rr}} + 2\mu {\dot \gamma _{rr}},$

式中${\dot \tau _{rr}}$为${\tau _{rr}}$对时间的一阶导数, γ_rr为周围介质的应变, $ {\dot{\gamma }}_{rr} $为γ_rr对时间的一阶导数, λ₁为介质的松弛时间, G为介质的弹性模量, μ为介质的黏度. 组织内微泡振动过程中, 求解介质连续方程可得组织应变率${\dot{\gamma }}_{rr}=-2\dfrac{{R}^{2}\dot{R}}{{r}^{3}}$, 应变${\gamma }_{rr}=-\dfrac{2}{3{r}^{3}}({R}^{3}- {R}_{0}^{3})$, 代入(4)式可得微泡界面处的应力${\tau _{rr}}\left| {_R} \right.$:

${\tau _{rr}}\left| {_R} \right. + {\lambda _1}{\dot \tau _{rr}}\left| {_R} \right. = - \frac{{4G}}{3}\left(1 - \frac{{R_0^3}}{{{R^3}}}\right) - 4\mu \frac{{\dot R}}{R}, $

式中, τ_rr|_R为R处的应力, R₀为微泡初始半径. 将(5)式除以r, 且从R到∞进行积分, 可得

$\begin{split}&\int_R^\infty {\frac{{{\tau _{rr}} + {\lambda _1}{{\dot \tau }_{rr}}}}{r}} {\rm{d}}r = \int_R^\infty {\frac{{2G{\gamma _{rr}} + 2\mu {{\dot \gamma }_{rr}}}}{r}} {\rm{d}}r \\=\;& \frac{1}{3}\left[ - \frac{{4G}}{3}\left(1 - \frac{{R_0^3}}{{{R^3}}}\right) - 4\mu \frac{{\dot R}}{R}\right].\end{split}$

引入变量$q = \displaystyle\int_R^\infty {\frac{{{\tau _{rr}}}}{r}} {\rm{d}}r$, 则(6)式可简化为^[18,19]

$q + {\lambda _1}\dot q + {\lambda _1}\frac{{\dot R{\tau _{rr}}\left| {_R} \right.}}{R} = \frac{1}{3}\left[ - \frac{{4G}}{3}\left(1 - \frac{{R_0^3}}{{{R^3}}}\right) - 4\mu \frac{{\dot R}}{R}\right], $

式中, $ \dot{q} $为q对时间的一阶导数, 由$\displaystyle\int_R^\infty \left(\frac{{\partial {\tau _{rr}}}}{{\partial r}}\right. + $$ \left.\dfrac{{3{\tau _{rr}}}}{r}\right) {\rm{d}}r = - {\tau _{rr}}\left| {_R} \right. + 3 q$, 进一步化简(1)式可得

$\begin{split}&R\ddot R\bigg(1 - \frac{{\dot R}}{c}\bigg) + \frac{3}{2}{\dot R^2}\bigg(1 - \frac{{\dot R}}{{3c}}\bigg) = \bigg(1 + \frac{{\dot R}}{c}\bigg)\\&\qquad\times\bigg(H - \frac{{{\tau _{rr}}| {_R - 3q} }}{\rho }\bigg)\\&\qquad + \frac{R}{c}\left[\dot H\left(1 - \frac{{\dot R}}{c}\right) - \frac{{{{\dot \tau }_{rr}}\left| {_R - 3\dot q} \right.}}{\rho }\right].\end{split}$

联立(2), (3), (5), (7)和(8)式可得描述黏弹性组织内空化动力学的Gilmore-Zener耦合模型.
2 -->

2.3.包膜的非线性黏弹性模型

包膜材料可增加微泡稳定性, 延长微泡在体内的循环时间, 但包膜黏弹特性会影响组织内包膜微泡动力学行为^[36,37]. 同时考虑脂类包膜微泡“仅压缩”和“剪切变稀”的特殊行为, 使用非线性修正的Kelvin-Voigt黏弹性模型描述脂类包膜材料的黏弹特性^[29-32]:

$S = {S_{\rm{e}}} + {S_{\rm{v}}} = \frac{{2\sigma (R)}}{R} + 4{\kappa _{\rm{s}}}\frac{{\dot R}}{{{R^2}}}, $

式中, 包膜黏弹性S包括弹性项S_e和黏性项S_v, 其中弹性项由包膜表面张力σ(R)随半径R变化的分段表达式描述^[29]:

$\sigma (R) =\! \left\{ \begin{aligned} &0 ,\qquad\qquad\qquad\qquad R \leqslant {R_{{\rm{buckling}}}}{\rm{ }},\\ &\chi \Bigg(\dfrac{{R^2}}{R_{{\rm{buckling}}}^2}- 1\Bigg) ,\;\;{R_{{\rm{buckling}}}} < R < {R_{{\rm{rupture}}}}{\rm{ }},\\ &{\sigma _{{\rm{tissue}}}},\qquad\qquad\qquad\;\; R \geqslant {R_{{\rm{rupture}}}}, \end{aligned}\right. $

式中, R_buckling为微泡半径的最小值, 微泡半径R低于该值, 微泡包膜将进行皱缩, 高于${R_{{\rm{rupture}}}} = $$ {R_{{\rm{buckling}}}}{\left(1\! +\! {{{\sigma _{{\rm{tissue}}}}}}/{\chi }\right)^{{1 / 2}}}$微泡将发生破裂, 在R_buckling和R_rupture两者之间时, 微泡包膜弹性存在线性关系; $ \chi $为包膜弹性模量; $ {\sigma }_{\mathrm{t}\mathrm{i}\mathrm{s}\mathrm{s}\mathrm{u}\mathrm{e}} $为组织的表面张力.
考虑脂类材料的流变特性, (9)式黏性项中包膜黏度$ {\kappa }_{\mathrm{s}} $被认为是包膜剪切速率($ \dot{R}/R $)的函数, 并对$ {\kappa }_{\mathrm{s}} $进行非线性修正^[31,32]:

${\kappa _{\rm{s}}} = {{{\kappa _0}}}/({{1 + \alpha {{{\rm{|}}\dot R{\rm{|}}} / R}}}), $

式中, $ {\kappa }_{0} $为包膜初始黏度, α为特征时间常数. 联合(9)式、(10)式和(11)式可得脂类微泡中包膜材料的非线性黏弹性模型.
2 -->

2.4.模型的初始条件、边界条件与数值求解

如图1所示, 处于初始平衡状态的微泡R = R₀, $ \dot{R}=0 $, 假设在振动过程中微泡内部气体满足理想气体定律, 故微泡内部压力p_i表示为

${p_i}(R) = \left({p_0} + \frac{{2{\sigma _0}}}{{{R_0}}}\right){\left(\frac{R}{{{R_0}}}\right)^{3\gamma }}, $

式中, σ₀为微泡的初始表面张力, γ为气体的多方指数.
对于微泡-组织界面的压力边界条件, 综合考虑了微泡包膜表面张力σ(R)、黏度$ {\kappa }_{\mathrm{s}} $和周围介质应力$ {\tau }_{rr} $(如图1所示), 可得在界面处组织内压力p_b为

${p_{\rm{b}}}(R) = {p_i}(R) - \frac{{2\sigma (R)}}{R} - 4{\kappa _{\rm{s}}}\frac{{\dot R}}{{{R^2}}} + {\tau _{rr}}\left| {_R} \right..$

联立(2)式、(3)式、(5)式、(7)—(13)式可得黏弹性组织内脂类包膜微泡声空化动力学模型. (5)式、(7)式和(8)式构成的非线性常微分方程组得不到解析解, 本文使用龙格-库塔算法进行迭代求解. 以下数值分析中, 超声声压为${p_{\rm{a}}} = - {p_{\rm{A}}}\sin (2\pi ft)$, 式中P_A为超声声压幅值, f为超声频率. 若无特殊说明时, 各物理参数取值如下^[18,32]: ρ = 1060 kg/m³, c₀ = 1540 m/s, p₀ = 1.013 × 10⁵ Pa, p_A = 1 MPa, f = 1 MHz, n = 7, $ {B=c}_{0}^{2}\rho /n-{p}_{0} $, G = 200 kPa, μ = 15 mPa·s, λ₁ = 3 × 10^–9 s, σ_tissue = 0.056 N/m, R₀ = 1 μm, χ = 0.4 N/m, κ₀ = 5 × 10^–9 kg/s, α = 0.5 μs, γ = 1.095.

2

-->

3.1.组织内脂类包膜微泡的声空化动力学行为

同时考虑微泡周围生物组织和微泡包膜的黏弹特性, 构建的理论模型(模型I)可描述临床常用脂类包膜微泡在生物组织内的声空化动力学行为. 该模型也可通过参数设置简化为生物组织内无包膜微泡(模型II)、牛顿流体水中脂类包膜微泡和无包膜微泡的动力学模型(模型III和IV), 表明该模型具有更广泛的适用性. 图2对比研究了脂类包膜微泡、无包膜微泡分别在生物组织和牛顿流体水中的声空化动力学行为, 阐释了组织和包膜黏弹性对微泡声空化动力学行为的影响, 并验证模型的适用性. 图2(a)和图2(b)分别给出了超声声压p_A = 0.1 MPa时, 微泡半径和振动速度随时间的变化曲线. 从图2(a)和图2(b)可以看出, 微泡在水中的振动幅度和振动速度都明显大于微泡在生物组织内的振动幅度和振动速度, 这表明生物组织黏弹性限制了微泡的声空化振动过程. 对于包膜材料黏弹性的影响, 对比发现包膜微泡振动呈现出明显的“仅压缩”行为, 且包膜微泡的振动幅度、振动速度都明显小于无包膜微泡, 这也说明包膜材料的黏弹性会减弱微泡的声空化振动过程. 图2(c)和图2(d)分别给出了p_A = 1 MPa时, 微泡半径和振动速度随时间的变化曲线. 同样地, 组织和包膜的黏弹性都会减弱微泡声空化的振动幅度和振动速度.
图 2 不同条件下微泡声空化动力学行为　(a), (b)和(c), (d)分别表示p_A为0.1 和1 MPa时, 微泡半径和振动速度随时间的变化曲线. I—IV分别表示组织内(I)脂类和(II)无包膜微泡、水中(III)脂类和(IV)无包膜微泡声空化动力学模型
Figure2. Bubble dynamics under different conditions. Panels (a), (b) and (c), (d) show the bubble radius and bubble wall velocity at p_A = 0.1 MPa and p_A = 1 MPa, respectively. The model I and II simulate a shelled microbubble and a clean microbubble oscillating in soft tissue, while the model III and IV in liquid water.

对比不同声压下微泡的动力学行为, 发现声压增大可显著地增大微泡振动幅度, 微泡缓慢膨胀而剧烈收缩, 微泡振动呈现出膨胀程度大于收缩程度的非线性特性. 在四种情况下, 微泡坍塌时最大速度的马赫数$ \dot{R}/c $都接近或超过1, 水中无包膜微泡坍塌速度马赫数最大($ \dot{R}/c=2.08 $), 如图2(d)所示. 相比于其他常用的空化动力学方程如Rayleigh-Plesset方程(适用范围为$ \dot{R}/c\ll 1 $^[18])和KM方程($\dot{R}/c \!<\! 1$^[18]), 本文中所用的Gilmore方程($\dot{R}/c < $$ 2.2$^[18])能够对较高声压下的微泡动力学行为(尤其是微泡急速的坍塌过程)进行更加准确地描述.
2 -->

3.2.组织内微泡声空化产生的应力分析

图3给出了组织内微泡声空化过程中, 微泡周围组织内应力随时间t和距离r的变化. 从图3(a)可清楚地看出, 组织内微泡在超声波作用下发生周期性膨胀与收缩的振动过程; 在微泡膨胀阶段, 微泡对周围组织造成挤压, 产生较大的挤压应力(负应力), 而在微泡收缩阶段, 微泡膨胀导致的形变组织会恢复到初始状态, 组织内产生较大的拉伸应力(正应力). 假设某一粒子(比如细胞)与微泡中心的初始间距为d₀, 如图1所示, 根据质量守恒原则, 在微泡声空化过程中该粒子与微泡中心间距d可以表示为$d\left(t\right)=\sqrt[\root{4}{3}]{{d}_{0}^{3}+{R\left(t\right)}^{3}-{R}_{0}^{3}}$^[11]. 图3(a)给出了在微泡声空化过程中3个粒子(d₀ = 2, 6 和11 μm)与微泡中心间距d随时间的变化曲线, 在微泡膨胀阶段d逐渐增大, 在收缩阶段d逐渐减小. 图3(b)给出了3个粒子所受应力随时间的变化曲线. 组织内距离微泡中心较近(2 μm)的粒子所受应力最大, 挤压应力最大值为0.42 MPa, 拉伸应力最大值可达1.52 MPa. 随着d₀增加, 组织内粒子所受到的应力急剧减小, 当d₀ = 11 μm时, 挤压应力最大值减小为0.03 MPa, 拉伸应力最大值减小为0.015 MPa.
图 3 组织内微泡声空化过程中, 微泡周围组织内机械应力的时间演化和空间分布　(a)微泡声空化动力学与组织内应力的时空分布, 实线、虚线和点线分别为d₀ = 2, 6和11 μm粒子的运动轨迹; (b) 3个粒子所受应力随时间的变化曲线; (c)不同时刻(图(a)中点线所示)组织内应力的空间分布, 蓝色和红色实线分别表示挤压应力和拉伸应力
Figure3. Temporal and spatial evolution of the mechanical stress distribution during acoustic cavitation in the tissue: (a) The bubble dynamics and spatiotemporal stress distribution, in which solid line, dashed line and dotted line correspond to the trajectories starting at 2, 6 and 11 μm from the bubble center, respectively; (b) the stress magnitudes as a function of time along each trajectory; (c) the stress distribution at the fixed times as indicated by vertical dotted lines in panel (a), with traces indicating the magnitudes of compressive stresses (solid blue) and tensile stresses (solid red).

图3(c)给出了不同时刻, 组织内应力的空间分布以及应力随距离r的变化关系. 从图3(c)可以看出, 较大的组织应力局部分布于微泡壁附近, 且随着r增大而急剧地减小, 但是不同时刻组织内应力衰减率不同, 其中最大拉伸应力(t = 0.726 μs)随r增大的衰减率最大. 实验研究表明组织内最终应力大于0.05 MPa时可造成组织机械损伤^[38], 若以0.05 MPa为组织损伤的应力阈值, 计算发现微泡声空化可对距离微泡中心小于6.3 μm范围内的组织造成机械损伤, 因此本文所构建的理论模型可用于预测与评价脂类包膜微泡声空化对周围组织的机械损伤.
按照组织弹性和黏性参数, 组织内应力可分为弹性和黏性分量^[11,13], 图4(a)给出了组织内微泡声空化过程中弹性应力$ {\tau }_{rr}^{\mathrm{E}} $的时间演化与空间分布. 从图4(a)可清晰地看出, 在微泡膨胀阶段弹性应力$ {\tau }_{rr}^{\mathrm{E}} $为挤压应力, 且随着微泡膨胀而逐渐增大, 当微泡膨胀至最大半径时$ {\tau }_{rr}^{\mathrm{E}} $到达最大值; 随着微泡收缩坍塌, $ {\tau }_{rr}^{\mathrm{E}} $迅速减小. 在气泡反弹的初始阶段弹性应力$ {\tau }_{rr}^{\mathrm{E}} $变为拉伸应力, 这是因为在此阶段微泡半径小于初始半径, 导致微泡周围组织处于拉伸状态. 图4(b)给出了声空化过程中组织应变γ_rr的时间演化与空间分布. 对比图4(a)和图4(b), 发现微泡声空化过程中组织内弹性应力$ {\tau }_{rr}^{\mathrm{E}} $与组织应变γ_rr具有相似的时间演化与空间分布规律, 结合(6)式分析表明组织内弹性应力$ {\tau }_{rr}^{\mathrm{E}} $主要由微泡声空化引起的组织应变γ_rr决定.
图 4 组织内微泡声空化过程中(a)弹性应力$ {\tau }_{rr}^{\mathrm{E}} $和(b)组织应变γ_rr的时间演化与空间分布, 以及(c)黏性应力$ {\tau }_{rr}^{\mathrm{v}} $和(d)组织应变率$ {\dot{\gamma }}_{rr} $的时间演化与空间分布
Figure4. Temporal and spatial evolution of the (a) elastic stress $ {\tau }_{rr}^{\mathrm{E}} $ and (b) strain γ_rr in the tissue; the temporal and spatial evolution of (c) the viscous stress $ {\tau }_{rr}^{\mathrm{v}} $ and (d) strain rate $ {\dot{\gamma }}_{rr} $, respectively.

图4(c)给出了组织内黏性应力$ {\tau }_{rr}^{\mathrm{v}} $的时间演化与空间分布. 在微泡膨胀阶段$ {\tau }_{rr}^{\mathrm{v}} $表现为挤压应力, 且随着微泡膨胀而逐渐减小, 在微泡膨胀至最大半径时$ {\tau }_{rr}^{\mathrm{v}} $达到最小值; 在微泡收缩阶段表现为拉伸应力, $ {\tau }_{rr}^{\mathrm{v}} $随着微泡收缩而逐渐增大, 当微泡收缩至最小值时$ {\tau }_{rr}^{\mathrm{v}} $达到最大值. 图4(d)给出了声空化过程中组织应变率$ {\dot{\gamma }}_{rr} $的时间演化与空间分布. 黏性应力$ {\tau }_{rr}^{\mathrm{v}} $与组织应变率$ {\dot{\gamma }}_{rr} $具有相似的时空分布规律, 这表明微泡声空化过程中组织内黏性应力$ {\tau }_{rr}^{\mathrm{v}} $与组织应变γ_rr无关, 而主要是由声空化引起的组织应变率$ {\dot{\gamma }}_{rr} $决定. 结合图3进一步分析可知, 组织内总应力由组织应变γ_rr和组织应变率$ {\dot{\gamma }}_{rr} $共同决定. 在微泡膨胀阶段和反弹初始阶段, 组织应变γ_rr起主要作用, 而在微泡收缩阶段, 组织应变率$ {\dot{\gamma }}_{rr} $起主要作用.
2 -->

3.3.关键参数对组织内声空化动力学和应力的影响

3 -->

3.3.1.微泡包膜材料的影响

图5给出了组织内无包膜和脂类包膜微泡在声压p_A = 0.1 MPa和p_A = 1 MPa下, 微泡的声空化动力学行为以及微泡周围组织内应力的时空分布. 在较低声压条件下, 对比图5(a)无包膜和图5(b)脂类包膜微泡的声空化动力学, 发现脂类包膜微泡的振动幅度更小且呈现出明显的“仅压缩”行为. 图5(c)给出了组织内最大挤压应力和最大拉伸应力随距离d₀的变化曲线, 可以明显看出组织内应力都随着距离增大而逐渐减小. 对比相同距离处的应力可发现包膜微泡声空化在周围组织内产生的应力明显小于无包膜微泡声空化产生的应力, 其中微泡膨胀阶段产生的挤压应力减弱程度更大, 这表明脂类微泡“仅压缩”行为限制了微泡膨胀, 减小了组织应变γ_rr, 进而减小了微泡声空化在周围组织产生的挤压应力.
图 5 微泡包膜对声空化动力学和组织内应力的影响, 包括p_A = 0.1 MPa时, 组织内(a)无包膜微泡和(b)脂类包膜微泡的声空化动力学以及微泡周围应力的时空分布; p_A = 1 MPa时, 组织内(d)无包膜微泡和(e)脂类包膜微泡的声空化动力学以及微泡周围应力的时空分布; (c)和(f)表示对应的最大挤压应力和最大拉伸应力随d₀的变化曲线
Figure5. Effects of bubble shell on the bubble dynamics and stress distributions in the tissues, including (a) the clean bubble and (b) shelled bubble at p_A = 0.1 MPa, (d) the clean bubble and (e) shelled bubble at p_A = 1 MPa; panels (c) and (f) show the corresponding maximum compressive stresses and tensile stresses, respectively.

图5(d)—(f)分别表示在p_A = 1 MPa时, 组织内(d)无包膜微泡、(e)脂类包膜微泡的声空化动力学与应力时空分布以及(f)最大应力随距离d₀的变化曲线. 相比于脂类包膜微泡, 无包膜微泡的声空化振动幅度、组织内最大挤压应力和最大拉伸应力都更大. 与低声压条件相比, 发现声压较高时微泡振动幅度和组织应力都更大, 其中微泡周围组织内最大拉伸应力增大约1000倍. 图5表明微泡包膜材料可减弱微泡声空化振动以及微泡空化在周围组织内产生的应力, 在声压较小时包膜材料的影响程度更大, 而在声压较大时包膜材料的影响较小, 几乎可以忽略. 这意味着在利用高强度超声进行超声治疗时可引入稳定性更高的包膜微泡/纳米相变液滴等人工空化核增强空化效应, 从而显著提高超声空化治疗效率.
3 -->

3.3.2.组织弹性模量和黏度的影响

图6(a)给出了不同弹性模量(G = 10, 100, 200, 500和1000 kPa)时, 组织内脂类包膜微泡半径随时间的变化曲线. 从图6(a)可以明显看出, 随着弹性模量G增大, 微泡的振动幅度和振动周期明显减小, 这说明弹性模量较大的组织(如病变的肝脏)将抑制微泡声空化振动过程. 图6(b)给出了不同d₀处(2, 6和11 μm)粒子所受最大挤压应力和最大拉伸应力随组织弹性模量G的变化曲线. 从图6(b)可以看出: 组织内最大挤压应力和最大拉伸应力都随着d₀增大而明显减小; 随着组织弹性模量G增大, 最大挤压应力和最大拉伸应力都逐渐减小, 这是因为组织的弹性模量增大限制了微泡声空化振动过程, 进而减弱了声空化对周围组织的机械力作用. 图6(c)给出了不同组织黏度(μ = 1, 5, 15, 50和100 mPa·s)时, 组织内脂类包膜微泡半径随时间的变化曲线. 随着组织黏度μ增大, 微泡的振动幅度和振动周期明显减小. 介质黏度增大使得微泡振动时需要克服介质中越来越大的分子间引力, 因此在黏度较大的介质中微泡比较难振动. 图6(d)给出了组织内不同d₀处(2, 6和11 μm)粒子所受最大挤压应力和最大拉伸应力随组织黏度μ的变化曲线. 由图6(d)可知: 随着d₀增大, 组织内最大挤压应力和最大拉伸应力都逐渐减小; 随着组织黏度μ增大, 最大挤压应力和最大拉伸应力先增大后减小, 在μ = 15 mPa·s时达到最大值.
图 6 (a)不同组织弹性模量(G = 10, 100, 200, 500和1000 kPa)下, 组织内微泡半径随时间的变化关系, (b)组织内不同d₀处(2, 6和11 μm)粒子所受最大挤压应力和最大拉伸应力随组织弹性模量G的变化关系; (c)不同组织黏度(μ = 1, 5, 15, 50和100 mPa·s)下, 组织内微泡半径随时间的变化关系, (d)组织内不同d₀处(2, 6和11 μm)粒子所受最大挤压应力和最大拉伸应力随组织黏度μ的变化关系
Figure6. (a) Bubble radius as a function of time, (b) maximum compressive and tensile stresses at different distances (d₀ = 2, 6 and 11 μm) in tissues under different shear moduli (G = 10, 100, 200, 500 and 1000 kPa); (c) the bubble radius as a function of time, (d) maximum compressive and tensile stresses at different distances (d₀ = 2, 6 and 11 μm) in tissues under different viscosities (μ = 1, 5, 15, 50 and 100 mPa·s).

3 -->

3.3.3.超声声压的影响

图7给出了超声声压变化(p_A = 0.1—5 MPa)对组织内微泡声空化动力学以及应力的影响. 如图7(a)所示, 随着激励超声声压增大, 组织内微泡振动幅度明显增大, 声压从0.1 MPa增至5 MPa时微泡最大半径R_max从1.016 μm增至24.45 μm. 图7(b)给出了R_max、微泡最小半径R_min和微泡膨胀-压缩比ΔR₊/ΔR_–(见插图)随声压变化的曲线. 当声压较小(p_A < 0.3 MPa)时, 微泡膨胀幅度明显小于压缩幅度(ΔR₊/ΔR_– < 1), 呈现出脂类包膜微泡典型的“仅压缩”行为; 随着声压增大, 微泡膨胀幅度明显大于压缩幅度(ΔR₊/ΔR_– > 1), 说明脂类包膜微泡典型的“仅压缩”行为仅在激励声压较小的情况下出现. 因为较高声压下的微泡剧烈膨胀可导致其脂类包膜发生破裂, 使得“仅压缩”行为消失; 同时, 较高声压下微泡非线性动力学行为更显著, 从而微泡膨胀程度远大于压缩程度. 该现象与水中脂类包膜微泡声空化动力学行为类似^[29]. 图7(c)和图7(d)分别给出了组织内不同d₀处粒子所受最大挤压应力和最大拉伸应力随激励超声声压的变化曲线. 从图7(c)和图7(d)可以清楚地看出: 随着超声声压增大, 最大挤压应力和最大拉伸应力都增大; 随着d₀增大, 粒子所受应力逐渐减小, 其中拉伸应力减小程度比挤压应力更大.
图 7 (a)不同超声声压(0.1—5 MPa)下, 微泡半径R随时间的变化关系; (b)微泡最大半径R_max、最小半径R_min和微泡膨胀-压缩比ΔR₊/ΔR_– (见插图)随超声声压的变化关系; 组织内不同d₀处(2, 6和11 μm)粒子所受(c)最大挤压应力和(d)最大拉伸应力随超声声压的变化关系
Figure7. (a) Bubble radius as a function of time at different ultrasound pressures (0.1?5 MPa); (b) the maximum bubble radius R_max, minimum bubble radius R_min and the ratio ΔR₊/ΔR_– of the positive and negative radius excursions as a function of ultrasound pressure; (c) maximum compressive stresses and (d) maximum tensile stresses in tissues at three distances (d₀ = 2, 6 and 11 μm) as a function of ultrasound pressure.

超声空化机械效应是聚焦超声进行局部释药、基因导入和新生血管/组织毁损等肿瘤治疗的重要物理机制. 考虑常用的脂类包膜微泡/纳米相变液滴作为人工空化核可显著增强空化治疗效率, 本文通过耦合空化动力学Gilmore方程、组织黏弹性Zener模型和脂类微泡包膜材料的非线性黏弹性模型, 建立了组织内脂类包膜微泡动力学模型, 该模型综合考虑了组织黏弹性与应力松弛、包膜非线性黏弹性且可准确描述高马赫数($ \dot{R}/c < 2.2 $)情况下的微泡动力学; 随后, 分析了组织内微泡的声空化动力学行为, 量化了微泡周围组织内机械应力的时空分布; 最后, 探究了包膜材料、组织黏弹性和超声声压等参数对微泡振动以及组织内机械应力的影响. 结果表明, 生物组织和包膜材料的黏弹性都会抑制微泡振动, 但由于组织黏弹性更大, 对微泡振动的抑制作用比包膜材料更大; 当超声声压较小时, 包膜微泡膨胀受限出现典型的“仅压缩”行为, 而随着声压增大, 微泡“仅压缩”行为消失, 出现膨胀幅度明显大于压缩幅度的非线性振动.
包膜微泡声空化可在周围组织内产生机械应力: 在微泡膨胀阶段, 应力为挤压应力(负应力), 而在微泡收缩阶段和反弹初始阶段, 应力为拉伸应力(正应力); 应力主要局部分布于微泡壁附近, 随着距离增大而显著减小, 其中拉伸应力衰减率明显大于挤压应力. 组织内机械应力可分为弹性应力和黏性应力: 弹性应力在微泡半径R大于初始半径R₀时为挤压应力, 在R < R₀时为拉伸应力, 它主要由组织应变决定. 黏性应力在微泡膨胀时为挤压应力, 微泡收缩时为拉伸应力, 它主要由组织应变率决定. 组织内包膜微泡声空化产生的应力比无包膜微泡更小, 但超声声压较大时, 包膜对应力的减小可忽略不计. 随着组织弹性模量增大, 组织内应力逐渐减小, 但随着组织黏度增大, 应力先增大后减小. 随着超声声压增大, 组织内微泡振动显著增强, 组织内应力明显增大. 本研究可为进一步阐明超声治疗中组织机械损伤的内在机制奠定理论基础. 本文以单个微泡作为研究对象, 后续将重点考虑微泡-微泡之间的相互作用, 进一步阐释组织内多泡声空化动力学和微泡周围组织内应力分布规律.