自驅動桿狀粒子在半柔性彈性環中的集體行為*

仲穎 施夏清

(蘇州大學物理科學與技術學院,軟凝聚態物理及交叉研究中心,蘇州 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].