曲率模量不同的兩組分膜泡形狀方程的數(shù)值解

周五斌，張劭光

（陜西師范大學(xué)物理學(xué)與信息技術(shù)學(xué)院，陜西西安710119）

曲率模量不同的兩組分膜泡形狀方程的數(shù)值解

周五斌，張劭光＊

（陜西師范大學(xué)物理學(xué)與信息技術(shù)學(xué)院，陜西西安710119）

采用變分法討論了旋轉(zhuǎn)對(duì)稱情況下曲率模量不相同的兩組分膜泡的歐拉－拉格朗日形狀方程及其邊界條件。通過(guò)雙向“打靶法”數(shù)值求解了兩組分膜泡在確定邊界條件下的形狀方程，計(jì)算了不同平均曲率模量比εκ和線張力系數(shù)λ下的平衡形狀。闡述了不同εκ和λ下平均曲率模量不相同的兩組分膜泡的形狀變化，導(dǎo)致這種變化的原因是膜泡兩組分的曲率能和線張力能相互競(jìng)爭(zhēng)的結(jié)果。計(jì)算結(jié)果說(shuō)明數(shù)值計(jì)算方法合理可行，此數(shù)值解可進(jìn)一步研究與實(shí)驗(yàn)相關(guān)的兩組分膜泡問(wèn)題。

兩組分膜泡；歐拉－拉格朗日方程；平均曲率模量

PACS：87.15.Aa

美國(guó)加州大學(xué)的Singer和Nicholson于1972年提出了關(guān)于生物膜結(jié)構(gòu)的流體鑲嵌模型。該模型認(rèn)為生物膜的基本結(jié)構(gòu)是由類脂類分子組成的雙層膜，同時(shí)蛋白質(zhì)分子和膽固醇則嵌在類脂雙層膜之中。類脂分子是雙親分子，有極性的頭端和疏水的尾端。極性的頭端與水接觸，而尾端因?yàn)槭杷畡t相互接觸，從而形成雙分子層膜。因此，人們猜測(cè)細(xì)胞膜的形狀是由雙層膜決定的。在此生物學(xué)模型的基礎(chǔ)上，為了解釋紅細(xì)胞的雙凹形狀，Canham和Helfrich獨(dú)立地提出了生物膜的曲率彈性模型，該模型認(rèn)為生物膜的形狀是在給定面積和體積的情況下，膜平均曲率能最小的形狀［1－2］。Deuling等對(duì)紅細(xì)胞的各種形狀在旋轉(zhuǎn)對(duì)稱下進(jìn)行了數(shù)值計(jì)算，得到了紅細(xì)胞的各種形態(tài)［3］。與此相關(guān)的實(shí)驗(yàn)工作也相繼展開(kāi)，人們發(fā)現(xiàn)不用考慮真實(shí)生物膜上磷脂分子的多樣性及蛋白質(zhì)骨架和膽固醇等，只用磷脂雙層亦能在試管中形成人造膜泡。理論及實(shí)踐相互促進(jìn)大大推動(dòng)了該理論的發(fā)展［4－5］。

脂筏模型認(rèn)為：真實(shí)的生物膜是由多種類脂分子和膽固醇（現(xiàn)不考慮蛋白質(zhì)）構(gòu)成，由此形成膜上不同成分的區(qū)域，扮演重要的生物學(xué)功能，脂筏的形成是由于不同類脂成分的相分離所導(dǎo)致。關(guān)于球形膜泡相分離的實(shí)驗(yàn)研究，已取得了很多進(jìn)展，人們發(fā)現(xiàn)由飽和磷脂分子、不飽和磷脂分子及膽固醇形成的巨型膜泡（Giant Vesicle），不同類脂分子之間的相互作用會(huì)導(dǎo)致相分離［6－7］，形成由飽和磷脂分子和膽固醇形成的Lo相（Liquid－Ordered Phase）子區(qū)域及由不飽和磷脂分子形成的 Ld相（Liquid－Disordered Phase）子區(qū)域。近來(lái)Yanagisawa等人［8］在球形的相分離的基礎(chǔ)上，加入山梨糖醇（Sorbitol），改變膜泡兩側(cè)的滲透壓，結(jié)果發(fā)現(xiàn)Lo區(qū)域發(fā)生向外或向內(nèi)吐出很多芽（Budding）的相分離，剩下的球形母泡將發(fā)生進(jìn)一步的相分離，變成三角海星形，最后變成凹盤形。凹盤形的上下兩面進(jìn)一步發(fā)生向外凸出，進(jìn)而發(fā)生發(fā)芽相變。他們還合成了旋轉(zhuǎn)對(duì)稱的長(zhǎng)橢球及扁橢球的兩組分膜泡［9］，發(fā)現(xiàn)這兩種形狀的相分離模式不同。Yanagisawa等人用橢球形的參數(shù)方程去近似表示得到的形狀，來(lái)估算兩種形狀的能量，與實(shí)驗(yàn)結(jié)果不完全符合。該形狀能否由相應(yīng)的歐拉 －拉格朗日方程來(lái)決定，如何理解這兩種相分離模式的轉(zhuǎn)變?nèi)允怯写鉀Q的問(wèn)題。

對(duì)于旋轉(zhuǎn)對(duì)稱性膜泡可以通過(guò)嚴(yán)格求解相應(yīng)的歐拉 －拉格朗日方程，來(lái)判斷該形狀是否存在，進(jìn)而給出相變的條件，黃聰?shù)热送ㄟ^(guò)這種方法得到了旋轉(zhuǎn)對(duì)稱情況下開(kāi)口膜泡的形狀［10］。Jülicher等人研究了旋轉(zhuǎn)對(duì)稱情況下兩組分膜泡的形狀，但他們只計(jì)算了兩個(gè)組分曲率模量相同時(shí)，線張力系數(shù)及滲透壓的作用［11］，對(duì)于實(shí)驗(yàn)上很重要的兩組分曲率模量不同所起的作用并沒(méi)有計(jì)算，而且對(duì)于實(shí)驗(yàn)觀察到的例如o－d－o（Order－Disorder－Order）這種形狀也未計(jì)算。本文采用雙向“打靶法”［12］數(shù)值求解了相應(yīng)的歐拉 －拉格朗日方程，得到了不同平均曲率抗彎模量比εκ和不同線張力系數(shù)λ下的一些兩組分膜泡的形狀，闡述了εκ和λ對(duì)兩組分膜泡形狀的影響。

1 兩組分膜泡的曲面彈性理論

1.1 模型

當(dāng)膜泡發(fā)生相分離形成兩組分共存的形狀時(shí)，膜泡由兩個(gè)組分α域和β域組成，其自由能由兩個(gè)組分的曲率能和兩個(gè)組分交界處的線張力能組成。另外，為了約束膜泡總面積和總體積，引進(jìn)兩個(gè)拉格朗日乘子∑和Δp，此時(shí)膜泡的總自由能［11］可表示為

上式中第一項(xiàng)是α域的曲率能，第二項(xiàng)是β域的曲率能，C1和C2表示膜泡任意一點(diǎn)的兩個(gè)主曲率，C和分別表示兩個(gè)組分的自發(fā)曲率和分別表示α域和β域的曲率模量，κ和κ分別表示α域和β域的高斯曲率模量。第三項(xiàng)是兩個(gè)組分交界處的線張力能，是沿著α域和β域的交界?α進(jìn)行的線積分，σ表示線張力系數(shù)。第四項(xiàng)和第五項(xiàng)是為了約束兩個(gè)組分的總面積而引入的，∑（α）和∑（β）表示兩個(gè)拉格朗日乘子。最后一項(xiàng)則是為了約束總體積引入的，Δp也表示拉格朗日乘子。

圖1 旋轉(zhuǎn)對(duì)稱的兩組分膜泡參數(shù)示意圖Fig.1 Schematic diagram of a two－domain vesicle

圖1是兩組分膜泡的參數(shù)示意圖，可視為理想的旋轉(zhuǎn)對(duì)稱曲面。實(shí)線代表α域，虛線代表β域，膜泡上任意一點(diǎn)的切線與水平方向的夾角為Ψ（s），s∈［0，s2］是輪廓線的弧長(zhǎng)，s＝s0＝0時(shí)表示膜泡南極點(diǎn)，s＝s2時(shí)表示膜泡的北極點(diǎn)，s＝s1時(shí)是兩組分的交界處。R（s）為任意一點(diǎn)到對(duì)稱軸的距離，Z則是旋轉(zhuǎn)軸，且滿足。

采用圖1所示的參數(shù)化方法，可以得到任意一點(diǎn)的兩個(gè)主曲率為C1＝Ψ·，C2＝siunΨ／R。且在圖（1）所示的參數(shù)化方法下dA＝2πRds，dV＝πR2sinΨds。則（1）式可表示為

其中R（s1）是s＝s1時(shí)的半徑，需要指出的是當(dāng)曲面的拓?fù)浣Y(jié)構(gòu)不變時(shí)，（2）式中涉及α域和β域的高斯曲率能的積分可用高斯 －波涅公式表示為

可見(jiàn)這一項(xiàng)只與兩域交界處的cosΨ值有關(guān)，而與兩域的形狀無(wú)關(guān)。因而它只影響兩域交界處的邊條件（參看后面（16）式），而不進(jìn)入相應(yīng)的歐拉 －拉格朗日方程（參看后面（11）式）。將（3）代入（2）式得

（4）式就是在旋轉(zhuǎn)對(duì)稱條件下兩組分膜泡的自由能表達(dá)式。

1.2 形狀方程及邊界條件的推導(dǎo)將（4）式化簡(jiǎn)可得

其中

在閉曲面的情況下膜泡的形狀方程可由δ＾F＝0得到，對(duì)（5）式做一階變分可以得到

其中i1和i2分別指i域的上邊界和下邊界，也就是說(shuō)α1＝s0，α2＝s1，β1＝s1，β1＝s2。而哈密頓量

因?yàn)棣摩贰ⅵ腞、δγ是任意的，為了使（8）式等于0，只能使系數(shù)等于0，即

將（6）式代入（10）式可得兩組分膜泡的形狀方程如下

其中Z、A、V分別表示任意一個(gè)組分的高度、面積、體積。邊界處將Ψ（s0）＝0，Ψ（s2）＝π，R（0）＝0，R（s2）＝0代入哈密頓量的表達(dá)式并使H（α）＝H（β）＝0可得γ（s0）＝γ（s2）＝0。于是就可以得到方程滿足的初始邊界條件

然而在兩個(gè)組分的交界處δΨ（s1）和δR（s1）都不等于0，將（8）式代入（7）式可得

化簡(jiǎn)上式可得

其中ε是一個(gè)無(wú)窮小量，將（6）式代入（14）、（15）式就可以得到在兩組分交界處的邊界條件

2 數(shù)值求解平衡形狀方程

為了方便編程計(jì)算，我們要將一些參數(shù)無(wú)量綱化，引進(jìn)一個(gè)無(wú)量綱參數(shù)t＝（s－s0）／（s2－s0），其中t∈［0，1］，s（t0）＝0，s（t1）＝s1，s（t2）＝s2，這樣R、Z、Ψ、s都成為t的函數(shù)。而膜泡的約化半徑為R0＝，約化體積為v＝V／（4πR／3），α部分的面積分?jǐn)?shù)為X＝A（α）／（A（α）+A（β）），兩個(gè)組分交界面處的線張力為λ＝σR0／κ（α），兩個(gè)組分曲率模量的比為εκ＝κ2／κ1。由于膜泡的體積不固定，在改變參數(shù)使膜泡形狀發(fā)生變化的時(shí)候，膜泡的體積可以自己調(diào)整，這就使得拉格朗日乘子ΔP＝0。另外，本文主要計(jì)算了εκ對(duì)兩組分膜泡形狀的影響，因此兩組分的自發(fā)曲率取，高斯曲率模量取

2.1 不同εκ下的膜泡形狀

圖2 不同εκ對(duì)應(yīng)的膜泡形狀Fig.2 The shape of vesicles with differentεκ

圖2是在約化線張力系數(shù)λ＝5.4，ΔP＝0，面積分?jǐn)?shù)X＝0.4時(shí)不同的εκ值下兩組分膜泡的形狀，實(shí)線對(duì)應(yīng)曲率模量為k1的組分的形狀，虛線對(duì)應(yīng)曲率模量為k2的組分的形狀。因?yàn)閷?shí)驗(yàn)上εκ的范圍通常為1.0～4.0，因此，本文只計(jì)算了εκ從1.0到4.0的結(jié)果。從圖中可以看出，隨著εκ的增加，膜泡的形狀發(fā)生了很大的變化，在εκ＝1.0的時(shí)候兩組分的交界處圓的半徑幾乎為0，在εκ增加到1.465的過(guò)程中，脖子慢慢打開(kāi)，兩組分交界處圓的半徑慢慢增加。到εκ＝1.465時(shí)，由于能量最小化膜泡會(huì)發(fā)生一個(gè)不連續(xù)相變，從一支解的穩(wěn)定部分到達(dá)另一支解的穩(wěn)定部分，然后隨著εκ的增加，兩組分交界處圓的半徑會(huì)越來(lái)越大，出現(xiàn)一系列上下不對(duì)稱的形狀。

2.2 不同λ下的膜泡形狀

圖3 不同λ對(duì)應(yīng)的膜泡形狀Fig.3 The shape of vesicles with differentλ

圖3是在εκ＝2.5，ΔP＝0，面積分?jǐn)?shù)X＝0.2時(shí)對(duì)應(yīng)不同的約化線張力系數(shù)λ值的兩組分膜泡的形狀，實(shí)線和虛線分別表示兩個(gè)組分的外形。從圖中可以看出，隨著λ的增加，兩組分交界處的圓的半徑一直減小，直到在λ＝9.774時(shí)兩組分交界處的半徑幾乎為0，再增加λ相當(dāng)于兩個(gè)單一組分的球形膜泡靠在一起，因此本文只計(jì)算了λ從0到9.774的結(jié)果。在增加λ的過(guò)程中，當(dāng)λ＝7.837時(shí)，因?yàn)橐鼓芰孔钚』?，兩組分膜泡會(huì)發(fā)生不連續(xù)相變，到達(dá)另一支解的穩(wěn)定部分。

3 結(jié)論

本文通過(guò)求解旋轉(zhuǎn)對(duì)稱情況下的兩組分膜泡的歐拉－拉格朗日方程，并用雙向“打靶法”得出了一些兩組分膜泡的平衡形狀，得到了不同εκ和不同λ下的平衡形狀，闡述了不同εκ和λ下兩組分膜泡的形狀變化，這種變化是膜泡兩組分的曲率能和線張力能相互競(jìng)爭(zhēng)引起的。說(shuō)明用雙向“打靶法”求解雙組分膜泡歐拉－拉格朗日方程的方法是可行的，為以后進(jìn)一步研究與實(shí)驗(yàn)相關(guān)的問(wèn)題提供了參考。本文是在無(wú)體積約束的情況下數(shù)值求解形狀方程，未來(lái)要研究加上體積約束之后的兩組分膜泡形狀方程的數(shù)值解，并且進(jìn)一步推導(dǎo)和求解三區(qū)域雙組分膜泡的平衡形狀方程。

［1］Canham P B.The minimum energy of bending as a possible explanation of the biconcave shape of the human red blood cell［J］.Journal of Theoretical Biology，1970，26：61－81.

［2］Helfrich W.Elastic properties of lipid bilayers：theory and possible experiments［J］.Zeitschrift Naturforschung.Teil：C，1973，28（11）：693－703.

［3］Deuling H J，Helfrich W.Red blood cell shapes as explained on the basis of curvature elasticity［J］.Biophysical Journal，1976，16：861－868.

［4］Seifert U.Congurations of fluid membranes and vesicles［J］.Advances in Physics，1997，46（1）：13－137.

［5］謝毓章，劉寄星，歐陽(yáng)鐘燦.生物膜泡曲面彈性理論［M］.上海：上海科學(xué)技術(shù)出版社，2003.

［6］Veatch S L，Keller S L.Organization in lipid membranes containing cholesterol［J］.Physical Review Letters，2002，89：268101－268104.

［7］Baumgart T，Hess S T，Webb W W.Imaging coexisting fluid domains in biomembrane models coupling curvature and line tension［J］.Nature，2003，425：821－824.

［8］Yanagisawa M，Imai M，Taniguchi T.Shape deformation of ternary vesicles coupled with phase separation［J］.Physical Review Letters，2008，100：148－152.

［9］Yanagisawa M，Imai M，Taniguchi T.Periodic modulation of tubular vesicles induced by phase separation［J］. Physical Review E，2010，82：051928.

［10］黃聰，張劭光.高斯曲率彈性模量對(duì)啞鈴形開(kāi)口膜泡形狀的影響［J］.陜西師范大學(xué)學(xué)報(bào)：自然科學(xué)版，2012，41（3）：31－35.

［11］Jülicher F，Lipowsky R.Domain－induced budding of vesicles［J］.Physical Review Letters，1993，70：2964－2967.

［12］William H Press，Saul A Teukolsky，William T Vetterling，et al.Numerical recipes in fortran 90［M］. London：Cambridge University Press，1996：1314－1328.

〔責(zé)任編輯 李 博〕

The numerical solutions of shape equation of two－domain vesicles with different bending rigidity

ZHOU Wubin，ZHANG Shaoguang＊
（School of Physics and Information Technology，Shaanxi Normal University，Xi′an 710119，Shaanxi，China）

The Euler－Lagrange shape equation and boundary conditions of two－domain vesicles which have different curvature modulus with rotational symmetry are investigated by variational method.Then the numerical solutions of shape equation of two－domain vesicles under definite boundary conditions are obtained through shooting method.The equilibrium shape at differentεκand different line tension coefficientλare confirmed.The shape change of two－domain vesicles which have different curvature modulus is find at the differentεκandλ.The reason of such phenomenon is bending energy and tension energy of two－domain vesicles competing with each other. The results show that the numerical method is reasonable for study the equilibrium shape of twodomain vesicles，which help for further study of experiment－related problems for two－domain vesicles with different curvature modulus.

two－domain vesicles；Euler－Lagrange equation；mean curvature modulus

O411.3

：A

1672－4291（2015）05－0043－05

10.15983／j.cnki.jsnu.2015.05.255

2015－01－21

國(guó)家自然科學(xué)基金（10374063）；中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金（GK261001071）

周五斌，男，碩士研究生，研究方向?yàn)槔碚撋镂锢怼－mail：zwb50＠163.com

＊通信作者：張劭光，男，副教授。E－mail：zhangsg＠snnu.edu.cn