胡晟^1,2, 赵勇^1,2, 胡海峰¹
1. 东北大学信息科学与工程学院, 辽宁 沈阳 110819;
2. 东北大学流程工业综合自动化国家重点实验室, 辽宁 沈阳 110819
收稿日期: 2015-03-03
基金项目: 国家自然科学基金资助项目(61273059, 61203206); 流程工业综合自动化国家重点实验室基础科研业务费资助项目(2013ZCX09).
作者简介: 胡晟(1984-),男,云南景洪人,东北大学博士研究生; 赵勇(1973-),男,辽宁沈阳人,东北大学教授,博士生导师.

摘要: 基于分子动力学提出了一种光诱导介电泳控制椭球粒子运动的数值模型.研究了光电芯片中椭球粒子承受的光诱导介电泳和斯托克斯阻力，采用Runge-Kutta方法计算不同长宽比粒子的自转速度.采用COMSOL有限元计算电场，借助Velocity-Verlet算法模拟了粒子受介电泳的运动轨迹.仿真结果表明，粒子长宽比越大，转速越快；在28, 30 μm位置处的椭球粒子，受正介电泳力向光斑运动，且沿着电场强度梯度方向行进，最高速度可达到312 μm/s.以上仿真的转动速度和运动轨迹都与实验保持了较好的一致性.
关键词: 光诱导介电泳光电镊子Runge-Kutta方法Velocity-Verlet方法椭球粒子
Numerical Model for Ellipsoid Experiencing Optically Induced Dielectrophoresis
HU Sheng^{1, 2}, ZHAO Yong^{1, 2}, HU Hai-feng^{1, 2}
1. School of Information Science & Engineering, Northeastern University, Shenyang 110819, China;
2.State Key Laboratory of Synthetical Automation for Process Industries, Northeastern University, Shenyang 110819, China
Corresponding author: ZHAO Yong, professor, E-mail: zhaoyong@ise.neu.edu.cn
Abstract: Based on the method of molecular dynamics, a numerical model was presented to simulate the behavior of such a ellipsoid experiencing optically induced dielectrophoresis (ODEP). The ODEP and Stokes force exerted on the ellipsoid were studied, and the Runge-Kutta method was used to calculate the self-rotation speed regarding different aspect ratio of ellipsoid. The electric field was solved by using COMSOL software, and the trajectory of the ellipsoid experiencing ODEP was simulated on the basis of the Velocity-Verlet algorithm. The simulation results indicated that the greater particle aspect ratio, the faster the speed. The ellipsoid at 28, 30 μm location was controlled by the positive ODEP, which could move towards the light pattern along the gradient of the electric field magnitude, and the maximum speed could reach 312 μm/s. The simulations results of the rotation speed and the trajectory are consistent with those in the experiment.
Key words: optically induced dielectrophoresisoptoelectronic tweezerRunge-Kutta methodVelocity-Verlet methodellipsoid 当前，细胞与蛋白质的富集、迁移、分离和输运成为生命科学和生物制药的一个热门研究领域^[1].传统的微粒控制有光学镊子、电泳、磁泳和微流泵^[2].相比这些微控制方法，介电泳技术具有无创性和操作便捷的优点^[3].随着介电泳技术的不断发展，如今已有光/电两门学科融合迸发的新颖技术——光诱导介电泳^{[4, 5]}.光诱导介电泳技术现已被应用于细胞富集、细胞计数器和细胞形态组织控制.另一方面，光诱导介电泳技术所需的光斑能量较低(一般约为100 nW)，光电薄膜的制作工艺也相对简单.本文提出光诱导介电泳力操控椭球细胞的数值模型，详细说明了椭球细胞受正介电泳力发生定向移动和自旋.仿真结果与文献[6, 7]的结果一致.
1 光诱导介电泳力电泳生物芯片的基本结构光诱导生物芯片主要采用一种光敏物质，其结构与太阳能电池类似，即俗称的三明治结构.光诱导生物芯片的上、下两块玻璃表面涂覆一层纳米量级厚度的氧化铟锡(ITO,indium tin oxide)材质，氧化铟锡玻璃除了有良好的透光性，还具备导电性.下表面的氧化铟锡玻璃采用等离子体增强化学气相沉淀法(PECVD,plasm enhanced chemical vapor deposition)镀上一层氢化非晶硅(a-Si:H)，该物质是一种光敏材料.a-Si:H光敏薄膜是由PECVD方法连续沉积50 nm厚的n+型氢化非晶硅(n+a-Si:H)层，1.5 μm厚的本征态氢化非晶硅(intrinsic a-Si:H)层和25 nm厚的碳化硅(SiC_x)层组成.光诱导介电泳生物芯片的主要结构如图 1所示.
图1(Fig. 1)

图1 光诱导介电泳生物芯片的结构示意图 Fig. 1 Structure schematic of optically induced dielectrophoresis biochip(a)—基本装置和粒子控制原理； (b)—截面图.

针对光诱导介电泳力控制悬浮细胞的工作原理，可以简单叙述如下.首先，用金属夹子将上、下两块ITO玻璃和信号发生器相连，信号发生器用于提供一定幅值和频率的电压信号.最后通过微米尺度的诱导光斑照射生物芯片表面，即可快速实现微粒子控制.空间微米光斑的产生有很多手段，目前国内外一般采用光学几何透镜、数字微镜器件(DMD,digital mirror display)，LED或LCD设备^[8].未有入射光斑照射时，氢化非晶硅具有较低的电导率，使大量的电场能量集中在a-Si:H层.根据电路串联理论可知，微粒溶液分获的电压十分微弱，因为该电场不足以引起粒子的极化运动.相反，当入射光照射到a-Si:H层后，氢化非晶硅物质内部会发生光生载流子现象，其宏观表征是电导率增大.较大的电势落于溶液层，为微粒电极化提供必要条件.由于入射光斑的尺度较小会使芯片内的电场分布具有非均匀特征，因此，微粒受非均匀电场引起的介电泳力作用发生定向移动.
2 椭球粒子介电泳基础理论2.1 椭球粒子的介电泳力光诱导介电泳力产生方式与传统金属介电泳力不同，但两者的物理本质一致.致力于形成一个强的非均匀电场迫使电中性微粒发生电偶极化，所以有学者也将诱导光斑形象地比喻为虚拟电极.光诱导介电泳力的表达式也与经典Pohl介电泳力一致：

式中：F_DEP，p和E分别为介电泳力、电偶极化强度和电场强度.介电泳力的矩阵形式为
式中，φ表示电势.椭球与规则球形颗粒的电偶极化强度不同.因为椭球粒子诱导电偶极矩方向一般不平行于空间电场方向，所以椭球粒子的介电泳力不能等效于球形表达式：，即p≠αE.另外，真实的生物细胞形状并非规则球体，介电泳理论也指出球形建模的细胞运动与实验观测结果具有一定的误差^[6].因此椭球的电偶极化强度应重新详细描述.式(3)首先说明椭球特征常数A_i和A，它们反应椭球的形态特征，式(3)可根据椭圆积分求解：

式中l_i是椭球在x，y，z三个分量的轴长(l_i=a,b或c).椭球的电偶极化强度可表示为

简化公式(4)，本文以k_i代替式(4)中的分母式，i代表空间矢量的x，y，z分量.
根据流体动力学，运动的物体移动过程中也承受阻力.椭球在溶液中的阻力被称为斯托克斯阻力，其表达式为

式中：η是液体的黏滞系数，本文仿真假定η为10^-3 Pa·s；ν_i是粒子在i轴的分速度.最后根据牛顿第二定律，即公式(6)，可求解椭球的运动轨迹和速度.为了避免公式的冗长，公式(6)仅表述椭球x分量的运动微分方程.y，z分量可根据轮换理论方便得到.

式中m为椭球的质量(m=4πρabc/3，ρ为仿真粒子的密度，传统生物细胞的密度近似于水溶液，因此本文的ρ=1 000 kg/m³).2.2 椭球粒子的扭矩在介电泳力的细胞控制实验里，细胞移动过程出现转动是一件比较常见的现象.究其原因是极化细胞的电偶极矩与空间电场指向具有一定差角.空间电场将迫使细胞发生自转，进而与电场方向平行一致.电偶极矩在电场中的扭矩可以由两者的矢量差积求得:

与上述介电泳力一样，扭矩T_DEP的矩阵形式改写为

同理椭球以z坐标轴旋转承受的斯托克斯阻力为

椭球在x轴处的转动惯量T_MI-x为

最终椭球粒子以x轴为法线在z-y平面上自转，其角度的常微分方程为

3 数值模拟和分析本文构建椭球在光电芯片横截面的二维数值模型.模型的电场分布采用COMSOL有限元软件求解，电场服从Laplace偏微分方程，虚拟电极边界条件上的电压值为5 V，属于第一类边界条件，虚拟电极电压服从高斯分布，高斯半高宽等于11.8 μm.而芯片两侧服从绝缘边界属于第二类边界条件.关于具体的边界条件设定，以及氢化非晶硅的光亮/光暗电导率可参考文献[9, 10].本文依据公式(11)模拟细胞的自旋角度.初始假定2种不同的椭球粒子：椭球粒子Ⅰ与椭球粒子Ⅱ，两椭球粒子a，b，c参数分别为[a 0.1a a](椭球Ⅰ),[0.9a 0.5a a](椭球Ⅱ).a为5 μm，两椭球粒子如图 2所示.
图2(Fig. 2)

图2 两类不同轴长椭球粒子的三维视图 Fig. 2 Three-dimensional view of two ellipsoids for different axial length

公式(11)中的θ(t)常微分方程采用经典的Runge-Kutta方法，该方法在若干点处求解等式的导数，以组合得到平均斜率,最终提高求解阶数，相比欧拉方法它具有较高的精度.
采用4阶Runge-Kutta方法求解，仿真结果如图 3所示.另外本文将长宽比定义为粒子长轴与短轴的比值，因为研究z轴方向的粒子旋转可得两椭球的长宽比分别为10和1.8.初始假设粒子与z轴的初始夹角138°，随着时间的增加，椭球长轴逐渐与电场线平行，两者夹角保持为0°或180°.从图 3a的仿真结果可发现，椭球Ⅰ比椭球Ⅱ更快与电场同向，且本仿真与文献[6]的椭球自转模拟结果相符.在此简化模型的复杂性，并使模拟粒子运动更加直观，模拟了椭球运动的二维模型.
图3(Fig. 3)

图3 两类椭球粒子自旋角度仿真 Fig. 3 Rotational simulation of two ellipsoids (a)—椭球随时间变化的角度； (b)—椭球粒子在x-y平面旋转示意图.

基于公式(6)的粒子受力描述，构建椭球粒子受光诱导介电泳力运动的数值模型.椭球粒子的运动模拟采用分子动力学Velocity-Verlet算法，该算法是一种用于求解牛顿运动方程的数值方法.优点在于数值稳定性比简单的欧拉方法高很多，并保持了物理系统中的时间可逆性与相空间体积元体积守恒的性质，具体相关编程可参考文献[11].因为粒子的转动速度非常快，所以粒子轨迹模拟时简化转动力矩的微分求解，设定时间步长Δt为40 ms.根据上述椭球转动的结果，此时间下椭球极轴已经与电场同向.不失一般性，假定椭球Ⅰ和椭球Ⅱ的相对介电常数ε_pⅠ为85，ε_pⅡ为20，溶液的相对介电常数ε_m为80.它们分别受到正、负介电泳力.图 4a为光电芯片的电场分布，光斑位置处的电场线较密，场强较大.因此椭球Ⅰ在正介电泳力作用下，在10.2 s左右运动至光斑位置，如图 4b所示.椭球Ⅱ受负介电泳力，被排斥到光电芯片边缘.两椭球粒子的运动轨迹都是切向于电场能量方向运动，最终实现了椭球运动轨迹的仿真.
图4(Fig. 4)

图4 椭球粒子运动仿真 Fig. 4 Simulation of ellipsoids movement (a)—基于COMSOL有限元方法仿真的电场分布； (b)—两类椭球粒子的运动轨迹模拟.

4 结 论1) 基于Runge-Kutta方法计算了两类椭球粒子随时间变化的自旋角度，长宽比越大的椭球其自转速度越快.
2) 根据有限元COMSOL的电场分布数值，采用Velocity-Verlet算法模拟了椭球在光电芯片内部受正、负介电泳力的运动轨迹，且仿真结果与光诱导介电泳实验观测具有较高的相似性.

参考文献

[1]	Pethig R.Dielectrophoresis:status of the theory,technology,and applications[J].Biomicrofluidics,2010,4(2):1-35.(1)
[2]	Ohta A T,Chiou P Y,Phan H L,et al.Optically controlled cell discrimination and trapping using optoelectronic tweezers[J].IEEE Journal of Selected Topics in Quantum Electronics,2007,13(2):235-243.(1)
[3]	Gagnon Z R.Cellular dielectrophoresis:applications to the characterization,manipulation,separation and patterning of cells[J].Electrophoresis,2011,32(18):2466-2487.(1)
[4]	Ohta A T,Chiou P Y,Han T H,et al.Dynamic cell and microparticle control via optoelectronic tweezers [J]. Journal of Microelectromechanical Systems,2007,16(3):491-499.(1)
[5]	Liang W F,Zhao Y L,Liu L Q,et al.Rapid and label-free separation of burkitt′s lymphoma cells from red blood cells by optically-induced electrokinetics[J].PLoS One,2014,9(3):1-9.(1)
[6]	Winter W T,Welland M E.Dielectrophoresis of non-spherical particles[J].Journal of Physics D:Applied Physics,2009,42(4):1-6.(3)
[7]	Yang S M,Yu T M,Huang H P,et al.Light-driven manipulation of picobubbles on a titanium oxide phthalocyanine-based optoelectronic chip [J].Applied Physics Letters,2011,98:153512-3.(1)
[8]	Hwang H,Park J K.Optoelectrofluidic platforms for chemistry and biology[J].Lab on a Chip,2011,11(1):33-47.(1)
[9]	倪中华,易红,朱树存.基于光诱导介电泳的微纳米生物粒子操纵平台关键技术[J].中国科学E辑:技术科学,2009,39(10):1635-1642. (Ni Zhong-hua,Yi Hong,Zhu Shu-cun,et al.Research on critical technology of micro/nano bio-particles manipulation platform based on light-induced dielectrophoresis[J].Science in China Series E:Technological Sciences,2009,39(10):1635-1642.)(1)
[10]	朱晓璐,尹芝峰,高志强.基于光诱导介电泳的微粒子过滤、输运、富集和聚焦的实验研究[J].中国科学E辑:技术科学,2011,41(3):334-342. (Zhu Xiao-lu,Yin Zhi-feng,Gao Zhi-qiang,et al.Experimental study on filtering,transporting,concentrating and focusing of microparticles based on optically induced dielectrophoresis[J].Science in China Series E:Technological Sciences,2011,41(3):334-342.)(1)
[11]	Satoh A.Introduction to practice of molecular simulation [M].Tokyo:Elsevier,2010:1-321.(1)

1

本文献在全文中的定位：
... 富集、迁移、分离和输运成为生命科学和生物制药的一个热门研究领域^[1].传统的微粒控制有光学镊子、电泳、磁泳和微流 ...



1

本文献在全文中的定位：
... 的微粒控制有光学镊子、电泳、磁泳和微流泵^[2].相比这些微控制方法，介电泳技术具有无创性和操作 ...



1

本文献在全文中的定位：
... 这些微控制方法，介电泳技术具有无创性和操作便捷的优点^[3].随着介电泳技术的不断发展，如今已有光/电两门学科融合迸发的新颖 ...



1

本文献在全文中的定位：
... 断发展，如今已有光/电两门学科融合迸发的新颖技术——光诱导介电泳^{[4, 5]}.光诱导介电泳技术现已被应用于细 ...



1

本文献在全文中的定位：
... 迸发的新颖技术——光诱导介电泳^{[4, 5]}.光诱导介电泳技术现已被应用于细胞富集、细胞计数器和细胞形态组 ...



3

本文献在全文中的定位：
... ，详细说明了椭球细胞受正介电泳力发生定向移动和自旋.仿真结果与文献[6, 7]的结果一致. ...

... 介电泳理论也指出球形建模的细胞运动与实验观测结果具有一定的误差^[6].因此椭球的电偶极化强度应重新详细描述.式(3)首先说明椭球特征常 ...

... 真结果可发现，椭球Ⅰ比椭球Ⅱ更快与电场同向，且本仿真与文献[6]的椭球自转模拟结果相符.在此简化模型的复杂性，并使模拟粒子运动更加 ...



1

本文献在全文中的定位：
... 力发生定向移动和自旋.仿真结果与文献[6, 7]的结果一 ...



1

本文献在全文中的定位：
... 般采用光学几何透镜、数字微镜器件(DMD,digital mirror display)，LED或LCD设备^[8].未有入射光斑照射时，氢化非晶硅具有较低的电导率，使大量的电场 ...



1

本文献在全文中的定位：
... .关于具体的边界条件设定，以及氢化非晶硅的光亮/光暗电导率可参考文献[9, 10].本文依据公式(11)模拟细胞的自旋角度 ...



1

本文献在全文中的定位：
... 化非晶硅的光亮/光暗电导率可参考文献[9, 10].本文依据公式(11)模拟细胞的自旋角度.初始假定2种不同的椭球粒子：椭球 ...



1

本文献在全文中的定位：
... 中的时间可逆性与相空间体积元体积守恒的性质，具体相关编程可参考文献[11].因为粒子的转动速度非常快，所以粒子轨迹模拟时简化转动力矩的微分求 ...