雙細(xì)胞光聲信號的有限元分析?

趙榮榮 楊 鵬 韓建寧

(中北大學(xué)信息與通信工程學(xué)院 太原 030051)

0 引言

光聲檢測作為一種高分辨率無標(biāo)簽的光學(xué)成像技術(shù)，結(jié)合了光學(xué)成像的高對比度和聲學(xué)成像的深度，可實(shí)現(xiàn)對處于自然狀態(tài)的生物粒子形態(tài)進(jìn)行無損傷測量，近年來發(fā)展迅速。光聲成像的本質(zhì)是光聲效應(yīng)，脈沖激光照射生物顆粒，生物顆粒由于吸收光能而溫度升高，體積發(fā)生膨脹和縮小，產(chǎn)生聲波，通過分析聲波，得到生物粒子的相關(guān)信息，細(xì)胞的頻域聲壓級曲線可用于評估細(xì)胞大小和形態(tài)等特征，作為疾病早期診斷的一種有前景的手段。紅細(xì)胞(Erythrocyte/red blood cell,RBC)中的血紅蛋白，憑借較高的光吸收系數(shù)，成為體內(nèi)光聲成像的理想內(nèi)源性造影劑[1]。根據(jù)紅細(xì)胞的頻域聲壓級曲線可以得到紅細(xì)胞的尺寸和形態(tài)等信息。細(xì)胞的形狀和功能之間有著密切的關(guān)系，因而細(xì)胞和細(xì)胞器形態(tài)的光聲測量，為了解細(xì)胞功能和疾病診斷提供重要的信息。

文獻(xiàn)[2]通過單細(xì)胞的頻域聲壓級曲線，可以實(shí)現(xiàn)對單個紅細(xì)胞的大小和形態(tài)的檢測，對正常細(xì)胞與異形細(xì)胞進(jìn)行定性鑒別。文獻(xiàn)[3]研究了紅細(xì)胞聚集體，用光聲信號的參數(shù)評估人體血液中紅細(xì)胞的聚集水平。然而，前者研究的對象是單個細(xì)胞，后者的研究對象是細(xì)胞聚集體，討論的內(nèi)容是光聲信號與細(xì)胞聚集程度的關(guān)系，并沒有關(guān)注細(xì)胞間光聲信號的相互影響。所以，一個細(xì)胞對另一個細(xì)胞在不同距離下的影響還沒有被關(guān)注。

人體內(nèi)正常紅細(xì)胞的存在狀態(tài)是游離的單個細(xì)胞，但細(xì)胞之間距離較近，信號會相互影響和干擾，且細(xì)胞間的距離是隨機(jī)的。在一項(xiàng)研究中[2]，提出了一種光聲快速定量單個紅細(xì)胞形態(tài)變化的方法，此方法比電阻抗和光散射方法更準(zhǔn)確，比血液涂片和光學(xué)干涉方法更快速簡便。這項(xiàng)研究表明，不同大小、形狀、取向和成分的紅細(xì)胞，在100 MHz 以上的光聲光譜中觀察到獨(dú)特的周期性變化的最小值和最大值，但是研究只提及光聲成像在測量大量樣本上的潛力，并沒有進(jìn)一步分析多細(xì)胞條件下的單細(xì)胞檢測情況。本文的目的是補(bǔ)充和拓展單細(xì)胞光聲檢測的內(nèi)容，討論距離對細(xì)胞光聲信號的影響，找到細(xì)胞間相互影響最小的距離。

1 分析方法

1.1 方法介紹

為了能夠在兩個細(xì)胞的光聲信號相互影響的條件下，通過調(diào)節(jié)細(xì)胞間的距離，獲得相互干擾最小的光聲信號，本文對單細(xì)胞和不同距離的雙細(xì)胞進(jìn)行有限元分析，研究建立了單細(xì)胞信號模型和不同距離的雙細(xì)胞信號模型。在理論模型中，一個紅細(xì)胞位于球形水環(huán)境的中心，稱其為主體紅細(xì)胞，另一個相同細(xì)胞位于其正上方，以400 MHz 的頻率計算單位波長，按照波長的倍數(shù)設(shè)置上方紅細(xì)胞與主體紅細(xì)胞的距離，得到了不同距離下主體紅細(xì)胞不同角度的頻域聲壓級曲線。另外，通過單細(xì)胞信號模型，計算出了沒有上方紅細(xì)胞干擾時，單個紅細(xì)胞不同角度的頻域聲壓級曲線。

判斷曲線相似性的方法主要有特征值法和相似性函數(shù)定義法。特征值法是對曲線特征參數(shù)進(jìn)行比較來度量兩個曲線的相似性，而相似性函數(shù)定義法是比較兩個曲線之間的距離來度量兩個曲線的相似性，相似性函數(shù)定義法在比較兩條曲線形狀的相似性上優(yōu)于特征值法。不同距離下紅細(xì)胞的頻域聲壓級曲線振幅相差不大，差異主要體現(xiàn)在曲線的斜率上，于是本文選擇了基于形狀判斷曲線相似性的離散frechet 距離的方法。基于離散frechet 距離的方法，以45?位置為例，通過frechet 數(shù)值度量細(xì)胞頻域聲壓級曲線的相似性，且假設(shè)frechet數(shù)值越小，相似性越高[4]。結(jié)果表明，frechet 數(shù)值最小的曲線與觀察得到的最相似的曲線是同一條曲線，即frechet 數(shù)值可以很好地度量頻域聲壓級曲線的相似性。

1.2 參數(shù)設(shè)置

本文利用COMSOL Multiphysics 實(shí)現(xiàn)了細(xì)胞的有限元分析，以檢測兩個紅細(xì)胞在不同距離下的光聲信號。研究建立了二維軸對稱模型，相比于三維模型，二維模型的計算量更小，且能產(chǎn)生有效的結(jié)果[5?7]。利用基本圖形的組合和分割，得到了雙凹形紅細(xì)胞幾何模型。同時，添加了壓力聲學(xué)瞬態(tài)物理場，建立了全局笛卡爾坐標(biāo)系。

正常紅細(xì)胞的形狀是雙凹形的圓餅狀，中間兩面微凹，邊緣較厚。紅細(xì)胞平均直徑約為7.8 μm，高為1～2 μm，體積為94 μm3。本文的理論模型是由Evans 等[8]開發(fā)的模型(雙凹形細(xì)胞位于10 μm 半徑球形系統(tǒng))擴(kuò)展得到。根據(jù)紅細(xì)胞的幾何物理參數(shù)，確定了紅細(xì)胞模型的幾何形狀，模型中的紅細(xì)胞直徑確定為7.82 μm，雙凹形細(xì)胞處于圓形的水環(huán)境中。為了進(jìn)行有限元分析，將細(xì)胞模型的幾何形狀劃分成自由剖分三角形網(wǎng)格，采用網(wǎng)格尺寸為標(biāo)準(zhǔn)。

表1 模型參數(shù)Table 1 Model parameters

在以往的研究中，已經(jīng)找到了紅細(xì)胞的最佳光吸收的激光束的波長，在血液樣本中，紅細(xì)胞可以被看作是入射光的主要吸收體和光聲信號的主要輻射源[9]。因此，本研究將細(xì)胞和水環(huán)境設(shè)置為理想狀態(tài)，只將紅細(xì)胞作為入射光輻射的吸收體和光聲信號輻射源，忽略外圍水對入射光的吸收和光聲信號的對外輻射，且紅細(xì)胞均勻吸收激光。激光束照射之前，細(xì)胞的初始壓力值為0 Pa，水的初始壓力值為0 Pa，激光束照射之后，細(xì)胞的初始壓力值為1 Pa，水的初始壓力值為0 Pa。紅細(xì)胞內(nèi)的質(zhì)量密度和聲速分別為1110 kg/m3和1650 m/s，與先前研究中測量的數(shù)據(jù)一致[10?11]。細(xì)胞外圍水環(huán)境的質(zhì)量密度和聲速分別為1000 kg/m3和1520 m/s。模型參數(shù)設(shè)置如表1 所示。時域的檢測范圍是0～10 ns，頻域的檢測范圍是0～1000 MHz，數(shù)據(jù)由COMSOL Multiphysics生成，數(shù)據(jù)的后處理由MATLAB R2014a完成。

1.3 單個紅細(xì)胞模型

在二維仿真模型中，以(0,0)點(diǎn)為圓心，以90 μm 為半徑，建立了水環(huán)境，紅細(xì)胞處于水環(huán)境中。如圖1(a)所示，雙凹形紅細(xì)胞處于圓形水環(huán)境的中心，在細(xì)胞的一側(cè)設(shè)置探測點(diǎn)。細(xì)胞模擬圖是細(xì)胞的側(cè)視圖，在側(cè)視圖中，建立直角坐標(biāo)系，設(shè)置水平的X軸和垂直的Y軸，X軸是細(xì)胞長直徑的方向，Y軸是細(xì)胞的厚度方向，且紅細(xì)胞的中心與圓形水環(huán)境的中心相互重合，x軸和y軸是主體紅細(xì)胞的對稱軸。紅細(xì)胞的尺寸和取向是影響單個紅細(xì)胞的光譜特征的因素，于是，選擇直徑為7.82 μm的雙凹形紅細(xì)胞平置于圓形水環(huán)境的中心。

圖1(b)是水環(huán)境中的細(xì)胞部分放大圖，以細(xì)胞中心為圓心，建立極坐標(biāo)系，以10 μm 為半徑，以?75?～75?為范圍，每隔15?布置一個探測點(diǎn)，檢測紅細(xì)胞的頻域聲壓級曲線。聲壓是大氣壓受到聲波擾動后產(chǎn)生的變化，在大氣壓強(qiáng)上的疊加一個聲波擾動引起的壓強(qiáng)變化，通過聲壓可以描述光聲信號的傳播情況。聲壓級定義為被測光聲信號的聲壓與基準(zhǔn)參考聲壓之比取以10為底的對數(shù)再乘以20，單位是分貝(dB)。

圖1 單細(xì)胞二維模型Fig.1 Single cell two-dimensional model

圖2 分別是單個細(xì)胞0?～75?、?15?～?75?的時域聲壓分布，通過觀察，由于曲線間重疊嚴(yán)重，很難得到曲線的相關(guān)規(guī)律，也難以得到細(xì)胞的形態(tài)特征。但是細(xì)胞的頻域聲壓級分布通過周期性變化的極大值和極小值來反映細(xì)胞的形態(tài)特征[12]。于是，對不同角度采集到的時域聲壓信號進(jìn)行傅里葉變換，得到頻域聲壓級分布[13]。

圖3 分別是單個紅細(xì)胞0?～75?、?15?～?75?的頻域聲壓級分布。探測點(diǎn)的角度不同，檢測到的頻域聲壓級曲線有著明顯的不同。觀察圖3(a)和圖3(b)，對比曲線15?與?15?、30?與?30?、45?與?45?、60?與?60?、75?與?75?，發(fā)現(xiàn)它們有著較高的相似度，這是因?yàn)閱蝹€紅細(xì)胞的形狀關(guān)于x軸對稱。但是，完全重合的部分大約是在0～500 MHz，超出500 MHz的范圍出現(xiàn)一些差異，這是由于聲波的散射。聲波的頻率越高，波長越短，個別分子或原子的運(yùn)動會影響聲波，且頻率越高，散射越強(qiáng)，導(dǎo)致在較高頻率時，對稱探測位置的仿真曲線不重合。

這里只分析從0?～75?的頻域聲壓級曲線分布情況。在0?方向，大約從200 MHz開始，出現(xiàn)周期性變化的最小值和最大值。在75?方向，整個頻譜只在800 MHz 左右出現(xiàn)了最小值。在15?、30?、45?、60?方向，大約分別從220 MHz、250 MHz、350 MHz、450 MHz開始，出現(xiàn)周期性變化的最小值和最大值。可以發(fā)現(xiàn)，從0?～75?，開始出現(xiàn)最小值的位置在x軸上不斷后移，在200 MHz～800 MHz 之間依次分布，且每條曲線的最小值之間的距離在不斷增大。在單個細(xì)胞的光聲頻譜中，0?位置可以觀察到最多的最小值點(diǎn)，75?位置觀察到最少的最小值點(diǎn)。同時，在這些曲線的第一個最小值(200 MHz)出現(xiàn)前，從0?～75?，曲線的振幅在不斷地增大。

圖2 單個紅細(xì)胞時域聲壓分布Fig.2 Time domain sound pressure distribution of a single RBC

圖3 單個紅細(xì)胞頻域聲壓級分布Fig.3 Sound pressure level distribution in frequency domain of a single RBC

1.4 不同距離的兩個紅細(xì)胞模型

在單個紅細(xì)胞模型的基礎(chǔ)上，以中心的紅細(xì)胞為主體紅細(xì)胞，添加另一個紅細(xì)胞。在本研究中，為了使細(xì)胞間的距離作為影響主體紅細(xì)胞光聲特征的唯一因素，需要將主體紅細(xì)胞的尺寸及取向和另一個紅細(xì)胞設(shè)置為相同狀態(tài)。主體紅細(xì)胞處于圓形水環(huán)境的圓心。另一個紅細(xì)胞處于主體紅細(xì)胞上方正對位置，其中心在y軸上。主體紅細(xì)胞在水環(huán)境中的位置固定，通過控制另一個細(xì)胞在y軸上位置，調(diào)節(jié)兩個細(xì)胞之間的距離，且兩個細(xì)胞的距離是兩個細(xì)胞中心之間的距離。如圖4(a)所示，是處于圓形水環(huán)境中兩個細(xì)胞距離為19 μm時的二維模型圖。

利用波的干涉原理，即兩個相同的波源，它們的兩列波疊加，在傳播區(qū)域會形成穩(wěn)定的加強(qiáng)區(qū)和穩(wěn)定的減弱區(qū)，加強(qiáng)區(qū)與減弱區(qū)相間分布，以此來調(diào)節(jié)兩個細(xì)胞之間的距離，使探測點(diǎn)的位置盡可能多地設(shè)置在加強(qiáng)區(qū)，使探測的光聲信號干擾最小。觀察圖3(a)和圖3(b) 可以發(fā)現(xiàn)，單個紅細(xì)胞的光聲頻譜在頻率是400 MHz 時，幾乎各個角度光譜振幅處于最大值或者最大值的附近，同時也表明，紅細(xì)胞的光聲頻譜信息在400 MHz 時最為豐富。細(xì)胞的光聲信號在37?C 水中傳播，水在37?C 的聲速是1520 m/s，根據(jù)頻率與波長的關(guān)系可知，400 MHz 的波長為3.8 μm。以3.8 μm 為單位，設(shè)置11組不同距離，分別是2.85μm(3/4波長)、3.8 μm、7.6 μm、11.4 μm、15.2 μm、19μm、22.8 μm、26.6 μm、30.4 μm、34.2 μm、38 μm。同樣，與單個紅細(xì)胞探測點(diǎn)的設(shè)置相同，以10 μm 為半徑，以?75?～75?為范圍，每隔15?布置一個探測點(diǎn)，如圖4(b)所示，描述了兩個細(xì)胞在距離為19 μm 時的相關(guān)設(shè)置。這樣，可以得到主體紅細(xì)胞在與另一個紅細(xì)胞不同距離時，各個角度的探測點(diǎn)的光聲信號。

圖4 雙細(xì)胞二維模型Fig.4 Two-cell two-dimensional model

2 分析結(jié)果

以主體紅細(xì)胞45?位置探測點(diǎn)的頻域聲壓級曲線為例，演示了在11 個不同的距離下，尋找與單個紅細(xì)胞最相似的曲線的過程，并進(jìn)行了效果對比。圖5(a)是主體紅細(xì)胞在45?位置下，距離為3/4λ ～4λ的頻域聲壓級曲線，將單個主體紅細(xì)胞與6 個不同距離下的頻域聲壓級曲線進(jìn)行對比。圖5(b)是主體紅細(xì)胞在45?位置下，距離為5λ ～10λ的頻域聲壓級曲線，將單個主體紅細(xì)胞與5 個不同距離下的頻域聲壓級曲線進(jìn)行對比。對比圖5(a)和圖5(b)可以發(fā)現(xiàn)，圖5(b)中有更多的曲線與紅色的單細(xì)胞曲線相似與重合，說明細(xì)胞間相互影響最小的距離，分布在5λ～10λ的可能性更高。

進(jìn)行頻域聲壓級曲線的定量分析時，主要尺度是曲線的振幅和斜率。根據(jù)圖5(a)和圖5(b) 可得，在26.6 μm距離下的曲線與單個紅細(xì)胞的曲線相似度最高，不論是振幅還是斜率都有著極高的相似度，兩個曲線幾乎重合，如圖6(a)所示，是主體紅細(xì)胞在45?位置下，距離為26.6 μm 的曲線與單個紅細(xì)胞45?位置的曲線進(jìn)行比較的結(jié)果。

圖5 45?位置不同距離與單細(xì)胞比較Fig.5 The different distances and at 45?position compared with single cells

圖6 45?位置最佳相似比較結(jié)果Fig.6 The best similarity comparison results of 45?position

圖6(b)是主體紅細(xì)胞在45?位置時，不同距離下的頻域聲壓級曲線與單個紅細(xì)胞的曲線進(jìn)行比較，得到的frechet 數(shù)值。此時，假設(shè)frechet 數(shù)值越小，相似度越高，通過觀察單個紅細(xì)胞在45?時的曲線與不同距離下主體紅細(xì)胞45?位置的曲線的相似程度，判斷假設(shè)是否正確。由圖6(b)可得，在26.6 μm 的距離下，frechet 數(shù)值最小。結(jié)果表明，frechet 數(shù)值最小的距離正是曲線相似度最高的距離，frechet數(shù)值可以很好地度量光聲曲線的相似度。

圖7 是主體紅細(xì)胞分別在0?、15?、30?、60?、75?位置時的最佳相似對比的結(jié)果。主體紅細(xì)胞的5 個位置執(zhí)行與上述45?位置相同的操作，即每個位置，將主體紅細(xì)胞11 個不同距離的頻域聲壓級曲線與單個紅細(xì)胞同位置的頻域聲壓級曲線對比，計算出11 個距離中frechet 數(shù)值最小的曲線，以曲線圖的方式呈現(xiàn)單個紅細(xì)胞同位置的頻域聲壓級曲線與frechet數(shù)值最小的曲線對比的結(jié)果。

圖7(a)是0?位置單個紅細(xì)胞與主體紅細(xì)胞34.2 μm 距離的對比結(jié)果；圖7(b) 是15?位置單個紅細(xì)胞與主體紅細(xì)胞30.4 μm 距離的對比結(jié)果；圖7(c)是30?位置單個紅細(xì)胞與主體紅細(xì)胞26.6μm距離的對比結(jié)果；圖7(d)是60?位置單個紅細(xì)胞與主體紅細(xì)胞19 μm 距離的對比結(jié)果；圖7(e)是75?位置單個紅細(xì)胞與主體紅細(xì)胞26.6 μm 距離的對比結(jié)果。除60?位置外，其他4 個位置曲線的重合度較高。

圖8 是主體紅細(xì)胞分別在?15?、?30?、?45?、?60?、?75?位置時的最佳相似對比的結(jié)果。主體紅細(xì)胞的5 個位置執(zhí)行與上述45?位置相同的操作，即每個位置，將主體紅細(xì)胞11 個不同距離的頻域聲壓級曲線與單個紅細(xì)胞同位置的頻域聲壓級曲線對比，計算出11 個距離中frechet 數(shù)值最小的曲線，以曲線圖的方式呈現(xiàn)單個紅細(xì)胞同位置的頻域聲壓級曲線與frechet 數(shù)值最小的曲線對比的結(jié)果。圖8(a)是?15?位置單個紅細(xì)胞與主體紅細(xì)胞30.4 μm 距離的對比結(jié)果；圖8(b)是?30?位置單個紅細(xì)胞與主體紅細(xì)胞30.4 μm 距離的對比結(jié)果；圖8(c)是?45?位置單個紅細(xì)胞與主體紅細(xì)胞38 μm 距離的對比結(jié)果；圖8(d)是?60?位置單個紅細(xì)胞與主體紅細(xì)胞38 μm 距離的對比結(jié)果；圖8(e)是?75?位置單個紅細(xì)胞與主體紅細(xì)胞26.6 μm 距離的對比結(jié)果。除?60?位置外，其他4個位置曲線的重合度較高。

在雙細(xì)胞模型中，主體紅細(xì)胞與另一個相同的細(xì)胞共有11 個不同的距離，每個距離下，有11 個不同角度的探測位置。將每個距離下的11 個角度的頻域聲壓級曲線設(shè)置為一組，計算11 個角度的頻域聲壓級曲線與同位置單個紅細(xì)胞頻域聲壓級曲線的frechet 數(shù)值，求取每組11 個角度的frechet 數(shù)值的均值。圖9 是主體紅細(xì)胞11 個距離下frechet數(shù)值的均值，隨著距離的增大，frechet均值整體呈現(xiàn)先下降后上升的趨勢，在26.6μm時達(dá)到最小值。即隨著距離的增大，主體紅細(xì)胞頻域聲壓級曲線與單個紅細(xì)胞同位置頻域聲壓級曲線的相似度先增大后減小，在26.6 μm 處相似度達(dá)到最高。由此可得，雙細(xì)胞模型中，兩個細(xì)胞在26.6 μm 時，光聲信號相互影響最小，且在30?、45?、75?、?75?位置達(dá)到最佳相似，有較高的相似度。

圖7 0?、15?、30?、60?、75?最佳相似對比Fig.7 The optimum similarity comparison of 0?,15?,30?,60?and 75?

圖8 ?15?、?30?、?45?、?60?、?75?最佳相似對比Fig.8 The optimum similarity comparison of ?15?, ?30?, ?45?, ?60?, ?75?

圖9 不同距離下frechet 數(shù)值的均值Fig.9 Mean value of frechet values at different distances

3 結(jié)論與討論

在先前的研究中，已經(jīng)實(shí)現(xiàn)了光聲快速定量單個紅細(xì)胞的形態(tài)，但是沒有進(jìn)一步討論在多細(xì)胞條件下如何檢測單細(xì)胞的光聲信號，即沒有考慮細(xì)胞信號之間的相互影響。而且，對光聲譜進(jìn)行定量分析時，采用的方法是分析光聲譜中某頻率的振幅和某段頻率的斜率，只提取了光聲譜的部分信息，適用于被測細(xì)胞的形狀有較大差異的情況。為此，本文討論了兩個紅細(xì)胞在不同距離下光聲信號之間的影響。對單個紅細(xì)胞和不同距離的兩個紅細(xì)胞進(jìn)行了有限元分析。以圓形水環(huán)境中心的紅細(xì)胞為主體紅細(xì)胞，在主體紅細(xì)胞右側(cè)的不同角度布置探測點(diǎn)，調(diào)整另一個紅細(xì)胞與主體紅細(xì)胞的距離，計算得到主體紅細(xì)胞不同距離不同角度下的頻域聲壓級曲線。引入frechet 距離來度量曲線的相似性，分析全部頻率段，量化曲線間的細(xì)微差異，將frechet 距離作為分析光聲曲線有力的工具。利用frechet數(shù)值來度量兩條光聲曲線的相似性，frechet 數(shù)值越小，兩條曲線越相似。其中，計算細(xì)胞間信號的相互影響，使用相似性函數(shù)定義法度量光聲曲線的相似性，在以往的研究中還沒有被討論過。

當(dāng)主體紅細(xì)胞的檢測角度一定時，在不同距離的主體紅細(xì)胞的光聲曲線中，存在frechet 數(shù)值最小的曲線，即這條曲線是此角度下的最佳相似曲線。同樣，觀察可得，當(dāng)兩個細(xì)胞距離為26.6 μm時，11 個角度的frechet 數(shù)值的均值最小，曲線相似度最高。

本文的實(shí)驗(yàn)原理是聲波干涉，這種方法不僅可以應(yīng)用于紅細(xì)胞，幾乎任何含有內(nèi)源性和外源性造影劑的細(xì)胞，都可以通過此方法找到細(xì)胞間光聲信號干擾最小的距離。由于實(shí)驗(yàn)條件的限制，本文研究的背景為體外檢測，首先討論了兩個相同紅細(xì)胞之間光聲信號的相互影響，再次，將研究多個不同細(xì)胞之間光聲信號的相互影響，最終實(shí)現(xiàn)體外同時檢測多個不同細(xì)胞的光聲信號的目標(biāo)，目的是補(bǔ)充和拓展光聲顯微技術(shù)中多細(xì)胞光聲檢測的內(nèi)容。