DNA超分子水凝膠的粗?；Ｅc模擬?

王曦1)2) 黎明1) 葉方富2)1) 周昕1)

1)(中國科學院大學物理科學學院,北京 100049)

2)(中國科學院物理研究所,北京 100190)

1 引 言

超分子水凝膠是一種利用超分子相互作用如氫鍵、疏水作用等進行交聯(lián),形成空間網絡的高分子聚合物[1,2].近年來,雙鏈DNA分子因如下原因成為頗受關注的水凝膠材料:較強的機械強度、良好的生物相容、容易通過序列設計對凝膠結構進行調控、DNA鏈之間的堿基互補配對具有足夠的強度和可逆性.純DNA水凝膠[3?5]和DNA橋接形成的各種凝膠[6,7]成為近年來凝膠領域的重要課題.

迄今為止,文獻中已經報道了多種DNA水凝膠[8,9].但對其結構與性能的關系,特別是對凝膠網絡的分子化學結構、分子鏈性質、拓撲孔洞結構、交聯(lián)密度與強度等結構參數及其與凝膠宏觀力學性質之間的關系的了解還十分缺乏,主要是缺少對相應體系的理論建模和實驗測量的比較研究[10].因此,構建DNA水凝膠的微觀模型并計算模擬其分子結構和宏觀性質的功能關系,對理解和設計DNA超分子水凝膠具有重要意義.

DNA水凝膠的宏觀性質主要依賴于其分子鏈段層次的結構與相互作用,其原子尺度的細節(jié)及水環(huán)境的影響在建立微觀模型時需要合理地近似和簡化.一般來說,對于DNA分子及其組成的材料,依據系統(tǒng)的大小及特征時間尺度,可采用不同的策略進行建模和模擬[11].當體系大小僅為百堿基對(base pair,bp)量級、所關注的過程在微秒量級及以下時(例如DNA分子自身的形變),可使用全原子模擬[12];當體系大小為千bp、時間尺度為百微秒量級(例如DNA分子與其他生物分子的相互作用),可使用保留一定細節(jié)的粗?；Ｐ?如MARTINI力場下的DNA模型[13];當體系大小超過萬bp、時間尺度為亞毫秒及以上量級,則需要使用粗?；潭雀叩哪Ｐ蚚14].通常模型中保留的分子細節(jié)越多,得到的結果越真實,但所需計算量也越大.對DNA水凝膠等由DNA分子模塊組裝而成的結構,一般需要忽略大部分DNA分子結構細節(jié),僅保留DNA分子模塊的形貌和模塊間的相互作用等特征以構建模型,進而計算模擬這些微觀特征與組裝結構及性能的關系.

本文以文獻[15]報道的一種純DNA自組裝水凝膠為例,建立粗?；Ｐ?模擬研究該DNA水凝膠的分子模塊、相互作用等微觀特征與凝膠形成、結構及性質的關系,得到了與實驗半定量或定性相符的結果.這個工作不僅有助于理解與優(yōu)化設計該類DNA水凝膠,也為更普遍地構建合理有效的粗粒化模型以研究各種DNA組裝體提供了經驗.

2 模型與方法

2.1 組裝模塊的粗?；Ｐ?/h3>

目前,Liu等[15,16]設計了可自組裝的模塊化純DNA水凝膠.該水凝膠通過預先設計序列的DNA形成小型多分支分子模塊,進而使這些分子模塊之間通過其末端堿基間的互補配對,自發(fā)組裝形成空間網絡結構.此類水凝膠可以快速自組裝,熱響應性和酶響應性良好,機械強度較大,在細胞培養(yǎng)、生物學3D打印等領域具有良好的應用前景.我們以此為例,建立粗粒化模型,模擬其微觀結構與性能的關系.

該DNA水凝膠由具有三條分支的Y形分子(Y-scaffold)和線狀的連接子(linker)構成[15].如圖1所示,Y形分子三條臂均為雙鏈DNA,每條臂的末端均有一段單鏈DNA(ssDNA)作為黏性末端.連接子為線狀雙鏈DNA(長約32 bp),兩端各有一段ssDNA作為黏性末端.Y形分子和連接子的末端ssDNA序列互補,兩者可以自發(fā)進行堿基互補配對,形成黏性連接.實驗中典型的參數是Y形分子每條雙鏈DNA臂長14 bp,連接子雙鏈DNA長約32 bp,所有ssDNA黏性末端長約8 bp.

由于該DNA水凝膠形成與構象改變對應的較大時間尺度(百微秒毫秒)與空間尺度(亞微米微米),常用的DNA粗?；Ｐ蚚17]保留細節(jié)較多,導致計算量過大而不適用.為此,我們建立了一個粗?；潭葮O高的模型.考慮到Y形分子的分支和連接子的雙鏈DNA段的長度遠小于雙鏈DNA的駐留長度(約為50 nm,150 bp),我們將其處理為硬桿.考慮到水凝膠網絡的主要相互作用為黏性末端間的連接,我們僅保留這些DNA分子的黏性末端作為粗?；Ｗ?其間有相互作用.水凝膠中的水分子被處理為隱性溶劑.最后得到的粗?；Ｐ腿鐖D1所示.

Y形分子粗?；癁橛扇齻€粒子(記為A)通過彈簧連接組成的等邊三角形,邊長為6σ.Y形分子本身可能具有一定程度的形變(例如,三條臂之間夾角的相對變化),這可以近似描述為三角形邊長的彈性伸縮.連接子粗粒化為兩個粒子(記為B)連接成線狀分子,考慮到DNA鏈的伸縮較小,可視為由硬桿或硬彈簧連接,鍵長為10σ.鍵長數值參考實驗中的設置,做了較小的近似.σ約為4 bp,更多粗?；瘏翟O置見后文.

圖1 (網刊彩色)粗粒化模型示意圖Fig.1.(color online)Schematics of the coarse-grained model

我們分別模擬了粗?；Ｗ娱g鍵連接為硬彈簧鍵和硬桿的情況,其中彈簧鍵的彈性系數為k=0.5ε/σ2,ε為能量參數.模擬表明這兩種情況對于感興趣的物理量,給出的結果基本一致(見附錄圖A1),因此下文中我們只展示彈簧鍵的結果.

組成Y形分子的單個粗粒化粒子A,其質量為Y分子單臂質量,約18對堿基質量;組成連接子的單個粗?；Ｗ覤,其質量為連接子的一半質量,約20對堿基質量.考慮到系統(tǒng)為過阻尼,慣性效應較小,可將A,B分子質量視為相等,均設為m0.

2.2 粗?；Ｗ又g的相互作用勢

粗?；Ｗ訛閟sDNA構成的黏性末端,粒子間的相互作用為堿基互補配對,具有如下性質:1)不同類型的黏性末端相互吸引配對,同一種黏性末端之間不相互吸引配對;2)一個黏性末端不能同時與多個黏性末端配對,即飽和性.

我們使用Lennard-Jones(LJ)勢能[18]作為粗?；Ｗ又g的相互作用,對于相互作用的粒子其勢能函數為

其中,rij為i粒子和j粒子間距,超過截斷半徑的分子間相互作用力忽略不計;εij和σij表示粒子間的作用勢能量和距離參數.粗粒化粒子A,B之間的相互作用勢的能量為ε,距離參數為σ,截斷距離為2.5σ.這個勢能可以在一定程度上模擬兩種黏性末端之間的配對和斷開.但是,現(xiàn)有的模型難以滿足黏性末端的飽和性,需要通過勢能設計來實現(xiàn).

此處的飽和性問題等價于一個A類粒子的力程范圍內不能有兩個B類粒子同時處于穩(wěn)態(tài)(反之亦然),如圖2中的插圖所示.

圖2 A粒子與其兩側對稱分布的兩個B粒子組成的三體系統(tǒng)的總勢能曲線,本圖表明兩個B粒子的過飽和態(tài)幾乎不會出現(xiàn),插圖為A粒子附近兩個B粒子處于過飽和狀態(tài)時的示意圖Fig.2.Total potential energy of the three-body system consisting of an A particle and two B particles which are placed symmetrically around A.This graph indicates that the over-saturated state can hardly occur.Small f i gure shows schematics of two B particles simultaneously interacting with an A particle and forming an over-saturated state.

為了實現(xiàn)飽和性,我們引入了同種粒子間的LJ排斥勢,通過較大范圍的排斥使得同種粒子難以同時存在于一個較小的區(qū)域.排斥勢也同時實現(xiàn)了同種黏性末端不相互吸引配對的特性.此處我們設置同種粒子之間相互作用勢的能量為ε,距離參數為4σ,截斷距離為勢能極小值點,約為4.49σ,同時將勢能函數進行向上平移修正,使其在截斷點處的勢能為0.

假設該三體系統(tǒng)能夠處于穩(wěn)態(tài),則根據對稱性,兩個B粒子與A粒子距離相等.總勢能為

UBB為排斥勢,隨rBB增大而減小,則θ=180°時Utotal最小,此時Utotal隨rAB的分布曲線如圖2所示,極小值點能量約為?0.03ε.而模擬中的最低溫度為0.04T0=0.04ε/kB,此時熱漲落能量大于穩(wěn)態(tài)勢阱深度.因此,在模擬溫度范圍內(T=0.04T0—0.30T0),該三體態(tài)非常不穩(wěn)定,極易受熱漲落破壞.經統(tǒng)計,在實際模擬過程中,過飽和的鍵占整個系統(tǒng)的比例小于0.01%,如附錄圖A2所示,證明我們的模型在飽和性方面是非常有效的.

2.3 粗?；瘏档恼鎸崋挝粚?/h3>

粗粒化模型的長度單位σ、質量單位m0等基本來源于實驗中的DNA參數,為了模型的簡化進行了部分微調.

粗粒化溫度單位由能量標定,即T0=ε/kB.例如文獻[15]實驗中使用的標準8 bp長度黏性末端,根據其DNA序列,由Santalucia近鄰熱力學參數[19],可以近似估計300 K溫度(即室溫Tr)下黏性末端配對的自由能ε約為10kBTr—16kBTr,實際數值受濃度、離子環(huán)境等因素影響而波動.即對應于此8 bp黏性末端,其真實溫度300 K對應于粗粒化溫度T=0.06T00.1T0.模型參數及其賦值列于表1.

表1 模型參數及其賦值Table 1.List of the model parameters and their values.

所有分子動力學模擬都采用LAMMPS軟件通過粗?；鲞M行模擬.模擬過程主要使用Langevin動力學[20].模擬體系采用周期性邊界條件,主要采用邊長為100σ的立方體盒子,部分模擬使用八倍、十六倍等體積的盒子.Y形分子與連接子數量比例為1:1.5以保持兩者黏性末端總數量相等,從而達到較高的配對效率,該比例與實驗相符[15].粗?；Ｗ訑祻?00至14400不等,由系統(tǒng)大小和粒子數密度決定,粗?；到y(tǒng)粒子數密度與實驗中的真實濃度對比見表2.模擬溫度范圍為0.04T0—0.30T0.模擬步長為0.005t,溫度阻尼系數為0.5t.單個模擬軌道時長為1億步至10億步.

表2 粒子數密度與真實Y形分子濃度對比(密度單位:個/(100σ)3)Table 2.Correspondence between the concentration of the coarse-grained particles and the concentration of Y-scaffold blocks used in experiments.

3 結果與討論

3.1 二維體系

二維情況下的模擬對于了解體系結構與性質具有一定參考價值.我們模擬了邊長為400σ的正方形二維盒子,粒子數為2880.

模擬發(fā)現(xiàn),在多個溫度模擬下,二維水凝膠體系呈現(xiàn)出兩種不同的狀態(tài).如附錄圖A3所示,低溫態(tài)(T=0.06)粒子分布不均勻,出現(xiàn)粒子聚集和較大空隙,處于聚集態(tài);高溫態(tài)(T=0.20)粒子分布較均勻,處于液態(tài).

圖3 (網刊彩色)平均勢能與Y分子平均配對數隨溫度的變化Fig.3.(color online)Average potential energy and pair number per Y-scaffold at different temperatures.

為了更精確地描述該系統(tǒng)的凝膠化程度,我們引入Y分子的平均配對數作為序參量以描述系統(tǒng)的交聯(lián)度.Y分子含三個黏性末端,其最高配對數為三,而連接子與Y分子的總黏性末端數相等,因此Y分子平均配對數越接近三,系統(tǒng)交聯(lián)度越高,其凝膠化程度越高.

低溫聚集態(tài)和高溫液態(tài)對應的勢能曲線和Y分子平均配對數曲線如圖3所示.隨著溫度上升,體系的平均勢能上升,Y分子平均配對數下降.

3.2 三維體系

如圖4所示,與二維情況類似,粗?；M體系在低溫(T≤0.06T0)下形成多孔海綿狀結構,幾乎所有分子都處于交聯(lián)網絡中.

圖4 (網刊彩色)三維水凝膠體系模擬結構圖Fig.4.(color online)Structure of simulated threedimensional hydrogel system.

3.2.1 勢能與交聯(lián)密度

如圖5所示,隨著溫度上升,體系的平均勢能上升,平均配對數下降,且低溫下近飽和,這與二維體系一致.說明系統(tǒng)存在兩個可能的狀態(tài),低溫下凝膠化程度高,高溫下程度低.升降溫曲線的重合表明系統(tǒng)隨溫度的轉變沒有通常結晶與熔化過程中的熱滯現(xiàn)象.

為了考察上述結果是否依賴于模擬體系的大小,我們對密度同樣為200μM的不同體積的系統(tǒng)進行了模擬,得到其平均配對數如圖6所示.可以看出,不同體積下平均配對數曲線基本重合,沒有顯示出有限尺度效應,說明我們的模擬結果的確可以反映出宏觀體系的特征.

圖5 (網刊彩色)勢能與平均配對數曲線與溫度的關系Fig.5.(color online)Average potential energy and pair number per Y-scaffold at different temperatures.

圖6 (網刊彩色)相同密度、不同體積的系統(tǒng)中平均配對數與溫度的關系Fig.6.(color online)Pair number per Y-scaffold of systems with different sizes and the same density.

3.2.2 擴散行為

為了進一步描述系統(tǒng)在兩種狀態(tài)下的性質,我們考察了粗粒化粒子擴散行為的轉變.均方位移(msd)是描述系統(tǒng)擴散性質的物理量.其定義為

一般msd滿足:

通常α＜1對應于次擴散(sub-diffusion),α＞1對應于超擴散(super-diffusion),而α=1對應于正常擴散.

我們對不同溫度下系統(tǒng)的msd進行計算,并在雙對數坐標下進行擬合,擬合得到直線斜率即為α.擬合時主要取msd曲線后半段,粒子擴散達到穩(wěn)態(tài)的區(qū)域,大致為后50%的數據點,在雙對數坐標下對應于大約t＞105之后的曲線.雙對數坐標下的msd曲線如圖7所示.

圖7 (網刊彩色)系統(tǒng)msd隨時間變化的曲線(雙對數坐標),圖中兩條虛線斜率為1Fig.7.(color online)The mean squared displacement of the system as a function of time t under different temperatures(in double logarithmic coordinates).The slope of each dashed line is 1.

從圖7可以看到,高溫(T≥0.10T0)下系統(tǒng)的msd曲線基本平行,斜率接近1,說明此時為常規(guī)擴散,表明系統(tǒng)為液態(tài);低溫(T≤0.06T0)下系統(tǒng)的msd曲線擬合直線斜率明顯小于1,表明系統(tǒng)處于凝膠態(tài),擴散受到限制.

3.2.3 Voronoi單胞與局域密度分布

為了進一步研究這個粗?；到y(tǒng)的微觀結構,我們以粗粒化粒子的位置為種子點計算了該系統(tǒng)的Voronoi單胞的體積分布.Voronoi圖是一種基于距離的空間劃分方法[21,22],該方法將空間劃分為若干區(qū)域,使得每個區(qū)域內的點距其種子點的距離比其他種子點更近,每個區(qū)域為一個Voronoi單胞.Voronoi單胞的體積大小反映了種子點的局域密度分布.

如圖8所示,隨著溫度下降,Voronoi單胞體積分布曲線的峰左移且升高.同時,在平均體積線右側(V=1100—2000位置),分布曲線升高.這表明在低溫下,系統(tǒng)中大量粒子集中于高密度態(tài),同時還有相當一部分粒子處于低密度態(tài),形成了海綿狀結構.而高溫下體積分布趨于平均線,大體積單胞也較少.表明高溫下系統(tǒng)結構較為平均,接近液態(tài).

圖8 (網刊彩色)Voronoi單胞體積分布(歸一化)Fig.8.(color online)Volume distribution of Voronoi cells

Voronoi單胞體積也可以用于描述凝膠中的孔洞尺寸,可以大致預測凝膠中的孔洞體積約為10000σ3—100000σ3(對應于實際尺寸約為30000—300000 nm3),將來或可與實驗測出的數值進行半定量對比驗證.

3.2.4 密度的影響

為了與實驗[15]進行對比,我們也研究了粒子密度對體系凝膠化程度的影響.在實驗中,DNA水凝膠模塊的濃度增加會使體系儲能模量增大,對應著凝膠化程度的增加.

圖9 (網刊彩色)不同模塊密度下平均配對數隨溫度的變化Fig.9.(color online)Pair number per Y-scaffold at different temperatures for systems with different densities.

我們的模擬結果如圖9所示,平均配對數與密度成正相關.相同溫度下,高密度系統(tǒng)其交聯(lián)程度更高,更容易凝膠化.溫度足夠低時,四條曲線基本歸于一點,系統(tǒng)的平均配對數基本飽和.這與實驗結果相符合.

3.3 微流變與力學性質

為了進一步研究系統(tǒng)的流變學性質,我們使用微流變學[23,24]方法對該系統(tǒng)進行研究.微流變學方法通過在系統(tǒng)中加入少量示蹤粒子,探測示蹤粒子的msd等參量,進而得到系統(tǒng)的一些性質.

此處我們使用的示蹤粒子為單個粗?；Ｗ覤,即組成連接子的粗?；Ｗ?分別在不同的溫度和不同大小的拉力下進行了模擬.

3.3.1 不同溫度下的被動微流變

被動微流變[25]中,對于示蹤粒子不施加外力.對于不同溫度下的系統(tǒng)進行多組模擬(每個條件下100條軌道),并對示蹤粒子進行了軌道平均.最終求出粒子的msd與時間的關系,如圖10所示.此外,我們仍按照前述方法,取后50%(即t＞105的穩(wěn)態(tài)區(qū)域)的數據點進行擬合,求出各線性區(qū)對應的指數α,得到其隨溫度的變化,如圖10中的插圖所示.

圖10 (網刊彩色)示蹤粒子的msd隨時間變化(雙對數坐標),圖中兩條虛線斜率為1,插圖為指數α與溫度的關系Fig.10.(color online)The mean squared displacement of the tracer as a function of time.The slope of each dashed line is 1.Small f i gure shows α for different temperatures

從圖10可以看到T≤0.08T0時曲線斜率α小于1,為次擴散,系統(tǒng)處于凝膠態(tài);T≥0.10T0時,曲線斜率α在1附近,為常規(guī)擴散,系統(tǒng)處于液態(tài).這表明系統(tǒng)的凝膠液態(tài)轉變溫度在0.06T0—0.10T0之間,與上文系統(tǒng)整體msd的結果大致符合.

3.3.2 不同拉力下的主動微流變

主動微流變[26]中,對于示蹤粒子施加外力以觀察其軌跡.示蹤粒子在拉力方向上的msd為

其中〈Δx2(t)〉=〈(x(t)?x(0))〉2為通常定義中的msd,表征粒子的擴散作用.但是此處的Δx(t)不僅包含粒子擴散,還包含拉力導致的相對t=0時的粒子移動.為了去除這一影響,需進行一個修正,減去〈Δx(t)〉2=〈(x(t)?x(0))〉2這一由外力拖拽造成的粒子偏移.

我們在示蹤粒子上加上一個x方向的牽引力,計算其牽引方向的msd,以探究其在非平衡狀態(tài)下的力學性質.此處我們先對每條軌道進行了窗口平均計算得出各自的msd,再進行軌道平均,得到如圖11所示的雙對數坐標下的msd.由圖11可以發(fā)現(xiàn),在較長的關聯(lián)時間下,隨著外力的增加,msd曲線表現(xiàn)出超線性行為.對長關聯(lián)時間下的區(qū)域即msd的后半段(t＞103.5)進行擬合,得到擴散指數α隨溫度的變化如圖11中的插圖所示.

圖11 (網刊彩色)示蹤粒子在力方向上的msd隨時間變化(雙對數坐標),插圖為指數α與外力的關系Fig.11.(color online)The mean squared displacement in the force direction of the tracer as a function of time.Small f i gure shows α for different forces.

由圖11可以看到,F≤ 0.005ε/σ時,α ＜1,示蹤粒子近似表現(xiàn)為擴散行為;F≥0.02ε/σ時,α＞1,示蹤粒子表現(xiàn)為超擴散.這表明示蹤粒子的外力增加到一個值(0.02ε/σ)時,會很快由擴散行為變?yōu)榻苿蛩龠\動.外力導致的示蹤粒子擴散超擴散轉變的臨界力在0.005ε/σ—0.02ε/σ之間.

3.4 與實驗的半定量對比

本文的粗?；Ｐ皖A言體系發(fā)生“凝膠-液態(tài)”轉變的溫度區(qū)間大致為

在黏性末端長度變化范圍較小的情況下,該模型對應的水凝膠的轉變溫度正比于黏性末端結合能.因此我們可以根據黏性末端的DNA序列預測水凝膠真實的“凝膠-溶液態(tài)”轉變溫度.

例如,對于文獻[15]實驗中采用的標準8 bp黏性末端,其結合能約為10kBTr—16kBTr,粗粒化溫度0.06T0—0.10T0對應的真實溫度約為300 K.而根據文獻[15]圖2(d)的結果,該DNA水凝膠的凝膠溶液態(tài)轉變溫度在25—50°C范圍,即298—323 K內.模擬與實驗基本符合.

文獻[15]的圖3比較了多組黏性末端的轉變溫度,表明黏性末端結合能越高,其轉變溫度越高.這與我們的結論一致.

轉變溫度與黏性末端結合能的關系可以對DNA水凝膠序列設計提供一定的指導.例如,對于模型對應的長度范圍內的DNA水凝膠,為了使其在300 K下成膠,其結合能不能低于10kBTr—12.5kBTr.

另外,根據我們的模擬,DNA水凝膠濃度越高,則相同溫度下其交聯(lián)程度越高,凝膠化程度越大.在文獻[15]圖S2中,水凝膠儲能模量隨濃度增加而增大,同樣表明凝膠化程度與濃度成正相關.本文模擬與實驗在定性上一致.

根據凝膠-液態(tài)轉變溫度范圍,我們可以推測出該DNA水凝膠成膠的交聯(lián)度條件大致為Y分子平均配對數大于2.7,即90%以上的黏性末端處于配對交聯(lián)狀態(tài).因此對于不同濃度下的水凝膠,可以通過交聯(lián)度預測其成膠溫度.這將來或可與實驗測出的數值進行半定量對比驗證.

4 結 論

我們設計了一種適用于DNA水凝膠的粗?；Ｐ?并針對一種DNA水凝膠進行了分子動力學模擬驗證.該模型對凝膠分子模塊進行了極大的簡化以減少計算量,并通過LJ排斥勢的設計實現(xiàn)了DNA配對的飽和性.

針對該DNA水凝膠的模擬結果發(fā)現(xiàn),該系統(tǒng)在不同溫度下存在較均勻的液體態(tài)和海綿狀的凝膠態(tài)兩種狀態(tài).兩種狀態(tài)間隨溫度改變發(fā)生連續(xù)轉變,且轉變受密度因素影響.我們深入研究了兩種狀態(tài)下系統(tǒng)的結構和擴散性質,并做出了結構上的預測.此外,通過微流變方法進一步得出系統(tǒng)的兩態(tài)轉變溫度約為0.06ε/kB—0.10ε/kB,并與文獻[15]中的一些結果進行對照,兩者基本一致.這一轉變溫度與黏性末端結合能的關系,可以對DNA水凝膠的設計和轉變溫度控制提供一定的指導.

另外,我們也預測該水凝膠成膠條件為黏性末端配對交聯(lián)比例大致超過90%.外力引起的粒子擴散超擴散轉變的臨界力范圍為0.005ε/σ—0.02ε/σ.

基于該粗?；Ｐ瓦M行模擬的結果與合作組的實驗結果基本相符,證明此模型對DNA水凝膠的模擬較為合理.對于更多類型的DNA水凝膠,或使用DNA鏈作為連接基團的膠體系統(tǒng),該模型由于其計算量小又滿足飽和性等優(yōu)勢,在模擬上具有潛在的應用前景.

附錄A

A1彈簧鍵和硬桿的對比

在模擬中,彈簧鍵使用諧振勢,為

其中r0為平衡位置距離,k為彈性系數.對于線狀連接子,彈簧鍵對應的是線狀雙鏈DNA的伸縮與彎曲造成的黏性末端之間距離的變化.本文所研究的雙鏈DNA長度遠小于其駐留長度,其伸縮和彎曲較小,可視為硬桿,故而連接子使用彈性系數較大的彈簧鍵,彈性系數k=0.5ε/σ2.對于Y形分子,彈簧鍵(鍵長為6σ)對應的是Y分子自身形變造成的黏性末端間的距離變化,主要來源為三條剛性臂間夾角在120°附近的較小變化.因此,同一個Y分子中的兩個A粒子之間的距離應該小于Y分子臂長的兩倍(約6.93σ),而且不會過于靠近.實際上,在我們的模型中,為了滿足黏性末端相互作用的飽和性,在組成Y分子的粗粒化粒子A之間具有較長距離的排斥勢(作用范圍4.49σ),這實際使得同一個Y分子中的兩個A粒子之間的距離也通常大于這個排斥勢范圍.因此彈簧鍵長變化范圍約為4.49σ—6.93σ,實際效果接近硬彈簧,因此我們采用較大的彈性系數,彈性系數k=0.5ε/σ2.

對于分子鍵為硬桿的模型,在Lammps軟件模擬中,為了實現(xiàn)鍵長不變,需要通過shake算法[27]對鍵長進行約束.該算法通過每一步對原子施加附加約束力,通過迭代計算使鍵長重置為平衡長度.

我們在彈簧和硬桿兩種模型下分別進行了模擬,得到在(100σ)3大小,總粒子數為1800的NVT系統(tǒng)中,Y分子平均配對數隨溫度分布曲線,為了進行對比,我們同時模擬了Y分子的彈簧鍵分別為硬彈簧(k=0.5ε/σ2)和軟彈簧(k=0.05ε/σ2)時的情況作為參考,模擬結果如圖A1所示.

圖A1 (網刊彩色)分子鍵為硬彈簧(k=0.5)、軟彈簧(k=0.05)和硬桿三種情況下Y分子平均配對數與溫度的關系(誤差棒為標準誤差)Fig.A1.(color online)Pair number per Y-scaffold at different temperatures,simulated by hard spring model(k=0.5),soft spring model(k=0.05)or stiffstick model

可以看到,硬桿和硬彈簧兩種分子鍵下,Y分子平均配對數基本重合,這表明硬彈簧和硬桿在系統(tǒng)交聯(lián)度等結構方面的影響基本相同.軟彈簧在低溫下Y分子平均配對數小于另外兩種模型.在實際模擬中,也多次出現(xiàn)彈簧鍵長過長的情況,偏離了實際系統(tǒng).可以認為彈性系數較小的軟彈簧模型不適合當前DNA水凝膠系統(tǒng),因此只采用硬彈簧和硬桿的結果.

A2過飽和概率統(tǒng)計

我們統(tǒng)計了系統(tǒng)在多個溫度下過飽和配對數占單種黏性末端總數的比例,每個溫度下統(tǒng)計2000幀,平均后得到過飽和比例如圖A2所示.

可以看到,在模擬的溫度范圍內,過飽和比例均小于0.01%.考慮到我們模擬的系統(tǒng)最大粒子數只有一萬至二萬,大部分模擬中粒子數不足一萬,可以認為我們的模型基本滿足了DNA鏈配對的飽和性.

圖A2 不同溫度下的過飽和比例(誤差棒為標準誤差)Fig.A2.The ratio of over-saturated pairs at different temperatures.

A3二維系統(tǒng)低溫與高溫態(tài)構象

模擬得到的二維系統(tǒng)構象如圖A3所示,可以看到低溫下大多數分子都處于連接狀態(tài),且單個Y分子所連接的連接子數目大多為2或3,出現(xiàn)分子聚集,系統(tǒng)形成密度不均勻的分布.高溫下,存在大量Y分子其相連的連接子數目小于2,系統(tǒng)形成密度較均勻的分布.

因為模型中不考慮DNA鏈的排斥體積,二維情況下可能存在少量鍵交疊的現(xiàn)象,而在三維情況下基本不會發(fā)生交疊.本文主要研究三維體系,因此不進一步考慮鍵交疊的影響.

[1]Foster J A,Steed J W 2010Angew.Chem.Int.Ed.49 6718

[2]Yu G,Yan X,Han C,Huang F 2013Chem.Soc.Rev.42 6697

[3]Topuz F,Okay O 2008Macromolecules41 8847

[4]Morán M C,Miguel M G,Lindman B 2010Soft Matter6 3143

[5]Um S H,Lee J B,Park N,Kwon S Y,Umbach C C,Luo D 2006Nat.Mater.5 797

[6]Angioletti-Uberti S,Mognetti B M,Frenkel D 2016PCCP18 6373

[7]Li C,Faulkner-Jones A,Dun A R,Jin J,Chen P,Xing Y,Yang Z,Li Z,Shu W,Liu D,Duncan R R 2015Angew.Chem.Int.Ed.54 3957

[8]Amiya T,Tanaka T 1987Macromolecules20 1162

[9]Topuz F,Okay O 2009Biomacromolecules10 2652

[10]Starr F W,Sciortino F 2006J.Phys.:Condens.Mater.18 L347

[11]Dans P D,Walther J,Gómez H,Orozco M 2016Curr.Opin.Struct.Biol.37 29

[12]Weiner S J,Kollman P A,Nguyen D T,Case D A 1986J.Comput.Chem.7 230

[13]Uusitalo J J,Ingólfsson H I,Akhshi P,Tieleman D P,Marrink S J 2015J.Chem.Theory Comput.11 3932

[14]Collepardo-Guevara R,Schlick T 2014Proc.Natl.Acad.Sci.USA111 8061

[15]Xing Y,Cheng E,Yang Y,Chen P,Zhang T,Sun Y,Yang Z,Liu D 2011Adv.Mater.23 1117

[16]Cheng E,Xing Y,Chen P,Yang Y,Sun Y,Zhou D,Xu L,Fan Q,Liu D 2009Angew.Chem.121 7796

[17]Ouldridge T E,Louis A A,Doye J P K 2010Phys.Rev.Lett.104 178101

[18]Lennard-Jones J E 1931Proc.Phys.Soc.43 461

[19]SantaLucia J 1998Proc.Natl.Acad.Sci.USA95 1460

[20]Schneider T,Stoll E 1978Phys.Rev.B17 1302

[21]Okabe A 1992Spatial Tessellations(New York:John Wiley&Sons)pp362–363

[22]Aurenhammer F 1991ACM Comput.Surv.23 345

[23]Winter D,Horbach J 2013J.Chem.Phys.138 12A512

[24]Mason T G,Weitz D 1995Phys.Rev.Lett.74 1250

[25]Mizuno D,Head D A 2008Macromolecules41 7194

[26]Choi S Q,Steltenkamp S,Zasadzinski J A,Squires T M 2011Nat.Commun.2 312

[27]Ryckaert J P,Ciccotti G,Berendsen H J C 1977J.Comput.Phys.23 327