非線性聲場下包膜微泡動力學與頻率響應分析*

索鼎杰 紀鎮祥 黃曉雲 靳杰 閆天翼?

1) (北京理工大學醫學技術學院,北京 100081)

2) (北京理工大學前沿交叉科學研究院,北京 100081)

包膜微泡在非線性聲場下的響應對高強度聚焦超聲治療具有重要意義.本文通過耦合Gilmore-Akulichev-Zener 模型與脂質包膜的非線性模型,采用KZK 方程仿真非線性聲場,分析在不同超聲頻率、不同聲壓以及不同包膜材料的黏彈性下的微泡動力學行為和微泡振蕩頻率,并進一步對比了實際測量聲場與KZK 方程仿真聲場下的微泡動力學行為和頻率響應.研究結果表明: 非線性聲場會導致微泡壁的瞬時運動速率減小,聲壓和頻率的改變對微泡動力學的影響與在線性聲場中類似;包膜材料的不同可以使振蕩頻率中的諧波分量發生改變,其中包膜材料的彈性對微泡的頻率響應影響較小,包膜的初始黏性和初始表面張力對微泡的振蕩頻率分布影響較大,當初始黏性越小時,二次分諧波的峰值越高,當初始表面張力越大時,主頻的峰值越高.本研究進一步闡明非線性超聲激勵包膜微泡的微泡動力學,為包膜微泡在非線性聲場下的頻率響應分析奠定理論基礎.

1 引言

聲空化[1-3]是指超聲波與傳播介質中的氣核之間的相互作用,在負聲壓的作用下使介質中氣核迅速擴張,出現空腔的現象;聲空化分為穩態空化和瞬態空化.穩態空化[4,5]是指氣泡穩定的振蕩,存在的時間較長,可用于聚焦超聲誘導的血腦屏障開放,通過調節穩態空化,可以降低腦損傷的風險,提高藥物的靶向遞送率;瞬態空化[6,7]是指聲壓達到慣性空化閾值后,氣泡劇烈的振蕩,存在的時間較短,瞬態空化可以產生沖擊波和微射流,可用于聚焦超聲組織摧毀[8,9]、血栓消融[10,11]、腫瘤治療[12,13]等.此外,空化效應還可能導致光學效應(聲致發光)[14,15]、化學效應[16,17]和其他熱效應[18,19].因此,對聲空化進行研究十分重要.

近年來,許多研究人員對單個微泡在超聲激勵下的動力學模型進行了研究和修正.Rayleigh-Plesset (RP)[20]方程通過將液體視為不可壓縮流體被廣泛運用于計算微泡動力學的基礎.Keller和Miksis[21]通過拓展RP 方程,將液體視為輕微可壓縮流體,使得方程在大氣泡和大振幅下更加適用.Doinikov[22]通過拓展Keller 和Miksis[21]的研究,將氣泡的平移運動方程和徑向運動方程進行耦合,研究了氣泡平移運動對微泡動力學的影響,隨后通過研究次級Bjerknes 力[23]和聲學流線的特性,表明次級Bjerknes 力對氣泡之間的相互作用起到了重要的作用,并且聲學流線的形成與次級Bjerknes 力密切相關,目前該模型已被廣泛應用.Yang 和Church[24]通過耦合Keller-Miksis (KM)方程與Kelvin-Voigt 方程,研究了軟組織中的聲空化閾值[25].Zilonova 等[26]通過耦合Gilmore-Akulichev 方程與Zener 黏彈性模型,使得方程在高馬赫數的情況下適用,并研究了高強度聚焦超聲在組織中的聲空化閾值.

為了降低聲空化閾值,通常使用微米涂層脂質包膜泡作為造影劑來提高空化治療時的安全性、穩定性和可控性.影響脂質包膜微泡聲空化的主要因素包括: 膜的黏性、膜的彈性、膜的初始表面張力等.Marmottant 等[27,28]提出了單個超聲造影劑氣泡的模型,該模型考慮了氣體微氣泡上脂質包膜的物理特性,目前已被廣泛地應用.于潔等[29]在此基礎上對脂質包膜材料的黏性參數進行了非線性修正,使得模型的應用更加廣泛.秦對等[30]通過耦合Gilmore-Akulichev-Zener (GAZ)模型與脂類包膜的非線性黏彈性模型,研究了微泡的聲空化動力學行為以及微泡周圍組織內機械應力.上述研究基本都是在線性波的驅動下進行數值分析,然而,高強度聚焦超聲在傳播時,會出現衍射、散射、衰減以及聲吸收的情況,導致波形畸變和諧波,產生非線性.通常用于模擬聚焦換能器發射超聲波的兩種常用非線性模型分別是KZK (Khokhlov-Zabolotskaya-Kuznetsov)[31]方程和Westervelt[32]方程.2019 年,Bakhtiari-Nejad 和Shahab[33]通過求解KZK 方程的數值解耦合KM 方程,研究了氣泡在非線性聲場下的徑向運動和平移運動.

本文主要采用KZK 方程進行高強度聚焦聲場仿真,耦合GAZ 方程與包膜的黏彈性模型,采用步長控制算法和Runge-Kutta5 (RK5)算法,探討在非線性聲場中不同超聲頻率、不同聲壓以及不同包膜材料的黏彈性下的微泡動力學行為和微泡振蕩頻率響應,進一步分析不同超聲參數和包膜材料對微泡振蕩頻率中諧波分量的影響,最后采用實際測量的非線性聲場數據與KZK 方程仿真的聲場數據進行對比,探究實際聲場和仿真聲場之間的差異.

2 理論模型

2.1 GAZ 模型

在數值仿真中,由于高強度聚焦超聲會導致微泡的大幅度振蕩,氣泡壁的運動速度高于1 Mach(1 Mach=340.3 m/s),GA 模型更適用高馬赫環境,并且Zener[26]模型相比于Kelvin-Voigt 模型能更好地適用于黏彈性組織,能較好地模擬微泡的徑向應力,因此本研究采用GAZ 模型:

其中,c0是流體中的聲速,R0是初始半徑,H是氣泡壁的焓變,μ是介質中的黏性,p0是靜態壓力,σ0是初始表面張力,ρ是組織密度,B和n是常數.B=-p0,p1是微泡內部的壓力,p是驅動聲壓,τrr|R是徑向應力,φ是組織的松弛時間,且q=.

2.2 包膜的非線性黏彈模型

包膜材料能提高微泡的穩定性,其黏彈特性能夠減小氣泡的壓縮和擴張.非線性修正的Kelvin-Voigt[28,29]模型能夠較好地考慮包膜微泡的黏彈性[34,35]:

其中,S1為包膜彈性項,S2為黏性項,且包膜表面張力σ(R) 的徑向變化如下:

其中,Rb為微泡半徑最小值,且Rb=,當R≤Rb時,微泡發生皺縮,Rr為微泡破裂半徑,且Rr=,且σr=σ0,當R≥Rr時,微泡發生破裂.其中k表示包膜微泡的黏度,λ為包膜微泡的彈性,且k=,k0為殼的初始黏性參數,α為特征時間常數.

2.3 驅動聲壓和數值條件

驅動聲壓采用KZK 方程的模擬,通過HIFUSimulator (Food and Drug Administration,Silver Spring,MD)[36]進行仿真:

其中p是在介質中的聲壓,z是軸向坐標(垂直于換能器表面),x和y為橫向坐標(平行于換能器表面),且延遲時間τ=t-z/c0,β為非線性系數,b為耗散系數.通過頻域的傅里葉級數解法展開可得

其中Cn為第n次諧波的振幅[37].

假設微泡內部的壓力服從范德瓦耳斯方程[38],則可以采用下式表示:

其中,h=R0/8.4 是微泡的范德瓦耳斯半徑,γ是多向性系數.

表1 參數名稱和值Table 1.Name and value of the parameters.

3 仿真結果分析與討論

3.1 非線性超聲的響應

超聲波在傳播過程中會出現衍射、散射、衰減現象,從而產生波形畸變,采用KZK 方程模擬超聲波在水中的非線性傳播.圖1 對比了線性超聲和非線性超聲在相同正壓下的微泡徑向運動行為和振蕩頻率響應(其中線性聲場采用驅動聲壓為:P=-Ptsin(2πft),Pt為振幅,f=1 MHz),初始半徑R0=1 μm,聲壓P=2.0 MPa.對比線性聲場和非線性聲場下的微泡動力學行為,可以發現超聲波傳播時的損耗導致負壓偏小,從而使得微泡的徑向運動減弱,氣泡壁的運動速度也相應減小,微泡的振蕩頻率峰值也隨之減小.圖2 對比了非線性聲場下不同聲壓和不同頻率下的微泡動力學行為和頻率響應,初始半徑R0=1 μm,當頻率固定在f=1.0 MHz,聲壓增大時,微泡的徑向運動顯著增強,從而使得振蕩響應頻率的峰值也增大.當聲壓固定在2.0 MPa 時,頻率增大時,微泡的徑向運動卻開始減小,這是由于微泡的振蕩主頻峰值減小和遠離共振頻率導致的.

圖1 線性聲場和非線性聲場下微泡動力學行為頻率響應 (a) 相同正壓下的線性超聲與非線性超聲波形圖;(b),(c) 微泡動力學行為的對比;(d) 微泡半徑變化的振蕩頻率響應,藍色線條表示驅動聲壓為KZK 方程仿真的非線性聲場,紅色線條表示線性聲場Fig.1.Dynamics behavior and frequency response of microbubbles under linear and nonlinear ultrasound fields: (a) Linear and nonlinear ultrasonic waveforms under the same positive pressure;(b),(c) comparison of bubble dynamics behaviors;(d) oscillation frequency response of the microbubble radius change,the blue lines represent the driving sound pressure as the nonlinear ultrasound field simulated by the KZK equation,and the red lines represent the linear ultrasound field.

圖2 (a),(b),(e) 不同聲壓下非線性超聲的微泡動力學行為和頻率響應的差異;(c),(d),(f) 不同頻率下非線性超聲的微泡動力學行為和頻率響應差異Fig.2.(a),(b),(e) Differences in the kinetic behavior and frequency response of microbubbles in nonlinear ultrasound at different acoustic pressures;(c),(d),(f) differences in microbubble kinetic behavior and frequency response of microbubbles in nonlinear ultrasound at different frequencies.

3.2 脂質包膜參數影響

圖3 對比了不同脂質包膜微泡的初始表面張力和初始黏性參數對微泡動力學行為和頻率響應的影響,初始條件設置如下:f=1 MHz,P=2.0 MPa,R0=1 μm.可以發現當膜的初始表面張力較大時,微泡的徑向運動增強,頻率響應中諧波的數量增加;當包膜的初始黏性減小時,微泡的徑向運動更加劇烈,頻率響應中二次分諧波的峰值升高.為了探究包膜材料對頻率響應的影響,采用同樣標準來判斷,即C=(P11max1表示最大響應頻率fmax1對于的頻譜值,P11max2表示第二響應頻率fmax2對于的頻譜值),若|C|越大,則微泡的頻譜峰值更高,若小于C＜0 且|C|越小,則表示分諧波得到了增強,若C＞0 且|C|越小,則表示次諧波得到了增強.圖4 對比了不同包膜參數對微泡動力學行為和頻率響應的影響,仿真條件都是:f=1 MHz,P=2.0 MPa,R0=1 μm.圖4(a),(b)分別表示包膜微泡的初始表面張力σ0在0—0.07 N/m 和彈性λ在0—1 N/m 時微泡動力學行為和頻率響應.可以發現包膜微泡的初始表面張力越大,微泡的徑向運動越劇烈,分諧波的頻率響應越大,但包膜微泡的彈性對微泡動力學行為和頻率響應幾乎沒有影響.圖4(c),(d)分別表示包膜的初始黏性k0在0.5×10-8—1.0×10-7kg/s 和初始彈性λ在0—1 N/m 時微泡動力學行為和頻率響應.可以發現當脂質包膜的初始黏性增大時,微泡的徑向運動減小,頻率響應中超諧波的分量增強;當包膜的初始黏性減小時,微泡的徑向運動增強,頻率響應也提高了分諧波的分量.圖4(e),(f)分別表示包膜的初始黏性k0在0.5×10-8—1.0×10-7kg/s和包膜微泡的初始表面張力σ0在0—0.07 mN/m時的微泡動力學行為和頻率響應,可以發現當包膜的初始黏性和初始表面張力一起改變時(包膜的初始黏性小且包膜微泡的初始表面張力越大),微泡振蕩的主頻和次諧波分量的峰值變高,因此在進行聚焦超聲消融時,可以通過改變脂質微泡的包膜材料,進一步調節信號的頻率分量,從而使分諧波或次諧波靠近共振頻率,提供聲空化效應.

圖3 不同包膜材料對微泡動力學行為和頻率響應影響 (a),(b) 包膜微泡的初始表面張力;(c),(d) 包膜的黏性參數Fig.3.Effects of different coating materials on bubble dynamic behavior and frequency response of encapsulated microbubbles:(a),(b) Initial surface tension of encapsulated microbubbles;(c),(d) viscosity parameters of encapsulated microbubbles.

圖4 不同包膜材料組合對微泡動力學行為和頻率響應影響 (a),(c),(e) 顏色表示微泡半徑變化;(b),(d),(f) 顏色表示C 的變化Fig.4.Effects of bubble shell on bubble dynamics and frequency response in the nonlinear ultrasound field: (a),(c),(e) Color bar of is the change of radius;(b),(d),(f) color bar is the change of C.

3.3 聲場測量結果與KZK 方程仿真對比

聲場實際測量的步驟為: 首先,信號發生器的輸入信號通過功率放大器進行放大并發送到匹配單元和聚焦換能器;然后,采用三軸定位系統將水聽器對準聚焦換能器的焦點位置,最后進行聲場測量.為了防止水中的氣泡發生空化,產生諧波,影響數值仿真,聚焦換能器和水聽器都在除氣水中進行測量.圖5(a)對比了真實的聚焦換能器在焦點處的超聲波形和KZK 方程仿真出的焦點處超聲波形,仿真參數如下:f=1 MHz,P=2.0 MPa,R0=1 μm.圖5(b),(c)是實際測量的聲場和KZK 方程仿真聲場下的微泡徑向運動對比,可以發現差距不大,但仿真下的氣泡壁運動速度與實際相差0.5 MHz.圖5(d)是微泡動力學行為的脈沖響應,可以發現真實的聲場在0.1 MHz 的頻譜響應較大,而在主頻1 MHz 的響應低于KZK 方程仿真的聲場.這可能是由于信號發生器發射的是脈沖類型,并且電壓在起始位置會發生損耗,導致聲壓的波形不一致,使得整個脈沖持續時間被當成是一個波形,導致0.1 MHz 的頻譜響應較大.

圖5 實際測量值與KZK 仿真對比圖 (a) 相同正壓下實際測量與KZK 方程仿真波形圖;(b),(c) 微泡動力學行為的對比;(d) 微泡半徑變化的振蕩頻率響應Fig.5.Actual measured values are compared with the KZK simulation: (a) Actual measurement and KZK equation simulation waveform diagram under the same positive pressure;(b),(c) comparison of the dynamic behavior of microbubbles;(d) oscillation frequency response of the microbubble radius change.

4 總結

微泡動力學分析是聚焦超聲應用于聲空化、開放血腦屏障、消融以及聲致發光中的重要理論基礎.考慮在高強度聚焦超聲治療中使用聚焦超聲換能器發射的超聲由于能量較高會在介質中傳播時產生非線性畸變,本文主要采用KZK 方程進行高強度聚焦聲場仿真,耦合Gilmore-Akulichev-Zener與包膜的黏彈性模型.相比之前的研究,我們在考慮超聲波傳播非線性效應的同時加入了脂質微泡模型,可以進行非線性聲場的脂質微泡的微泡動力學仿真和微泡振蕩頻率響應,使微泡動力學仿真和頻率響應下的環境更加趨近于真實聲場.研究結果表明: 當包膜材料不變時,超聲波傳播時的非線性會導致微泡的擴張減弱,但會產生更多的諧波分量,包膜微泡的振蕩頻率分量中諧波含量增大.當超聲聲壓增大時,微泡的徑向運動增強,各頻率的分量也增高;當超聲頻率靠近振蕩頻率時,微泡的徑向運動增強,頻率分量中主頻的峰值增高.

當改變包膜材料時,包膜微泡的初始黏性越小,振蕩頻率中分諧波的分量越大,包膜微泡的初始表面張力越大,振蕩頻率中出現的諧波分量越多,分諧波的分量越多.通過調節不同的包膜參數材料組合,可以發現包膜微泡的徑向運動行為與振蕩頻率的改變比較一致,當諧波靠近共振頻率時,微泡的聲空化效應得以增強;振蕩主頻分量越高,微泡的聲空化效應也會相應增強.本研究進一步為非線性超聲消融奠定了理論基礎.由于本文以單個微泡作為對象,且未將微泡內部與外部的熱傳遞和熱效應考慮進去,后續將重點考慮多微泡的徑向運動和平移運動,并結合微泡內部與外部的熱效應,進一步分析多微泡在非線性聲場激勵下的微泡動力學行為和熱能傳遞.