袁亮
Yuan Liang
新疆大學(xué)機(jī)械工程學(xué)院,新疆烏魯木齊市 830047
School of Mechanical Engineering, Xinjiang University, Urumqi, Xinjiang, 830047
生物細(xì)胞的功能主要取決于細(xì)胞內(nèi)物質(zhì)積極和定向運(yùn)動1。這些物質(zhì)非常多樣,包括細(xì)胞內(nèi)蛋白質(zhì)復(fù)合物,細(xì)胞骨架聚合物,核糖核蛋白顆粒,和膜細(xì)胞。由于熒光顯微鏡功能強(qiáng)大,可以對活細(xì)胞內(nèi)移動結(jié)構(gòu)直接觀察2,這就使得它成為我們研究這些物質(zhì)的運(yùn)動的強(qiáng)有力工具。為了提取所產(chǎn)生的圖像序列的動力學(xué)數(shù)據(jù),我們采用實(shí)時或延時視頻錄像以在線或離線的方式跟蹤大量隨時間變化的運(yùn)動目標(biāo)。
神經(jīng)細(xì)胞內(nèi)運(yùn)動的重要例子是軸突上的運(yùn)動,這是沿神經(jīng)細(xì)胞軸突上的亞細(xì)胞和分子結(jié)構(gòu)運(yùn)動3,4。軸突的一個有趣的特點(diǎn)是:他們是非常細(xì)長的圓柱形突起,通常直徑上只有幾個微米,但長度可達(dá)幾個毫米、厘米、甚至米。因此,軸突通常是將一個復(fù)雜的三維運(yùn)動模型減少到一維模型,其中大部分沿軸突的長軸向前或向后移動,如圖1所示。
圖1 神經(jīng)絲蛋白質(zhì)在培養(yǎng)的老鼠皮質(zhì)神經(jīng)元軸突中運(yùn)動的例子。
神經(jīng)絲蛋白質(zhì)在軸突上運(yùn)動,測量直徑約10納米,長幾微米5,是軸突細(xì)胞骨架的主要組成部分。它對齊于平行軸突的長軸做間歇的雙向移動,這種移動是停頓時間長短不同的短時間向前或向后運(yùn)動6,7。該微絲移動迅速平均速率約0.5微米/秒,但其整體運(yùn)動速度卻要慢得多,因?yàn)樗麄円ù罅康臅r間暫停7。由于運(yùn)動的持續(xù)時間,以及移動和暫停狀態(tài)之間頻率的轉(zhuǎn)變都是高度可變的,因此其運(yùn)動可視為隨機(jī)過程。為了分析運(yùn)動動力性能,我們必須在連續(xù)幀視頻中跟蹤神經(jīng)絲蛋白質(zhì)。然而,目前這一過程完全由手動完成,這大大增加了勞動強(qiáng)度和人為錯誤影響8。因此,建立一個完全自動跟蹤方法是必要的,這樣可以提高跟蹤效率和消除人工操作引起的主觀性和潛在的可變性。
視頻跟蹤已應(yīng)用于各種不同類型的細(xì)胞內(nèi)運(yùn)動,例如9,10,11。傳統(tǒng)的細(xì)胞內(nèi)跟蹤方法分為兩個主要步驟:檢測和定位,對應(yīng)。在檢測和定位的步驟中,目標(biāo)的位置是通過局部圖像降噪和分割來確定的。關(guān)鍵是對應(yīng)步驟,就是要在連續(xù)幀的視頻之間建立同一物體的對應(yīng)關(guān)系。由于若干因素的影響,包括快速運(yùn)動,物體變形,圖像噪聲,對象合并和分裂等等,這使得跟蹤細(xì)胞內(nèi)運(yùn)動的目標(biāo)變得較為困難。Jaqaman等人12,13提出了一種跟蹤算法,被稱為線性分配法,這種方法可以在對應(yīng)步驟中解決對象消失,合并和分裂的問題。Yang等人14基于Kalman濾波方法設(shè)計(jì)了可靠的跟蹤大規(guī)模密集反平行粒子運(yùn)動的方法。然而,他們的方法對低信噪比圖像無法工作15,16。對于跟蹤有噪聲圖像序列中的物體,Smal17等人提出了一種使用粒子濾波的魯棒算法。他們的方法考慮到衍射極限對象的光學(xué)模糊,用一個點(diǎn)擴(kuò)展函數(shù)模型計(jì)算觀測模型的可能性。然而,這種方法使用通用粒子濾波,增加了計(jì)算負(fù)擔(dān)。這是因?yàn)樵谕ㄓ昧W訛V波算法中,為保證數(shù)據(jù)的可靠性,必須使用大量粒子,而每個觀察的粒子都需要一個二維的模型計(jì)算。當(dāng)粒子的數(shù)量增加時,總的計(jì)算時間就會快速增加。
在本文中,我們提出了一種新的空間約束的粒子濾波跟蹤方法。這種方法最鮮明的特點(diǎn)是在神經(jīng)絲蛋白質(zhì)的跟蹤運(yùn)動中利用軸突的幾何形狀。簡單地說,我們利用了軸突狹長的突起結(jié)構(gòu)并且軸突上運(yùn)動的物體必須在軸突內(nèi)運(yùn)動。因此,在空間上對先驗(yàn)動態(tài)狀態(tài)能進(jìn)行約束,從而明顯降低了在下一個視頻中可能性狀態(tài)分布的方差。這樣,在不影響跟蹤精度的情況下可以減少一些粒子的使用量,從而提高了跟蹤效率和準(zhǔn)確性。
正如圖1所示,培養(yǎng)的神經(jīng)元軸突跟蹤軌跡為一條光滑曲線。為了建立這個模型的形狀,我們使用三次樣條插值獲得近似的分段多項(xiàng)式軸突路徑。利率曲線被節(jié)點(diǎn)分為多條線段,每條線段是由不同的多項(xiàng)式樣條插值函數(shù)建模而成。[u; v]是任何像素在圖像坐標(biāo)系統(tǒng)的位置。對于給定節(jié)點(diǎn)的水平坐標(biāo)u,每一段曲線段的垂直坐標(biāo)v由三次樣條插值函數(shù)給定:
這里i代表節(jié)點(diǎn)下標(biāo)(1,2, …, n),c0i,c1i,c2i和 c3i是每個部分的多項(xiàng)式系數(shù)并可以使用實(shí)際的像素坐標(biāo)內(nèi)曲線計(jì)算。這樣,假如曲線被分為n - 1段,4n個參數(shù)能被確定。為了獲得這4n個未知參數(shù),要求有4n個條件:位置,連續(xù)性,一階導(dǎo)數(shù)的連續(xù)性,二階導(dǎo)數(shù)在節(jié)點(diǎn)兩邊界條件。一個自然樣條,即在其中的兩端不封閉,邊界條件是讓二階導(dǎo)數(shù)在兩個端點(diǎn)處等于零。圖2顯示仿照這種方式的軸突實(shí)例。由于神經(jīng)絲蛋白質(zhì)的運(yùn)動受到軸突的限制,我們使用這個路徑定義軸突的中軸線。為了得到神經(jīng)絲蛋白質(zhì)路徑,我們在整個視頻序列中將最大像素強(qiáng)度投影到同一平面上(圖2b),然后我們手動地分出節(jié)點(diǎn)(圖2c)。
圖2 用三次樣條差值建立軸突模型。 (a)運(yùn)動神經(jīng)絲蛋白質(zhì),(b)是整個視頻序列每一圖像的最大強(qiáng)度投影,(c)是通過重疊在最大強(qiáng)度投影圖像上的三次樣條插值(白線)獲得的曲線擬合結(jié)果。
在跟蹤應(yīng)用中,動態(tài)狀態(tài)的元素選擇依據(jù)是對象的運(yùn)動特性。要跟蹤一個非對稱的移動物質(zhì),如神經(jīng)絲蛋白質(zhì),所需的元素包括位置,速度和方向。為方便起見,我們將神經(jīng)絲蛋白質(zhì)的外觀模型設(shè)定為一個邊界框長度(l)和寬度(w)的矩形的方塊。那么此框就表示濾波算法中的粒子。然后,我們根據(jù)移動神經(jīng)絲的長度和軸突的彎曲程度分別為每個視頻序列選擇邊界框的尺寸。對于直的軸突一個窄框是足夠的,但是對于有急彎的軸突,我們必須增加邊界框的寬度,來保證能在所有幀中采樣到完整的神經(jīng)絲蛋白質(zhì)。在時間t的邊界框的位置、速度和方向可以由動態(tài)定義,這里表示在圖像中坐標(biāo)系方塊中心的位置,θt表示方塊的轉(zhuǎn)角是對應(yīng)的速度。動態(tài)狀態(tài)可表示為如下形式:
這里Δt是兩禎之間的時間間隔。Ν表示正態(tài)分布,σu,σv和σθ分別是分布和θt的方差。我們這里使用正態(tài)分布是由于重要性密度函數(shù)具有高斯形式。在方程(2)中的均值和方差可以根據(jù)經(jīng)驗(yàn)來選擇。應(yīng)該指出的是,我們除了位置和速度之外還使用轉(zhuǎn)角來確定的動態(tài)狀態(tài),這與許多跟蹤方法不同。方向變量也許用于跟蹤軸突運(yùn)輸中的球狀物質(zhì)時不那么重要,但是當(dāng)用于跟蹤細(xì)胞骨架聚合物如神經(jīng)絲蛋白質(zhì)時卻非常重要,因?yàn)檫@些物質(zhì)是長的不對稱結(jié)構(gòu)。
使用三個決定動態(tài)狀態(tài)的向量元素序列(u, v,和θ)的結(jié)果是產(chǎn)生大量的粒子,這將會降低粒子濾波算法的工作效率。但是,粒子的數(shù)量可以通過約束多個向量元素來降低。
由于神經(jīng)蛋白質(zhì)必須在軸突路徑上運(yùn)動,所以粒子的方向在某個特定位置上必須與軸突的位置保持一致。這也就是,在時刻t的方向θt應(yīng)當(dāng)是軸突方程v(ut)的切角方向。相應(yīng)地,動態(tài)狀態(tài)模型方程(2)可以改寫為:
并且
我們稱含有方向約束的粒子濾波算法為方向約束粒子濾波(orientationally constrained particle filtering(OCPF)).
除了對神經(jīng)絲蛋白質(zhì)方向進(jìn)行約束外,軸突也限制了蛋白質(zhì)位置。為了對位置約束建模,我們僅限粒子濾波算法中的粒子分布在三次樣條曲線擬合程序所定義的軸突上??紤]到三次樣條曲線插值的不準(zhǔn)確,我們用一個寬度為2d?的很窄的條來代表軸突,這里窄帶的中心線就是軸突的三次樣條曲線(圖3所示)。d?是由曲線的擬合誤差決定的。假如在時刻t一個粒子到中心線的距離(dt)小于d?,這個粒子就被接受;如果dt>d?,神經(jīng)經(jīng)蛋白質(zhì)在這個位置的可能性就為零,這個粒子就是被拒絕。這種運(yùn)動神經(jīng)絲蛋白質(zhì)的位置約束就可以建模為:
現(xiàn)在整個先驗(yàn)動態(tài)狀態(tài)方程就能改寫為:
聯(lián)合位置約束方程(6) 和方向約束方程(3)和(4),我們用接受-拒絕機(jī)制來產(chǎn)生粒子,具體如下三個步驟:
Step 1: 根據(jù)(3)和(4)產(chǎn)生粒子tx′;
Step 2: 計(jì)算距離(td);
Step 3: 假如 dt<d*,接受這個粒子xt′為xt;否則,返回Step 1.
現(xiàn)在粒子可以從位置約束動態(tài)模型和方向約束動態(tài)模型中產(chǎn)生。我們稱這個方法為空間約束粒子濾波(spatially constrained particle filtering (SCPF)).
圖3 t時刻粒子在軸突邊界內(nèi)的分布概略圖。實(shí)線代表三次樣條插值描述的軸突中心線,虛線代表軸突邊界。黑色矩形框代表在定義的軸突約束下的軸突內(nèi)粒子的分布。粗黑色矩形框表示在這個最大似然估計(jì)例子下的粒子,它對應(yīng)于神經(jīng)絲蛋白質(zhì)的位置。
實(shí)驗(yàn)設(shè)置如下:所有的實(shí)驗(yàn)序列都是從小鼠大腦皮質(zhì)分離的神經(jīng)元中記錄下來的,該項(xiàng)工作由美國俄亥俄州立大學(xué)Anthony Brown教授實(shí)驗(yàn)室中完成。從顯微鏡測得的原始圖像為512×512像素,倍率為0.131μm/像素。首先,我們將進(jìn)行跟蹤初始化和參數(shù)描述,然后我們將介紹的跟蹤實(shí)驗(yàn)的結(jié)果。
所有神經(jīng)絲蛋白質(zhì)都大約10納米寬。經(jīng)驗(yàn)上,我們認(rèn)為用10個像素的(大約650納米)矩形框?qū)挾葁足夠在我們的圖像中包圍一個神經(jīng)絲蛋白質(zhì)。移動神經(jīng)絲的長度是易變的,但平均是5-10微米長。這樣,我們設(shè)定矩形框的寬度為10個像素,長度隨跟蹤神經(jīng)絲長度的不同而相應(yīng)變化。從中心到兩邊的距離d?都選為5個像素。在水平和垂直方向的粒子分布的高斯噪聲方差(uσ和vσ)大致設(shè)定為25像素,粒子方向的方差(θσ)設(shè)定為0.5 rad。在觀測似然性的計(jì)算中,參數(shù)λ設(shè)定為20。這種跟蹤算法使用Python在2.1GHz的英特爾酷睿2雙核計(jì)算機(jī)上運(yùn)行。為了計(jì)算神經(jīng)絲蛋白質(zhì)的實(shí)際運(yùn)行速度,我們使用了MetaMorph軟件手動跟蹤神經(jīng)絲蛋白質(zhì)的運(yùn)動。這個速度用作參考值來計(jì)算自動粒子跟蹤算法中的誤差。
我們的實(shí)驗(yàn)包括通用粒子濾波跟蹤,方向約束粒子濾波跟蹤,空間約束粒子濾波跟蹤,并對其性能進(jìn)行了評估和比較。圖4 (a) 和(b) 顯示出采用通用粒子濾波算法中的100和200個粒子的跟蹤結(jié)果。從圖中可以看出無論在位置還是在方向上都出現(xiàn)了大量的跟蹤誤差,即便200個粒子的跟蹤結(jié)果會比100個粒子略為好一些。這樣,100個和200個粒子的通用粒子濾波算法無法用于跟蹤這些錄像中的神經(jīng)絲蛋白質(zhì)。圖4 (c) 和(d) 顯示了采用50個和100個粒子的方向約束粒子濾波算法的跟蹤結(jié)果。正如我們在圖4 (c)中所看到的,50個粒子的方向約束粒子濾波算法在某些禎中仍然存在著跟蹤誤差。圖4(d) 顯示出了100個粒子的跟蹤精度有所提高,但仍然不能滿足所有禎。圖4 (e) 顯示了僅僅50個粒子的空間約束粒子濾波算法表現(xiàn)出了最佳的跟蹤結(jié)果。
圖4 三種不同粒子濾波算法跟蹤結(jié)果的比較。這些數(shù)據(jù)對應(yīng)于視頻中的第42,50,52,55,61,和77幀。(a) 100個粒子的通用粒子濾波算法(GPF) 。(b) 200個粒子的通用粒子濾波算法(GPF) 。(c) 50個粒子的方向約束的粒子濾波算法(OCPF)。(d) 100個粒子的方向約束粒子濾波算法(OCPF) 。(e) 50個粒子的空間粒子濾波算法(SCPF) 。在(a),(b),(c),(d)和(e)上排表示粒子在跟蹤時的分布,下排表示跟蹤結(jié)果。
圖5 三種不同粒子濾波算法估計(jì)的神經(jīng)絲蛋白質(zhì)運(yùn)動速度和實(shí)際的運(yùn)動速度的比較。 (a)200個粒子的通用粒子濾波算法(GPF)。(b)50個粒子的方向約束的粒子濾波算法(OCPF)。(c)100個粒子的方向約束粒子濾波算法(OCPF)。(d)50個粒子的空間粒子濾波算法(SCPF)。注意:使用 GPF和OCPF獲得的估計(jì)速度和實(shí)際速度有相當(dāng)大的差異,而使用SCPF獲得的估計(jì)速度與實(shí)際速度很相近。
圖6 三種不同粒子濾波跟蹤算法跟蹤誤差的比較。(a)200個粒子的通用粒子濾波算法(GPF)。(b) 50個粒子的方向約束的粒子濾波算法(OCPF)。(c)100個粒子的方向約束粒子濾波算法(OCPF)。(d)50個粒子的空間粒子濾波算法(SCPF)。
為了進(jìn)一步分析空間約束粒子濾波算法的跟蹤效果,每個跟蹤實(shí)驗(yàn)我們使用通用粒子濾波算法、方向約束粒子濾波算法和空間粒子濾波算法分別重復(fù)執(zhí)行了100次,并計(jì)算每個時間間隔神經(jīng)絲蛋白質(zhì)速度跟蹤誤差的平均值和標(biāo)準(zhǔn)偏差。跟蹤誤差定義為實(shí)際速度(通過有經(jīng)驗(yàn)專家手工跟蹤獲得)和估計(jì)速度(用粒子濾波算法法是計(jì)算有效粒子的數(shù)量(Neff),具體公式如下:自動跟蹤獲得)的比值。圖5顯示了用200個粒子的通用粒子濾波算法,50個和100個粒子的方向約束粒子濾波算法以及50個粒子的空間約束粒子濾波算法的跟蹤速度和手動跟蹤的速度相比較。圖6顯示出了相應(yīng)的均值跟蹤誤差。200個粒子的通用粒子濾波在每一幀中都顯示出了大量的跟蹤誤差;50個和100個粒子的方向約束粒子濾波算法也不能得到令人滿意的跟蹤結(jié)果。空間約束粒子濾波算法產(chǎn)生的跟蹤誤差最小,產(chǎn)生的速度也最接近于手動跟蹤結(jié)果(如圖5(d)和圖6(d))。這些結(jié)果表明通用粒子濾波算法無法用于跟蹤延時視頻圖像序列中的運(yùn)動神經(jīng)絲蛋白質(zhì)。跟蹤效果能通過方向約束來提高,最好是同時使用方向和位置這兩個約束。
圖7 不同跟蹤算法的有效粒子數(shù)量和計(jì)算時間的比較。(a) 有效粒子數(shù)量。星號,空心圓和空心三角分別代表使用200個粒子的GPF,100個粒子的OCPF,和50個粒子的 SCPF所獲得的數(shù)據(jù)。(b) 計(jì)算時間。
為了進(jìn)一步檢驗(yàn)通用粒子濾波算法,方向約束粒子濾波算法和空間約束粒子濾波算法的運(yùn)行效率,我們比較了他們的運(yùn)行時間。在所有粒子濾波算法中,計(jì)算時間原則上由粒子的數(shù)量和根據(jù)觀測模型的粒子估計(jì)來決定。本論文中,我們在不同算法中使用了同一觀測模型,因此計(jì)算時間基本上與粒子使用的數(shù)量成正比。這樣50個粒子的空間約束粒子濾波算法的計(jì)算時間是100個粒子的方向約束粒子濾波算法計(jì)算時間的一半,是200個粒子的通用粒子濾波算法計(jì)算時間的四分之一(圖7(b))。因此,與通用粒子濾波算法和方向約束粒子濾波算法相比,空間約束粒子濾波算法對在軸突上跟蹤神經(jīng)絲蛋白質(zhì)運(yùn)動更準(zhǔn)確和效率更高。
軸突神經(jīng)絲蛋白質(zhì)運(yùn)動行為的非線性和不穩(wěn)定性對自動跟蹤算法提出了挑戰(zhàn)。由于粒子濾波在處理非線性和不確定性運(yùn)動的優(yōu)越性能,它很自然的被用來解決這個問題。然而,由于通用粒子濾波必須使用大量的粒子來實(shí)現(xiàn)有效的跟蹤,這使得它的計(jì)算強(qiáng)度較高。在本文中,我們提出了一個新的空間約束的粒子濾波方法跟蹤神經(jīng)絲蛋白質(zhì)。這種方法利用神經(jīng)絲蛋白質(zhì)的運(yùn)動范圍被限制在一個狹長的軸突內(nèi)這一特點(diǎn),對神經(jīng)絲蛋白質(zhì)可能的運(yùn)動方向和位置進(jìn)行了動力學(xué)限制。為了將這兩種約束用到粒子濾波算法中,我們利用三次樣條插值建立軸突路徑的模型,使粒子的位置限制在一個狹窄的路徑上,位置的方向限制在軸突路徑曲率的切線方向。同時使用這兩個約束增加了有效粒子數(shù),從而減少了動態(tài)模型的不確定性。為了進(jìn)行測試,我們比較了空間約束的粒子濾波與通用粒子濾波,以及其中只考慮方向約束的粒子濾波。結(jié)果表明,空間約束的粒子濾波方法在效率上大大高于通用或定向約束的粒子濾波方法。因此有必要采取位置和方向的兩個限制實(shí)現(xiàn)最優(yōu)跟蹤效率。
1 D. Bray. Cell Movements: From Molecules to Motility [M]. 2nd ed. New York: Garland Science, 2000.
2 R. D. Goldman and D. L. Spector. Live Cell Imaging: A Laboratory Manual [M].2nd ed. New York: Cold Spring Harbor Laboratory, 2009.
3 A. Brown. Slow axonal transport, in New Encyclopedia of Neuroscience [M]. L. R.Squire, Ed. Oxford, U.K.: Academic, 2009,1–9.
4 A. Brown. Axonal transport of membranous and non-membranous cargoes:A unified perspective [J]. The Journal of Cell Biology, 2003, 160(6):817–821.
5 R. Perrot, R. Berges, A. Bocquet, and J.Eyer. Review of the multiple aspects of neurofilament functions, and their possible contribution to neurodegeneration [J].Cellular and Molecular Neurobiology, 2008,38(1): 27–65.
6 A. Brown. Slow axonal transport: Stop and go traffic in the axon [J]. Nature Reviews Molecullar Cell Biology, 2000,1(2): 153–156.
7 A. Brown. Axonal transport. In“Neuroscience in the 21st century”, ed. D.W.Pfaff [M], Springer, New York, 2013,255-308.
8 A. Brown and P. Jung. A critical reevaluation of the stationary axonalcytoskeleton hypothesis[J]. Cytoskeleton,2013, 70(1):1-11.
9 E. Meijering, I. Smal, and G. Danuser.Tracking in molecular bioimaging [J]. IEEE Signal Process Magzine, 2006, 23(3):46–53.
10 I. F. Sbalzarini and P. Koumoutsakos.Feature point tracking and trajectory analysis for video imaging in cell biology[J]. Journal of Structural Biology, 2005,151(2): 182–195.
11 J. Zhu and L. Yuan. Neurofilament tracking by detection in fluorescence microscopy images [C], in proceeding of international conference on image processing, 2013, 3123-3127.
12 K. Jaqaman, D. Loerke, M. Mettlen, H.Kuwata, S. Grinstein, S. Schmid, and G.Danuser. Robust single-particle tracking in live-cell time-lapse sequences [J]. Nature Methods, 2008, 5(8): 1212–1221.
13 K. Jaqaman and G. Danuser.Computational image analysis of cellular dynamics: A case study based on particle tracking [J]. Cold Spring Harb Protocal,2009, 2009(12): pdb.top65.
14 G. Yang, A. Matov, and G. Danuser.Reliable tracking of large scale dense antiparallel particle motion for fluorescence live cell imaging [C]. Proceding of IEEE Conference on Computer Vision and Pattern Recognition, San Diego, CA,2005: 9–17.
15 M. K. Cheezum,W. F.Walker, andW. H.Guilford. Quantitative comprison of algorithms for tracking single fluorescent particles [J]. Biophysical Jurnal, 2001, 81(4):2378–2388,.
16 B. C. Carter, G. T. Shubeita, and S.P.Gross. Tracking single particles:A user-friendly quantitative evaluation [J].Physical Biology, 2005, 2(1): 60–72,.
17 I. Smal,K.Draegestein, N.Galjart,W.Niessen, and E.Meijering. Particle filtering for multiple object tracking in dynamic fluorescence microscopy images:Application to microtubule growth analysis[J]. IEEE Transacitons on Medical Imaging,2008, 27(6): 789–803.