木質(zhì)素過(guò)氧化酶對(duì)雙酚A類化合物的分子動(dòng)力學(xué)模擬

摘要：為了解雙酚A類化合物（BPAs）與木質(zhì)素過(guò)氧化物酶（LiP）結(jié)合位點(diǎn)的結(jié)合構(gòu)象及其相互作用機(jī)理，采用SYBYL-X1.3中分子對(duì)接的方法對(duì)BPA、BPF、BPAF、BPQ、BPB、BPC、BPS、BPE、BPAP等具有突出特點(diǎn)的雙酚A類化合物與LiP進(jìn)行對(duì)接研究，并利用GROMACS軟件進(jìn)行分子動(dòng)力學(xué)模擬。通過(guò)對(duì)對(duì)接最穩(wěn)定化合物BPQ和常用化合物BPA的模擬參數(shù)具體分析，找到了LiP的活性位點(diǎn)為HIS39、MET172、ALA175、HIS176、ILE42、ALA179、PHE193、ARG43、ILE235和PHE46；找到了該酶被識(shí)別的氫鍵結(jié)合位點(diǎn)中的氨基酸殘基，BPQ通過(guò)氫鍵連接了PRO83.O（id141）、ASP183.OD2（id470）和ASN182. HD21（id1221），而B(niǎo)PA通過(guò)氫鍵連接了PRO83.O（id141）、ARG43.HH21（id818）。氨基酸殘基對(duì)LiP與雙酚A類小分子化合物的降解起到重要作用。

關(guān)鍵詞：木質(zhì)素過(guò)氧化物酶；雙酚A；分子對(duì)接；分子動(dòng)力學(xué)模擬

中圖分類號(hào)：Q554 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：0439-8114（2016）01-0220-06

DOI：10.14088/j.cnki.issn0439-8114.2016.01.058

雙酚A類化合物是一種類雌激素，在塑料制品生產(chǎn)過(guò)程中可作為阻燃劑、抗氧化劑、硬化防止劑、聚合抑制劑和增塑劑等[1]，在水杯、奶瓶、食品包裝、飲料容器、醫(yī)療設(shè)備上被廣泛使用，且越來(lái)越多地應(yīng)用在新產(chǎn)品的開(kāi)發(fā)上。由于其難降解，該類化合物會(huì)緩慢地釋放到周圍的環(huán)境中，并在環(huán)境中積累，通過(guò)空氣、水、土壤以及食品污染物（塑料垃圾）等途徑引起人體畸形、癌變、內(nèi)分泌失調(diào)、新陳代謝紊亂[2，3]，成為了繼“溫室效應(yīng)”，“臭氧層空洞”后第三大威脅人類健康的熱點(diǎn)。因此，尋找并用正確的方法處理環(huán)境中的雙酚A類化合物已迫不及待。

作為木質(zhì)素生物降解過(guò)程當(dāng)中的酶類，木質(zhì)素過(guò)氧化物酶（Lignin Peroxidase，LiP）以過(guò)氧化氫為電子受體并以鐵撲啉作為輔基，其催化反應(yīng)遵循傳統(tǒng)過(guò)氧化或者氫氧化催化機(jī)理[4]。1988年Tien等[5]首次從白腐真菌屬的黃孢原毛平革的限制性培養(yǎng)基中發(fā)現(xiàn)了LiP并且認(rèn)定該酶為氧化還原酶。研究發(fā)現(xiàn)，LiP的作用底物比較廣泛，主要有酚類和非酚類芳香化合物質(zhì)，也包括苯并芘等多環(huán)芳烴[6]。

分子對(duì)接是研究小分子與大分子物質(zhì)相互作用模式、生物大分子之間的相互識(shí)別的常用方法之一[7]。分子動(dòng)力學(xué)模擬是利用計(jì)算機(jī)軟件對(duì)生物或化學(xué)反應(yīng)的過(guò)程進(jìn)行模擬（包括反應(yīng)的條件，如溫度、水環(huán)境、氣壓、酸堿性等），還可以用來(lái)預(yù)測(cè)蛋白質(zhì)的結(jié)構(gòu)、研究酶促反應(yīng)的過(guò)程與機(jī)制以及蛋白質(zhì)（酶）與小分子（配體）的相互作用機(jī)理等，為某些酶促反應(yīng)、藥物設(shè)計(jì)與藥理作用研究提供理論依據(jù)[8，9]。林錦霞等[10]用分子對(duì)接的方法研究了短小芽孢桿菌木聚糖酶與底物木聚糖的作用，也有學(xué)者利用分子對(duì)接和分子動(dòng)力學(xué)模擬的方法研究蛋白質(zhì)配體和受體的作用機(jī)理[11]。目前，有關(guān)LiP催化降解雙酚A類化合物的研究鮮見(jiàn)報(bào)道。因此，本研究選取了幾種常見(jiàn)的雙酚A類化合物，利用分子對(duì)接軟件和分子動(dòng)力學(xué)模擬軟件，將LiP與幾種小分子雙酚A類化合物進(jìn)行對(duì)接，尋找酶的活性位點(diǎn)，并對(duì)對(duì)接效果最好的小分子物質(zhì)進(jìn)行動(dòng)力學(xué)模擬，了解其相關(guān)作用機(jī)制，為L(zhǎng)iP催化代謝雙酚A類化合物提供理論指導(dǎo)。

1 材料與方法

1.1 數(shù)據(jù)來(lái)源

試驗(yàn)中選取的9種化合物包括：2，2-雙（4-羥基苯基）丙烷（簡(jiǎn)稱BPA）、4，4-二羥基二苯基甲烷（簡(jiǎn)稱BPF）、2，2′-二（4-羥基苯基）-六氟丙烷（簡(jiǎn)稱BPAF）、2，2-雙（2-羥基-5-聯(lián)苯基）丙烷（簡(jiǎn)稱BPQ）、2，2-二（4-羥基苯基）丁烷（簡(jiǎn)稱BPB）、4，4′-（1-甲基亞乙基）雙（2-甲基苯甲醚）（簡(jiǎn)稱BPC）、4，4-二羥基二苯砜（簡(jiǎn)稱BPS）、4，4′-亞乙基雙苯酚（簡(jiǎn)稱BPE）、4，4-（1-苯乙基）雙酚（簡(jiǎn)稱BPAP），其共同特點(diǎn)就是具有兩個(gè)或多個(gè)芳環(huán)并共同連接到一個(gè)中心C上，各種化合物的結(jié)構(gòu)式如圖1所示。

1.2 計(jì)算及模擬方法

LiP在RCSB Protein Data Bank—RCSB PDB中有多種構(gòu)型，本試驗(yàn)參照馮義平等[12]關(guān)于LiP催化去除水中雌激素的研究，選取酶晶體結(jié)構(gòu)1B82，以此物質(zhì)作為分子對(duì)接和動(dòng)力學(xué)模擬的蛋白質(zhì)受體。

LiP與雙酚A類化合物的分子對(duì)接：分子結(jié)構(gòu)進(jìn)行能量最小化的優(yōu)化（Tripos力場(chǎng)進(jìn)行），采用最陡下降法-共軛梯度法[13]，設(shè)置能量收斂值為0.005 kJ/（mol·？魡），優(yōu)化次數(shù)設(shè)置為1 000，最后進(jìn)行加電荷并負(fù)載Gasteiger-Hückel電荷[14]。為了研究的簡(jiǎn)化，除去LiP（1B82.pdb）中兩條完全對(duì)稱的B鏈中的一條，為蛋白質(zhì)加載電荷加氫，固定所有的氫原子。采用AMBER7 FF99力場(chǎng)對(duì)蛋白質(zhì)進(jìn)行修復(fù)和能量最小化優(yōu)化。對(duì)接過(guò)程以CH4、C=O、N-H為分子探針，指定原型分子的探測(cè)方式為配體，閾值默認(rèn)值為0.5，其他參數(shù)為默認(rèn)值[15]，進(jìn)行對(duì)接。對(duì)接之后找出對(duì)接效果最好或有代表性的物質(zhì)進(jìn)行分子動(dòng)力學(xué)模擬。

SYBYL-X1.3軟件中Total_Score為總的Surflex-Dock得分，其值愈大說(shuō)明對(duì)接復(fù)合物越穩(wěn)定[16]?；谶x擇小分子的原則，本試驗(yàn)中主要考慮以下兩點(diǎn)：一是選取對(duì)接效果最好的進(jìn)行模擬；二是從實(shí)際情況出發(fā)，雙酚A普遍存在于自然環(huán)境中，而且雙酚A類物質(zhì)含量大大低于BPA，如何去除BPA的危害也將是今后研究的重點(diǎn)。本試驗(yàn)?zāi)M將選取BPA和另外一種對(duì)接效果最好的小分子分別進(jìn)行分子動(dòng)力學(xué)模擬。

分子動(dòng)力學(xué)模擬：體系模擬在Linux系統(tǒng)中的GROMACS[17]軟件上進(jìn)行。采用GROMOS96 43A1力場(chǎng)條件，定義晶胞和添加溶劑。模擬系統(tǒng)使用立方體水盒子，蛋白質(zhì)與水盒子壁的最小距離為1 nm，并且加入了Na+，用以維持系統(tǒng)環(huán)境的性質(zhì)為中性。用最陡下降法進(jìn)行能量最小化，設(shè)置小分子化合物與蛋白質(zhì)分子為同一整體，并進(jìn)行系統(tǒng)平衡[18]。系統(tǒng)平衡包括正則系統(tǒng)NVT平衡和等溫等壓NPT平衡，兩步平衡以保證完全平衡每個(gè)系統(tǒng)的性質(zhì)[19]。最后，在水環(huán)境恒溫300 K的系統(tǒng)平衡條件下，對(duì)試驗(yàn)進(jìn)行了10 000步數(shù)動(dòng)力學(xué)模擬。

2 結(jié)果與分析

2.1 分子對(duì)接結(jié)果分析

將LiP與雙酚A類小分子化合物進(jìn)行分子對(duì)接，得到的分子對(duì)接結(jié)果如圖2所示，其結(jié)果數(shù)據(jù)見(jiàn)表1。由表1可知，BPQ的總值得分最高，說(shuō)明其與Lip對(duì)接最牢固，與LiP的對(duì)接穩(wěn)定性比其他雙酚A類化合物好。由于BPA是雙酚A類化合物的典型代表，而B(niǎo)PQ與LiP的對(duì)接得分最高，對(duì)接效果最好，所以選擇BPA和BPQ進(jìn)行模擬分析。

過(guò)氧化酶結(jié)合位點(diǎn)的結(jié)構(gòu)特性對(duì)于其與底物的結(jié)合有重要意義。通過(guò)SYBYL軟件，可以找出小分子化合物與LiP的結(jié)合位點(diǎn)，圖3中BPQ、BPA活性位點(diǎn)處的氨基酸殘基均為HIS39、MET172、ALA175、HIS176、ILE42、ALA179、PHE193、ARG43、ILE235和PHE46。在距離LiP活性位點(diǎn)5 ？魡范圍內(nèi)的蛋白質(zhì)結(jié)合位點(diǎn)的氨基酸殘基有HIS39、PRO145、

GLU146、PRO147、GLU40、MET172、LEU173、SER174、ALA175、HIS176、ILE42、VAL178、ALA179、ALA180、ARG43、PHE193、ILE235、LEU44、SER237、PHE265、VAL45、PHE269、LEU272、PHE46和HIS47。

氫鍵在受體和配體的相互結(jié)合過(guò)程中發(fā)揮著重要作用，研究結(jié)合位點(diǎn)上的氨基酸與配體鍵的距離可以反映受體和配體之間是否發(fā)生相互作用。利用SYBYL軟件還可以識(shí)別氫鍵連接的位點(diǎn)，即找到并生成氫鍵的供體和受體。圖4顯示小分子化合物BPQ、BPA通過(guò)氫鍵連接的氨基酸殘基位點(diǎn)和小分子在酶腔體中的位置。從計(jì)算結(jié)果得知，當(dāng)BPQ與LiP進(jìn)行對(duì)接時(shí)，氫鍵與酶的3個(gè)氨基酸殘基相連，它們分別是PRO83.O（id141）、ASP183.OD2（id470）和ASN182.HD21（id1221）；當(dāng)BPA與LiP進(jìn)行對(duì)接時(shí)，氫鍵與酶的2個(gè)氨基酸殘基相連，分別是PRO83.O（id141）和ARG43.HH21（id818）。將二者進(jìn)行比較可知，BPQ、BPA與LiP對(duì)接之間的氫鍵有一個(gè)共同的受體，那就是氨基酸殘基PRO83.O（id141），說(shuō)明該氨基酸殘基在雙酚A類化合物與LiP的反應(yīng)過(guò)程中起著重要作用。在BPQ、BPA與LiP對(duì)接的過(guò)程中，不僅存在著氫鍵的作用，而且還存在著疏水力的作用。LiP活性位點(diǎn)附近親水面和疏水面如圖5所示。LiP活性位點(diǎn)附近的親水面和疏水面分別作為氫的受體和供體的兩個(gè)原子，如果二者在3.5 ？魡距離范圍內(nèi)，則可能存在親水作用，如果蛋白質(zhì)中的某個(gè)疏水集團(tuán)抑制劑分子在6.5 ？魡距離內(nèi)，則二者之間可能存在某種疏水相互作用。由表2中的作用原子鍵長(zhǎng)可知，BPA和BPQ分別在蛋白質(zhì)活性位點(diǎn)的空腔內(nèi)形成親水相互作用。

2.2 動(dòng)力學(xué)模擬結(jié)果分析

分子動(dòng)力學(xué)模擬是在水環(huán)境中完成，體系的溫度為300 K，壓力為1.018個(gè)大氣壓[20]。溫度、壓力及勢(shì)能均達(dá)到收斂，說(shuō)明分子動(dòng)力學(xué)模擬正常。模擬過(guò)程中將配位蛋白質(zhì)（LiP）相對(duì)于初始結(jié)構(gòu)的均方根偏差RMSDs值作為衡量體系是否穩(wěn)定的指標(biāo)，即檢查能量項(xiàng)的收斂情況[21]。在本試驗(yàn)中蛋白質(zhì)復(fù)合物相對(duì)于初始結(jié)構(gòu)的均方根偏差如圖6所示。在10 000 ps的分子動(dòng)力學(xué)模擬過(guò)程中，小分子BPA和BPQ對(duì)蛋白質(zhì)的影響不同，未結(jié)合任何小分子的蛋白質(zhì)在3 000 ps后基本趨于穩(wěn)定，而結(jié)合了BPA的蛋白質(zhì)復(fù)合體的RMSDs值波動(dòng)明顯，且在8 000 ps時(shí)呈上升趨勢(shì)，這說(shuō)明小分子BPA的結(jié)合增加了蛋白質(zhì)柔性，而結(jié)合了BPQ的蛋白質(zhì)復(fù)合體RMSDs表現(xiàn)與未配位蛋白R(shí)MSDs一致，說(shuō)明了BPQ與酶結(jié)合更穩(wěn)定，BPA比BPQ對(duì)LiP酶構(gòu)象的誘導(dǎo)更顯著。

酶和底物結(jié)合以后，不僅LiP的構(gòu)象會(huì)改變，小分子的構(gòu)象也會(huì)發(fā)生變化，這種變化使得底物更接近于過(guò)渡態(tài)，可以降低活化能。除了分析蛋白質(zhì)結(jié)構(gòu)的變化外，本試驗(yàn)同時(shí)分析了小分子配體BPQ、BPA在分子動(dòng)力學(xué)模擬過(guò)程中結(jié)構(gòu)的變化。圖7為小分子配體相對(duì)于初始結(jié)構(gòu)的均方根偏差（RMSDs）。BPQ的RMSDs值在0～3 000 ps時(shí)波動(dòng)明顯，3 000 ps之后則在0.275 nm周圍穩(wěn)定波動(dòng)，說(shuō)明BPQ在結(jié)合位點(diǎn)附近經(jīng)歷構(gòu)象的改變；而B(niǎo)PA的RMSDs值在3 000～8 000 ps起伏變化較大，說(shuō)明在酶催化反應(yīng)中底物BPA經(jīng)歷了較多的構(gòu)象選擇；通過(guò)對(duì)比兩種小分子配體RMSDs值可知，BPA比BPQ的上下波動(dòng)更明顯，說(shuō)明BPA不容易與LiP形成穩(wěn)定結(jié)構(gòu)，在酶空腔中較難尋找到最適構(gòu)象。

回旋半徑用于度量結(jié)構(gòu)的緊密度，其未配位蛋白質(zhì)（LiP）骨架碳原子及其配位后蛋白碳骨架回旋半徑（Rg）[22]如圖8所示。BPQ體系和BPA體系的Rg不同程度變大，說(shuō)明體系發(fā)生了膨脹。其中兩個(gè)系統(tǒng)的Rg值變化趨勢(shì)基本是相同的，說(shuō)明未配位蛋白質(zhì)、蛋白質(zhì)與BPQ結(jié)合的復(fù)合物構(gòu)象改變的影響是適度的；但是可以明顯看出，酶在配位后，BPA的蛋白質(zhì)的回旋半徑變化范圍比BPQ的大，證明BPQ對(duì)酶的結(jié)合更穩(wěn)定，容易形成過(guò)渡態(tài)從而加快酶對(duì)BPA催化，而B(niǎo)PA對(duì)LiP有更多誘導(dǎo)。

三次蛋白質(zhì)（LiP）骨架碳原子的基于殘基分析的均方根波動(dòng)（RMSFs）[23]分析如圖9顯示。由圖9可以看出，無(wú)配體和有配體的復(fù)合物結(jié)構(gòu)對(duì)部分殘基的柔性的改變是一樣。結(jié)果表明，蛋白質(zhì)與小分子反應(yīng)時(shí)氨基酸殘基柔性發(fā)生改變，即小分子與氨基酸殘基發(fā)生碰撞位點(diǎn)的變化。

BPA和BPQ經(jīng)過(guò)10 000 ps的分子動(dòng)力學(xué)模擬后，在6.5 ？魡之內(nèi)與蛋白活性位點(diǎn)的相互作用力主要是疏水作用，其情況如圖10所示（睫毛鍵表示疏水力）。由圖10可以看出，6.5 ？魡距離內(nèi)，BPQ與活性位點(diǎn)Val184、Glu146、His82、Ile85、Phe193、Arg43、Glu40、His39、Ala180、Asp183、Asn182、Val181、Pro145、Pro83、Pro147、His47總共16個(gè)氨基酸殘基存在疏水作用；同樣，BPA與Met172、Leu169、Val158、Glu146、His176、Ile154、Ile235、Ser237、Pre145、Pre147、His47、Leu173、Phe46共13個(gè)氨基酸殘基也存在疏水作用。雙酚A類小分子與LiP之間的非鍵相互作用有利于底物的代謝。

3 小結(jié)

通過(guò)分子對(duì)接和分子動(dòng)力學(xué)模擬進(jìn)行LiP與BPQ、BPA的相互作用研究，結(jié)果表明LiP催化降解BPAs的反應(yīng)能夠自發(fā)進(jìn)行，小分子與酶之間通過(guò)氫鍵連接，BPQ、BPA與LiP的結(jié)合位點(diǎn)處的氨基酸殘基為HIS39、MET172、ALA175、HIS176、ILE42、ALA179、PHE193、ARG43、ILE235和PHE46。對(duì)接時(shí)小分子在蛋白質(zhì)空腔內(nèi)是親水性質(zhì)且呈現(xiàn)彎曲結(jié)構(gòu)，在10 000 ps的分子動(dòng)力學(xué)模擬過(guò)程后，BPQ和BPA在LiP活性位點(diǎn)疏水力發(fā)揮極大作用。由此可見(jiàn)，雙酚A類化合物中疏水作用激活了LiP。另一方面，BPQ和BPA經(jīng)歷了較大的結(jié)構(gòu)變化，底物分子BPA的結(jié)構(gòu)變化更為明顯，BPA引起的LiP結(jié)構(gòu)的變化也比BPQ要大，可以看出LiP可能被小分子BPA的誘導(dǎo)比BPQ顯著，并在此基礎(chǔ)上蛋白質(zhì)發(fā)揮了更進(jìn)一步的催化代謝作用。因此，推測(cè)LiP與BPQ的結(jié)合能力比BPA強(qiáng)，但是BPA分子結(jié)構(gòu)以及酶體積和形狀的變化程度更大，BPA比BPQ在酶活性中心也會(huì)經(jīng)歷更多的催化代謝途徑。

參考文獻(xiàn)：

[1] EUN J K，DONGHO L，BONG C C，et al. Association between urinary levels of bisphenol-A and estrogen metabolism in Korean adults[J]. Science of the Total Environment， 2014，470：1401-1407.

[2] 郭永梅.雙酚A的危害及相關(guān)限制法規(guī)[J].現(xiàn)代食品科技， 2012，28（5）：549-551.

[3] 國(guó)先芬，王貽鑫，耿存珍，等.環(huán)境激素雙酚A的毒理學(xué)研究進(jìn)展[J].環(huán)境科學(xué)與技術(shù)，2013，36（2）：86-93.

[4] CHOINOWSKI T，BIODIG W，WINTERHALTER K H.The crystal structure of lignin peroxidase at 1.70？魡 resolution reveals a hydroxyl group on the Tryptophan 171：A novel radical site formed during the redox cycle[J]. Molecular Biosystems，1999， 286：809-827.

[5] TIEN M，KIRK T K. Lignin peroxidase of Phanerochaete chrysosporium[J].Methods Enzyoml，1988，161：238-249.

[6] 張朝暉，邵 紅，周曉云，等. 三相鼓泡塔生物反應(yīng)器培養(yǎng)黃孢原毛平革菌合成木素過(guò)氧化物酶系[J]. 應(yīng)用與環(huán)境生物學(xué)報(bào)，2003， 9（3）：288-292.

[7] 姚 遠(yuǎn).幾類重要蛋白質(zhì)的結(jié)構(gòu)、功能及催化反應(yīng)機(jī)理的理論研究[D].長(zhǎng)春：吉林大學(xué)， 2008.

[8] HOU T J，ZHU L L，XU X J，et al. Diffusion of benzene in MCM-22 zeolite：A molecular dynamics study[J]. Acta Physico Chimica Sinica，2000，16（8）：701-707.

[9] WANG J L，SUN M，MIAO F M，et al. Computer chemical simulation-a new method of molecular conformational analysis[J]. Chinese Journal of Structural Chemistry，2000，19：281-287.

[10] 林錦霞，張燎原，張光亞，等.計(jì)算機(jī)模擬短小芽孢桿菌木聚糖酶與底物木聚糖的對(duì)接[J].生物工程學(xué)報(bào)，2007，23（4）：715-718.

[11] COLOSI L M，HUANG Q，WEBER W J，et al. Quantitative structure activity relationship based quantification of the impacts of enzyme-substrate binding on rates of peroxidase-mediated reactions of estrogenic phenolic chemicals[J].Journal Of The American Oil Chemists Society，2006，128：4041-4047.

[12] 馮義平，毛 亮，董仕鵬，等.過(guò)氧化物酶催化去除水體中酚類內(nèi)分泌干擾物的研究進(jìn)展[J].環(huán)境化學(xué)，2013，32（7）：1218-1225.

[13] 呂家楨，陸小華，周 健，等.化學(xué)工程中的分子動(dòng)力學(xué)模擬[J]. 化工學(xué)報(bào)，1998，49（1）：64-70.

[14] LU Y K，WANG J，HU Y，et al. Molecular modeling studies of vascular endothelial growth factor receptor tyrosine kinase inhibitors combining molecular docking and 3D-QSAR methods[J]. Structural Chemistry，2013，32：679-694.

[15] CAI W S，CHRISTOPHE C.Frontiers in high-performance， large-scale molecular dynamics.35 years of molecular-dynamics simulations of biological systems[J]. Acta Chimica Sinica，2013， 71（2）：159-168.

[16] CHEN C，LI W Z. Molecular dynamics simulation of hydrogen bonding characteristics in aqueous glycerol solutions[J].Acta Physico-Chimica Sinica，2009， 25（3）：507-512.

[17] LUTTMANN E，ENSIGN D L，VAIDYANATHAN V. Accelerating molecular dynamic simulation on the cell processor and playstation[J]. Journal of Computational Chemistry，2009，30（2）：268-274.

[18] BERENDSEN H，POSTMA J，GUNSTEREN W F.Molecular dynamics with coupling to an external bath[J]. The Journal of Chemical Physics，1984，81：36-84.

[19] WANG O，DONINI C，REYES P. Biomolecular simulations： Recent developments in force fields， simulations of enzyme catalysis， protein-ligand， protein-protein， and protein-nucleic acid noncovalent interactions[J]. Annual Review of Biophysics and Biomolecular Structure，2001，30：211-243.

[20] JORGENSEN W L，CHANDRASEKHAR J，MADURA J D， et al. Comparison of simple potential functions for simulating liquid water[J]. Physical Chemistry Chemical Physics，1983， 79（2）：926-935.

[21] 李利華，趙蔡斌，閔鎖田，等.基于配體-受體理論的計(jì)算機(jī)輔助藥物分子設(shè)計(jì)方法及應(yīng)用[J].西北藥學(xué)雜志，2007，22（5）：282-285.

[22] NIU J N，QIANG Y H.Molecular dynamics simulation on structure of water molecules in a kaolinite-water system[J]. Acta Physica Sinica， 2009，25（6）：1167-1172.

[23] MA G Z，JIANG Y J，YU Q S，et al.Molecular dynamics simulation of the norfloxacin-dna complex[J].Acta Chimica Sinica，2007，65（4）：284-288.