JCIM | Benchmark教你CADD如何選擇半經驗量子力學方法？

背景

在過去的幾十年中，信息技術的飛速發展為科學研究開辟了許多新領域。其中之一是計算機輔助藥物設計，它利用了現在強大的計算能力來預估基于靶標相互作用的各種藥物的結合特性。其中的一種常用方法是基于解Schrodinger方程的量子力學方法。量子力學方法（QM）方法在理論上無疑是精確的，但當藥物設計中經常遇到較大的體系（如蛋白質），QM需要更大的計算資源才能在合理的時間內完成計算，因此面對較大體系時無法同時兼顧準確度及計算時間。引入經驗參數的半經驗量子力學（SQM）方法是一種可以在時間和準確性之間進行折中，以換取近似方法。例如SQM方法可用于評分，以估算配體對蛋白質的親和力。但是，此類方法的參數測試通?；趥鹘y化學模型的數據集，盡管該模型系統包含基礎及不同類型的測試集，但它們可能無法準確真實地代表藥物設計中的蛋白質-配體相互作用。對于非共價相互作用，這些數據集通常僅包含簡單的有機分子。即從實際的蛋白質-配體復合物的結構中建立了一個很小的非共價復合物數據集，但是其中的相互作用是由非常小的片段代表的。為了使模型系統與藥物設計應用更加接近，數據集必須變得更大，更多樣化。由此提出了兩個新型數據集，它們來源于片段化的蛋白質-配體復合物。

方法

用于PLF547數據集構建的17種蛋白質配體復合物基于PDB數據庫中的以下晶體結構：10GS，2CET，2FVD，2OBF，2P4Y，2VOT，2VW5，2XB8、2YKI，2ZX6、3G0W，3GCU ，3GNW，3JVS，3NOX，3PE2和4GID。從參考文獻中獲得了用于QM計算的結構（包含氫原子）。這些配合物中的配體尺寸范圍為37至95個原子，其凈電荷為-1、0或+1。它們都含有芳香族雜環。五個配體包含二價的四面體硫原子，其中四個和三個分別包含F和Cl原子。為了生成與配體相互作用的蛋白質片段，C和Cα之間的鍵被切割并被H原子以0.767?的距離封端（保證其不帶電），將主鏈片段化為代表主鏈段的N個甲基甲酰胺殘基。Cα和Cβ之間的鍵斷裂（脯氨酸和甘氨酸除外），并在0.769?處加氫原子以將側鏈與主鏈分開。片段化通過Cuby4框架的片段化界面實現自動化。片段化后，選擇側鏈和主鏈的片段，該片段具有至少一個距配體4?的原子。此外，由于空間問題（非物理上緊密接觸），去除了對應于2YKI的酪氨酸和2VW5的賴氨酸的兩個片段。生成了構成PLF547數據集的547個蛋白片段-配體復合物。為了測試溶劑化能的方法，省去了另外8個碎片。這些片段包括以下內容（列出的數據集中的pdb文件的名稱）：2P4Y_27_tyr，2VOT_12_glu，2XB8_20_arg，3JVS_15_asp，3JVS_18_glu，3JVS_21_glu，3NOX_17_asp，3PE2_15_arg。除了每個配體的這些片段外，與配體（片段的并集）相距4?之內的蛋白質殘基的相同選擇也用作完整活性位點的模型。在這17種蛋白質-配體活性位點復合物中，有2種被排除在外， 其余15種配合物及其周圍的蛋白質構成了我們數據庫PLA15數據集的活性位點模型部分。用于構建PLA15數據集的蛋白由259到584個原子組成，可提供有關實際應用中SQM方法性能的獨特信息。

PLF547數據集根據相互作用的片段（蛋白質片段和配體）是否帶電荷，相互作用之間的距離以及蛋白質衍生部分的性質而分為幾組?；谒槠Ｐ碗姾傻幕鶊F是離子-離子（例如天冬氨酸側鏈與非零凈電荷的配體相互作用），離子-中性（相互作用的配偶體之一具有非零凈電荷，另一個為中性）或中性。我們還將系統根據相互作用之間的距離進行分類，分為“短”（片段與配體之間的最短接觸低于范德華半徑總和的90％），“平衡”（介于90％和110％之間） ），“長”（從110％到130％）和“遠處”（片段模型，其中相互作用的伙伴的最接近原子大于范德華半徑之和的130％）。根據蛋白質部分的基團為“骨架”（蛋白質骨架片段，以N甲基甲酰胺模型表示），“芳香族”（組氨酸，色氨酸，苯丙氨酸和酪氨酸的芳香族側鏈），“非極性”（丙氨酸的非極性側鏈） ，亮氨酸，異亮氨酸，纈氨酸和蛋氨酸），“極性”（絲氨酸，蘇氨酸，天冬酰胺，谷氨酰胺，半胱氨酸和脯氨酸的側鏈），“陰離子”（天冬氨酸和谷氨酸的側鏈）和“陽離子”（側鏈）賴氨酸，精氨酸和質子化組氨酸的鏈）。一些結果以箱形圖的形式表示數據集中誤差的分布。用中線劃分的框包含50％的數據點（中心兩個四分位數）。其寬度是相應中央四分位數寬度的1.5倍。超出此誤差范圍的其余復合物（離群值）用點表示。

對于PLF547數據集中片段復合物相互作用能的基準檢驗方法，我們選擇了基于顯式相關MP2-F12計算的復合方案，并采用ΔCCSD（T）方法校正。一個片段的總基準相互作用能為：

由于體系的大?。ㄆ渲幸恍┏^100個原子），MP2-F12方法以cc-pVDZ-F12基組進行計算。研究表明使用較小的基組計算的結果要優于傳統的MP2/CBS外推法。高階校正項ΔCCSD（T）通過CCSD（T）的相關性彌補了MP2/CBS對相互作用能計算的誤差。ΔCCSD（T） 部分為DLPNO-CCSD（T）/ aug-cc-pVDZ相互作用能與MP2/cc-pVDZ相互作用能之差。MP2和MP2-F12的計算是在TURBOMOLE 7.3中進行的，DLPNO-CCSD（T）的計算是在ORCA 4.0.1中進行的。

配體與完整的活性位點模型（PLA15數據集）之間的相互作用無法使用上述方案直接計算，因為這些系統包含500多個原子。因此，采用了多體效應通過DFT-D3級別（B3LYP-D3/ DZVP-DFT）計算的片段方案，其中i遍歷了活性位點的所有片段，因此所有成對的貢獻都在MP2-F12 + DLPNO-CCSD（T）級別上進行處理，整個系統的DFT-D3計算僅考慮了非可加性。后者主要受多體極化效應的影響，應該在DFT級別上對此進行充分說明。由于碎片相互作用能的基準是溶劑，因此我們使用基于BP/def2-TZVPD計算的COSMO-RS溶劑模型。實驗通過將Turbomole與COSMOTHERM X17程序結合進行，計算公式如下：

在數據集的片段-配體復合物部分，對QM和DFT方法進行了測試：包括MP2/ aug-cc-pVDZ，B3LYP-D3/def2-QZVP，BLYP-D3/def2-QZVP和BLYP-D3/DZVP-DFT。在D3校正中，使用了Becke- Johnson阻尼。對于以DZVP-DFT為基礎的計算，我們使用了重新參數化進行矯正。計算在TURBOMOLE 7.0程序下進行。在SQM方法中，選擇了以下方法進行測試：AM1，PM6，具有D3H4校正的PM6，PM7，具有D3H4和D3H5校正的DFTB和GFN2-xTB。我們在半經驗方法中添加了HF-3c進行比較，AM1，PM6和PM7已在MOPAC中進行了計算，DFTB +軟件中的DFTB，GFN2-xTB使用了方法作者提供的XTB代碼在ORCA進行了基于HF-3c的計算。并在Cuby4下進行數據集的計算；對于溶劑化，將HF/6-31G * /SMD方法作為另一種高級方法，并將以下SQM方法作為測試對象：PM6/COSMO，PM6/COSMO2，PM7/COSMO和PM7/COSMO2，DFTB3/SMD和DFTB3/PCM。

結果

圖1.相互作用能的相對誤差分布

圖片來源JCIM

本文的主要結果是對PLA15數據集中的活性位點模型的近似方法進行了評估。AM1和PM6相比帶有校正的SQM方法的性能差很多，平均相對誤差分別為55.9％和31.5％。PM7相對誤差為-24.1％，PM6-D3H4和DFTB3-D3H4的誤差分別為-9.4％和-11.1％。其中最好的方法GFN2-xTB誤差為（8.1％），它提供了SQM最小的平均誤差，但誤差分布較大。第二好的方法DFTB3-D3H5（-8.4％）具有最小誤差分布，并在較小的基組上達到了DFT-D3的精度（B3LYP-D3/DZVP-DFT，平均誤差為-7.2％）。僅考慮誤差的隨機部分（相對于平均值的平均絕對偏差）時，DFTB3-D3H5在SQM方法中表現最佳，其次是DFTB3-D3H4，然后是PM6-D3H4和GFN2-xTB。

表1.測試方法的RMSE（以kcal/mol為單位）

圖片來源JCIM

總相互作用能準確度和不同電荷下能量的準確度：

對QM方法，MP2和B3LYP-D3/def2-QZVP較準確的結果。相比之下，BLYP-D3/def2-QZVP準確性有所降低，最明顯的是對帶電物質進行計算的結果?結果表明對大多數體系的精度很好，但對于復雜的系統具有明顯的誤差。其中BLYP-D3/DZVP-DFT的準確度更低，?所有這些QM方法都為中性片段復合物提供了準確結果。但是比SQM計算昂貴。PLF547數據集（體系包含從44到114個原子）上采用BLYP-D3/DZVP-DFT計算一個體系的平均CPU時間是540 s，使用HF-3c的CPU時間是170 s（此處使用的是較小的基組，但是HF比GGA DFT更昂貴），而使用PM6或PM7則不到一秒鐘。

對于SQM方法，PM6-D3H4，DFTB3-D3H4和DFTB3 D3H5對于中性配合物計算最為準確，GFN2-xTB在所有SQM方法中均具有最低的RMSE，但總體而言效果不佳。它僅勝過AM1，PM6和HF-3c方法。沒有進行任何校正的PM6嚴重低估了整體的相互作用能，但PM7在高估了相互作用，這可以歸因于高估了長程色散作用，DFTB3-D3H5誤差最少，其次是DFTB3-D3H4，PM6-D3H4和GFN2-xTB。

圖2.相對于基準方法，片段數據集上測試方法的誤差分布

圖片來源JCIM

表2.按照相互作用片段間距離排序的測試方法的RMSE

圖片來源JCIM

?

表3.按照蛋白質類型歸類的測試方法的RMSE

圖片來源JCIM

不同相互作用距離的準確性比較：

當根據相互作用片段的最近距離對PLF547配合物進行分類時，對于某些SQM方法，RMSE不會隨距離增大而適當減小。相反，對于QM DFT和SQM DFTB3方法，范德華半徑總和的130％以上的復合物的RMSE大于平衡距離附近的復合物。PM6-D3H4由于沒有離域誤差，因此在更長和更遠的距離上比其他任何半經驗方法都具有優勢，并且在較短距離上的性能與DFTB3-D3H5相當。對于包括GFN2-xTB在內的DFTB3方法，很難在更長的距離上進行比較，因為它們也受到上述自交互誤差的影響。PM7的精度不會受到自相互作用誤差的影響，但隨著分子間的分離，PM7的精度不會提高，就像PM6-D3H4一樣，其色散校正在較大的間距時被高估了。

蛋白質部分類型的錯誤：

蛋白質片段按與配體相互作用的蛋白質種類進行分類。這些基團是主鏈（具有N-甲基甲酰胺模型），芳香族（His，Trp，Phe或Tyr的側鏈），非極性（Ala，Leu，Ile，Val或Met的側鏈），極性（N的側鏈）。Ser，Thr，Asn，Gln，Cys或Pro），陰離子（Asp或Glu的側鏈）和陽離子（Lys，Arg或質子化的His的側鏈）。在非極性復合物中，所有對非共價相互作用進行了校正的SQM方法都表現良好。區別是，對于SQM方法，涉及芳環的配合物的描述更困難。則相互作用（包括被高估的硫）的強烈影響導致HF-3c對于非極性配合物有較高RMSE，對于一些配體中氯原子的相互作用描述也并不準確，在極性和離子基團中有較大的差異。對于HF-3c方法，最大誤差來自對陰離子的描述，并且在大多數組中，該結果優于或等同于PM6-D3H4和DFTB3方法。但HF-3c方法對于陽離子-蛋白片段復合物的描述相對較好（見表3）。當蛋白質的骨架NH基團與配體（通常為氫鍵）形成緊密接觸時，PM7往往會高估相互作用。在HB104氫鍵數據集中也觀察到此錯誤。

表4.對于測試方法，片段配體配合物形成時能量變化的RMSE

圖片來源JCIM

除DFTB3/PCM會產生較大的誤差，尤其是對于中性系統,所有SQM方法的準確性均相當。相對于我們的基準，某些SQM方法甚至提供了比昂貴得多的HF/6-31G */SMD更精確的結果。COSMO-RS被認為是這類系統的最佳溶劑化方法之一，但仍未達到相互作用能計算（方法，基準計算）中基準方法達到的精度水平。因此，在這種情況下不能直接說明SQM優于QM方法。但具有COSMO2性能的PM6和PM7提供的精度最接近COSMO-RS。首先是DFTB3 / SMD，HF / 6-31G */SMD，然后是帶有COSMO的PM6和PM7。DFTB3 / PCM的準確性最低。當分析單個片段復合物時，對于大多數弱相互作用系統，DFTB3 / PCM低估了相互作用（始終提供更高的正值）。這可能歸因于排斥能的過高估計。對于較強的相互作用時，通常是在帶電基團或離子對之間，會過高的估計相互作用能。在這種情況下排斥力大于吸引力，誤差不能消除，因此PCM對于中性和帶電系統都會產生相對較大的誤差。對于弱相互作用的中性體系，根據系統不同誤差不同。相對于各自的基準，通過這些模型得出的溶劑化相互作用能估計值比通過測試的SQM方法進行的相互作用能計算得出的絕對誤差更大。這意味著在進一步開發藥物設計方法時，最好集中精力于對溶劑作用的改進。大多數情況下是由于訓練集中存在帶點體系而導致電荷體系的描述有所改善。通過中性體系的RMSE對比，COSMO2相對于COSMO略微改善了對非極性物質的計算結果，對于中性極性物質的描述精度稍差。

在為PLF547數據集選擇的基組（B-P/def2-TZVPD，COSMO-RS）的合理水平下計算活性位點溶劑化相互作用能是不可行的。重要的是量化大型模型中溶劑化能變化的整體誤差，以便將其與相互作用能的誤差進行比較。因此，此處討論的PLA15結果可通過各個片段中的溶劑化能差之和來近似。所有SQM模型都低估了能量，而HF/6-31G */SMD則高估了能量。DFTB3/SMD是性能最好的方法，而所有COSMO方法都具有相似的準確性。對比這些方法的RMSE，PM6/COSMO2和PM7/COSMO2是最準確的，但是除DFTB3/PCM和HF/6-31G */SMD以外，所有方法都具有較好的準確度。根據相對百分比和RMSE，PM6 / COSMO2隨（參考）能量的絕對大小而變化，而DFTB3 / SMD以基準參考的一小部分表示時，會提供更好的結果。DFTB3 / SMD（與DFTB3 / PCM不同）為中性配體配合物提供了更準確的估計值，因為沒有離子對相互作用，因此能量的絕對大小趨向于更小。從圖3中也可以明顯看出，DFTB3/PCM精度不足。在PM6和PM7中從COSMO過渡到COSMO2并不能僅在基準數據集中改善結果。對此之前的測試對碳酸酐酶II抑制劑進行評分，結果表明COSMO2系統地改善了得分與實驗結合自由能的相關性。?

圖3.蛋白質-配體復合的溶劑化能量變化的誤差

圖片來源JCIM

分析與總結

1. 這項工作引入了兩個數據集PLF547，PLA15，用于測試半經驗方法（SQM），重點是它們在計算機輔助藥物設計中的應用。PLF547數據集包含與配體相互作用的547個蛋白片段（氨基酸側鏈和主鏈片段），而PLA15集包含15個完整的活性位點AA5Q配體模型。相互作用能量基于顯式相關的MP2-F12和DLPNO-CCSD（T）校正的組合進行計算，可解決基組的收斂和高階相關效應帶來的誤差。

2. 使用這些數據集評估了該領域中使用較為廣泛的幾種半經驗QM方法的準確性?；跀祿?，PM6-D3H4方法提供了最好的結果，其次是DFTB3-D3H5。比起這些方法的實際應用，更重要的是它們在PLA15數據集的配體-活性位點相互作用的計算結果。其中DFTB3-D3H5是最成功的方法，其精度在較小的基組中接近DFT-D3計算的準確度，且誤差分布遠小于其他測試方法。

3.?對于結合時溶劑化能的計算，DFTB3/SMD，PM6/COSMO2和PM7/COSMO2溶劑化模型最準確。但溶劑模型的精確度仍然比氣相相互作用能的精確度差得多。因此，溶劑化限制了在藥物設計中的應用，應在以上方面做出改進。

參考文獻

Benchmarking of Semiempirical Quantum-Mechanical Methods on Systems Relevant to Computer-Aided Drug Design.?Kristian Kr?íz? and Jan R?eza?c?. J. Chem. Inf. Model. 2020. DOI: 10.1021/acs.jcim.9b01171