基于检测的荧光显微图像中的神经丝蛋白质跟踪*

袁 亮，朱俊达

(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