MD/QM/CSM方法計算β-環(huán)糊精及其衍生物對二甲四氯的包合機(jī)理

周玉梅，周寶晶，聶雪玫，葉仁龍，貢雪東，朱衛(wèi)華，肖鶴鳴

（南京理工大學(xué)化工學(xué)院，分子與材料計算研究所，江蘇 南京 210094）

MD/QM/CSM方法計算β-環(huán)糊精及其衍生物對二甲四氯的包合機(jī)理

周玉梅，周寶晶，聶雪玫，葉仁龍，貢雪東，朱衛(wèi)華，肖鶴鳴

（南京理工大學(xué)化工學(xué)院，分子與材料計算研究所，江蘇 南京 210094）

利用β-環(huán)糊精（β-CD）及其衍生物控制二甲四氯農(nóng)藥分子的釋放是目前研究的熱點。本文使用一種最近發(fā)展的分子動力學(xué)/量子力學(xué)/連續(xù)介質(zhì)溶劑模型（MD/QM/CSM）方法計算β-CD及其甲基、羥丙基衍生物對二甲四氯（MPCA）農(nóng)藥分子的包合機(jī)理。揭示了溶劑效應(yīng)在決定對這3種CD相對包合能力時發(fā)揮了主導(dǎo)作用，而真空包合自由能順序與溶劑中的正好相反。算得的3個體系包合能力強(qiáng)弱排序為：DM-β-CD>HP-β-CD>β-CD，與實驗測量結(jié)果一致，且線性相關(guān)度達(dá)到R=0.99。結(jié)果表明：MD/QM/CSM方法對于計算不同主體的超分子復(fù)合物體系的相對穩(wěn)定性是可靠的。

計算化學(xué); 熱力學(xué)性質(zhì); 隱溶劑模型; 包合自由能; 環(huán)糊精

環(huán)糊精（cyclodextrin，簡稱 CD）包合農(nóng)藥分子是一類重要的主客體超分子。CD包合農(nóng)藥后可以改變農(nóng)藥分子的理化性質(zhì)，在農(nóng)藥的穩(wěn)定性、生物活性、毒性、對疏水性農(nóng)藥的增溶、光解催化作用和在手性農(nóng)藥分離等方面有著廣泛的應(yīng)用，而且CD還在農(nóng)藥污染治理和農(nóng)藥殘留檢測方面有著廣泛應(yīng)用[1-3]。但是天然CD在應(yīng)用中還有一定的局限性[4-5]，如缺少顯示電子轉(zhuǎn)移、光致變色等功能性基團(tuán)，在紫外、熒光等光譜中是惰性的，無法借助各種光學(xué)儀器對主客體間的相互作用進(jìn)行研究；在水中的溶解度較小，也限制了它的應(yīng)用。對β-CD進(jìn)行適當(dāng)?shù)男揎椇铣尚滦偷男阅軆?yōu)異的β-CD衍生物變得至關(guān)重要，常用的衍生物包括甲基化CD和羥丙基化CD等。

2-甲基-4-氯苯氧乙酸又稱二甲四氯（2-methyl-4-chlorophenoxy acetic acid，簡稱MCPA，圖1）屬于苯氧乙酸類除草劑，選擇性高，內(nèi)吸性好。二甲四氯作為除草劑一般以鹽或者酯的形式存在，因其成本低、速度快、無殘留、對后茬作物安全等優(yōu)勢，被廣泛使用于水稻、小麥、甘蔗、玉米等水旱田間的闊葉雜草和莎草的防除[6-7]。為了提高其水溶性，采用 CD及其衍生物對其進(jìn)行包合[8-9]。2014年Fernando等[9]研究表明，CD及其衍生物與MPCA以1∶1的化學(xué)計量包合時，天然β-CD、二甲基-β-環(huán)糊精（dimethyl-β-cyclodextrin，簡稱DM-β-CD）和羥丙基-β-環(huán)糊精（hydroxypropyl-β-cyclodextrin，簡稱HP-β-CD）都可以增加MPCA的水溶性，且對MCPA的增溶效果：DM-β-CD>HP-β-CD>β-CD。這些研究者們還對DM-β-CD/MCPA和β-CD/MCPA兩種包合物做了分子動力學(xué)模擬，但未深入計算包合機(jī)理。

圖1 二甲四氯（MCPA）的分子結(jié)構(gòu)

主客體包合能力的理論計算無論從學(xué)術(shù)還是應(yīng)用角度看都有很重要的意義。用于預(yù)測主客體包合能力的理論計算方法有十幾種[10-13]，但是這些理論計算方法的精確性依然有待進(jìn)一步改進(jìn)。與實驗值相比，計算結(jié)果要么是均方根（RMS）誤差較大，要么是相關(guān)線性度不高。從計算結(jié)果的線性相關(guān)度角度考慮，DFT/PCM方法較其他方法有顯著的優(yōu)勢。不同于其他方法需要對大量包合構(gòu)象進(jìn)行統(tǒng)計熱力學(xué)平均，DFT/PCM方法對主客體包合能力計算時只采用一個構(gòu)象。這個構(gòu)象可以是由分子對接得到的，分子對接由于受到力場的限制，對接獲得的構(gòu)象也會有一定的局限性，從而影響了量子力學(xué)（QM）計算的精確性。以DFT/PCM方法為依托，將取得代表性構(gòu)象的方法進(jìn)行改進(jìn)，設(shè)計了一個新的理論模型來計算環(huán)糊精及其衍生物包合客體分子的相對能力，即分子動力學(xué)/量子力學(xué)/連續(xù)介質(zhì)溶劑模型（MD/QM/CSM）方法。建立的模型不僅計算真空中的能量，還考察了體系的溶劑效應(yīng)，使得CD/客體分子包合自由能的計算更為完善。該模型應(yīng)用于CD包合酰胺類農(nóng)藥分子和硝基化合物，可以對復(fù)合體系的包合自由能進(jìn)行正確排序，并且獲得較高的線性度。

這里將該方法用來研究 3種 CD，即 β-CD、DM-β-CD、HP-β-CD對MPCA農(nóng)藥分子的包合能力。由于主體分子比客體分子大很多，研究不同的CD對客體分子的包合效果影響更具有挑戰(zhàn)性。本文的目標(biāo)是：①檢驗當(dāng)主體 CD分子不同時，MD/QM/CSM方法的可靠性；②在原子分子水平上揭示CD包合二甲四氯的包合機(jī)理。

1 理論與方法

1.1 MD/QM/CSM方法

圖2 MD/QM/CSM方法計算流程

圖2 是MD/QM/CSM方法的計算流程。首先，通過MD計算獲得主客體分子及其復(fù)合物在溶液中的大量熱力學(xué)構(gòu)象。對MD模擬平衡后的運(yùn)動軌跡進(jìn)行聚類分析，然后對最大兩個聚類的所有勢能極小值構(gòu)象進(jìn)行進(jìn)一步優(yōu)化，從中提取出勢能最低的構(gòu)象。將由此產(chǎn)生的最穩(wěn)定構(gòu)象稱為代表性包合構(gòu)象，并將其用于包合自由能計算。基于熱力學(xué)循環(huán)，包合自由能可以分為兩部分：真空中自由能的變化和溶劑效應(yīng)引起的自由能變化。為了減小計算量，采用半經(jīng)驗QM方法計算氣相包合過程的焓變、熵變。溶劑效應(yīng)采用單點 DFT連續(xù)介質(zhì)溶劑模型（CSM）方法計算。最后獲得溶液中總的包合自由能，用來判斷不同CD對MPCA包合能力的強(qiáng)弱。

1.2 結(jié)構(gòu)模型

主體分子選擇 3種，即 β-CD、DM-β-CD和HP-β-CD。β-CD的最初結(jié)構(gòu)是從蛋白質(zhì)數(shù)據(jù)庫獲取（PDB code：1DMB），DM-β-CD的初始結(jié)構(gòu)是在β-CD的基礎(chǔ)上在2位和6位上全部引入甲基，取代度為2[11]；HP-β-CD是不同取代方式的混合物。一般認(rèn)為β-CD的2位上的—OH最容易發(fā)生親電取代[14-15]。本文的計算中，HP-β-CD結(jié)構(gòu)是在 β-CD的2位上全部引入—CH2CHOHCH3，取代度為1。

將CD的質(zhì)心設(shè)為坐標(biāo)原點，Z軸正方向是從CD的大口端指向小口端的方向。MCPA分子進(jìn)入CD時，其主軸平行于Z軸方向，可能采用不同的初始取向（圖3）。為了對MPCA最佳包合取向進(jìn)行預(yù)判，先用AutoDock軟件對主客體進(jìn)行分子對接，然后將對接后數(shù)目最多，能量最低處的構(gòu)象取向作為包合體系的初始構(gòu)象的取向。

圖3 MCPA進(jìn)入環(huán)糊精空腔的初始取向示意圖

1.3 分子動力學(xué)模擬

MD模擬使用AMBER 12程序進(jìn)行。MCPA分子力場參數(shù)是采用 generalized Amber force field（gaff）力場和由Gaussian09在HF/6-31G（d）水平上計算獲得的 RESP電荷。3種 CD均采用q4md-CD分子力場參數(shù)。q4md分子力場是結(jié)合GLYCAM04和Amber99SB力場對環(huán)糊精衍生物進(jìn)行描述。

將包合物以及主客體小分子分別置于TIP3P型水盒子中。體系先用共軛梯度法和最速下降法將能量最小化，然后再進(jìn)行MD模擬，采用NPT系統(tǒng)。首先固定體系，采用Langevin升溫方法經(jīng)400ps將體系升溫到300K；然后在無約束恒溫、1atm壓強(qiáng)下模擬 2ns。與氫原子相連的鍵用 SHAKEMD 算法限制，非鍵相互作用的截斷距離設(shè)為 9?。長程靜電相互作用采用Particle Mesh Ewald方法計算。整個模擬過程都采用周期性邊界條件，模擬的積分步長為 2fs。

MD模擬完成后，對軌跡中的構(gòu)象進(jìn)行聚類分析。以主客體分子的所有原子作為對構(gòu)象計算RMS的基準(zhǔn)，調(diào)節(jié) RMS閾值，使得構(gòu)象的類數(shù)為 7左右。

對聚類結(jié)果進(jìn)行分析，將前兩類中所有勢能局部最小的構(gòu)象提取出來。然后在采用周期性邊界條件下，使用共軛梯度和最速下降法，將它們能量最小化。產(chǎn)生的全局能量最低的構(gòu)象即為代表性構(gòu)象，被用作包合自由能量子力學(xué)計算的初始結(jié)構(gòu)。

1.4 基于熱力學(xué)循環(huán)的自由能計算

采用熱力學(xué)循環(huán)將溶液中復(fù)合物體系的包合自由能變化分成兩部分，如式（1）。

真空中包合自由能包括焓變和熵變兩部分，如式（2）。

由于半經(jīng)驗量子力學(xué)方法 PM3 較適合 CD體系且計算量適中[16]，真空中自由能變化采用該方法。

溶劑化能的計算由式（3）給出。

2 結(jié)果與討論

2.1 分子對接

AutoDock對接的結(jié)果表明：MCPA分子進(jìn)入3種CD的取向是不一樣的，見圖4。β-CD/MCPA和DM-β-CD/MCPA 體系采用圖 3取向Ⅰ，而HP-β-CD/MCPA體系采用圖3取向Ⅱ。構(gòu)建MD模擬的初始構(gòu)象時，將客體小分子放在CD腔外略靠大口端處。小分子取向和AutoDock給出的結(jié)果保持一致。

2.2 MD軌跡分析

MD模擬3種CD包合體系的主客體間質(zhì)心距離隨時間的變化曲線如圖5所示。在T≈500ps前后，客體分子進(jìn)入CD空腔，之后一直被包合在腔內(nèi)。從圖 5可以看出，MD運(yùn)動趨于穩(wěn)定之后，β-CD/ MCPA和DM-β-CD/MCPA這兩種包合體系的主客體小分子間質(zhì)心距離都是在 1.0左右上下波動，HP-β-CD/MCPA包合體系的主客體小分子間質(zhì)心距離在2.0左右上下波動。MCPA分子進(jìn)入3種CD深度波動大小順序為β-CD

圖4 AutoDock中3種CD與MCPA分子對接結(jié)果

圖5 包合體系的主客體間質(zhì)心距離隨時間的變化

2.3 代表性包合構(gòu)象分析

對MD軌跡進(jìn)行聚類之后，提取出3種CD包合MPCA的代表性構(gòu)象（圖6）。圖6中，左圖是CD的大口端朝上，右圖是從CD的小口端看向大口端的視圖。從圖6可以看出，MCPA進(jìn)入β-CD和DM-β-CD，角度會發(fā)生傾斜，不是垂直進(jìn)入，但基本是從氯原子端進(jìn)入。β-CD的剛性比較好，空腔變化不明顯。對于DM-β-CD/MCPA體系，DM-β-CD的6位上的甲基取代基會向空腔外伸展，DM-β-CD的空腔在變扁的同時小口端又變寬。MCPA分子進(jìn)入HP-β-CD的角度會發(fā)生很大的變化，從初始羧基端進(jìn)入變成了MCPA分子橫在腔內(nèi)，大口端2位上的羥丙基顯著向外伸展，這就導(dǎo)致了形成包合物后HP-β-CD的空腔變扁。與圖4相比較，HP-β-CD變化很大，這是因為羥丙基和水溶液之前較強(qiáng)的相互作用導(dǎo)致的。

2.4 溶液中包合自由能計算

表1為3種CD對MCPA的包合自由能計算的結(jié)果。計算結(jié)果 ΔGcal為正，原因可能在于：①真空中的熵過大，研究表明真空的熵值與溶劑中的熵值有所差別[17]；②在使用CSM計算溶劑效應(yīng)時使用不同的原子半徑、不同的算法所得出的結(jié)果都會有差別[12]。所以計算得出的包合自由能的絕對值并不一定可靠，相對值才具有價值。總的包合自由能變化的相對大小順序為DM-β-CDHP-β-CD>β-CD，這與實驗結(jié)果一致，見表 1。對理論計算與實驗包合自由能值進(jìn)行線性擬合（圖7），兩者之間存在良好的線性相關(guān)，且線性相關(guān)系數(shù)R達(dá)到0.99。這表明用MD/QM/CSM方法預(yù)測不同改性CD對客體分子包合能力的相對大小是可行的。

圖6 從MD計算和聚類分析產(chǎn)生的最優(yōu)包合構(gòu)象

表1 3種CD/MCPA體系真空中的熱力學(xué)量、溶劑效應(yīng)及總的包合自由能變化 單位：kcal·mol?1

圖7 實驗值與理論計算值關(guān)系

進(jìn)一步對理論計算真空包合自由能與實驗包合自由能值進(jìn)行了線性擬合（圖7），結(jié)果表明真空包合自由能的順序和溶液中的正好相反。在真空中，HP-β-CD的變形最大，相應(yīng)的熵增加最為顯著，但主客體的相互吸引也最強(qiáng)。而DM-β-CD與MPCA的相互吸引最弱，但熵增加最小。考慮到溶劑效應(yīng)后，總的自由能變化才能夠與實驗相吻合。所以僅有真空包合自由能并不能正確排序，只有綜合考慮了真空包合自由能和溶劑效應(yīng)才能可靠地預(yù)測出不同CD的相對包合能力。最近對CD復(fù)合物體系的3個理論研究工作也得出了同樣的結(jié)論[18-19]。此外，溶劑效應(yīng)在決定對這3種CD相對包合能力時發(fā)揮了主導(dǎo)作用。

3 結(jié) 論

本文利用最近發(fā)展的MD/QM/CSM方法計算3種 CD對 MCPA農(nóng)藥分子的包合機(jī)理，得出以下結(jié)論。

（1）將MD/QM/CSM方法用于不同種CD對農(nóng)藥分子的包合研究，得到與實驗一致的結(jié)果，這表明此方法對于計算不同主體的超分子復(fù)合物體系的相對穩(wěn)定性是可靠的。

（2）MD/QM/CSM方法不僅可以應(yīng)用于CD在環(huán)境修復(fù)、醫(yī)藥行業(yè)、手性拆分等方面的研究，還可以應(yīng)用于其他超分子化學(xué)體系，如環(huán)番、杯芳烴、葫蘆脲等超分子體系。

（3）MD/QM/CSM方法還可以指導(dǎo)改性CD新材料的設(shè)計、合成和應(yīng)用，通過理論計算為實驗提供理論指導(dǎo)，可以減小難度、減少工作量、降低成本。

[1]Flaherty R J，Nshime B，DeLaMarre M，et al. Cyclodextrins as complexation and extraction agents for pesticides from contaminated soil[J]. Chemosphere，2013，91（7）：912-920.

[2]Bian H，Chen J，Cai X，et al. Inclusion complex of butachlor with β-cyclodextrin：Characterization，solubility，and speciation-dependent adsorption[J]. Journal of Agricultural and Food Chemistry，2009，57（16）：7453-7458.

[3]王兆守，劉麗花，邵宗澤，等. 農(nóng)藥殘留檢測新方法研究進(jìn)展[J]. 化工進(jìn)展，2008，27（9）：1370-1374.

[4]Del Valle E M. Cyclodextrins and their uses：A review[J]. Process Biochemistry，2004，39（9）：1033-1046.

[5]楊建平，巨曉潔，褚良銀. 兩親性環(huán)糊精的合成研究進(jìn)展[J]. 化工進(jìn)展，2008，27（1）：21-25.

[6]陳嵐，權(quán)宇珩，李志勇. 臭氧降解土壤 2，4-D 污染物的動力學(xué)分析[J]. 化工學(xué)報，2014，65（8）：3255-3260.

[7]鄧?yán)W(xué)，黃啟谷. 農(nóng)藥用表面活性劑研究進(jìn)展[J]. 化工進(jìn)展，2009，28（12）：2199-2204.

[8]Garrido E M，Santos M，Silva P，et al. Host-guest complexes of phenoxy alkyl acid herbicides and cyclodextrins. MCPA and β-cyclodextrin[J]. Journal of Environmental Science and Health，Part B，2012，47（9）：869-875.

[9]Jorge G，F(xiàn)ernando C，Manuel M F，et al. Microencapsulation of herbicide MCPA with native β-cyclodextrin and its methyl and hydroxypropyl derivatives： An experimental and theoretical investigation[J]. Journal of Molecular Structure，2014（1061）：76-81.

[10]王騰，蔡文生，邵學(xué)廣. 環(huán)糊精的理論計算[J]. 化學(xué)進(jìn)展，2010，22（5）：803-811.

[11]宋樂新，王海名，楊燕. 鄰-甲氧基苯酚和α-，β-環(huán)糊精包合現(xiàn)象的理論與實驗研究[J]. 化學(xué)學(xué)報，2007，65（16）：1593-1599.

[12]Muddana H S，Varnado C D，Bielawski C W，et al. Blind prediction of host-guest binding affinities：A new SAMPL3 challenge[J]. Journal of Computer-Aided Molecular Design，2012，26（5）：475-487.

[13]張慶東，李玉星，王武昌. 四氫呋喃-甲烷水合物穩(wěn)定性的分子動力學(xué)模擬[J]. 化工進(jìn)展，2014，33（s1）：98-105.

[14]Zhao M，Hao A，Li J，et al. New cyclomaltoheptaose (β-cyclodextrin) derivative 2-O-(2-hydroxybutyl) cyclomaltoheptaose：Preparation and its application for the separation of enantiomers of drugs by capillary electrophoresis[J]. Carbohydrate Research，2005，340（8）：1563-1565.

[15]Yong C W，Washington C，Smith W. Structural behaviour of 2-hydroxypropyl-β-cyclodextrin in water：Molecular dynamics simulation studies[J]. Pharmaceutical Research，2008，25（5）：1092-1099.

[16]Sayede A，F(xiàn)erreira M，Bricout H，et al. Interaction of water-soluble triphenylphosphines with β-cyclodextrin：A quantum chemistry study[J]. Journal of Physical Organic Chemistry，2011，24（12）：1129-1135.

[17]Kua J，Krizner H E，De Haan D O. Thermodynamics and kinetics of imidazole formation from glyoxal，methylamine，and formaldehyde：A computational study[J]. The Journal of Physical Chemistry A，2011，115（9）：1667-1675.

[18]Zhang H，Tan T，Hetényi C，et al. Quantification of solvent contribution to the stability of noncovalent complexes[J]. Journal of Chemical Theory and Computation，2013，9（10）：4542-4551.

[19]Wickstrom L，He P，Gallicchio E，et al. Large scale affinity calculations of cyclodextrin host-guest complexes：Understanding the role of reorganization in the molecular recognition process[J]. Journal of Chemical Theory and Computation，2013，9（7）：3136-3150.

A theoretical study on the microencapsulation of herbicide MCPA with native β-cyclodextrin and its derivatives by a molecular dynamics/quantum mechanics/continuum solvent model approach

ZHOU Yumei，ZHOU Baojing，NIE Xuemei，YE Renlong，GONG Xuedong，ZHU Weihua，XIAO Heming
（Computational Institute for Molecules and Materials，School of Chemical Engineering，Nanjing University of Science and Technology，Nanjing 210094，Jiangsu，China）

Microencapsulation is an effective technology to reduce the environmental losses of pesticides during their use. This paper reports the complexation mechanism of chlorophenoxy herbicide MCPA with native β-cyclodextrin (CD) and its methyl and hydroxypropyl derivatives by a molecular dynamics/quantum mechanics/continuum solvent model (MD/QM/CSM) approach developed recently. The results reveal that the solvation effect plays a dominant role in determining the relative stabilities of the three CD/MPCA inclusion complexes. However，the results in vacuum environment are opposite to those in solvent media. The calculated relative binding affinities of the three systems in aqueous solution decrease in the order：DM-β-CD>HP-β-CD>β-CD，which agrees with the experimental results，with a correlation coefficient R=0.99. This demonstrates that the MD/QM/CSM method is able to determine the relative stability of host-guest systems with different host molecule.

computational chemistry; thermodynamic properties; continuum solvent modeling; binding affinity; cyclodextrin

TQ 028.8

A

1000-6613（2015）12-4185-06

10.16085/j.issn.1000-6613.2015.12.009

2015-03-30；修改稿日期：2015-04-26。

中央高校基本科研業(yè)務(wù)費(fèi)專項項目（30915011314）。

周玉梅（1989—），女，碩士研究生。E-mail meiyuzhou@qq.com。聯(lián)系人：周寶晶，副教授，碩士生導(dǎo)師，研究方向為理論計算與應(yīng)用化學(xué)。E-mail bzhou@njust.edu.cn。