基于檢測的熒光顯微圖像中的神經絲蛋白質跟蹤*

袁 亮，朱俊達

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

基于檢測的熒光顯微圖像中的神經絲蛋白質跟蹤*

袁 亮1，朱俊達2

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

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

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

1 引言

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

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

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

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

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

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

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

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

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

本文的主要工作如下：

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

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

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

2 我們的跟蹤方法

2.1 軸突模型和神經絲檢測

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

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

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

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

(1)

(2)

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

(3)

(4)

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

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

(5)

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

應當指出的是，在我們的問題中標簽集是二進制的，上面的公式也滿足二元勢函數子模塊的要求[14]，即:

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

2.2 端點細化

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

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

(6)

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

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

(7)

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

2.3 關聯檢測

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

3 實驗結果

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

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

Figure 4 Neurofilament fluorescence microscopy image sequences and the tracking results

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

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

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

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

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

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

4 結束語

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

[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基金項目：國家自然科學基金資助項目(31460248,61262059)；教育部留學回國人員科研啟動基金資助項目；新疆優秀青年科技創新人才培養項目(2013721016);新疆大學博士啟動基金資助項目

TP391.42

A

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

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

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

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