当前，细胞与蛋白质的富集、迁移、分离和输运成为生命科学和生物制药的一个热门研究领域^[1].传统的微粒控制有光学镊子、电泳、磁泳和微流泵^[2].相比这些微控制方法，介电泳技术具有无创性和操作便捷的优点^[3].随着介电泳技术的不断发展，如今已有光/电两门学科融合迸发的新颖技术——光诱导介电泳^{[4, 5]}.光诱导介电泳技术现已被应用于细胞富集、细胞计数器和细胞形态组织控制.另一方面，光诱导介电泳技术所需的光斑能量较低(一般约为100 nW)，光电薄膜的制作工艺也相对简单.本文提出光诱导介电泳力操控椭球细胞的数值模型，详细说明了椭球细胞受正介电泳力发生定向移动和自旋.仿真结果与文献[6, 7]的结果一致.

光诱导生物芯片主要采用一种光敏物质，其结构与太阳能电池类似，即俗称的三明治结构.光诱导生物芯片的上、下两块玻璃表面涂覆一层纳米量级厚度的氧化铟锡(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层后，氢化非晶硅物质内部会发生光生载流子现象，其宏观表征是电导率增大.较大的电势落于溶液层，为微粒电极化提供必要条件.由于入射光斑的尺度较小会使芯片内的电场分布具有非均匀特征，因此，微粒受非均匀电场引起的介电泳力作用发生定向移动.

光诱导介电泳力产生方式与传统金属介电泳力不同，但两者的物理本质一致.致力于形成一个强的非均匀电场迫使电中性微粒发生电偶极化，所以有学者也将诱导光斑形象地比喻为虚拟电极.光诱导介电泳力的表达式也与经典Pohl介电泳力一致：

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

椭球的电偶极化强度可表示为

简化公式(4)，本文以k_i代替式(4)中的分母式，i代表空间矢量的x，y，z分量.

根据流体动力学，运动的物体移动过程中也承受阻力.椭球在溶液中的阻力被称为斯托克斯阻力，其表达式为

最后根据牛顿第二定律，即公式(6)，可求解椭球的运动轨迹和速度.为了避免公式的冗长，公式(6)仅表述椭球x分量的运动微分方程.y，z分量可根据轮换理论方便得到.

在介电泳力的细胞控制实验里，细胞移动过程出现转动是一件比较常见的现象.究其原因是极化细胞的电偶极矩与空间电场指向具有一定差角.空间电场将迫使细胞发生自转，进而与电场方向平行一致.电偶极矩在电场中的扭矩可以由两者的矢量差积求得:

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

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

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

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

本文构建椭球在光电芯片横截面的二维数值模型.模型的电场分布采用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)—两类椭球粒子的运动轨迹模拟.

1) 基于Runge-Kutta方法计算了两类椭球粒子随时间变化的自旋角度，长宽比越大的椭球其自转速度越快.

2) 根据有限元COMSOL的电场分布数值，采用Velocity-Verlet算法模拟了椭球在光电芯片内部受正、负介电泳力的运动轨迹，且仿真结果与光诱导介电泳实验观测具有较高的相似性.