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

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

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

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

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

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

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

計算化學; 熱力學性質; 隱溶劑模型; 包合自由能; 環糊精

環糊精（cyclodextrin，簡稱 CD）包合農藥分子是一類重要的主客體超分子。CD包合農藥后可以改變農藥分子的理化性質，在農藥的穩定性、生物活性、毒性、對疏水性農藥的增溶、光解催化作用和在手性農藥分離等方面有著廣泛的應用，而且CD還在農藥污染治理和農藥殘留檢測方面有著廣泛應用[1-3]。但是天然CD在應用中還有一定的局限性[4-5]，如缺少顯示電子轉移、光致變色等功能性基團，在紫外、熒光等光譜中是惰性的，無法借助各種光學儀器對主客體間的相互作用進行研究；在水中的溶解度較小，也限制了它的應用。對β-CD進行適當的修飾合成新型的性能優異的β-CD衍生物變得至關重要，常用的衍生物包括甲基化CD和羥丙基化CD等。

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

圖1 二甲四氯（MCPA）的分子結構

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

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

1 理論與方法

1.1 MD/QM/CSM方法

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

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

1.2 結構模型

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

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

圖3 MCPA進入環糊精空腔的初始取向示意圖

1.3 分子動力學模擬

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

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

MD模擬完成后，對軌跡中的構象進行聚類分析。以主客體分子的所有原子作為對構象計算RMS的基準，調節 RMS閾值，使得構象的類數為 7左右。

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

1.4 基于熱力學循環的自由能計算

采用熱力學循環將溶液中復合物體系的包合自由能變化分成兩部分，如式（1）。

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

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

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

2 結果與討論

2.1 分子對接

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

2.2 MD軌跡分析

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

圖4 AutoDock中3種CD與MCPA分子對接結果

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

2.3 代表性包合構象分析

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

2.4 溶液中包合自由能計算

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

圖6 從MD計算和聚類分析產生的最優包合構象

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

圖7 實驗值與理論計算值關系

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

3 結 論

本文利用最近發展的MD/QM/CSM方法計算3種 CD對 MCPA農藥分子的包合機理，得出以下結論。

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

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

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

[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]王兆守，劉麗花，邵宗澤，等. 農藥殘留檢測新方法研究進展[J]. 化工進展，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]楊建平，巨曉潔，褚良銀. 兩親性環糊精的合成研究進展[J]. 化工進展，2008，27（1）：21-25.

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

[7]鄧坤學，黃啟谷. 農藥用表面活性劑研究進展[J]. 化工進展，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，Fernando 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]王騰，蔡文生，邵學廣. 環糊精的理論計算[J]. 化學進展，2010，22（5）：803-811.

[11]宋樂新，王海名，楊燕. 鄰-甲氧基苯酚和α-，β-環糊精包合現象的理論與實驗研究[J]. 化學學報，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]張慶東，李玉星，王武昌. 四氫呋喃-甲烷水合物穩定性的分子動力學模擬[J]. 化工進展，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，Ferreira 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。

中央高校基本科研業務費專項項目（30915011314）。

周玉梅（1989—），女，碩士研究生。E-mail meiyuzhou@qq.com。聯系人：周寶晶，副教授，碩士生導師，研究方向為理論計算與應用化學。E-mail bzhou@njust.edu.cn。