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

1)(中國(guó)科學(xué)院大學(xué)物理科學(xué)學(xué)院,北京 100049)

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

1 引 言

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

迄今為止,文獻(xiàn)中已經(jīng)報(bào)道了多種DNA水凝膠[8,9].但對(duì)其結(jié)構(gòu)與性能的關(guān)系,特別是對(duì)凝膠網(wǎng)絡(luò)的分子化學(xué)結(jié)構(gòu)、分子鏈性質(zhì)、拓?fù)淇锥唇Y(jié)構(gòu)、交聯(lián)密度與強(qiáng)度等結(jié)構(gòu)參數(shù)及其與凝膠宏觀力學(xué)性質(zhì)之間的關(guān)系的了解還十分缺乏,主要是缺少對(duì)相應(yīng)體系的理論建模和實(shí)驗(yàn)測(cè)量的比較研究[10].因此,構(gòu)建DNA水凝膠的微觀模型并計(jì)算模擬其分子結(jié)構(gòu)和宏觀性質(zhì)的功能關(guān)系,對(duì)理解和設(shè)計(jì)DNA超分子水凝膠具有重要意義.

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

本文以文獻(xiàn)[15]報(bào)道的一種純DNA自組裝水凝膠為例,建立粗?；Ｐ?模擬研究該DNA水凝膠的分子模塊、相互作用等微觀特征與凝膠形成、結(jié)構(gòu)及性質(zhì)的關(guān)系,得到了與實(shí)驗(yàn)半定量或定性相符的結(jié)果.這個(gè)工作不僅有助于理解與優(yōu)化設(shè)計(jì)該類DNA水凝膠,也為更普遍地構(gòu)建合理有效的粗?；Ｐ鸵匝芯扛鞣NDNA組裝體提供了經(jīng)驗(yàn).

2 模型與方法

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

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

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

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

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

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

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

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

2.2 粗?；Ｗ又g的相互作用勢(shì)

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

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

其中,rij為i粒子和j粒子間距,超過截?cái)喟霃降姆肿娱g相互作用力忽略不計(jì);εij和σij表示粒子間的作用勢(shì)能量和距離參數(shù).粗?；Ｗ覣,B之間的相互作用勢(shì)的能量為ε,距離參數(shù)為σ,截?cái)嗑嚯x為2.5σ.這個(gè)勢(shì)能可以在一定程度上模擬兩種黏性末端之間的配對(duì)和斷開.但是,現(xiàn)有的模型難以滿足黏性末端的飽和性,需要通過勢(shì)能設(shè)計(jì)來實(shí)現(xiàn).

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

圖2 A粒子與其兩側(cè)對(duì)稱分布的兩個(gè)B粒子組成的三體系統(tǒng)的總勢(shì)能曲線,本圖表明兩個(gè)B粒子的過飽和態(tài)幾乎不會(huì)出現(xiàn),插圖為A粒子附近兩個(gè)B粒子處于過飽和狀態(tài)時(shí)的示意圖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.

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

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

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

2.3 粗?；瘏?shù)的真實(shí)單位對(duì)應(yīng)

粗粒化模型的長(zhǎng)度單位σ、質(zhì)量單位m0等基本來源于實(shí)驗(yàn)中的DNA參數(shù),為了模型的簡(jiǎn)化進(jìn)行了部分微調(diào).

粗?；瘻囟葐挝挥赡芰繕?biāo)定,即T0=ε/kB.例如文獻(xiàn)[15]實(shí)驗(yàn)中使用的標(biāo)準(zhǔn)8 bp長(zhǎng)度黏性末端,根據(jù)其DNA序列,由Santalucia近鄰熱力學(xué)參數(shù)[19],可以近似估計(jì)300 K溫度(即室溫Tr)下黏性末端配對(duì)的自由能ε約為10kBTr—16kBTr,實(shí)際數(shù)值受濃度、離子環(huán)境等因素影響而波動(dòng).即對(duì)應(yīng)于此8 bp黏性末端,其真實(shí)溫度300 K對(duì)應(yīng)于粗?；瘻囟萒=0.06T00.1T0.模型參數(shù)及其賦值列于表1.

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

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

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

3 結(jié)果與討論

3.1 二維體系

二維情況下的模擬對(duì)于了解體系結(jié)構(gòu)與性質(zhì)具有一定參考價(jià)值.我們模擬了邊長(zhǎng)為400σ的正方形二維盒子,粒子數(shù)為2880.

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

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

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

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

3.2 三維體系

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

圖4 (網(wǎng)刊彩色)三維水凝膠體系模擬結(jié)構(gòu)圖Fig.4.(color online)Structure of simulated threedimensional hydrogel system.

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

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

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

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

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

3.2.2 擴(kuò)散行為

為了進(jìn)一步描述系統(tǒng)在兩種狀態(tài)下的性質(zhì),我們考察了粗?；Ｗ訑U(kuò)散行為的轉(zhuǎn)變.均方位移(msd)是描述系統(tǒng)擴(kuò)散性質(zhì)的物理量.其定義為

一般msd滿足:

通常α＜1對(duì)應(yīng)于次擴(kuò)散(sub-diffusion),α＞1對(duì)應(yīng)于超擴(kuò)散(super-diffusion),而α=1對(duì)應(yīng)于正常擴(kuò)散.

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

圖7 (網(wǎng)刊彩色)系統(tǒng)msd隨時(shí)間變化的曲線(雙對(duì)數(shù)坐標(biāo)),圖中兩條虛線斜率為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,說明此時(shí)為常規(guī)擴(kuò)散,表明系統(tǒng)為液態(tài);低溫(T≤0.06T0)下系統(tǒng)的msd曲線擬合直線斜率明顯小于1,表明系統(tǒng)處于凝膠態(tài),擴(kuò)散受到限制.

3.2.3 Voronoi單胞與局域密度分布

為了進(jìn)一步研究這個(gè)粗?；到y(tǒng)的微觀結(jié)構(gòu),我們以粗?；Ｗ拥奈恢脼榉N子點(diǎn)計(jì)算了該系統(tǒng)的Voronoi單胞的體積分布.Voronoi圖是一種基于距離的空間劃分方法[21,22],該方法將空間劃分為若干區(qū)域,使得每個(gè)區(qū)域內(nèi)的點(diǎn)距其種子點(diǎn)的距離比其他種子點(diǎn)更近,每個(gè)區(qū)域?yàn)橐粋€(gè)Voronoi單胞.Voronoi單胞的體積大小反映了種子點(diǎn)的局域密度分布.

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

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

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

3.2.4 密度的影響

為了與實(shí)驗(yàn)[15]進(jìn)行對(duì)比,我們也研究了粒子密度對(duì)體系凝膠化程度的影響.在實(shí)驗(yàn)中,DNA水凝膠模塊的濃度增加會(huì)使體系儲(chǔ)能模量增大,對(duì)應(yīng)著凝膠化程度的增加.

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

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

3.3 微流變與力學(xué)性質(zhì)

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

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

3.3.1 不同溫度下的被動(dòng)微流變

被動(dòng)微流變[25]中,對(duì)于示蹤粒子不施加外力.對(duì)于不同溫度下的系統(tǒng)進(jìn)行多組模擬(每個(gè)條件下100條軌道),并對(duì)示蹤粒子進(jìn)行了軌道平均.最終求出粒子的msd與時(shí)間的關(guān)系,如圖10所示.此外,我們?nèi)园凑涨笆龇椒?取后50%(即t＞105的穩(wěn)態(tài)區(qū)域)的數(shù)據(jù)點(diǎn)進(jìn)行擬合,求出各線性區(qū)對(duì)應(yīng)的指數(shù)α,得到其隨溫度的變化,如圖10中的插圖所示.

圖10 (網(wǎng)刊彩色)示蹤粒子的msd隨時(shí)間變化(雙對(duì)數(shù)坐標(biāo)),圖中兩條虛線斜率為1,插圖為指數(shù)α與溫度的關(guān)系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時(shí)曲線斜率α小于1,為次擴(kuò)散,系統(tǒng)處于凝膠態(tài);T≥0.10T0時(shí),曲線斜率α在1附近,為常規(guī)擴(kuò)散,系統(tǒng)處于液態(tài).這表明系統(tǒng)的凝膠液態(tài)轉(zhuǎn)變溫度在0.06T0—0.10T0之間,與上文系統(tǒng)整體msd的結(jié)果大致符合.

3.3.2 不同拉力下的主動(dòng)微流變

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

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

我們?cè)谑聚櫫Ｗ由霞由弦粋€(gè)x方向的牽引力,計(jì)算其牽引方向的msd,以探究其在非平衡狀態(tài)下的力學(xué)性質(zhì).此處我們先對(duì)每條軌道進(jìn)行了窗口平均計(jì)算得出各自的msd,再進(jìn)行軌道平均,得到如圖11所示的雙對(duì)數(shù)坐標(biāo)下的msd.由圖11可以發(fā)現(xiàn),在較長(zhǎng)的關(guān)聯(lián)時(shí)間下,隨著外力的增加,msd曲線表現(xiàn)出超線性行為.對(duì)長(zhǎng)關(guān)聯(lián)時(shí)間下的區(qū)域即msd的后半段(t＞103.5)進(jìn)行擬合,得到擴(kuò)散指數(shù)α隨溫度的變化如圖11中的插圖所示.

圖11 (網(wǎng)刊彩色)示蹤粒子在力方向上的msd隨時(shí)間變化(雙對(duì)數(shù)坐標(biāo)),插圖為指數(shù)α與外力的關(guān)系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ε/σ時(shí),α ＜1,示蹤粒子近似表現(xiàn)為擴(kuò)散行為;F≥0.02ε/σ時(shí),α＞1,示蹤粒子表現(xiàn)為超擴(kuò)散.這表明示蹤粒子的外力增加到一個(gè)值(0.02ε/σ)時(shí),會(huì)很快由擴(kuò)散行為變?yōu)榻苿蛩龠\(yùn)動(dòng).外力導(dǎo)致的示蹤粒子擴(kuò)散超擴(kuò)散轉(zhuǎn)變的臨界力在0.005ε/σ—0.02ε/σ之間.

3.4 與實(shí)驗(yàn)的半定量對(duì)比

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

在黏性末端長(zhǎng)度變化范圍較小的情況下,該模型對(duì)應(yīng)的水凝膠的轉(zhuǎn)變溫度正比于黏性末端結(jié)合能.因此我們可以根據(jù)黏性末端的DNA序列預(yù)測(cè)水凝膠真實(shí)的“凝膠-溶液態(tài)”轉(zhuǎn)變溫度.

例如,對(duì)于文獻(xiàn)[15]實(shí)驗(yàn)中采用的標(biāo)準(zhǔn)8 bp黏性末端,其結(jié)合能約為10kBTr—16kBTr,粗?；瘻囟?.06T0—0.10T0對(duì)應(yīng)的真實(shí)溫度約為300 K.而根據(jù)文獻(xiàn)[15]圖2(d)的結(jié)果,該DNA水凝膠的凝膠溶液態(tài)轉(zhuǎn)變溫度在25—50°C范圍,即298—323 K內(nèi).模擬與實(shí)驗(yàn)基本符合.

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

轉(zhuǎn)變溫度與黏性末端結(jié)合能的關(guān)系可以對(duì)DNA水凝膠序列設(shè)計(jì)提供一定的指導(dǎo).例如,對(duì)于模型對(duì)應(yīng)的長(zhǎng)度范圍內(nèi)的DNA水凝膠,為了使其在300 K下成膠,其結(jié)合能不能低于10kBTr—12.5kBTr.

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

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

4 結(jié) 論

我們?cè)O(shè)計(jì)了一種適用于DNA水凝膠的粗粒化模型,并針對(duì)一種DNA水凝膠進(jìn)行了分子動(dòng)力學(xué)模擬驗(yàn)證.該模型對(duì)凝膠分子模塊進(jìn)行了極大的簡(jiǎn)化以減少計(jì)算量,并通過LJ排斥勢(shì)的設(shè)計(jì)實(shí)現(xiàn)了DNA配對(duì)的飽和性.

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

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

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

附錄A

A1彈簧鍵和硬桿的對(duì)比

在模擬中,彈簧鍵使用諧振勢(shì),為

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

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

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

圖A1 (網(wǎng)刊彩色)分子鍵為硬彈簧(k=0.5)、軟彈簧(k=0.05)和硬桿三種情況下Y分子平均配對(duì)數(shù)與溫度的關(guān)系(誤差棒為標(biāo)準(zhǔn)誤差)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分子平均配對(duì)數(shù)基本重合,這表明硬彈簧和硬桿在系統(tǒng)交聯(lián)度等結(jié)構(gòu)方面的影響基本相同.軟彈簧在低溫下Y分子平均配對(duì)數(shù)小于另外兩種模型.在實(shí)際模擬中,也多次出現(xiàn)彈簧鍵長(zhǎng)過長(zhǎng)的情況,偏離了實(shí)際系統(tǒng).可以認(rèn)為彈性系數(shù)較小的軟彈簧模型不適合當(dāng)前DNA水凝膠系統(tǒng),因此只采用硬彈簧和硬桿的結(jié)果.

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

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

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

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

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

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

因?yàn)槟Ｐ椭胁豢紤]DNA鏈的排斥體積,二維情況下可能存在少量鍵交疊的現(xiàn)象,而在三維情況下基本不會(huì)發(fā)生交疊.本文主要研究三維體系,因此不進(jìn)一步考慮鍵交疊的影響.

[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