邱 星，楊紫貝，夏詩琪，王海燕

(北京師范大學(xué) 物理學(xué)系，北京 100875)

布朗運(yùn)動自發(fā)現(xiàn)以來被很多科學(xué)家深入研究，布朗運(yùn)動的理論在物理、化學(xué)甚至經(jīng)濟(jì)領(lǐng)域得到了廣泛的應(yīng)用[1-3]．但生物體中的主動運(yùn)動，比如分子馬達(dá)在軌道上的定向運(yùn)動，細(xì)胞的遷移和微生物的定向運(yùn)動等，這些運(yùn)動主體可以吸收能量，使系統(tǒng)遠(yuǎn)離熱平衡狀態(tài)[4-6]，基于熱漲落的布朗運(yùn)動理論(本文稱之為被動布朗運(yùn)動理論)很難解釋這些現(xiàn)象．因此，為了揭示生物體中主動輸運(yùn)的物理機(jī)制，研究人員提出具有能量子庫的活性布朗粒子模型[7-12]．研究表明，具有能量子庫的活性布朗粒子模型不但可以解釋生物體內(nèi)定向運(yùn)動的物理機(jī)制，而且可以加深人們對非平衡統(tǒng)計物理的理解.

活性布朗粒子模型最初由德國洪堡大學(xué)Werner Ebeling團(tuán)隊提出[10,11]，其重要特性就是摩擦系數(shù)是粒子速度的非線性函數(shù)，發(fā)現(xiàn)活性布朗粒子從環(huán)境中吸收能量，儲存于能量子庫，并能夠?qū)⒛芰哭D(zhuǎn)化為定向運(yùn)動．Yu. M. Romanovsky等進(jìn)一步研究了內(nèi)部振子與能量子庫耦合的活性布朗馬達(dá)模型，并將該模型應(yīng)用于轉(zhuǎn)動馬達(dá)的研究中[8]．A. Fiasconaro等研究了非對稱棘輪勢中的活性布朗粒子在泊松噪聲驅(qū)動下運(yùn)動的定向性和能量轉(zhuǎn)化效率對隨機(jī)能量關(guān)聯(lián)時間的依賴關(guān)系[9]． 浙江大學(xué)鄧茂林[13-15]等學(xué)者從理論上研究了活性布朗粒子的動力學(xué)，得到在拋物型勢能和高斯白噪聲驅(qū)動下的穩(wěn)態(tài)解析解，發(fā)現(xiàn)在四維相空間上，系統(tǒng)處于遠(yuǎn)離平衡態(tài)時的擴(kuò)散極限環(huán)．由于活性布朗粒子能夠從環(huán)境中吸收能量，活性布朗粒子系統(tǒng)遠(yuǎn)離熱力學(xué)平衡狀態(tài)，呈現(xiàn)出復(fù)雜的動力學(xué)行為和反常擴(kuò)散特征,還存在很多富有挑戰(zhàn)性的問題[16].

本文基于德國洪堡大學(xué)Werner Ebeling研究小組提出的活性布朗粒子模型[10，11]，分析了內(nèi)部能源子庫的活性布朗粒子模型的阻尼特性，通過與被動布朗粒子比較，分別在自由場、諧振子勢和非對稱周期勢中利用了Matlab軟件模擬活性布朗粒子的隨機(jī)動力學(xué)行為，發(fā)現(xiàn)非對稱周期勢壘的高度是活性布朗粒子定向運(yùn)動的關(guān)鍵因素之一．本研究可以為統(tǒng)計物理教學(xué)提供一個有趣的實例，同時為大學(xué)生利用物理規(guī)律解決實際問題提供可以借鑒的思路和方法.

1 模型與方法

1.1 活性布朗粒子模型

為研究微小生物體的運(yùn)動特性，本文不考慮微小的生物體內(nèi)部結(jié)構(gòu)，把微小的生物體抽象為具有內(nèi)部能源子庫的布朗粒子．內(nèi)部能源子庫總能量用e(t)表示，q(r)是單位時間在r處流入庫的能量，內(nèi)部能源庫的平衡方程可以表示為

(1)

其中，c是生物體內(nèi)部損耗系數(shù)，代表了生物體內(nèi)部能量損耗；d(v)是能量轉(zhuǎn)化為動能的比率，可以表達(dá)為

d(v)=d2v2

(2)

d2是轉(zhuǎn)換系數(shù)，且d2≥0．v是粒子運(yùn)動速度．考慮整個系統(tǒng)的能量，E0是系統(tǒng)的機(jī)械能，系統(tǒng)機(jī)械能平衡方程:

(3)

其中，γ0是阻尼系數(shù)，U(r)是勢能函數(shù)，E0(t)是系統(tǒng)的機(jī)械能．在系統(tǒng)中，機(jī)械能的變化來自兩個方面：能量子庫輸入的能量：d2ev2；由于摩擦損耗的能量：γ0v2．考慮隨機(jī)力的影響，活性布朗粒子的運(yùn)動方程為

(4)

D=kTγ0

(5)

k是玻耳茲曼常數(shù)，T是環(huán)境溫度．與描述被動布朗粒子的郎之萬方程相比，式(4)增加了一項d2ev，這一項表明內(nèi)部能源庫在速度v方向上提供了加速度．通過把速度v的系數(shù)合并，我們重新定義摩擦系數(shù)γ(v):

γ(v)=γ0-d2e(t)

(6)

(7)

在平衡時，子庫的能量e0是流入能量通量和系統(tǒng)損耗能量的系數(shù)的比值，損耗主要來自系統(tǒng)內(nèi)部損耗和輸出轉(zhuǎn)化為動能．將式(7)帶入式(6)，我們得到準(zhǔn)靜態(tài)過程中的摩擦系數(shù)與粒子速度的關(guān)系，如圖1 所示.

圖1 摩擦系數(shù)和速度(γ0=2.0,c=0.2,q0=1.0,d2=1.0)

(8)

1.2 模擬方法

本文采用的數(shù)值解法是改進(jìn)-歐拉法，又稱預(yù)估-修正算法:

(9)

un表示函數(shù)在運(yùn)算進(jìn)行n步時的函數(shù)值，h是每一步的步長，xn是運(yùn)算進(jìn)行到第n步時自變量的值，f(xn,un)是un對xn的導(dǎo)數(shù)．歐拉法的思想是函數(shù)在n+1步的值由第n步的值和第n步函數(shù)的微分決定，預(yù)估-修正算法則通過計算第n步的導(dǎo)數(shù)，得到歐拉法中un+1的值，再計算出n+1的導(dǎo)數(shù)f[xn+1,un+hf(xn,un)]，兩個導(dǎo)數(shù)的平均值才是第n步計算需要的導(dǎo)數(shù)值，這種方法的優(yōu)勢是可以提高計算的精度.

在具體的模擬中，本文所有的數(shù)據(jù)均是模擬1 000個粒子經(jīng)歷100 s的運(yùn)動，這是因為改進(jìn)歐拉法數(shù)值解方程選取了較短的步長(h=0.001)．為了便于模擬，對式(4)進(jìn)行無量綱處理得到

(10)

2 結(jié)果與討論

2.1 活性布朗粒子在自由場中的運(yùn)動

在圖2中，我們模擬了活性布朗粒子和被動布朗粒子位移的分布，結(jié)果表明活性布朗粒子的位移呈現(xiàn)雙峰分布．在正負(fù)半軸的雙峰都是正態(tài)分布，峰的中心位置關(guān)于原點對稱，出現(xiàn)粒子數(shù)的幾率也近似相等，說明在自由場中，活性布朗粒子沒有定向運(yùn)動．為了與被動布朗運(yùn)動比較，在圖3中，我們模擬了被動布朗粒子位移的分布，發(fā)現(xiàn)被動布朗粒子在原點附近呈正態(tài)分布，也不存在定向運(yùn)動.

圖2 活性布朗粒子位移分布[x=x0，y=y0代表位移在(x0，x0+0.25)范圍內(nèi)的粒子數(shù)為y0. t=100，γ=2.0，D0=0.01，d2=1.0，q0=3.0，c=0.1]

圖3 被動布朗粒子位置分布(x=x0，y=y0代表著位移在[x0，x0+0.2)范圍內(nèi)的粒子數(shù)為y0. t=100，γ=2.0，D0=0.01)

無論是活性還是被動的布朗粒子，它們在自由場中都未發(fā)生定向運(yùn)動，說明粒子的運(yùn)動屬于自由擴(kuò)散．在圖4中，我們分別計算了活性布朗粒子和被動布朗粒子的均方位移隨時間的變化．從均方位移的數(shù)值上看，活性布朗粒子均方位移幾乎是被動布朗粒子均方位移的104倍，表明活性布朗粒子擴(kuò)散更主動，隨機(jī)運(yùn)動的距離更遠(yuǎn)．被動布朗粒子的均方位移與時間成正比，這符合愛因斯坦的理論[1]．而活性布朗粒子的均方位移隨時間變化偏離了線性關(guān)系，類似于二次曲線，屬于反常擴(kuò)散行為[16].

圖4 活性和被動布朗粒子均方位移((a)活性布朗粒子的-t圖像，(b)被動布朗粒子的-t圖像. 其中：γ=1.0，D0=0.01，d2=1.0，q0=3.0，c=0.1)

2.1 活性布朗粒子在諧振子勢中的運(yùn)動

通常用諧振子勢表示單分子實驗中光阱對粒子的作用，本文選取二維諧振子勢:

(11)

這個勢能產(chǎn)生與位移相反的力：F=-k(x+y).

在諧振子勢中，首先在一定的時間內(nèi)分別模擬了活性布朗粒子和被動布朗粒子的運(yùn)動，如圖5 所示，我們觀察到在諧振子勢中粒子的自由擴(kuò)散受到了限制，對比圖5(a)和(b)還可以看出一定時間內(nèi)活性布朗粒子擴(kuò)散區(qū)域更大．然后，我們研究了勁度系數(shù)k對活性布朗粒子運(yùn)動的影響， 如圖6所示．我們可以看出：當(dāng)k變小時，粒子的均方位移整體變大，圖像上下震動的幅度變大，粒子達(dá)到平衡所需的時間變長；而k增大時，粒子的均方位移整體變小，圖像上下震動的幅度變小，粒子達(dá)到平衡所需的時間變短．同時，我們發(fā)現(xiàn)在諧振子勢中粒子的平均位移接近零，這說明在諧振子勢中，粒子運(yùn)動并沒有產(chǎn)生定向運(yùn)動.

圖5 粒子在二維諧振子勢中的運(yùn)動((a)活性布朗粒子；(b)被動布朗粒子. 其中：γ=0.2，D0=0.01，d2=1.0，q0=3.0，c=0.1，二者模擬時間均為100 s)

圖6 諧振子勢中活性布朗粒子均方位移和位移平均圖像((a)均方位移；(b)位移平均. 參數(shù)設(shè)置： γ=1，D0=0.01，d2=2.0，q0=3.0，c=0.1)

2.3 活性布朗粒子在非對稱周期勢中的運(yùn)動

以上研究表明，活性布朗粒子在自由場和對稱的諧振子勢中不存在定向運(yùn)動，因此，我們考慮非對稱周期勢對粒子運(yùn)動影響，選取組合正弦函數(shù)來模擬非對稱周期勢:

(12)

x0是調(diào)配周期函數(shù),使x=0時，V為最小值，在式(12)中，經(jīng)過計算選取x0=-3.239．V0用來調(diào)節(jié)勢壘高度．在這個周期勢中，勢能變化較慢一側(cè)為a側(cè)，另一側(cè)則為b側(cè)．首先，在相同的勢壘高度和模擬時間條件下，我們計算了活性布朗粒子在非對稱周期勢中的位移分布，同時比較了被動布朗的情況，如圖7所示.

圖7 非對稱周期勢中布朗粒子的運(yùn)動((a)活性布朗粒子位移分布；(b)被動布朗粒子最終位移分布．活性布朗粒子參數(shù)設(shè)置：γ=1，D0=0.01，d2=2.0，q0=2.5，c=0.1，V0=1/2．被動布朗粒子參數(shù)：二者模擬時間均為100 s)

從圖7我們可以觀察到活性布朗粒子的位移分布不再對稱，出現(xiàn)了明顯的定向運(yùn)動，而被動布朗粒子的位移分布相對初始位置是對稱的，說明沒有出現(xiàn)定向運(yùn)動．然后，我們研究了勢壘高度對活性布朗粒子平均速度的影響，如圖8所示．根據(jù)F.Schweitzer等人[10]提出了活性布朗粒子維持“上坡”運(yùn)動的限制條件:

(13)

圖8 (a) d2和V(x)的示意圖,q0=1.0，γ=1.0，V0=1/8．(b) 不同V0下活性布朗粒子的平均速度-t時間圖像，參數(shù)：γ=1.0，D0=0.01，d2=1.0，q0=1.0，c=0.1

3 結(jié)論

細(xì)胞內(nèi)活性物質(zhì)的自推進(jìn)運(yùn)動一直是生物和醫(yī)學(xué)研究的熱點問題，而活性布朗粒子模型是研究自推進(jìn)運(yùn)動基礎(chǔ)的模型之一．本文基于經(jīng)典的活性布朗運(yùn)動模型，首先對方程進(jìn)行無量綱處理，然后采用預(yù)估-修正算法數(shù)值解方程，最后用Matlab作圖分析運(yùn)動情況.

在自由場中，我們的計算表明活性布朗粒子不存在定向運(yùn)動，而且均方位移隨時間非線性變化．其次，在諧振子勢中，我們發(fā)現(xiàn)粒子的自由擴(kuò)散受到了限制，但在相同時間內(nèi)，活性布朗粒子的擴(kuò)散區(qū)域比被動布朗粒子大．最后，在非對稱周期勢中，我們計算了活性布朗粒子的位置分布，同時計算不同勢壘高度下，活性布朗粒子平均速度隨時間的變化，結(jié)果表明活性布朗粒子能否將儲存子庫的能量轉(zhuǎn)化為定向運(yùn)動依賴于勢壘高度．本研究有助于理解為生物中主動輸運(yùn)的物理機(jī)制，同時可以加深我們對非平衡統(tǒng)計物理的認(rèn)識.