基于檢測的熒光顯微圖像中的神經(jīng)絲蛋白質(zhì)跟蹤*

袁 亮，朱俊達

(1.新疆大學(xué)機械工程學(xué)院，新疆 烏魯木齊 830047;2.澳門大學(xué)電子與計算機工程系,澳門 519000)

基于檢測的熒光顯微圖像中的神經(jīng)絲蛋白質(zhì)跟蹤*

袁 亮1，朱俊達2

(1.新疆大學(xué)機械工程學(xué)院，新疆 烏魯木齊 830047;2.澳門大學(xué)電子與計算機工程系,澳門 519000)

神經(jīng)絲蛋白質(zhì)是沿神經(jīng)軸突運輸?shù)墓δ苄缘鞍拙酆衔?，神?jīng)絲蛋白質(zhì)的運動研究對于如神經(jīng)退行性疾病的診斷等應(yīng)用是非常重要的。傳統(tǒng)的方法在很大程度上依賴于在熒光顯微鏡圖像下手工標(biāo)記神經(jīng)絲蛋白質(zhì)。描述了一種基于檢測跟蹤的自動神經(jīng)絲運動分析方法，用這種方法提取出沿軸突運動的神經(jīng)絲軌跡并描繪成一條參數(shù)化曲線。首先,將軸突分解成塊，然后利用馬爾可夫隨機場圖形標(biāo)簽來確定包含運動神經(jīng)絲的軸突塊，最后將神經(jīng)絲的首端和尾端位置細(xì)化到亞像素精度。神經(jīng)絲運動的實際延時熒光圖像序列實驗表明了所提方法的有效性和可靠性。

馬爾科夫隨機場(MRF);圖形切割;熒光顯微鏡

1 引言

軸突是一個拉長的細(xì)胞結(jié)構(gòu)，該結(jié)構(gòu)將一個神經(jīng)細(xì)胞延伸至其他的神經(jīng)細(xì)胞或器官來完成信息傳遞。它是一個細(xì)長的圓柱形結(jié)構(gòu)，往往直徑只有幾微米，但長度可達毫米、厘米甚至米。神經(jīng)絲是非常復(fù)雜的蛋白質(zhì)聚合物，直徑約10納米，長約幾十微米。它在軸突上的運動對軸突的生長和存活至關(guān)重要[1,2]。神經(jīng)絲平行于軸突的長軸以雙向方式隨機運動，如圖1所示。圖1a～圖1f為一個神經(jīng)絲蛋白質(zhì)從左向右運動的軌跡，圖1g～圖1k為視頻3中一個神經(jīng)絲蛋白質(zhì)從左向右運動的軌跡，并且在第j幀出現(xiàn)短時間的暫停。比例尺長度等于5微米。研究神經(jīng)絲的運動在神經(jīng)學(xué)領(lǐng)域中具有重要意義。

Figure 1 Two examples of neurofilament movement in axons of cultured cortical neurons

由于通過熒光顯微鏡可以直接觀察活細(xì)胞中的移動結(jié)構(gòu)，所以它被普遍應(yīng)用在神經(jīng)絲的運動研究中[3，4]。為了分析神經(jīng)絲的移動動力性能，我們必須在連續(xù)幀的熒光顯微鏡圖像上跟蹤神經(jīng)絲。然而，在目前許多的神經(jīng)學(xué)研究中，大多采用手工標(biāo)記方法來跟蹤神經(jīng)絲的運動，這種手動跟蹤不但勞動強度大，而且容易出現(xiàn)人為錯誤。總的來說，對神經(jīng)絲及細(xì)胞內(nèi)其它物質(zhì)的自動跟蹤方法會給生物醫(yī)學(xué)研究帶來越來越多利益，相關(guān)的工作內(nèi)容在接下來的一節(jié)中進行探討。

視覺跟蹤技術(shù)已被廣泛應(yīng)用于細(xì)胞內(nèi)各種物質(zhì)的運動跟蹤。 Smal I等人[5]將粒子濾波應(yīng)用于微管(Microtubules)的生長動態(tài)熒光顯微鏡圖像分析。在前期工作中[6]，我們研究了軸突中神經(jīng)絲運動的空間約束粒子濾波跟蹤方法。軸突的路徑被加到粒子濾波的動態(tài)模型中，作為先驗信息來限制神經(jīng)絲運動的方向和位置，不但顯著降低了粒子濾波方法的計算成本，同時也提高了精度。 Sargin M E[7]和Koulgi P[8]等人在隱馬爾可夫模型HMM的基礎(chǔ)上描述了不同的方法來跟蹤微管的伸長、縮短和滑移。雖然我們所提出的方法也是基于概率推理的圖形化模型，但是與這些現(xiàn)有工作的不同點在于：

(1)在馬爾可夫隨機場的基礎(chǔ)上，我們有效地利用圖形切割技術(shù)來解決問題;

(2)我們將軸突分解成塊來約束軸突，從而大大降低了研究的復(fù)雜性;

(3)方法中的端點細(xì)化提高了神經(jīng)絲的跟蹤精度。

在近幾年中，另一個值得注意的發(fā)展是使用kymographs方法來跟蹤圖像序列中的細(xì)胞內(nèi)物質(zhì)[9～13]。這種方法是將物體的像素強度投影到其運動路徑上，從而建立一個線性強度輪廓。在視頻的每幀中重復(fù)這一過程，然后將所得的線性強度輪廓堆疊就可以創(chuàng)建一個二維圖像，其中一個軸代表沿軸突的距離，而另一個軸代表時間。使用這種方法，沿著軸突移動的物體在kymograph圖中創(chuàng)建斜軌跡。這時物體的移動速度就可以通過測量軌跡的斜率獲得。然而，當(dāng)物體的移動速度相對于視頻幀的速率過快時，這種方法可能會出現(xiàn)問題，因為在神經(jīng)絲運動視頻中，視頻幀之間必須使用相對長的時間間隔(4或5 s)以盡量減少這些微弱的衍射極限結(jié)構(gòu)的圖像漂白現(xiàn)象，但是這種長時間間隔會使kymograph圖中產(chǎn)生不連續(xù)的物體運動軌跡。

本文提出了一種新的基于檢測跟蹤的神經(jīng)絲運動分析方法。首先，將圖像中的軸突分解成多個塊，然后利用馬爾可夫隨機場MRF(Markov Random Field)對神經(jīng)絲進行檢測。為了精確定位神經(jīng)絲的首尾端，必須進一步提高檢測精度。于是我們將檢測到的神經(jīng)絲對應(yīng)到相關(guān)聯(lián)的連續(xù)幀中，形成連續(xù)的軌跡。馬爾可夫隨機場的圖像處理技術(shù)被廣泛應(yīng)用于計算機視覺領(lǐng)域中，本論文首次嘗試將該技術(shù)應(yīng)用于神經(jīng)絲運動分析。

本文的主要工作如下：

(1)充分利用軸突的幾何形狀來分析神經(jīng)絲的運動。這樣，只有軸突周圍的圖像區(qū)域才涉及計算。

(2)檢測神經(jīng)絲時，首先進行塊的神經(jīng)絲檢測，這樣不但保證其空間上的連續(xù)性。同時，相比于對每個像素操作，大大降低了計算時間。這種檢測可以轉(zhuǎn)化成二進制圖形標(biāo)記問題，并完全可以使用基于多項式插值的圖形切割技術(shù)來解決。

(3)端點細(xì)化主要是把首尾端的精度從塊提高到亞像素級別，并提高處理熒光成像的光學(xué)模糊的準(zhǔn)確性。

2 我們的跟蹤方法

2.1 軸突模型和神經(jīng)絲檢測

我們首先將軸突分解成多個同樣大小的塊。神經(jīng)細(xì)胞的軸突沿一條平滑的曲線運動，如圖2a所示。在本文中，和我們以前做的工作[6]一樣，我們使用三次樣條插值法獲得近似于軸突路徑的分段多項式曲線。然后，沿軸突生成具有固定寬度和高度的矩形塊，如圖2b所示。

Figure 2 Example of axon/neurofilament fluorescence microscopy image and the detection result based on block decomposition

將軸突分解成多個塊之后，神經(jīng)絲的檢測就是找到其中包含神經(jīng)絲的軸突塊。我們把含有神經(jīng)絲的塊標(biāo)記為1，否則標(biāo)記為0。將沿著軸突的所有塊的標(biāo)簽組態(tài)表示為x={x0,x1,…,xn-1}，其中xi從二進制標(biāo)簽集L={0,1}中取值。根據(jù)圖像觀察結(jié)果y={y0,y1,…,yn-1}，那么一個標(biāo)簽組態(tài)x的后驗概率為P(x|y)，從而具有最大后驗概率標(biāo)準(zhǔn)的最佳標(biāo)簽組態(tài)就可以推斷為X*=arg maxP(x|y)。

考慮到一維結(jié)構(gòu)的軸突塊和相鄰塊標(biāo)簽之間的空間相干性，我們認(rèn)為一個特定塊的標(biāo)簽取決于它的本地圖像觀察結(jié)果yi以及直接相鄰的標(biāo)簽。變量之間的相關(guān)性滿足馬爾可夫概率，相應(yīng)的無向概率圖是一個馬爾可夫隨機場。在這種情況下，可以寫成以下形式的后驗概率：

(1)

(2)

一元勢函數(shù)對本地觀察結(jié)果進行了建模，而二元勢函數(shù)則考慮相鄰塊之間的相互作用。通過獨立考慮每個數(shù)據(jù)塊的塊標(biāo)簽，一元勢函數(shù)取為塊標(biāo)簽概率對數(shù)的相反數(shù)：

(3)

(4)

其中θ是通過試驗得到的參數(shù)向量。圖像強度的平均值和標(biāo)準(zhǔn)偏差，以及每個塊的圖像強度對比組成了分類的特征向量，而且正、負(fù)樣本塊，都是從訓(xùn)練分類器的手動標(biāo)記圖像中選擇。

考慮到存在的空間連續(xù)性問題，我們建議相鄰塊在二元勢函數(shù)下采取相同的標(biāo)簽，除非圖形觀測有很大差異。我們事先定義了一個依賴于相似性的平滑度，如式(5)所示：

(5)

其中,hi和hj的強度直方圖是由圖像觀察結(jié)果yi和yj得到的。b(·)是兩個直方圖之間的Bhattacharyya系數(shù)，λ參數(shù)用于控制平滑度。當(dāng)兩個直方圖完全不同時，b(·)的值為0，相同時則為1。因此，這種情況不利于對外觀相似的相鄰塊分配不同的標(biāo)簽。

應(yīng)當(dāng)指出的是，在我們的問題中標(biāo)簽集是二進制的，上面的公式也滿足二元勢函數(shù)子模塊的要求[14]，即:

只要在一個有向圖中找到最小S-T切割(最大流-最小切割)，公式(2)中定義的能量E(x；y)可以很容易地進行最小化，這個定向圖的節(jié)點為xis，邊緣容量取決于一元和二元勢函數(shù)。由于版面的限制，我們不詳細(xì)討論s-t線圖的構(gòu)建和最小s-t切割問題的解決，由讀者自行參考文獻[15]。

2.2 端點細(xì)化

通過對軸突塊進行檢測可以確定神經(jīng)絲的大致位置。然而，位置的精確度受到軸突塊尺寸的限制。對于神經(jīng)絲的精確跟蹤來說，這樣的首尾端位置精度還不夠。此外，由于神經(jīng)絲微弱的衍射極限結(jié)構(gòu)的圖像漂白，在捕獲的熒光顯微鏡圖像中，首尾端的位置可能會出現(xiàn)模糊。在這兩種情況下，需要一個額外的步驟來完善端點的位置。

假設(shè)從軸突到神經(jīng)絲的圖像強度變化是一個平穩(wěn)的過渡，我們使用廣義Logistic曲線[16]，也被稱為理查茲曲線，對沿軸突方向的過渡邊界附近進行強度變化建模：

(6)

其中，A、K、B、ν、Q和M是數(shù)據(jù)擬合的曲線參數(shù)。

理查茲曲線的形式在每個側(cè)面曲線允許非對稱增長率，并適合我們研究的問題的強度轉(zhuǎn)變規(guī)律。從檢測步驟中挑選出邊界塊，用kernel法對沿軸突的強度進行平滑化計算，以減輕不規(guī)則的強度波動和影像噪聲。擬合過程如圖3所示。其強度變化率最大的點，即擬合曲線的拐點，被選擇作為神經(jīng)絲的端點。對于Richards的曲線，拐點的計算公式如式(7)所示：

(7)

Figure 3 Endpoint refinement. Solid and dashed boxes represent the neurofilament and background blocks, respectively圖3 端點細(xì)化。實線和虛線方框分別代表神經(jīng)絲和背景塊(僅示出一端)

2.3 關(guān)聯(lián)檢測

在各個幀中檢測到的神經(jīng)絲需要與在連續(xù)幀中檢測到的一一對應(yīng)并關(guān)聯(lián)在一起，以便獲得連續(xù)的軌跡。這是多目標(biāo)跟蹤技術(shù)中常用的數(shù)據(jù)關(guān)聯(lián)步驟。由于神經(jīng)絲的圖像外觀上很相似，我們必須依靠神經(jīng)絲的長度和位置信息來進行關(guān)聯(lián)。一個重要的貪婪(Greedy)關(guān)聯(lián)步驟是將一個神經(jīng)絲的檢測圖像和下一幀中距離最小、長度最相似的檢測圖像相關(guān)聯(lián)，由于受到約束，任何兩個神經(jīng)絲的運動不能彼此重疊。 基于神經(jīng)絲運動的1-D性質(zhì)，這樣的關(guān)聯(lián)方案滿足我們的應(yīng)用需求。

3 實驗結(jié)果

我們測試了基于顯微熒光圖像序列的跟蹤方法的性能，本節(jié)中將討論其中三個序列的實驗結(jié)果。實驗設(shè)置如下：所有的實驗序列都是從小鼠大腦皮質(zhì)分離的神經(jīng)元中記錄下來的，該項工作由美國俄亥俄州立大學(xué)AnthonyBrown教授在實驗室中完成。從顯微鏡測得的原始輸出圖像為512×512像素，倍率為0.131μm/像素，并且通過裁剪去除了不相關(guān)區(qū)域。我們的方法是，先由基于圖形工具的Boykov算法[15]來完成最低的s-t切割算法，然后通過Python/numpy和最低的s-t切割算法來實現(xiàn)。視頻序列1和視頻序列2的軸突塊的大小是10個像素，視頻序列3是5個像素，平滑度的先驗因素λ為2。跟蹤實驗使用Python語言在2.1GHz的英特爾酷睿2雙核計算機上運行。

圖4顯示出視頻序列1和視頻序列2的圖像幀和跟蹤結(jié)果。視頻序列1中，只有中間的神經(jīng)絲在運動，所以我們只繪制了這個圖像來進行闡述。在圖中首尾端用十字線標(biāo)記，兩個序列的定性跟蹤結(jié)果十分準(zhǔn)確。

Figure 4 Neurofilament fluorescence microscopy image sequences and the tracking results

為了定量評價跟蹤性能，我們將用我們的方法得到的神經(jīng)絲速度和有經(jīng)驗的用戶使用MetaMorph軟件人工標(biāo)記獲得的速度進行對比。得知：人工標(biāo)記的方法適用于視頻序列2。用我們的方法計算神經(jīng)絲速度的步驟如下：首先計算出在連續(xù)幀之間的神經(jīng)絲首端移動的像素距離，然后基于0.131μm/s的放大系數(shù)和5s的幀間隔條件轉(zhuǎn)化得到所求速度。計算結(jié)果表明該速度非常接近手動標(biāo)記獲得的速度。如圖5所示。

Figure 5 Comparison of the neurofilament velocities between manual labeled method and our method for sequence 2

為了進一步驗證基于檢測的跟蹤方法在跟蹤精確度和效率方面的優(yōu)勢，我們將已發(fā)表的基于粒子濾波的神經(jīng)絲的跟蹤結(jié)果(具體算法和實驗參看文獻[6])與基于檢測的跟蹤方法進行對比，圖6顯示了粒子濾波的跟蹤結(jié)果[6]。從圖6和圖5的對比可以看出，基于檢測的跟蹤方法在跟蹤的精確度上明顯優(yōu)于粒子濾波的方法，尤其在速度變化較大的地方。在效率方面，上述基于粒子濾波的方法完成對圖像序列2的跟蹤大致需要140s，而基于檢測的方法完成同樣的跟蹤大致只需要100s。因此，我們可以看出，基于檢測跟蹤的方法在自動跟蹤神經(jīng)絲中，無論在準(zhǔn)確性還是效率方面都優(yōu)于基于粒子濾波的方法。

Figure 6 Comparison of the neurofilament velocities between manual labeled method and particle filtering for sequence 2

對一個具有挑戰(zhàn)性的快速運動的短神經(jīng)絲圖像序列，如圖7所示,我們用kymograph方法和我們自己的方法進行了對比分析。我們發(fā)現(xiàn)，由于短神經(jīng)絲的快速運動，kymograph軌跡圖中有相當(dāng)多的不連續(xù)，如圖7a所示。這對于應(yīng)用kymograph方法的跟蹤技術(shù)提出了一個挑戰(zhàn)，因為它們需要連續(xù)強度來恢復(fù)軌跡。相比較而言，對于相同的圖像序列，用我們的方法獲得的結(jié)果更令人滿意，如圖7b所示。

Figure 7 Comparison of the results of kymograph and our method on sequence 3

4 結(jié)束語

在本文中，我們提出了一種新的基于檢測跟蹤的神經(jīng)絲運動分析方法。這種方法利用馬爾可夫隨機場建模方法和圖形切割技術(shù)有效地解決了檢測問題。將軸突的限制、軸突塊的檢測和端點細(xì)化相結(jié)合，有效地保證了效率和準(zhǔn)確性。熒光顯微鏡圖像中的神經(jīng)絲運動實驗的結(jié)果表明，我們提出的方法十分理想。

[1]BrayD.Cellmovements:Frommoleculestomotility[M]. 2edition.NewYork:GarlandScience, 2000.

[2]BrownA.Axonaltransportofmembranousandnonmembranouscargoes[J].JournalofCellBiology, 2003, 160(6):817-821.

[3]BrownA.Live-cellimagingofslowaxonaltransportinculturedneurons[J].MethodsinCellBiology, 2003, 71:305-323.

[4]GoldmanRD,SpectorDL.Livecellimaging:Alaboratorymanual[M]. 2edition.NewYork:ColdSpringHarborLaboratoryPress, 2009.

[5]SmalI,DraegesteinK,GaljartN,etal.Particlefilteringformultipleobjecttrackingindynamicfluorescencemicroscopyimages:Applicationtomicrotubulegrowthanalysis[J].IEEETransactionsonMedicalImaging, 2008, 27(6):789-804.

[6]YuanL,ZhengYF,ZhuJ,etal.Objecttrackingwithparticlefilteringinfluorescencemicroscopyimages:Applicationtothemotionofneurofilamentsinaxons[J].IEEETransactionsonMedicalImaging, 2012, 31(1):117-130.

[7]SarginME,AltinokA,RoseK,etal.Deformabletrellis:Opencontourtrackinginbioimagesequences[C]∥ProcofIEEEInternationalConferenceonAcoustics,Speech,andSignalProcessing, 2008:561-564.

[8]KoulgiP,SarginME,RoseK,etal.Graphicalmodel-basedtrackingofcurvilinearstructuresinbio-imagesequences[C]∥ProcofInternationalConferenceonPatternRecognition, 2010:2596-2599.

[9]ZhangK,OsakadaY,XieW,etal.Automatedimageanalysisfortrackingcargotransportinaxons[J].MicroscopyResearchandTechnique, 2010, 74(7):605-613.

[10]WelzelO,BoeningD,StroebelA,etal.Determinationofaxonaltransportvelocitiesviaimagecross-andautocorrelation[J].EuropeanBiophysics, 2009, 38(7):883-889.

[11]StepanovaT,SmalI,vanHarenJ,etal.History-dependentcatastrophesregulateaxonalmicrotubulebehavior[J].CurrentBiology, 2010, 20(11):1023-1028.

[12]SmalI,GrigorievI,AkhmanovaA,etal.Microtubuledynamicsanalysisusingkymographsandvariable-rateparticlefilters[J].IEEETransactionsonImageProcessing, 2010, 19(7):1861-1876.

[13]RacineV,SachseM,SalameroJ,etal.Visualizationandquantificationofvesicletraffickingonathree-dimensionalcytoskeletonnetworkinlivingcells[J].JournalofMicroscopy, 2007, 225(3):214-228.

[14]KolmogorovV,ZabihR.Whatenergyfunctionscanbeminimizedviagraphcuts? [J].IEEETransactionsonPatternAnalysisandMachineIntelligence, 2004, 26(2):147-159.

[15]BoykovY,KolmogorovV.Anexperimentalcomparisonofmin-cut/max-flowalgorithmsforenergyminimizationinvision[J].IEEETransactionsonPatternAnalysisandMachineIntelligence, 2004, 26(9):1124-1137.

[16]RichardsFJ.Aflexiblegrowthfunctionforempiricaluse[J].JournalofExperimentalBotany, 1959, 10(2):290-301.

YUAN Liang,born in 1972,PhD,associate professor,his research interests include computer vision, and image processing.

Neurofilaments tracking by detection in fluorescence microscopy images

YUAN Liang1,ZHU Jun-da2

(1.School of Mechanical Engineering,Xinjiang University,Urumqi 830047;2.Department of Electrical and Computer Engineering,University of Macau,Macau 519000,China)

Neurofilaments are functional protein polymers that are transported along the axonal processes of nerve cells. Studying the movement of neurofilaments is important in applications such as diagnosing neurodegenerative diseases. Traditional methods largely rely on manual labeling of the neurofilaments in fluorescence microscopy images. An automated method for analyzing the motion of neurofilaments based on tracking-by-detection is proposed. The axon along which the neurofilaments move is extracted and represented by a parametric curve. Firstly, the axon is decomposed into blocks. Secondly, the blocks containing the moving neurofilaments are determined by graph labeling using Markov Random Field. Finally, the leading and trailing locations of a neurofilament are refined to sub-pixel accuracy. Experiments on real time-lapse fluorescence image sequences of neurofilament movement demonstrate the efficacy and efficiency of our proposed method.

Markov Random Field (MRF);graph cut;fluorescence microscopy

1007-130X(2015)01-0119-06

2013-04-24;

2013-09-26基金項目：國家自然科學(xué)基金資助項目(31460248,61262059)；教育部留學(xué)回國人員科研啟動基金資助項目；新疆優(yōu)秀青年科技創(chuàng)新人才培養(yǎng)項目(2013721016);新疆大學(xué)博士啟動基金資助項目

TP391.42

A

10.3969/j.issn.1007-130X.2015.01.018

袁亮(1972-),男,新疆人，博士，副教授，研究方向為計算機視覺與圖像處理。E-mail:ylhap@163.com

通信地址：830047 新疆烏魯木齊市延安路1230號新疆大學(xué)機械工程學(xué)院

Address:School of Mechanical Engineering,Xinjiang University,1230 Yan’an Rd,Urumqi 830047,Xinjiang,P.R.China