自由能方法及應用（三）煉金術自由能計算方法在藥物發現中的歷史

理論基礎

盡管自由能微擾（FEP）方法直到1980年代才引起計算機輔助藥物設計（CADD）領域的關注，其理論基礎還要回溯到1930年代的Peierls的工作，隨后在對非極性氣體的性質進行理論探索時，Zwanzig推導了FEP的主方程（參見等式1）：

FEP方程的內在美在于它可以通過對狀態0的勢能引導的配置進行采樣然后計算采樣配置中兩個狀態的勢能，從而計算兩個系綜狀態（狀態0和狀態1）之間的自由能差。在等式1中，??0代表在狀態0下的玻爾茲曼平均值，T是溫度，kB是玻爾茲曼常數，U1和U0分別是狀態1和狀態0的勢能。在這里，我們使用亥姆霍茲（Helmholtz）自由能差（ΔA）代替吉布斯（Gibbs）自由能差（ΔG），但是在ΔPV值通常較小的凝聚相模擬中，兩者的差異可忽略不計。從表面上看，使用該主方程似乎很簡單，但FEP計算的收斂性并不容易實現，尤其是對于兩個端態集合的平衡分布幾乎沒有重疊的體系。Kirkwood在1935年描述的使用耦合參數的想法極大地有利于FEP計算。簡而言之，就是經由耦合參數λ引入一系列中間狀態。隨著λ從0遞增到1，中間狀態的勢能（Uλ）從U0變為U1。然后計算相鄰狀態之間的ΔAs，所有ΔAs的總和得出總ΔA。用于將中間狀態定義為λ函數的電勢函數的典型形式是參考狀態U0的電勢與U1的電勢的線性組合（參見等式2）。?

另一種有效估計自由能差的方法是Bennett在1976年開發的Bennett acceptance ratio（BAR），通過類比FEP方程，Bennett推導出了等式3，其中引入了加權函數（w）。?

通過最小化自由能差的期望平方誤差，Bennett獲得了方程4的最佳加權函數，進而可以通過對ΔA進行迭代試驗來求解。

與Zwanzig的方程式不同，BAR分析需要在兩種狀態下進行采樣配置。2008年，多狀態的BAR（MBAR）方法得到了報道，該方法通過組合來自多個狀態而不是僅來自兩個狀態的模擬數據來估算自由能。另一種常規方法，即熱力學積分（TI），是基于Kirkwood在液體理論方面的工作得出的，其通過對狀態λ相對于λ的玻爾茲曼平均勢能導數進行積分來計算ΔA（參見等式5）。

技術支持

為了進行煉金術自由能（AFE）計算，基于上述理論基礎，我們需要生成代表性結構并評估相應的勢能。因此，模擬仿真技術和力場的發展成為AFEM計算的關鍵因素。一種主要的采樣技術是蒙特卡洛（MC）模擬。簡而言之，MC模擬從給定的配置開始，然后以最終采樣的配置恢復實際整體的真實分布的方式接受或拒絕下一個建議的配置。另一種主要的采樣技術是分子動力學（MD）模擬。MD模擬涉及沿時間軸的采樣配置。1964年，Rahman在液態氬方面發表了具有里程碑意義的MD模擬工作，并首次使用了Lennard-Jones勢描述了氬原子。隨著適用于生物分子的力場的發展，1977年，McCammon，Gelin和Karplus在Levitt和Warshel的Lifson實驗室開發的模擬程序的基礎上，首次對蛋白質進行了MD模擬。同一年Berendsen開發了“SHAKE”算法，該算法約束了高頻重原子形成的氫鍵，這一貢獻使得在可用計算機資源允許的情況下將MD模擬中使用的時間步長增大2倍來進行更加充分的采樣成為可能。此外，維持溫度恒定和壓力恒定的算法的實現使MD模擬能夠采樣物理上實際的配置。生物分子力場的發展是基于Lifson實驗室在1960年代后期所做的開創性工作發展而來，勢能函數的形式參見公式6。其中前三項分別表示鍵長、鍵角和二面角對勢能的貢獻，后兩項分別表示非鍵相互作用的范德華相互作用和成對靜電相互作用。在1980年代，基于Lifson力場模型開發了其他早期力場和程序，包括Karplus等人開發的CHARMM、Kollman等人開發的AMBER、Jorgensen等人開發的OPLS以及Berendsen和van Gunsteren等人開發的GROMOS。

進行AFE計算的另一項重要的技術支持是開發強大而準確的水模型。藥物分子的絕對結合自由能和相對結合自由能都受溶劑影響。因此，精確建模描述水和小分子的溶劑化至關重要。隱性水模型，如Poisson–Boltzmann或Generalized Born模型，可用于模擬溶劑對藥物結合的影響，但是許多水分子結合相關受體或水分子-藥物相互作用是特定的，因此需要使用顯式的水分子模型以獲得更高質量的結果。例如，活性位點處的水分子的置換在自由能預測中可以發揮顯著作用。根據使用了多少個額外的虛擬帶電原子（即沒有vdW勢項），水模型可以分為三點、四點、五點和六點模型。Jorgensen是開發用于分子模擬的實用而精確的水模型的先驅。他開發了TIP（transferable intermolecular potentials）模型來創建用于水、酒精和乙醚模擬的溶劑模型，并在此基礎上創建了TIPS2，TIP3P和TIP4P水模型，然后將這些水模型與一系列實驗結構和熱力學性質進行了比較，這些性質為將來對水模型進行評估奠定了基礎。同樣在1981年，Berendsen和他的同事們設計了SPC（simple point charge）水模型，并通過自我能量校正改進了SPC模型，創建了SPC/E水模型。目前TIP3P和SPC/E水模型得到了最廣泛的應用。盡管如此，開發更精確的水模型仍然是一個充滿挑戰的研究領域。

AFE早期研究

Postma、Berendsen和Haak在1982年報道了關于在水中形成空腔的FEP計算。這項工作幾年后，Jorgensen和Ravimohan于1985年計算了稀溶液中甲醇和乙烷的相對水合自由能，其結果與實驗基本吻合。這預示著FEP計算在CADD中的潛力，因為相對溶劑化自由能在確定一個共同受體位點上兩個配體的相對結合自由能方面起著主要作用。除了將FEP應用于相對水合自由能計算外，Jorgensen及其同事還應用FEP為SN1反應（（CH3）3CCl的解離）建立了平均力勢（PMF）分布，從而開啟了溶劑化對反應過程影響的研究方向。Jorgensen的工作鼓舞了Kollman及其同事將FEP應用于其AMBER MD模擬程序，并且他們運用FEP/MD方法研究了包括有機分子和蛋白質抑制劑體系在內的各種系統。這些早期研究顯示了AFEM在不同體系上的應用潛力，但使AFEM廣泛應用于CADD的主要推動力是實現了熱力學循環的吉布斯或亥姆霍茲自由能等狀態函數的冪運算，從而簡化了相對結合自由能（RBFE）計算的計算。圖1所示的上下部過程分別表示配體L1和配體L2與受體位點R的絕對結合自由能（ABFE）。通過連接圖1的上部和下部過程，我們發現ΔΔA=ΔA1–ΔA2可以通過兩個垂直過程表示的煉金術轉化獲得ΔΔA=ΔA3–ΔA4。因此，兩個配體的RBFE可以通過以下兩個煉金術過程的AFE計算獲得：R + L1→R + L2和RL1→RL2。以這種方式轉換問題會將兩個計算量大的計算轉換為兩個計算量大的AFE計算。除了使這些計算更容易處理之外，熱力學循環概念還提供了誤差消除功能，從而提高了計算出的自由能的精度。

圖1 熱力學循環：R是受體，L1和L2是兩個不同的配體，RL1代表與R結合的L1，RL2代表與R結合的L2。

圖片源自JCIM.

AFE在CADD中的現代研究

CADD中的第一個現代AFE研究是Wong和McCammon關于“酶和抑制劑的動力學與設計”的研究，并首次證明了“熱力學循環擾動法”在CADD中的適用性。1986年，通過在AMBER程序中實施的FEP方法，Kollman及其同事研究了一對與嗜熱菌素結合的抑制劑（即膦酰胺和膦酸酯），計算得出的ΔΔA值為4.21±0.54 kcal/mol，與4.10 kcal/mol的實驗值非常接近。此外，在McCammon和科爾曼Kollman組的兩項研究中，抑制劑的相對溶劑化自由能對計算出的ΔΔA做出了重要貢獻，這進一步說明了準確描述有機分子與水分子之間相互作用的重要性。隨后，Kollman及其同事將研究擴展到除了配體結合之外，還進行了位點特異性誘變對酶催化的影響。經過1980年代的初步研究后，CADD中的AFE研究在1990年代得到了進一步的發展，但是盡管MC模擬在宿主-客體系統的早期FEP研究和相對溶劑自由能的計算中是成功的，由于圍繞采樣效率的問題，其在蛋白質模擬中的使用受到限制。在一個蛋白質區域中，一個小幅度的主鏈角度的變化可能會導致遠端部位大幅度的移動，因此，總體接受率將很小。Jorgensen于1997年在CADD中進行了前兩次現代MC/FEP研究與實驗取得了很好的一致性，但兩項研究均使蛋白質主鏈保持固定。MC/FEP在2002年實現了突破，當時Ulmschneider和Jorgensen在MCPRO程序中為MC模擬實現了具有柔性鍵長和二面角的協調旋轉算法，這種實現允許在MC模擬中進行蛋白質主鏈采樣并使得MC/FEP能夠開始應用于蛋白質-配體體系。

AFE當前研究進展

在過去的20年中，隨著計算機性能的增強和對計算研究的興趣，人們已花費大量精力來改進CADD中的AFEM。我們無法總結所有貢獻，但將重點介紹一些重要的研究進展，包括采樣技術的改進以及力場的改進和自動化。除了所做的改進之外，我們還討論了一些仍需解決的問題，以提高AFEM的效率。

為了在AFE計算中達到收斂性和準確性，對所有相關構象進行采樣非常重要。但是，由于不同構象之間的勢壘很大，因此在典型AFE計算的時間尺度上，采樣通常是不完整的。因此，AFE計算結果通常取決于初始結構。目前已經有多種增強采樣方法得到了開發，其中一種主要方法是將AFEM與副本交換方法（REM）耦合。溫度REM涉及并行運行具有不同溫度的MD或MC模擬的副本，最近開發了更有效的副本交換方法并將其與FEP結合使用，即采用溶質退火的副本交換（REST1和REST2）方法。對于標準FEP計算未能捕獲構象變化的系統，已證明FEP/REST2可解決采樣問題并能夠精確再現實驗性ΔΔAs。另一種增強采樣的方法是使用圖形處理單元（GPU）。最近，FEP、TI和MBAR都已經在GPU上得到了實現，并且伴隨著2個數量級的加速，可以實現更多采樣。隨著計算機和GPU的性能不斷增強，大型蛋白質-配體數據集上的AFE計算可以快速且常規地進行。

當前第二項重大進展是在力場方面。用于類藥物有機分子的第一代力場是在1990年代末和2000年代初開發的，包括Merck大分子力場、OPLS、GAFF和CHARMM中的CGenFF。使用基于GPU的FEP/REST2方法和OPLS2.1力場的FEP+程序（Schr?dinger, Inc.）針對八個蛋白-配體靶標進行了追溯驗證。在八個系統上觀察到的平均誤差約為1 kcal/mol，并且據報道還為兩個前瞻性藥物設計項目可靠地預測了真正的陽性化合物。最近，GAFF v1.8力場與AMBER蛋白力場FF14SB相結合使用AMBER GPU TI實現針對相同的數據集進行了驗證。與實驗相比，所計算的ΔΔAs的RMSE（0.62–1.83 kcal/mol）與通過GPU FEP/REST2和OPLS2.1力場進行計算獲得的RMSE（0.93–1.41 kcal/mol）相當。通過重新擬合肽的二面角參數以及配體電荷模型的改進，Harder及其合作者基于OPLS2.1力場開發了OPLS3力場。最近，對配體二面角類型的進一步廣泛優化和對配體部分電荷分配的優化導致了OPLS3e力場得以報道。使用GPU FEP計算對前面八個蛋白質-配體數據集測試了OPLS3和OPLS3e力場，兩者均顯示出明顯改善，RMSE分別為0.70–1.25 kcal/mol和0.63–1.17 kcal/mol。

GAFF2是基于GAFF開發的，其包括的范德華相互作用、鍵長、鍵角和二面角參數的重新參數化有望進一步改善力場。最近，分別使用Flair中的MBAR和GROMACS中的非平衡TI對GAFF2結合AMBER蛋白力場（FF14SB和FF99SB-ILDN）在大型蛋白質數據集上進行了測試，對于AMBER FF14SB/GAFF2和AMBER FF99SB-ILDN/GAFF2，前面提到的八個蛋白質-配體測試集的RMSE分別為0.29-1.68 kcal/mol和0.90-1.56 kcal/mol。所有測試均表明，當前的類藥物分子力場與不同的AFEM結合使用的性能對于回顧性研究結果是令人滿意的。

當前力場使用原子類型來描述不同化學環境中的不同原子，然后根據原子類型分配不同的非鍵合和鍵合參數。CADD中AFEM的一個主要問題是，在藥物開發項目中，類藥物配體的化學空間很大，以至于現有的力場可能無法覆蓋所有的配體原子類型。因此，需要進一步的原子類型和參數化工作來解決這一缺陷。最近，一種新的力場格式，即SMIRKS Native Open Force Field（SMIRNOFF）格式，被用于開發開放力場模型。使用SMIRNOFF格式時，力場參數是通過直接化學感知分配的，而不是基于原子類型的。新格式具有更大的靈活性和簡便性，并且還可以避免基于原子類型的力場出現某些問題。另一種力場，即QUantum mechanical BEspoke（QUBE）力場也可以規避傳統力場的參數傳遞性問題。QUBE力場從QM電子密度得出系統特定的非鍵合參數。目前它已被用于研究蛋白質-配體結合自由能以及蛋白質動力學，并顯示出令人鼓舞的結果。

最后也是十分重要的一項進展是自動化。此前設置不同類型的自由能計算可能十分不方便、耗時且容易出錯，因此使這些方法更易于使用的關鍵步驟是簡化設置，運行模擬然后分析結果。目前已經開發了許多工具可以幫助自動執行AFE計算的各個階段，使用這些工具和程序，可以以更快、更可靠的方式執行大規模計算的設置和分析，這將大大有利于AFEM在制藥領域中的應用。但是，盡管可以使用多方位檢查來評估AFE計算的可靠性，包括向前和向后擾動的滯后、周期閉合錯誤以及相鄰的λ窗口之間的重疊，親自實踐的經驗仍然扮演著十分重要角色。因此還需要用戶在相關領域的專業知識和培訓，例如對程序適用范圍的理解，如果一切自動化，用戶將失去對設置模擬詳細信息的控制權，這可能會導致意外或意外的結果。

總結

總而言之，在過去30多年的時間里，我們在發展AFEMs方面的經驗有起有落，但可以相信通過最終的努力，我們最終可能會擁有一種能夠可靠地加速藥物開發的工具。最后的努力包括解決在工業環境中進行AFEM的一些剩余的關鍵的挑戰。Merck KGaA在最近關于AFEM的實際使用的研究中，列出了在藥物設計過程中使用AFEM時遇到的一系列問題。這些問題包括由于構象變化引起的蛋白質結構不確定性，配體結合模式的不確定性，X射線結構中缺失的成分，測定條件和模擬條件之間的差異，金屬和輔因子以及蛋白質和配體的質子化狀態的不確定等。此外，對配體互變異構的正確處理仍然是一個挑戰，還有對于AFE計算，涉及電荷變化、骨架變化、芳環到脂肪鏈等的擾動也帶來了一系列挑戰。隨著計算能力的不斷提高和新型分子表示形式的出現，未來幾年AFE計算技術的發展必將更加完善。

參考文獻

Lin Frank Song and Kenneth M. Merz, Jr.* Evolution of Alchemical Free Energy Methods in Drug Discovery. Journal of Chemical Information and Modeling. 2020. DOI: 10.1021/acs.jcim.0c00547