自驱动杆状粒子在半柔性弹性环中的集体行为*

仲颖 施夏清

(苏州大学物理科学与技术学院,软凝聚态物理及交叉研究中心,苏州 215006)

在生物体系的活性系统中,杆状粒子在弹性半柔性边界中的受限行为极为常见.本文研究了二维情况下,自驱动杆状粒子受限在半柔性弹性环中的集体行为.改变系统的粒子数及噪声强度,系统显示明显的自驱吸附有序态、无序态及中间的过渡态.通过表征弹性环内部粒子的径向极性大小和空间分布的非球度性对这些状态进行了刻画.进一步对弹性环中心附近粒子密度的分析,发现环中心气态粒子分布存在一个与边界高密度区域共存的饱和平台,出现类似吸附转变的粒子分布.在过渡区间,体系内存在较大的涨落会导致弹性环出现异常形变.非对称的粒子分布对弹性环整体的迁移具有重要贡献,系统在过渡区间能获得相对较强的定向迁移.

1 引 言

活性物质通过消耗自身携带的或从环境中吸收的能量来实现各种力学运动,是一类典型的具有多级结构的非平衡体系,展现出丰富的集体动力学行为.不同于平衡态,活性物质常常具有很强的密度与取向序的耦合,产生极为复杂的时空演化行为[1−8].在一些大尺度的模型和理论研究中,活性物质在取向上的对称性破缺,会导致系统在空间上的密度不均匀分布,比如自驱动粒子在周期性边界中会自发形成一些大的相分离结构[9−12].除此之外,活性物质也会自发演化出特异性的形态,比如在细胞骨架纤维的实验中,可以看到肌动蛋白丝和微管,会形成一些涡旋、波形的图案分布[13].这其中粒子的形状会对系统的集体行为产生较大影响,从最简单的球形[9,14,15],到哑铃型双球[16,17]、蠕虫形粒子[18]以及杆状粒子[11,12,19−21].

生命系统中常常存在非常明显的空间受限条件.迁移的细胞个体内部就有复杂而有序的物质结构.细胞膜下覆盖的肌动蛋白微丝与分子马达的动态结构对于细胞迁移具有至关重要的作用.比如在肌动蛋白贡献细胞迁移[22,23]的过程中,具有活性物质特性的蛋白纤维在细胞膜附近生长并对其施加压力[24−26],所有过程都是在细胞膜内进行的.所以在活性物质,尤其是相关生命系统中,边界束缚条件或者界面作用是个值得理论研究注意的问题.

自驱动粒子和界面的相互作用会对内部自驱动粒子的集体运动行为产生很大影响[27].在类似于通道的开放的受限系统中,自驱动粒子会在边界附近形成大量集聚,并根据通道的具体特性形成一些对应分布[17,28,29].在封闭系统内部,自驱动粒子同样会在界面位置形成非常稳定的堆积分布[30,31].目前,对更有趣的软受限条件研究得比较少.已有的一些研究发现,在该条件下系统大尺度的粒子涨落会诱发柔性边界的明显形变,反过来又会影响粒子的均匀分布[32−34].

本文研究了大量自驱动杆状粒子在柔性边界中的动力学行为.简单的球形自驱动粒子模型,往往不能很好地体现活性物质中个体的形状各向异性.杆状粒子在非平衡受限时会有丰富的堆积行为[35].我们的模型使用了弹性杆体系,不仅可以实现柔性边界的弹性伸缩,还可以在保证杆特性的前提下,对杆的算法进行优化.我们用弹性杆链接而成的柔性环作为约束边界,统计了大数目自驱动杆状粒子受限在该柔性环中的分布情况和整体运动.在该柔性环提供的弹性约束下,其内部的自驱动杆状粒子同样也倾向于在环边界上聚集分布,但在中心位置会保持一个相对稳定的气态密度分布.大部分情况下这些自驱动杆状粒子都具有明显的对称性分布,但在有序无序转变区域,粒子大尺度的涨落与膜的大尺度形变两者相互影响,会形成各向异性的粒子分布,同时伴随较强的整体定向迁移.可以看到,这部分系统整体的运动以及形变,一定条件下发生明显极化,并伴随较大尺度的迁移,与一些细胞迁移运动行为类似,因而该系统的研究对调控这类集体运动状态行为有一定的参考作用.

2 模拟模型和方法

模拟一个二维平面下的系统,如图1(a)所示,该系统由两部分组成:其一是Nr个原长为Lr的自驱动杆状粒子;另一部分是一条柔性环,由Nl个原长为Ll的杆首尾链接而成,柔性环将所有自驱动杆状粒子约束在其内部.

图1 (a)系统组成的示意图,颜色代表杆身的取向;(b)杆间碰撞受力示意图Fig.1.(a)The schematic diagram of this system,and the rods are colored according to their angle with respect to the radial direction;(b)the interaction between rods.

在该模型中的自驱动杆和环的结构单元都是采用相同构造的杆,这种杆沿其径向两端是直径等宽的半圆形结构,同时所有的杆都具有固定的杆宽r0(r0=1),如图 1(b)所示.这样设定的杆可以比较方便地计算杆间最短距离,进而方便判断近邻杆间的碰撞及计算相互作用.

所有的这些杆状粒子,杆身受力不能弯曲.但沿杆身方向是弹性的,且只能沿径向弹性伸缩.如果使用固定杆长的硬杆链接而成的环,其中近邻的杆在受到碰撞时,需要很高精度的迭代保证最近邻杆间的运动刚性.而杆身弹性的设定可以有效简化这方面的计算,通过杆自身的弹性约束保持环中杆间的连续.此外我们选取较大的杆身弹性系数k,对于内部的自驱动杆状粒子,可以视为是固定杆长的硬杆.

为方便表述,约定所有杆两端的圆心位置分别为该杆的正负两端(正负端形状完全对称),杆的取向 θ 则定义为由其负端到正端的指向.每根杆的位置信息可由其质心位置 ri,取向 θi和实际的杆长 li确定.相应的杆的正负两端的位置分别为特别地,对于环中的相邻杆,规定两者相连处的异号端点位置是始终重合的,即环上杆 i的“ + ”端端点位置,与它顺时针方向的相邻杆 i +1 的“-”端端点位置重合,满足的关系[36].

杆之间的相互作用是截断的简谐弹簧势

其中r是两杆之间的最短距离,r0为固定的杆宽,φ0为 势能的强度大小.通过在 r=r0处的截断,体系为纯弹簧排斥势.在体积排斥效应下,发生碰撞的杆状粒子有平行排列倾向.

体系中所有的杆侧面都是光滑的,即碰撞时不考虑杆之间的滑动摩擦.在这种杆的模型中,碰撞时杆所受的排斥力垂直于杆身或作用在端点位置,如图1(b)所示.杆间的碰撞所受的相互作用可以等效到端点位置处.实际碰撞位置的排斥受力和端点上的等效受力,其合力和力矩存在如下的等量关系:

其中 Fi,j是碰撞时i杆身上所受的与j杆相互作用的排斥力,分别对应杆正负两端位置的等效受力,λ和λ0分别是碰撞受力位置和杆质心,沿杆取向 θ 在杆身上的相对位置,λ 在负端取0正端取1,杆质心的相对位置为 λ0=0.5 .方程(2)联立可解得碰撞时杆正负两端的等效受力分别为:

直接作用在端点处的碰撞受力同样符合上面的结果.

环内的自驱动粒子除了杆间的碰撞排斥作用外,沿着其取向 θ 还受恒定大小的自驱动力 sθ 持续牵引,从而具有自我推进的能力.

其中 μ 是迁移率,这里使用各向同性的迁移率.j∈ Ω 是 与当前杆发生碰撞的近邻杆.ξi(t)是高斯白噪声,定义噪声强度为 η .杆两端位置上的高斯白噪声可以给杆提供一个有效的角度上的扰动.

对于柔性环上的杆,由于不存在自驱动力其动力学方程为

这里也忽略了噪声的作用.特别地,柔性环上的相邻两杆之间不存在空间排斥作用,即 j ∈Ω′是除去环上相邻两杆后与当前杆发生碰撞的近邻杆.由于柔性环上相邻两杆端点重合,环上杆i正端端点的实际位移与它顺时针方向的相邻杆 i +1 的负端端点实际位移相同,即

在本文模拟中,将原长 Lr=2 自驱动杆受限于半柔性环内.半柔性环由 Nl=200,原长为Ll=1的杆首尾链接而成.将系统整体放在足够大的二维周期性边界内(Lx=Ly=200).选取足够小的时间步长 d t=0.0005τ,保证数值稳定性及必要的精度,其中 τ=1 是模拟的单位时间量程.根据动力学方程(4)式和(5)式对粒子的位置进行更新.经过足够长时间演化后,可对系统所处的稳定分布进行分析.本文主要研究系统在自驱动杆的粒子数Nr和噪声强度 η 构成的二维参数空间中的统计动力学行为.

3 模拟结果和讨论

封闭空间中的自驱动粒子倾向于在边界附近聚集[30−34].不同于球状粒子,本文模型由杆状粒子组成.由图2中的快照可以看出,这样的自驱动杆在柔性环附近形成比较规则的极化液晶态排列.由于杆粒子间向列型的相互作用,整个柔性环内的自驱动杆状粒子基本是中心对称分布的,此时用系统平均取向模的大小表示整体的极性程度,柔性环内反向粒子的取向相互抵消,系统的极性序接近于零.但从快照上可以看出,内部自驱动杆状粒子无序分布和有序聚集在环边界上是两种不同的分布,而这两种分布的极性序都接近零无法很好地区分开.

图2 三种典型分布的快照,自驱动杆粒子数Nr均为1500,噪声大小η分别为0.10,0.20和0.50,依次对应(a)自驱吸附有序态、(b)过渡态和(c)无序态.粒子颜色代表取向,同图 1Fig.2.The snapshots of three regions with fixed particle number Nr=1500 for different noise levels,and,respectively,with(a)η=0.10,self-propelled particle absorbed ordered region,(b)η=0.20 transient region,and(c)η=0.50 disordered phase.The color represents the radial direction as Fig.1.

为分析柔性环内自驱动杆状粒子角度上的分布,定义一个径向极性序参量

其中θi是杆i的取向,φi是杆i到环心位置相对位移的方向.该极性序是粒子取向与其相对位移方向两者夹角余弦值的平均,反应了内部自驱动杆状粒子在沿环质心向外方向上的取向有序程度.极性序SP趋近0时,内部自驱动杆的取向是各向均匀的;而当极性序 SP趋近1时,这些粒子基本都是背离环质心指向环外.对于单个粒子,其夹角的余弦值可以为负,但由于柔性环边界会聚集内部的自驱动粒子,系统整体平均后的极性序 SP基本都是正的,且极性序 SP的值越大,表明在该参数点下,内部自驱动粒子在柔性环边界的聚集程度越高.

3.1 极性序相图

本文主要研究体系的密度和噪声对系统形态的影响.系统改变自驱动杆的粒子数 Nr和杆端的噪声强度 η,测量各参数空间点的径向极性序 SP的值,可以得到如图 3(a)所示的相图.明显地,根据极性序 SP值的大小,相图中从最左边极性序 SP接近1的有序相区域,逐渐过渡到右侧极性序 SP接近0的无序区.有序区主要集中在粒子数 Nr较大,噪声强度 η 较低的区域,对应的快照如图2(a)所示.大部分的自驱动杆状粒子都指向环外方向,集中排列在柔性环上,且可以构成完整的层状分布.同时剩余的粒子在中心区域形成角度和位置都比较均匀地分布.无序区域则主要分布在粒子数Nr较小或噪声强度 η 较大的区域,如图 2(c)所示,其内部自驱动粒子的取向是无序的,均匀分布在环内.过渡区间主要分布在这两相之间的区域,如图2(b)所示,外层的自驱动杆无法形成完整的层状稳定排布,而分别集中成反向的两个集团或异向的多个集团.外侧的柔性环也因此有明显的变形,中心区域同样存在一定密度的比较均匀的无序气态分布.

3.2 粒子分布特性

由三个相区不同的粒子分布可知,内部自驱动粒子除了角度分布上具有向外的极性取向,粒子本身的空间位置分布也存在各向异性的情况.为分析粒子位置分布的不均匀性,根据所有自驱动杆的位置信息定义体系分布非球度 Δ .

首先由所有杆的质心位置计算惯量张量Q,其元素分别是

〈···〉t表示系统稳定状态各时刻的时间平均.显然非球度 Δ=0 对应于粒子在各个方向上均匀分布,非球度 Δ 的值越大,表明粒子分布的各向异性越明显,当非球度 Δ=1 时粒子基本分布在一条直线上.

图3 改变噪声强度 η和弹性环中自驱动杆粒子数 Nr 得到的相图(a)比较径向极性序参 Sp 大小得到的热力图;(b)比较非球度 Δ 大小得到的热力图,其中转变区域具有极大值Fig.3.Phase diagrams for self-propelled rods in elastic-ring with varying the noise strength η and the number of selfpropelled rods Nr,and the order parameter corresponding to(a)the radial polarity SP and(b)the asphericity Δ .We have maximal asphericity Δ in the transition region.

通过计算各参数点的非球度 Δ,可以得到如图3(b)的热度图,在有序相和无序相区域,系统的非球度 Δ 都是接近0,这两个相区内粒子的位置分布都是各项同性的.无序区由于粒子取向的无序性,各向同性的粒子位置均匀分布在柔性环内;而有序区柔性环中心虽然也存在类似无序区的粒子均匀分布,但它的极性序 SP主要由聚集在柔性环边界层状分布粒子贡献.这部分粒子在边界高密度堆积形成稳定的层状液晶分布,使系统整体处于类似汽液共存的动态平衡中.而过渡区具有相对较大的非球度 Δ,这是由于粒子数限制以及噪声影响,外层的自驱动粒子很难在环边界处形成完整的层状结构.在这个区间外层粒子由于角度上的偏离和数量上的减少,以及转变区间涨落的增强,对整体的径向极性序 SP贡献会有所降低.

内部自驱动粒子的位置分布除了在过渡相区会有明显的各向异性外,也容易在柔性环边界位置聚集,形成类似吸附相分离的密度分布差异.我们将柔性环边界附近聚集形成的比较高密度分布的粒子排布区域划分为高密度态.而中间比较无序分布的区域,划分为低密度区.由图 2可以看到,在有序区和过渡区(图2(a)和图2(b)),内外侧粒子分布存在明显的密度差异,而在无序区则不明显.通过比较两种密度态的密度差异,来表征粒子分布的相分离程度和粒子在环边界的聚集情况.

由于柔性环的形状易变,以及杆状粒子自身的各向异性,直接用杆的质心位置计算粒子密度不太合适.我们根据杆的质心位置先得到每个粒子的泰森多边形,计算各粒子相应的占据面积 Si,该面积对应每个粒子相对自由的运动范围.然后定义粒子的数密度,

其中ω 为符合高密度区或低密度区的粒子序号集合,|ω|为集合ω中的元素个数.相应的 ψin=表示系统的中心粒子数密度,ωin是沿环边界等比例收缩后中心气态区域内的粒子集合.是系统的外层粒子数密度,ωout即杆质心在沿环一圈宽度为杆长 Lr的环状区间内的粒子集合.由于内部自驱动粒子在环边界聚集排列,或低粒子数 Nr时与柔性环碰撞,环边界附近 ωout区域长时间尺度上稳定有粒子存在,而靠近环边界粒子的泰森多边形面积会明显小于内部粒子,最后得到的会 稍大于

为比较高低密度两相的分离程度,可以定义约化密度差为

约化密度差P的值越大说明系统内外两侧粒子分布的密度差异越大,P值越接近0则对应系统内部的粒子分布越均匀.通过计算系统的约化密度差可以得到如图4(a)所示的热度图,低噪声有序区域明显对应的区域约化密度差P值较高,高低密度两相分离明显.高噪声无序相区域的P值很低,接近均匀分布.这与图2快照中的粒子分布是一致的.可以看到约化密度差P与极性序 SP有类似的分布,随着噪声强度 η 的减弱,系统从无序相进入有序相的过程中,柔性环附近内部粒子的堆积程度也相应增强.极性序 SP主要由环附近堆积的粒子贡献,而自驱动粒子在环边界的堆积程度由约化密度差P刻画,所以极性序 SP的大小一定程度上与约化密度差P正相关.

图4 (a)改变噪声强度 η和自驱动杆粒子数 Nr,比较约化密度差P得到的热力图;(b)不同噪声强度 η 下,弹性环中心附近粒子数密度随自驱动杆粒子数 Nr 的变化趋势Fig.4.(a)Phase diagram of the reduced density difference P for self-propelled rods with varying the noise strength η and the number of self-propelled rods Nr;(b)density of central particles,ψ in,versus the particle number N r for different noise strength η .

注意到约化密度差的相图在高噪声强度 η和低粒子数 Nr区域,P的值会有反向的变化.此处的约化密度差P随粒子数 Nr的增长而减小.单独比较中心粒子数密度 ψin的变化,由图4(b)可以看到,不同噪声下随粒子数 Nr的 增加 ψin值最终都会有一个稳定的平台出现,对应一个饱和密度的中心低密度态的存在.该饱和密度与噪声强度有关,噪声越强饱和密度越高.而 ψin在达到饱和密度前随粒子数 Nr线性增加.在相图右侧高噪声区域,粒子数Nr比较小的情况下,自驱动杆状粒子很难在环边界附近形成稳定的集聚,此时外层区域 ωout内的粒子是与环边界碰撞的少数自驱动杆,因而外层粒子数密度 ψout比较稳定.此时中心粒子数密度 ψin未达到饱和,增加粒子数Nr会减小内部粒子的泰森多边形面积,即中心粒子数密度 ψin会随粒子数Nr增加而迅速增大.整体的约化密度差P主要受中心粒子数密度ψin的影响而减小.而当中心粒子数密度ψin达到饱和后,继续增加粒子数 Nr,中心粒子数密度 ψin基本保持稳定.柔性环附近的粒子受挤压,导致外层粒子数密度 ψout小幅增大.此时约化密度差P的值会相应有所增大.所以约化密度差相图中P值变化的极值位置,应该对应系统中心粒子数密度 ψin刚达到饱和密度的参数点,比较发现这两者基本是符合的.

在图 4(b)中,低噪声如 η=0.10 时,系统的中心粒子数密度 ψin会先有个小幅回落才能稳定在饱和密度.这是因为低噪声区域,粒子数 Nr较小时,内部的杆状粒子在自驱动作用下容易聚集在环边界位置,但未能形成完整的层状排列.噪声和粒子间的碰撞都参与到两侧的粒子交换中.当粒子数Nr增大到形成完整序列时,只有内侧层状排列的粒子在噪声作用下,与环中心区域的低密度态的粒子发生交换.此时的中心粒子数密度 ψin会有所降低.随噪声强度的增加,粒子碰撞贡献的交换逐渐减少,该密度回落过程也逐渐减弱,ψin达到饱和平台的曲线也愈加平滑.而当噪声较大时,如噪声强度 η=0.4,0.5对应的ψin曲线几乎重合,且没有明显水平的饱和平台.此时由于噪声太大,粒子很难在环边界形成有效的堆积,主要依靠粒子间排斥将柔性环撑开.且随粒子数增大,柔性环拉伸后对内部粒子的压力也相应增大,中心粒子数密度 ψin会随粒子数变化而有一个小幅的增长.

3.3 整体迁移

由于内部自驱动粒子分布的不均匀性,系统整体会有一定迁移运动,为分析系统整体的动力学行为,我们测量了柔性环质心的均方位移,

图5 弹性环及杆状粒子质心均方位移随时间的变化(a)粒子数Nr 为 1500 时,噪声大小η为0.10,0.20和0.50所在三个区区域 的比较;(b)粒 子数 Nr =1000,η 为0.25,0.30和0.50下无序态时的对比Fig.5.Mean-squared displacement(MSD)for the center of mass of particle and elastic ring:(a)Noise levels η=0.10,η=0.20,and η=0.50 for Nr=1500;(b)noise levels for η=0.25,η=0.30,and η=0.50 with particle number Nr=1000 in the disordered regime.

其中 R(t)是柔性环质心的位置.如图5(a)所示,我们测量了足够长时间尺度下,三个相区各自的均方位移及其斜率(补充材料movie1.mov,movie2.mov,movie3.mov,分别对应η=0.10,0.20,0.50),其中有序相和过渡相的均方位移斜率均接近2.无序相在长时间尺度时的均方位移斜率则接近1.系统在有序相和过渡相区域几乎都是在做整体的定向迁移运动,而在无序区则接近随机游走.由于中心低密度态的粒子数比较少,且在角度上均匀分布,其对系统整体的运动影响很弱,系统的整体移动主要由沿环附近堆积排列的粒子贡献.当有序相和过渡相在形成稳定分布后,其沿环附近堆积的粒子不会有很大的变化.由于很难实现完全的对称,未被完全相互抵消的自驱动贡献,会导致系统整体沿某一方向定向迁移.从图5(a)可以看出,由于有序相具有更高的对称性,多余的自驱动贡献比过渡相要低很多,其均方位移的尺度也要低一个数量级左右.而无序态的粒子分布和角度上的分布都很均匀,所以其瞬时的多余自驱动贡献在长时间尺度上是类似高斯白噪声的分布,其整体的运动也接近随机游走.

另外还比较了无序相区域不同噪声下的均方位移曲线,如图5(b)所示,随着噪声强度 η 的增加,系统自身对称性更高,其均方位移尺度的数量级会相应降低,同时其长时间的扩散系数也相应减小,接近无规的随机运动.在无序相中,系统没有稳定的外层堆积分布,其整体的迁移主要来自噪声和极性涨落的竞争.当噪声较弱时系统的整体运动趋向于定向迁移,而噪声很强时系统会有类似随机游走的运动行为.值得指出的是质心的运动有两种可能,一种是整个弹性环形变导致质心运动而整体可能并没有发生明显的迁移,另一种是系统整体的迁移.在目前的体系中,系统的形变带来的影响远远小于整体迁移的作用,形变导致的质心迁移对均方位移的贡献是可以忽略的.

4 结 论

本文研究了在二维条件下,自驱动杆状粒子受限在一条可伸缩的柔性环中的动力学行为.通过随机动力学模拟,发现根据内部自驱动粒子整体的径向极性,可以从低噪声的有序区过渡到高噪声低密度时的无序区.同时粒子的空间分布也会产生变化.最主要的区别在于柔性环附近内部粒子的堆积方式: 有序区内部自驱动粒子可以在环附近形成完整的层状排布,而无序区则由于高噪声影响,自驱动粒子无法在相应区域形成稳定的高密堆积.特别地,我们发现无序和有序转变的过渡区间环上有稳定的粒子堆积,但由于粒子数限制无法形成完整的层状排列,分散的堆积集团会导致柔性环明显形变.

除此之外,根据约化密度的异常转变,发现中心粒子数密度存在一个与噪声相关的饱和密度,该密度对应环边界附近位置开始形成稳定的粒子堆积.

系统整体的运动与这些边界集聚的粒子分布有关.在有序区和过渡区域中,环心的运动接近定向迁移.但由于有序相环边界位置粒子的层状分布更均匀,其整体定向迁移的强度会非常小.而当噪声增大,环边界附近稳定的粒子堆积逐渐减少,系统整体的运动也逐渐趋向无规则的随机游走.在我们的模型中,柔性环在过渡区会有明显形变,且此时系统整体的定向迁移最为显著.研究这类受限条件下的活性物质体系的动力学行为,对探讨细胞形变迁移方面的机制具有一定的参考意义[37].目前在一些人工受限体系,比如将分子马达与杆状的微丝、微管系统受限到液滴表面,研究其中活性液晶态的动力学演化也正受到越来越多人的关注[38,39].