工業(yè)酶研究中的計(jì)算化學(xué)方法

劉海燕

工業(yè)酶研究中的計(jì)算化學(xué)方法

劉海燕

中國(guó)科學(xué)技術(shù)大學(xué) 生命科學(xué)學(xué)院，安徽 合肥 230026

文中介紹用于工業(yè)酶研究特別是用于指導(dǎo)酶工程的主要計(jì)算化學(xué)方法，包括分子力學(xué)力場(chǎng)和分子動(dòng)力學(xué)模擬、量子力學(xué)以及量子力學(xué)/分子力學(xué)結(jié)合模型、連續(xù)介質(zhì)靜電模型以及分子對(duì)接等。文中從以下兩個(gè)角度分別概要地介紹這些方法：一是方法本身的基本概念、原始計(jì)算結(jié)果、適用條件和優(yōu)缺點(diǎn)等；二是通過(guò)計(jì)算得到有價(jià)值的信息，指導(dǎo)突變體和突變庫(kù)設(shè)計(jì)。

分子力學(xué)，分子動(dòng)力學(xué)，量子力學(xué)/分子力學(xué)，連續(xù)介質(zhì)模型，酶工程

酶的工業(yè)應(yīng)用已有百年歷史，酶催化具有高效、高特異性和高選擇性、環(huán)境友好等優(yōu)點(diǎn)，在食品、農(nóng)業(yè)、醫(yī)藥、化工等不同行業(yè)廣泛應(yīng)用[1-2]。由于工業(yè)應(yīng)用環(huán)境遠(yuǎn)遠(yuǎn)不同于酶在自然界所處的生命環(huán)境，天然酶的性質(zhì)、催化功能等和其應(yīng)用環(huán)境通常不匹配，或者不是最優(yōu)的。此時(shí)，需要借助酶工程對(duì)酶的天然氨基酸序列進(jìn)行改造，以改進(jìn)其性能[3]。最常用的酶工程策略是構(gòu)建突變庫(kù)進(jìn)行篩選，即實(shí)驗(yàn)室定向進(jìn)化[4]。有效的定向進(jìn)化的必要前提之一是：相對(duì)于潛在有益突變體在庫(kù)中所占的比例而言，接受篩選的突變庫(kù)的庫(kù)容 (即庫(kù)中包含的突變體數(shù)) 要足夠大。突變體庫(kù)的庫(kù)容常常受到篩選方法、可用資源等客觀條件的限制。此時(shí)要解決的關(guān)鍵問(wèn)題是如何提高突變庫(kù)中有效突變體的占比。對(duì)酶序列、結(jié)構(gòu)與重要特性的關(guān)系的深入了解，有助于發(fā)現(xiàn)突變熱點(diǎn)，限定突變范圍，實(shí)現(xiàn)優(yōu)質(zhì)突變庫(kù)設(shè)計(jì)。計(jì)算化學(xué)方法是獲取這方面認(rèn)識(shí)的重要手段。有研究表明，相對(duì)于隨機(jī)突變庫(kù)，基于計(jì)算設(shè)計(jì)的蛋白質(zhì)突變庫(kù)可以將有效突變體占比提高多個(gè)數(shù)量級(jí)[5]。對(duì)一些困難的酶工程或蛋白質(zhì)工程課題而言，計(jì)算可帶來(lái)的大幅度改進(jìn)可能足以決定課題的最終成敗，而不再僅僅局限于效率的提升。實(shí)際上，計(jì)算化學(xué)、計(jì)算生物學(xué)方法已成功實(shí)現(xiàn)從頭設(shè)計(jì)具有天然酶不具備的催化功能的人工酶。由于本專輯中已有其他綜述專門討論對(duì)氨基酸序列進(jìn)行自動(dòng)優(yōu)化設(shè)計(jì)的方法，本文將主要介紹對(duì)給定氨基酸序列的酶進(jìn)行模擬和分析的計(jì)算方法。當(dāng)然，研究者可以用這些方法分別研究野生型和突變體，再對(duì)結(jié)果進(jìn)行比較。

長(zhǎng)期以來(lái)，對(duì)蛋白質(zhì)特別是酶的研究一直是計(jì)算化學(xué)研究的重要前沿[6-8]。主要方法包括基于經(jīng)典分子力學(xué)力場(chǎng)的分子動(dòng)力學(xué)模擬 (經(jīng)典MD)[9]、量子力學(xué) (QM)[10]以及量子力學(xué)/分子力學(xué) (QM/MM) 結(jié)合方法[8,11-12]、分子間復(fù)合物預(yù)測(cè)即分子對(duì)接 (Docking)[13]、定量分析靜電和溶劑效應(yīng)的可極化連續(xù)介質(zhì)模型 (PCM) 如泊松-玻爾茲曼模型 (PB)[14]，以及一些基于幾何性質(zhì)的模型等。本文將從以下兩個(gè)角度分別概要地介紹這些方法：一是關(guān)于方法本身，包括基本原理、原始計(jì)算結(jié)果、適用條件和 (潛在) 優(yōu)缺點(diǎn)等；二是如何用這些方法獲得工程相關(guān)信息，如更深入理解催化相關(guān)機(jī)制、對(duì)不同突變體相對(duì)于野生型的性質(zhì)或功能變化進(jìn)行理論預(yù)測(cè)或解釋，從而指導(dǎo)定向進(jìn)化所需的優(yōu)質(zhì)突變庫(kù)設(shè)計(jì)，或者基于對(duì)原始計(jì)算結(jié)果的分析建議特定的突變位點(diǎn)和突變類型等。

1 基于經(jīng)典分子力學(xué)力場(chǎng) (MM) 的分子動(dòng)力學(xué)模擬 (MD)

1.1 方法簡(jiǎn)介

我們暫不考慮酶催化所涉及的化學(xué)變化，僅考慮由于分子熱運(yùn)動(dòng)引起的酶構(gòu)象變化、酶與反應(yīng)物 (或產(chǎn)物) 之間非共價(jià)復(fù)合物的形成與解離等過(guò)程。在這些過(guò)程中，分子的電子狀態(tài)沒(méi)有發(fā)生變化 (例如沒(méi)有共價(jià)鍵的斷裂或生成)，適用分子力學(xué)力場(chǎng)模型。所謂分子力學(xué)力場(chǎng)，是用經(jīng)驗(yàn)的數(shù)學(xué)函數(shù)來(lái)表示分子體系勢(shì)能對(duì)幾何構(gòu)型 (即組成分子體系的所有原子的空間坐標(biāo)) 的依賴關(guān)系(圖1A)。換言之，如果我們用代表所有原子的空間坐標(biāo)，V()代表分子力場(chǎng)勢(shì)能，則當(dāng)分子從一種構(gòu)象1變到另一種構(gòu)象2后，其勢(shì)能變化：

?V=V(2)–V(1)。

依據(jù)熱力學(xué)理論，分子中的原子一直在進(jìn)行熱運(yùn)動(dòng)，即在隨時(shí)間不斷變化；此外，當(dāng)我們進(jìn)行實(shí)驗(yàn)觀察時(shí)，樣品總是由大量分子組成 (單分子實(shí)驗(yàn)例外)，不同分子處于不同構(gòu)象態(tài)。所以，從動(dòng)力學(xué)角度，我們需要考慮構(gòu)象隨時(shí)間的變化；從熱力學(xué)角度，我們需要考慮分子具有不同構(gòu)象的概率分布。分子動(dòng)力學(xué)模擬 (MD) 是考察這兩個(gè)方面性質(zhì)的最直接的模型(圖1B)。在進(jìn)行MD模擬時(shí)，我們從某個(gè)初始構(gòu)象出發(fā)，在每個(gè)時(shí)間點(diǎn)根據(jù)當(dāng)前構(gòu)象和勢(shì)能函數(shù)計(jì)算出作用在每個(gè)原子上的力 (力是勢(shì)能函數(shù)對(duì)原子坐標(biāo)的負(fù)導(dǎo)數(shù))，數(shù)值積分牛頓運(yùn)動(dòng)方程，得到下一個(gè)時(shí)間點(diǎn)的構(gòu)象；重復(fù)該過(guò)程，即得到構(gòu)象隨時(shí)間演化的軌跡。其間，可采用特殊算法模擬環(huán)境因素 (如溫度、壓強(qiáng)等) 對(duì)分子運(yùn)動(dòng)的影響。根據(jù)熱力學(xué)原理，當(dāng)時(shí)間間隔足夠長(zhǎng)，同一分子在不同時(shí)間點(diǎn)的構(gòu)象與熱力學(xué)平衡態(tài)下不同分子的構(gòu)象這二者的概率分布是相同的 (即時(shí)間平均等價(jià)于系綜平均)。因此，如果MD模擬的時(shí)間足夠長(zhǎng)，模擬獲得的構(gòu)象集合可以作為在特定熱力學(xué)平衡態(tài)下分子構(gòu)象分布的樣本。基于此原理，我們可以基于MD得到的時(shí)間軌跡分析體系熱力學(xué)平衡態(tài)下的任意可觀察性質(zhì)。

MD提供了一種以原子分辨率全面解析構(gòu)象變化動(dòng)力學(xué)變化過(guò)程和重要微觀量熱力學(xué)分布的強(qiáng)大計(jì)算工具，這對(duì)闡釋酶等復(fù)雜生物大分子機(jī)器的設(shè)計(jì)原理和工作機(jī)制尤為重要。由于目前大分子結(jié)構(gòu)的實(shí)驗(yàn)解析方法還只能提供時(shí)空平均的靜態(tài)結(jié)構(gòu)，MD模擬在相關(guān)研究中有不可替代的功能。在此前提下，MD工具本身還在不斷改進(jìn)完善之中。在方法學(xué)上，MD的主要局限來(lái)源于兩個(gè)方面：一是分子力場(chǎng)模型的精度問(wèn)題；二是模擬時(shí)間有限難以實(shí)現(xiàn)對(duì)構(gòu)象空間的充分采樣問(wèn)題。對(duì)于第一個(gè)問(wèn)題，近年來(lái)分子力場(chǎng)已有較大改進(jìn)，對(duì)生物大分子特別是蛋白質(zhì)體系構(gòu)象平衡的熱力學(xué)描述準(zhǔn)確度提高了，成功模擬了多種蛋白分子從頭折疊到天然結(jié)構(gòu)的過(guò)程[15-16]。在模擬時(shí)間方面，由于計(jì)算機(jī)軟硬件發(fā)展，目前可以用常規(guī)計(jì)算硬件 (如課題組用多核服務(wù)器) 完成對(duì)通常大小的體系 (如數(shù)百個(gè)殘基的酶分子在水溶液中的體系) 微秒量級(jí)的模擬。在此時(shí)間尺度下，可觀察到結(jié)構(gòu)域或環(huán)區(qū)開(kāi)合等過(guò)程。如果可用計(jì)算資源更多，對(duì)底物結(jié)合/解離等過(guò)程的直接模擬也可以實(shí)現(xiàn)。如果要研究時(shí)間尺度在模擬可及范圍之外的過(guò)程 (如別構(gòu)蛋白的大尺度功能變化等)，可以采用增強(qiáng)采樣方法[17]，但使用者對(duì)MD理論要有較為深入的了解。

目前絕大部分MD模擬應(yīng)用涵蓋的時(shí)間尺度為納秒到微秒，對(duì)構(gòu)象空間的采樣多被限制在初始結(jié)構(gòu)附近 (對(duì)單結(jié)構(gòu)域蛋白而言，通常是均方根位移在3–4 ?幅度內(nèi)的結(jié)構(gòu)漲落)。故而需要使用合理的初始結(jié)構(gòu)作為MD的輸入，模擬結(jié)果才有意義。多數(shù)情況下人們使用實(shí)驗(yàn)測(cè)定的晶體結(jié)構(gòu)或者基于同源蛋白比較建模得到的結(jié)構(gòu)作為MD的初始結(jié)構(gòu)。在模擬酶和底物復(fù)合物時(shí)，常常需要根據(jù)空酶或酶與其他分子復(fù)合物的結(jié)構(gòu)對(duì)復(fù)合物初始結(jié)構(gòu)進(jìn)行建模，方法包括使用分子對(duì)接或用底物直接替換晶體結(jié)構(gòu)中的其他小分子 (如抑制劑)。進(jìn)行MD模擬前還需要構(gòu)建能刻畫體系中所有化學(xué)單元的分子力場(chǎng)。當(dāng)被模擬的體系中包括底物小分子時(shí)，常會(huì)遇到MD軟件包中提供的標(biāo)準(zhǔn)分子力場(chǎng)不覆蓋該底物小分子的情況。此時(shí)可使用能自動(dòng)生成小分子力場(chǎng)的工具軟件[18-19]。在使用軟件自動(dòng)生成的力場(chǎng)進(jìn)行長(zhǎng)時(shí)間MD模擬之前，應(yīng)人工核對(duì)力場(chǎng)文件并用其進(jìn)行短時(shí)間模擬試算。

1.2 在酶研究中的應(yīng)用

從MD模擬得到的信息可以不同方式應(yīng)用于指導(dǎo)酶工程改造[20]。例如，通過(guò)比較常溫和高溫MD模擬，可以預(yù)測(cè)酶分子哪些區(qū)域的結(jié)構(gòu)穩(wěn)定性對(duì)環(huán)境溫度可能最為敏感。在這些區(qū)域引入脯氨酸點(diǎn)突變、二硫鍵等，有可能提升酶的耐熱性[21-24]。提高穩(wěn)定性的另一策略是設(shè)計(jì)能形成更多的表面氫鍵和鹽鍵的突變體[25-26]。在實(shí)驗(yàn)驗(yàn)證這類突變體之前，可平行模擬野生型和突變體，理論上評(píng)估突變是否可能達(dá)到預(yù)期效果[27-28]。除溫度外，MD還可以用于分析環(huán)境pH、溶劑等的變化對(duì)蛋白構(gòu)象及其穩(wěn)定性的影響[29-30]。

除穩(wěn)定性外，MD還被應(yīng)用于預(yù)測(cè)有可能顯著影響與底物結(jié)合/產(chǎn)物釋放相關(guān)的構(gòu)象動(dòng)力學(xué)的熱點(diǎn)殘基，為設(shè)計(jì)能改變底物選擇性、反應(yīng)選擇性、產(chǎn)物釋放速率等的突變或突變庫(kù)提供依據(jù)[31-32]。用MD研究底物/反應(yīng)選擇性的方法之一是比較(初始結(jié)構(gòu)) 不同的酶-底物復(fù)合物的模擬結(jié)果，預(yù)測(cè)親和力更高 (或反應(yīng)活性更高) 的底物或結(jié)構(gòu)狀態(tài)。計(jì)算親和力 (或反應(yīng)活性) 的嚴(yán)格定量方法是自由能計(jì)算[33-34]。由于自由能計(jì)算計(jì)算量大，目前大多數(shù)應(yīng)用中使用了定性方法預(yù)測(cè)：對(duì)相對(duì)親和力的定性判別可以小分子–大分子復(fù)合物結(jié)構(gòu)的穩(wěn)定性、平均分子間相互作用能等為依據(jù)；對(duì)反應(yīng)活性的定性判別依據(jù)則包括催化和反應(yīng)官能團(tuán)之間相對(duì)幾何構(gòu)型分布等[35]。這類定性判別結(jié)果可用作設(shè)計(jì)定向進(jìn)化序列文庫(kù)的依據(jù)。此外，MD模擬也可以用于分析底物結(jié)合/產(chǎn)物解離孔道周邊的熱點(diǎn)殘基[36-37]。這類應(yīng)用涉及對(duì)小分子從蛋白質(zhì)解離的解離路徑的模擬，如存在模擬時(shí)間尺度不足的困難，可使用增強(qiáng)采樣技術(shù)來(lái)克服[38-39]。

2 量子力學(xué) (QM) 以及量子力學(xué)/分子力學(xué) (QM/MM) 結(jié)合模型

2.1 方法簡(jiǎn)介

要模擬酶催化中的化學(xué)步驟，如共價(jià)鍵生成和斷裂、電子轉(zhuǎn)移、不同電子態(tài)之間的躍遷等，需要使用量子力學(xué) (QM) 模型。目前計(jì)算化學(xué)中常用的QM模型分為從頭算 ()、密度泛函理論(DFT)、半經(jīng)驗(yàn)方法等幾種類型[40]。其中，半經(jīng)驗(yàn)方法計(jì)算代價(jià)最低。但它們是非第一性原理方法，計(jì)算結(jié)果的可靠性對(duì)特定體系和問(wèn)題的依賴性高。從頭算和DFT方法都屬于第一性原理方法，具普適性。實(shí)際DFT模型中可能比從頭算包含了更多具經(jīng)驗(yàn)性的理論近似，但DFT能以非常高的計(jì)算效率來(lái)處理電子相關(guān)能。此外，對(duì)很多化學(xué)反應(yīng)問(wèn)題，最好的DFT模型對(duì)反應(yīng)程中能量變化等關(guān)鍵參數(shù)的計(jì)算誤差已能小至約1 kcal/mol左右，其結(jié)果已經(jīng)足以作為判別特定催化機(jī)制、反應(yīng)路徑化學(xué)上是否合理的依據(jù)。

給定分子的幾何構(gòu)型后，可用QM計(jì)算其能量。這被稱為單點(diǎn)計(jì)算 (即只處理幾何構(gòu)型空間中的一個(gè)點(diǎn))。QM模型更多被用于分子幾何構(gòu)型的優(yōu)化，即從一個(gè)初始構(gòu)型出發(fā)經(jīng)連續(xù)改變后找到一個(gè)局部穩(wěn)定結(jié)構(gòu) (能量比相鄰結(jié)構(gòu)低)，也可以用于發(fā)現(xiàn)連接反應(yīng)物和產(chǎn)物的最低能量路徑和路徑上的過(guò)渡態(tài)。這類計(jì)算中要考慮和比較不同的幾何構(gòu)型，一般要進(jìn)行數(shù)十到上千次的單點(diǎn)計(jì)算，計(jì)算量較大。節(jié)省計(jì)算量的一種常用策略是先使用效率高、精度有限的QM模型進(jìn)行大范圍的反應(yīng)路徑搜索優(yōu)化，在搜尋到的最低能量構(gòu)型/路徑附近再使用更高精度的模型完成構(gòu)型優(yōu)化，或進(jìn)行單點(diǎn)計(jì)算。

目前，將第一性原理QM方法應(yīng)用于整個(gè)酶分子計(jì)算量大，基本僅限于單點(diǎn)計(jì)算，還缺乏實(shí)用性。對(duì)大分子常用的是QM/MM模型 (圖2)[11]。在此模型中，分子體系被劃分為至少兩部分：直接參與化學(xué)反應(yīng)的部分用QM模型處理，其余部分用分子力學(xué) (MM) 處理。有多種不同的策略可以處理QM-MM間的邊界和相互作用[41]。在第一性原理QM/MM模型中，QM計(jì)算花費(fèi)的代價(jià)遠(yuǎn)超MM。因此，對(duì)QM區(qū)域大多使用構(gòu)型優(yōu)化的方法來(lái)預(yù)測(cè)或模擬其幾何結(jié)構(gòu)，對(duì)MM部分可使用分子動(dòng)力學(xué)模擬采樣[42]。這意味著計(jì)算結(jié)果可能對(duì)體系QM區(qū)的初始結(jié)構(gòu)比較敏感。此時(shí)需要用不同初始結(jié)構(gòu)模型進(jìn)行計(jì)算以獲得可靠結(jié)果。如果QM部分使用半經(jīng)驗(yàn)方法[43]或經(jīng)驗(yàn)價(jià)鍵理論[44-45]，則可能通過(guò)較長(zhǎng)時(shí)間的QM/MM MD采樣更充分地探索構(gòu)象空間，降低初始結(jié)構(gòu)的影響。

圖2 量子力學(xué)(QM)/分子力學(xué) (MM)模型

2.2 方法應(yīng)用

QM模型[10]和QM/MM模型[41]都被廣泛應(yīng)用于酶催化反應(yīng)化學(xué)機(jī)制的理論預(yù)測(cè)和檢驗(yàn)。其結(jié)果可以幫助我們判別哪些關(guān)鍵殘基參加了化學(xué)反應(yīng)過(guò)程，找到反應(yīng)的限速步驟，建立反應(yīng)中間物和過(guò)渡態(tài)的結(jié)構(gòu)模型，分析它們?nèi)绾闻c酶環(huán)境相互作用等。相對(duì)于QM團(tuán)簇模型，QM/MM模型能更真實(shí)地模擬化學(xué)反應(yīng)中心所處的酶環(huán)境。QM/MM已被廣泛應(yīng)用于從理論上預(yù)測(cè)/檢驗(yàn)酶催化的化學(xué)機(jī)制，分析預(yù)測(cè)環(huán)境氨基酸殘基對(duì)催化過(guò)程的可能影響[46]。原則上，這些結(jié)果可用于指導(dǎo)以提升催化活力、改變特異性或選擇性為目標(biāo)的定向進(jìn)化突變庫(kù)設(shè)計(jì)。更有挑戰(zhàn)性的研究，是基于從QM或QM/MM預(yù)測(cè)的過(guò)渡態(tài)結(jié)構(gòu)模型從頭設(shè)計(jì)新的活性中心，得到全新的人工酶[47]。

3 靜電連續(xù)介質(zhì)模型

3.1 方法原理

酶催化幾乎總是在特定溶液環(huán)境中完成的。溶劑效應(yīng)對(duì)酶的性質(zhì)有至關(guān)重要的影響。計(jì)算化學(xué)處理溶劑效應(yīng)的模型分為兩類：一類是顯溶劑模型，例如，在分子力學(xué)力場(chǎng)或QM模型中，每個(gè)溶劑分子和其中的每個(gè)原子是被顯式地包含在模型中的；另一類是隱溶劑模型或連續(xù)介質(zhì)模型[48]，模型中不包含溶劑分子和原子，而是用所謂的“溶劑平均場(chǎng)”來(lái)處理溶劑效應(yīng)。顯溶劑模型的優(yōu)點(diǎn)是能夠以完全一致的方式處理溶質(zhì)和溶劑，真實(shí)地模擬溶質(zhì)-溶劑間氫鍵、鹽鍵等特異性相互作用。其缺點(diǎn)是溶劑分子數(shù)目多，耗費(fèi)的計(jì)算量大。此外，溶劑隨機(jī)漲落對(duì)體系總能量的貢獻(xiàn)很大，必須進(jìn)行長(zhǎng)時(shí)間的模擬采樣平均才能消除漲落的影響。隱溶劑模型刻畫的是溶劑的平均效應(yīng)，溶劑的熱力學(xué)漲落已經(jīng)被平均了。

為簡(jiǎn)化處理，在隱溶劑模型中我們通常把非極性溶劑效應(yīng) (疏水效應(yīng)) 與極性溶劑效應(yīng)分開(kāi)。經(jīng)驗(yàn)表明，非極性溶質(zhì)的溶劑化自由能與其溶劑可接近表面積 (SASA) 成正比。因此，對(duì)這部分常使用SASA溶劑化模型。該模型中的參數(shù)包括計(jì)算SASA所需的原子半徑、溶劑分子半徑(水分子為1.4 ?)，以及溶劑化自由能正比于SASA的比例常數(shù)。這些參數(shù)一般是通過(guò)擬合小分子溶劑化自由能的實(shí)驗(yàn)值來(lái)確定的。

考慮極性溶劑效應(yīng)最常用的模型把被溶劑占據(jù)的區(qū)域視為具有特定介電常數(shù) (水為78.4) 的連續(xù)介質(zhì)，溶質(zhì)區(qū)域則視為被低介電常數(shù) (常用值為2–8) 或真空 (介電常數(shù)為1) 介質(zhì)占據(jù)(圖3A)。連續(xù)介質(zhì)會(huì)被溶質(zhì)區(qū)的電荷分布所產(chǎn)生的靜電場(chǎng)極化，由此產(chǎn)生的極化電荷分布又會(huì)在溶質(zhì)區(qū)產(chǎn)生靜電場(chǎng)，作用于溶質(zhì)電荷上。極化電荷產(chǎn)生的電場(chǎng)被稱作反應(yīng)場(chǎng)。因此，靜電連續(xù)介質(zhì)模型又被稱為反應(yīng)場(chǎng)模型。在溶劑區(qū)中無(wú)游離離子的連續(xù)介質(zhì)模型中，空間靜電勢(shì)與空間電荷分布之間的關(guān)系滿足泊松方程。對(duì)包含游離離子的溶液環(huán)境而言，離子的空間分布會(huì)受到空間靜電勢(shì)的影響。考慮這一因素，空間靜電勢(shì)與空間電荷分布之間的關(guān)系滿足泊松-玻爾茲曼方程 (PB方程)。PB方程是關(guān)于三維空間中靜電勢(shì)分布與電荷和電介質(zhì)分布之間關(guān)系的偏微分方程，可以用數(shù)值方法求解。求解酶等大分子體系PB方程的最常用數(shù)值方法是有限差分法 (FD)，合稱為FDPB模型(圖3B)[14]。用FDPB可以基于溶質(zhì)的空間電荷分布計(jì)算三維空間的靜電勢(shì)，進(jìn)而計(jì)算靜電自由能等其他性質(zhì)。在小分子體系的QM計(jì)算中，反應(yīng)場(chǎng)往往用分子表面的面電荷分布產(chǎn)生的電場(chǎng)來(lái)等效替代，相應(yīng)模型稱為可極化連續(xù)介質(zhì) (PCM) 模型。

MM/PBSA模型[49]是一種把顯溶劑模型和隱溶劑模型結(jié)合起來(lái)的方法：用顯溶劑分子力場(chǎng)模型進(jìn)行分子動(dòng)力學(xué)模擬，對(duì)溶質(zhì)構(gòu)象空間進(jìn)行抽樣；再對(duì)MD得到的溶質(zhì)構(gòu)象用隱溶劑模型進(jìn)行溶劑化自由能計(jì)算，用溶質(zhì)的分子力場(chǎng)能量 (MM) 和溶劑化自由能 (PBSA) 之和來(lái)評(píng)估溶質(zhì)構(gòu)象的穩(wěn)定性。該方法的好處是：用PBSA代替MM的溶劑部分，消除溶劑漲落；對(duì)顯溶劑MD得到的溶質(zhì)構(gòu)象集合進(jìn)行PBSA計(jì)算，避免PBSA結(jié)果對(duì)單一溶質(zhì)構(gòu)象的依賴性。

3.2 方法應(yīng)用

連續(xù)介質(zhì)模型的重要應(yīng)用之一是研究酶分子中帶電氨基酸側(cè)鏈基團(tuán)的質(zhì)子化狀態(tài)。軟件PROPKA通過(guò)求解PB方程計(jì)算不同質(zhì)子化狀態(tài)的靜電自由能來(lái)預(yù)測(cè)各可解離基團(tuán)的pKa[50]。酶分子的表面靜電勢(shì)分布是影響酶的底物選擇性的重要因素。給定酶分子的空間結(jié)構(gòu)和質(zhì)子化狀態(tài)，用FDPB方法可以計(jì)算酶分子的表面靜電勢(shì)分布，也可以預(yù)測(cè)氨基酸殘基突變或環(huán)境pH改變、離子濃度變化等對(duì)表面靜電勢(shì)的影響[14]。

在用QM簇模型研究酶催化的化學(xué)步驟時(shí)，常常需要使用PCM模型來(lái)模擬環(huán)境對(duì)反應(yīng)區(qū)的靜電影響。如果反應(yīng)過(guò)程涉及電荷分布的重大變化，不使用連續(xù)介質(zhì)的真空QM計(jì)算結(jié)果是不合理的，甚至可能導(dǎo)致錯(cuò)誤的定性結(jié)論。在QM/MM模型中，反應(yīng)中心周圍一般包含用MM方式處理的顯溶劑分子，一般不需要考慮連續(xù)介質(zhì)反應(yīng)場(chǎng)。但如果反應(yīng)前后體系的凈電荷發(fā)生變化 (如氧化還原電勢(shì)計(jì)算)，則很可能需要考慮體系邊界以外的溶液環(huán)境對(duì)反應(yīng)自由能的貢獻(xiàn)，此時(shí)可使用連續(xù)介質(zhì)模型對(duì)QM/MM結(jié)果進(jìn)行修正。

作為一種兼顧效率和精度的方法，MM/PBSA可以用于分析蛋白質(zhì)-蛋白質(zhì)、蛋白質(zhì)-小分子復(fù)合物的親和力[49]。為實(shí)現(xiàn)誤差抵消，常規(guī)的做法是：對(duì)復(fù)合物進(jìn)行顯溶劑的分子動(dòng)力學(xué)模擬，得到構(gòu)象集合；對(duì)每一個(gè)復(fù)合物構(gòu)象，分別計(jì)算復(fù)合物整體、組成復(fù)合物的每個(gè)單體的MM/PBSA能量；用整體的MM/PBSA能量與單體MM/PBSA能量之差對(duì)全部構(gòu)象的平均值來(lái)近似估計(jì)結(jié)合自由能。該方法可用于分析影響底物親和力的熱點(diǎn)殘基，也可以用于預(yù)測(cè)突變體的底物選擇性變化。

4 其他方法

4.1 分子對(duì)接

分子對(duì)接 (Docking) 泛指基于單體的結(jié)構(gòu)預(yù)測(cè)復(fù)合物結(jié)構(gòu) (和親和力) 的計(jì)算過(guò)程。小分子與蛋白質(zhì)對(duì)接是基于結(jié)構(gòu)進(jìn)行藥物虛擬篩選的核心工具，為此已開(kāi)發(fā)了多種算法[13]。這些算法和模型同樣可應(yīng)用于底物-酶復(fù)合物的對(duì)接。藥物虛擬篩選需要考慮大量不同小分子，出于計(jì)算效率考慮，分子對(duì)接計(jì)算中往往不考慮受體的結(jié)構(gòu)變化 (或僅僅考慮側(cè)鏈的結(jié)構(gòu)變化)。與虛擬篩選不同，在底物-酶對(duì)接研究中往往只需要考慮一種或數(shù)種不同底物，原則上可以更充分考慮酶的結(jié)構(gòu)變化。實(shí)現(xiàn)這一點(diǎn)的最直接方法是通過(guò)MD等構(gòu)象采樣方法得到多樣的酶結(jié)構(gòu)，分別與底物進(jìn)行對(duì)接。在底物-酶對(duì)接中，往往還可以利用底物與催化官能團(tuán)的相對(duì)空間排布來(lái)篩選/評(píng)價(jià)對(duì)接結(jié)果。

4.2 基于幾何結(jié)構(gòu)預(yù)測(cè)小分子孔道

大量實(shí)驗(yàn)研究發(fā)現(xiàn)，一些遠(yuǎn)離活性中心的突變對(duì)酶的催化性能可產(chǎn)生很大影響。其中一些位點(diǎn)可能是通過(guò)改變底物結(jié)合/產(chǎn)物釋放孔道發(fā)揮作用，孔道大小、孔道周邊殘基的物理化學(xué)性質(zhì)等可以改變底物/產(chǎn)物的通過(guò)速率，影響底物選擇性等。用孔道預(yù)測(cè)的方法可以找到相關(guān)的熱點(diǎn)殘基，為定向進(jìn)化庫(kù)設(shè)計(jì)提供依據(jù)。目前有數(shù)種基于幾何結(jié)構(gòu)的方法可以預(yù)測(cè)蛋白質(zhì)表面凹坑、內(nèi)部空腔、連接不同區(qū)域的孔道等[51-53]。這些方法使用靜態(tài)空間結(jié)構(gòu)作為輸入，多采用幾何學(xué)和圖論方法實(shí)現(xiàn)預(yù)測(cè)，計(jì)算效率高。

4.3 活性中心比較方法

目前，蛋白質(zhì)三維結(jié)構(gòu)數(shù)據(jù)庫(kù) (PDB) 中已積累了大量不同結(jié)構(gòu)類型、不同家族的酶的三維結(jié)構(gòu)數(shù)據(jù)。如果我們比較不同的酶，會(huì)發(fā)現(xiàn)其中一些盡管整體結(jié)構(gòu)序列不相似，但活性中心相似程度較高 (典型例子如絲氨酸蛋白酶共有的催化三聯(lián)體活性中心)?；钚灾行慕Y(jié)構(gòu)比較方法[54-55]能夠用于自動(dòng)檢索與當(dāng)前酶活性中心相似的其他酶的活性中心。把多個(gè)相似的活性中心在三維空間中疊合到一起，分析不同活性中心之間的異同能夠?yàn)橥蛔兾稽c(diǎn)選擇提供有價(jià)值的信息。

5 總結(jié)

為表述清晰，我們對(duì)上述方法的介紹是分類進(jìn)行的。在實(shí)際應(yīng)用中，不同類型的方法互相并不排斥。它們可以多種方式結(jié)合起來(lái)使用，以更好地回答我們所關(guān)心的問(wèn)題。例如，在酶-底物復(fù)合物模擬中，分子對(duì)接可以用來(lái)獲得模擬的初始構(gòu)象；MD模擬得到的構(gòu)象集合，可以用于孔道預(yù)測(cè)分析、分子對(duì)接、QM/MM模擬等；可以建立刻畫QM或QM/MM模型得到的過(guò)渡態(tài)的MM模型，用于長(zhǎng)時(shí)間的經(jīng)典MD模擬，以分析構(gòu)象漲落對(duì)化學(xué)過(guò)程的影響，或用于模擬大量突變體，實(shí)現(xiàn)基于MD模擬的突變體虛擬篩選；我們已經(jīng)提到的MM/PBSA方法，是MD與連續(xù)介質(zhì)模型的結(jié)合等。

用計(jì)算化學(xué)方法對(duì)蛋白質(zhì)等生物大分子體系的研究已有40多年的歷史。這些方法在自身不斷發(fā)展的同時(shí)，在工業(yè)酶研究中也越來(lái)越廣泛地得到使用。我國(guó)在計(jì)算化學(xué)與工業(yè)酶工程這兩個(gè)學(xué)科方向的研究隊(duì)伍均日益擴(kuò)大，研究能力快速提升。隨著這兩個(gè)學(xué)科方向的交叉結(jié)合日漸緊密，計(jì)算化學(xué)在酶工程中的應(yīng)用將會(huì)不斷拓寬和深入。蛋白質(zhì)工程、定向進(jìn)化等技術(shù)曾經(jīng)對(duì)工業(yè)酶研究產(chǎn)生了巨大的影響。計(jì)算方法的未來(lái)發(fā)展，特別是新酶設(shè)計(jì)方法的突破，將有望為進(jìn)入合成生物學(xué)時(shí)代的工業(yè)酶研究帶來(lái)新的技術(shù)突破。

[1] Choi JM, Han SS, Kim HS. Industrial applications of enzyme biocatalysis: current status and future aspects. Biotechnol Adv, 2015, 33(7): 1443–1454.

[2] Singh R, Kumar M, Mittal A, et al. Microbial enzymes: industrial progress in 21st century. 3 Biotech, 2016, 6(2): 174.

[3] Rastetter WH. Enzyme engineering. Appl Biochem Biotechnol, 1983, 8(5): 423–436.

[4] Zeymer C, Hilvert D. Directed evolution of protein catalysts. Annu Rev Biochem, 2018, 87: 131–157.

[5] Sun MG, Seo MH, Nim S, et al. Protein engineering by highly parallel screening of computationally designed variants. Sci Adv, 2016, 2(7): e1600692.

[6] Garcia-Viloca M, Gao JL, Karplus M, et al. How enzymes work: analysis by modern rate theory and computer simulations. Science, 2004, 303(5655): 186–195.

[7] Caddell Haatveit K, Garcia-Borràs M, Houk KN. Computational protocol to understand P450 mechanisms and design of efficient and selective biocatalysts. Front Chem, 2018, 6: 663.

[8] Hu H, Yang WT. Free energies of chemical reactions in solution and in enzymes withquantum mechanics/molecular mechanics methods. Annu Rev Phys Chem, 2008, 59: 573–601.

[9] Karplus M. Molecular dynamics simulations of biomolecules. Acc Chem Res, 2002, 35(6): 321–323.

[10] Siegbahn PE, Himo F. Recent developments of the quantum chemical cluster approach for modeling enzyme reactions. J Biol Inorg Chem, 2009, 14(5): 643–651.

[11] Warshel A, Levitt M. Theoretical studies of enzymic reactions: dielectric, electrostatic and steric stabilization of the carbonium ion in the reaction of lysozyme. J Mol Biol, 1976, 103(2): 227–249.

[12] Gao JL, Truhlar DG. Quantum mechanical methods for enzyme kinetics. Annu Rev Phys Chem, 2002, 53: 467–505.

[13] Warren GL, Andrews CW, Capelli AM, et al. A critical assessment of docking programs and scoring functions. J Med Chem, 2006, 49(20): 5912–5931.

[14] Gilson MK, Honig BH. Energetics of charge-charge interactions in proteins. Prot-Struct Funct Genet, 1988, 3(1): 32–52.

[15] Robustelli P, Piana S, Shaw DE. Developing a molecular dynamics force field for both folded and disordered protein states. Proc Natl Acad Sci USA, 2018, 115(21): E4758–E4766.

[16] Shaw DE, Maragakis P, Lindorff-Larsen K, et al. Atomic-level characterization of the structural dynamics of proteins. Science, 2010, 330(6002): 341–346.

[17] Bernardi RC, Melo MCR, Schulten K. Enhanced sampling techniques in molecular dynamics simulations of biological systems. Biochim Biophys Acta, 2015, 1850(5): 872–877.

[18] Malde AK, Zuo L, Breeze M, et al. An automated force field topology builder (ATB) and repository: version 1.0. J Chem Theory Comput, 2011, 7(12): 4026–4037.

[19] Vanommeslaeghe K, MacKerell AD Jr. Automation of the CHARMM General Force Field (CGenFF) I: bond perception and atom typing. J Chem Inf Model, 2012, 52(12): 3144–3154.

[20] Childers MC, Daggett V. Insights from molecular dynamics simulations for computational protein design. Mol Syst Des Eng, 2017, 2(1): 9–33.

[21] Joo JC, Pack SP, Kim YH, et al. Thermostabilization ofxylanase: computational optimization of unstable residues based on thermal fluctuation analysis. J Biotechnol, 2011, 151(1): 56–65.

[22] Watanabe M, Matsuzawa T, Yaoi K. Rational protein design for thermostabilization of glycoside hydrolases based on structural analysis. Appl Microbiol Biotechnol, 2018, 102(20): 8677–8684.

[23] Yu HR, Huang H. Engineering proteins for thermostability through rigidifying flexible sites. Biotechnol Adv, 2014, 32(2): 308–315.

[24] Pikkemaat MG, Linssen ABM, Berendsen HJC, et al. Molecular dynamics simulations as a tool for improving protein stability. Protein Eng, 2002, 15(3): 185–192.

[25] Perl D, Mueller U, Heinemann U, et al. Two exposed amino acid residues confer thermostability on a cold shock protein. Nat Struct Biol, 2000, 7(5): 380–383.

[26] Ladenstein R, Antranikian G. Proteins from hyperthermophiles: stability and enzymatic catalysis close to the boiling point of water//Antranikian G, Ed. Biotechnology of Extremophiles. Berlin, Heidelberg: Springer, 1998, 61: 37–85.

[27] Sakuraba S, Kono H. Spotting the difference in molecular dynamics simulations of biomolecules. J Chem Phys, 2016, 145(7): 074116.

[28] Tang L, Liu HY. A comparative molecular dynamics study of thermophilic and mesophilic ribonuclease HI enzymes. J Biomol Struct Dyn, 2007, 24(4): 379–392.

[29] Suplatov D, Panin N, Kirilin E, et al. Computational design of a pH stable enzyme: understanding molecular mechanism of penicillin acylase’s adaptation to alkaline conditions. PLoS ONE, 2014, 9(6): e100643.

[30] Gu W, Wang TT, Zhu J, et al. Molecular dynamics simulation of the unfolding of the human prion protein domain under low pH and high temperature conditions. Biophys Chem, 2003, 104(1): 79–94.

[31] Ying HX, Wang J, Shi T, et al. Engineering of lysine cyclodeaminase conformational dynamics for relieving substrate and product inhibitions in the biosynthesis of L-pipecolic acid. Catal Sci Technol, 2019, 9(2): 398–405.

[32] Dodani SC, Kiss G, Cahn JKB, et al. Discovery of a regioselectivity switch in nitrating P450s guided by molecular dynamics simulations and Markov models. Nat Chem, 2016, 8(5): 419–425.

[33] Mikulskis P, Genheden S, Ryde U. A large-scale test of free-energy simulation estimates of protein-ligand binding affinities. J Chem Inf Model, 2014, 54(10): 2794–2806.

[34] Liu HY, Mark AE, van Gunsteren WF. Estimating the relative free energy of different molecular states with respect to a single reference state. J Phys Chem, 1996, 100(22): 9485–9494.

[35] Petrovic D, Bokel A, Allan M, et al. Simulation-guided design of cytochrome P450 for chemo- and regioselective macrocyclic oxidation. J Chem Inf Model, 2018, 58(4): 848–858.

[36] Kingsley LJ, Lill MA. Substrate tunnels in enzymes: structure-function relationships and computational methodology. Proteins, 2015, 83(4): 599–611.

[37] Li GY, Yao PY, Gong R, et al. Simultaneous engineering of an enzyme’s entrance tunnel and active site: the case of monoamine oxidase MAO-N. Chem Sci, 2017, 8(5): 4093–4099.

[38] Hu Y, Liu HY. Case study on temperature-accelerated molecular dynamics simulation of ligand dissociation: inducer dissociation from the Lac repressor protein. J Phys Chem A, 2014, 118(39): 9272–9279.

[39] Lüdemann SK, Lounnas V, Wade RC. How do substrates enter and products exit the buried active site of cytochrome P450cam? 1. Random expulsion molecular dynamics investigation of ligand access channels and mechanisms. J Mol Biol, 2000, 303(5): 797–811.

[40] McDouall JJW. Computational Quantum Chemistry—Molecular Structure and Properties. Cambridge, UK: RSC Publishing, 2013.

[41] Senn HM, Thiel W. QM/MM methods for biomolecular systems. Angew Chem Int Ed, 2009, 48(7): 1198–1229.

[42] Zhang YK, Liu HY, Yang WT. Free energy calculation on enzyme reactions with an efficient iterative procedure to determine minimum energy paths on a combinedQM/MM potential energy surface. J Chem Phys, 2000, 112(8): 3483–3492.

[43] Liu HY, Müller-Plathe F, van Gunsteren WF. A combined quantum/classical molecular dynamics study of the catalytic mechanism of HIV protease. J Mol Biol, 1996, 261(3): 454–469.

[44] Hwang JK, King G, Creighton S, et al. Simulation of free energy relationships and dynamics of SN2 reactions in aqueous solution. J Am Chem Soc, 1988, 110(16): 5297–5311.

[45] Amrein BA, Steffen-Munsberg F, Szeler I, et al.: computer-aided directed evolution of enzymes. IUCrJ, 2017, 4(1): 50–64.

[46] Liu HY, Zhang YK, Yang WT. How is the active site of enolase organized to catalyze two different reaction steps? J Am Chem Soc, 2000, 122(28): 6560–6570.

[47] Kiss G, ?elebi-?l?üm N, Moretti R, et al. Computational enzyme design. Angew Chem Int Ed, 2013, 52(22): 5700–5725.

[48] Roux B, Simonson T. Implicit solvent models. Biophys Chem, 1999, 78(1/2): 1–20.

[49] Genheden S, Ryde U. The MM/PBSA and MM/GBSA methods to estimate ligand-binding affinities. Expert Opin Drug Discov, 2015, 10(5): 449–461.

[50] Olsson MH, Sondergaard CR, Rostkowski M, et al. PROPKA3: consistent treatment of internal and surface residues in empirical papredictions. J Chem Theory Comput, 2011, 7(2): 525–537.

[51] Pet?ek M, Ko?inová P, Ko?a J, et al. MOLE: a voronoi diagram-based explorer of molecular channels, pores, and tunnels. Structure, 2007, 15(11): 1357–1363.

[52] Yaffe E, Fishelovitch D, Wolfson HJ, et al. MolAxis: efficient and accurate identification of channels in macromolecules. Proteins, 2008, 73(1): 72–86.

[53] Chovancova E, Pavelka A, Benes P, et al. CAVER 3.0: a tool for the analysis of transport pathways in dynamic protein structures. PLoS Comput Biol, 2012, 8(10): e1002708.

[54] Konc J, Miller BT, Stular T, et al. ProBiS-CHARMMing: web interface for prediction and optimization of ligands in protein binding sites. J Chem Inf Model, 2015, 55(11): 2308–2314.

[55] Ren JY, Xie L, Li WW, et al. SMAP-WS: a parallel web service for structural proteome-wide ligand-binding site comparison. Nucleic Acids Res, 2010, 38(Web Server issue): W441–W444.

Computational chemistry approaches in studies on industrial enzymes

Haiyan Liu

School of Life Sciences, University of Science and Technology of China, Hefei 230026, Anhui, China

We review major computational chemistry techniques applied in industrial enzyme studies, especially approaches intended for guiding enzyme engineering. These include molecular mechanics force field and molecular dynamics simulation, quantum mechanical and combined quantum mechanical/molecular mechanical approaches, electrostatic continuum models, molecular docking, etc. These approaches are essentially introduced from the following two angles for viewing: one is about the methods themselves, including the basic concepts, the primary computational results, and potential advantages and limitations; the other is about obtaining valuable information from the respective calculations to guide the design of mutants and mutant libraries.

molecular mechanics, molecular dynamics, quantum mechanics/molecular mechanics, electrostatic continuum model, enzyme engineering

10.13345/j.cjb.190293

劉海燕 中國(guó)科學(xué)技術(shù)大學(xué)生命科學(xué)學(xué)院教授。于中國(guó)科學(xué)技術(shù)大學(xué)獲學(xué)士(1990年) 和博士 (1996年) 學(xué)位。曾在瑞士蘇黎世高工物理化學(xué)實(shí)驗(yàn)室、美國(guó)杜克大學(xué)化學(xué)系、北卡羅林納大學(xué)教堂山分校生物化學(xué)與生物物理系等機(jī)構(gòu)學(xué)習(xí)和從事博士后研究。2001入選中國(guó)科學(xué)院“百人計(jì)劃”。主要研究方向?yàn)榈鞍踪|(zhì)設(shè)計(jì)及其在合成生物學(xué)中的應(yīng)用、蛋白質(zhì)結(jié)構(gòu)功能的計(jì)算生物學(xué)與生物信息學(xué)。

劉海燕. 工業(yè)酶研究中的計(jì)算化學(xué)方法. 生物工程學(xué)報(bào), 2019, 35(10): 1819–1828.

Liu HY. Computational chemistry approaches in studies on industrial enzymes. Chin J Biotech, 2019, 35(10): 1819–1828.

July 3, 2019;

August 19, 2019

Supported by: National Natural Science Foundation of China (No.21773220).

Haiyan Liu. Tel/Fax: +86-551-63607451; E-mail: hyliu@ustc.edu.cn

國(guó)家自然科學(xué)基金 (No. 21773220) 資助。

(本文責(zé)編 陳宏宇)