馮晨博,馬維強,程 潤,王 駿

（南京大學物理學院，南京，210093）

蛋白質(zhì)是地球生命體系中重要的功能高分子，在新陳代謝、物質(zhì)運輸、形成細胞骨架、催化反應、免疫反應等過程中都發(fā)揮著積極的作用.在蛋白質(zhì)結(jié)構(gòu)形成過程和功能運動中，分子內(nèi)相互作用是決定這些行為的核心物理要素，而在各種分子內(nèi)相互作用中，疏水相互作用有十分重要的作用和意義.據(jù)估計，在單域球蛋白中，疏水相互作用占總體折疊自由能的70%～80%，這一現(xiàn)象啟發(fā)了多種簡化模型，如HP 模型等，這些模型為蛋白質(zhì)折疊物理理論的研究提供了基本物理基礎(chǔ)和典型物理模型.然而，疏水相互作用不同于經(jīng)典相互作用等微觀相互作用，它來源于蛋白質(zhì)分子與溶劑水分子的相互作用、水分子之間的相互作用以及水分子部分的熵效應，要把這些效應唯象定義為蛋白質(zhì)自身自由度的函數(shù)，常用的方法是通過對所有溶劑坐標積分來定義平均力勢.但事實上，在微觀時間和空間尺度上，還存在來自動力學的更多復雜性.在這種情況下，蛋白質(zhì)體系中的疏水相互作用一定是蛋白質(zhì)原子坐標的復雜多體函數(shù)，這使得描述疏水相互作用具有內(nèi)秉的復雜性，成為一件非平凡的工作，定量刻畫蛋白質(zhì)-溶液之間相互影響的自由能函數(shù)一直是蛋白質(zhì)物理研究的重要內(nèi)容［1-2］，常見的近似模型包括連續(xù)化溶劑（Continuum Solvent）模型［3］和原子溶劑化參數(shù)方法［4］（Atomic Solvation Parameters Approach）.前者用均勻可極化的連續(xù)介質(zhì)代替添加進系統(tǒng)中的水分子來降低計算耗時，后者假設(shè)原子的自由能與其暴露在溶劑中的表面積成正比，將蛋白質(zhì)的自由能寫為：

其中，系數(shù)σi與原子種類有關(guān)，可以根據(jù)實驗結(jié)果擬合得到.Ai為每個原子的溶劑可及表面積（Solvent-Accessible Surface Area，SASA）［5］，具體是指探針溶劑分子滾過待計算分子的范德華表面時，其球心可到達的表面，等效于將待計算分子中的每個原子半徑擴大一個探針分子半徑后構(gòu)成的范德華表面，而探針的大小反映了溶劑的微觀結(jié)構(gòu)特征.目前，研究溶劑可及表面積已成為分析蛋白質(zhì)體系疏水相互作用的典型手段，在蛋白質(zhì)與水的相互作用中，探針水分子半徑通常取1.4 ?.目前已有多種相關(guān)手段可以計算SASA，但這些方法在計算成本和精確度方面仍難以平衡，相比于力場中其他種類的相互作用，一定程度上限制了蛋白質(zhì)模型的簡化.因此，找到高精度和高效率的疏水相互作用（或SASA）的計算方法受到了蛋白質(zhì)模型研究的廣泛關(guān)注.

近年來，人工智能技術(shù)發(fā)展迅猛，其不僅在圖像/視頻識別［6-9］、自然語言分析［10-13］、游戲博弈［14-17］等經(jīng)典領(lǐng)域發(fā)揮十分重要的作用，還在疾病診斷［18-19］、藥物設(shè)計［20-21］、結(jié)構(gòu)生物學［22-23］、材料計算［24-25］乃至高能物理［26-27］等眾多領(lǐng)域也發(fā)揮了十分積極的作用.這歸因于深度學習（深度神經(jīng)網(wǎng)絡(luò)）方法在因果推斷、模型表示等方面的突出彈性，也為研究復雜多體關(guān)聯(lián)提供新的手段和工具.本文嘗試運用深度神經(jīng)網(wǎng)絡(luò)的方法，基于深度勢能架構(gòu)，以解析計算結(jié)果為標記，對SASA這一物理量進行學習，實現(xiàn)高精度的SASA 預測，相應的計算速度也顯著高于解析計算.這一模型實現(xiàn)了對蛋白質(zhì)體系多體相互作用的重構(gòu)，給出了一種準確且高效的計算手段，為提升蛋白質(zhì)模擬效率提供支撐.

1 方法

監(jiān)督學習是一種機器學習范式，是從成對的樣本和與其相關(guān)聯(lián)的標簽中嘗試學習一個函數(shù)，將輸入的特征向量映射為輸出，從而基于任意的可能輸入來推斷問題的答案.對于蛋白質(zhì)SASA的預測問題，本文借鑒了基于深度神經(jīng)網(wǎng)絡(luò)的分子動力學模擬方案——深度勢能分子動力學（Deep Potential Molecular Dynamics，DPMD）［28］.

1.1 模型為了降低問題的復雜度，將蛋白質(zhì)的總SASA 拆分為單個原子的SASA 之和，而單原子SASA 的預測依賴于原子的自身性質(zhì)和其近鄰環(huán)境.對于每個待計算的中心原子，首先提取與其相交的所有近鄰原子，并求出近鄰原子與中心原子的相對坐標，保證系統(tǒng)的平移對稱性.然后，將近鄰原子按照與中心原子的距離排序，保證系統(tǒng)的置換對稱性.最后，使用兩個最近鄰原子建立坐標軸的參考向量（如圖1 所示，圖中原子j和k為原子i的最近鄰和次近鄰原子，l為任意待旋轉(zhuǎn)的近鄰原子.為局域坐標系的單位向量），并通過旋轉(zhuǎn)矩陣調(diào)整所有近鄰原子的空間取向：

圖1 原子i 的局域坐標系示意圖Fig.1 The schematic diagram of the local coordinate related to the Atom i

其中，Rij=Ri-Rj，Rik=Ri-Rk，e[R]的作用是將R轉(zhuǎn)化為單位向量.至此，具有平移、旋轉(zhuǎn)和置換對稱性的局域坐標系得以建立.在對原始坐標做整體平移、旋轉(zhuǎn)或置換任意原子編號后，局域坐標系下的構(gòu)型不會發(fā)生改變，對應了這些操作下不變的SASA.這壓縮了原始數(shù)據(jù)的冗余，降低了勢能面的維度.

最后，將旋轉(zhuǎn)后的近鄰坐標除以其與中心原子的距離，以分開輸入角度和距離信息.據(jù)此，每個近鄰原子與三條局域坐標軸夾角的余弦值、與中心原子的距離倒數(shù)以及自身半徑，共計5 維向量作為近鄰的描述符被輸入網(wǎng)絡(luò).為了保證輸入向量定長，只有離中心原子最近的固定個數(shù)的近鄰信息得到保留.圖2 展示了不超過一定近鄰數(shù)的原子SASA 之和占全部原子SASA 的比例（SASA ratio）與近鄰數(shù)M的關(guān)系.這一比例在近鄰數(shù)小于40 時快速增長，近鄰數(shù)為56 時對應比例已超過99%（如圖中紅色虛線所示），即包含更多近鄰的原子SASA 之和不超過數(shù)據(jù)集中所有樣本SASA 總和的1%，可以忽略，這反映了空間堆積的飽和性，此時即使選取更大的近鄰數(shù)閾值，對SASA 的影響也很小.故選取截斷近鄰數(shù)Mcut=56，可以兼顧信息的完整性和效率.最后，為了便于神經(jīng)網(wǎng)絡(luò)的訓練，所有輸入數(shù)據(jù)都被標準化，即減去平均值后再除以標準差.

圖2 一定近鄰數(shù)以內(nèi)的原子SASA 之和所占比例與近鄰數(shù)M 的關(guān)系Fig.2 The ratio of SASA with an assigned number of neighbors

1.2 神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)和訓練過程得到的近鄰信息被輸入一個全連接的前饋網(wǎng)絡(luò)，數(shù)據(jù)從輸入層到輸出層單向傳播.該神經(jīng)網(wǎng)絡(luò)共有七個隱層，每層節(jié)點數(shù)分別為(1000，500，240，120，60，30，10)，每層節(jié)點與下一層節(jié)點之間存在連接和偏置權(quán)重矩陣.數(shù)據(jù)經(jīng)過權(quán)重矩陣的線性變換后，還會被施以非線性變換，以增強網(wǎng)絡(luò)的非線性擬合能力.每層的非線性激活函數(shù)為Sigmoid.圖3 直觀展示了SASA 的預測流程和神經(jīng)網(wǎng)絡(luò)的架構(gòu).圖3a 展示了完整的計算過程，Ri為蛋白質(zhì)中原子坐標，N為蛋白質(zhì)原子個數(shù)，{Rij}為原子i的近鄰坐標集合，j為截斷近鄰數(shù)，{Dij}為原子i的環(huán)境的描述符集合，Si為原子i的暴露面積的預測值，SASA 為蛋白質(zhì)總SASA 的預測值.圖3b更加詳細地展示了每個通道中Si的計算過程.

圖3 SASA 預測方法和神經(jīng)網(wǎng)絡(luò)架構(gòu)示意圖Fig.3 Schametic diagram of SASA prediction workflow (a) and the related neural network architecture (b)

網(wǎng)絡(luò)的損失函數(shù)由兩部分組成：

losssingle為輸出值與標簽的均方誤差（Mean-Square Error，MSE），losssum考察批次的整體誤差.訓練前期psingle占據(jù)主導，使權(quán)重從隨機值快速擬合；訓練后期psum逐漸占據(jù)主導，幫助穩(wěn)定整體誤差.只選用MSE作為損失函數(shù)會使網(wǎng)絡(luò)的訓練效果下降，且批次的整體誤差會變得更不穩(wěn)定.

訓練過程中使用Adam 優(yōu)化器調(diào)整權(quán)重，每次從數(shù)據(jù)集中隨機抽取1000 個數(shù)據(jù)組成一個批次（batch）輸入網(wǎng)絡(luò)，訓練100 個周期.每個周期調(diào)整權(quán)重的次數(shù)為訓練集的樣本數(shù)除以每個批次中的樣本數(shù)，即N(train dataset)/N(batch)，平均每個樣本在每個周期輸入網(wǎng)絡(luò)一次.此處取N(batch)=1000，約為蛋白質(zhì)包含的典型原子數(shù).學習率的初始值為1×10-3，隨訓練周期衰減，前50 個周期內(nèi)每周期衰減0.95 倍，后50 個周期調(diào)整為0.85 倍，最終值為2.68×10-8.為了緩解模型可能的過擬合，網(wǎng)絡(luò)應用L2 權(quán)重正則化，為網(wǎng)絡(luò)的損失函數(shù)添加了一項與權(quán)重值大小有關(guān)的懲罰：

使網(wǎng)絡(luò)傾向于選擇參數(shù)值分布的熵更小的簡單模型，從而限制模型的復雜度.

1.3 數(shù)據(jù)集蛋白質(zhì)數(shù)據(jù)來自SCOPe（Structural Classification of Proteins Extended）數(shù)據(jù)集［29-30］，選用其中一個子集，其包含的蛋白質(zhì)相互之間的序列相似度低于95%，既排除了相似度過高的同源序列，保障了相關(guān)數(shù)據(jù)的獨立性，也避免數(shù)據(jù)集中出現(xiàn)過多的重復構(gòu)型.SCOPe95 數(shù)據(jù)集共包含30201 個蛋白質(zhì)，每個蛋白質(zhì)包含的原子個數(shù)跨越三個數(shù)量級，大多數(shù)蛋白質(zhì)的原子數(shù)約為1000 個.由于氫原子的體積過小，對SASA的影響可以忽略，所以在處理數(shù)據(jù)時所有氫原子都被去除.剩余原子共4120 萬個，其中訓練集、驗證集和測試集分別占80%，10%和10%.

數(shù)據(jù)集的標簽來源于ARVO 工具計算的解析值［1］，使用相對簡單的球極平面投影法［31］.計算時每個原子半徑都增加了1.4 ?，即一個水分子的半徑.大約40%的原子暴露面積為0.圖4 展示了單原子SASA 的分布情況.相應的統(tǒng)計信息：SASA 最小為0，最大為90.3 ?2，中位數(shù)為0.49 ?2，均值為6.95 ?2，標準差為12.5 ?2.因此，暴露面積的中位數(shù)僅為0.49 ?2.它們完全被埋在蛋白質(zhì)內(nèi)部，即使周圍仍有空隙，水分子仍然無法與之接觸.考慮到加入水分子半徑后的單個原子總表面積約120 ?2，則蛋白質(zhì)體系中約有6%的面積可與溶劑接觸.

圖4 蛋白質(zhì)中20 萬原子SASA 的直方圖統(tǒng)計Fig.4 Histogram of atomic SASA in proteins

2 結(jié)果與討論

2.1 模型訓練使用Nvidia RTX 2070 訓練約16 h 后得到的學習曲線如圖5 所示，圖中對應訓練集的結(jié)果用紅線表示，對應驗證集的結(jié)果用藍線表示，包括訓練集和驗證集的平均絕對誤差（圖5a）和整體誤差losssum（圖5b）.測試集中，隨著訓練的進行，最終MAE和整體誤差分別達到0.29 ?2和14.8 ?2.由于驗證集的整體誤差抖動幅度較大，圖5b 應用了滑動平均，滑動窗口為10 個周期，驗證集MAE下降到0.294 ?2，表明網(wǎng)絡(luò)沒有出現(xiàn)過擬合，且曲線與訓練集MAE貼合緊密.對于整體誤差，訓練集呈平穩(wěn)下降趨勢，但驗證集在前期偶爾出現(xiàn)強烈的反彈，在后期隨著psum的增大逐漸趨于平穩(wěn)，與訓練集相差不大，收斂于14.2 ?2.結(jié)合訓練集和驗證集的結(jié)果可以證明，在嘗試調(diào)整網(wǎng)絡(luò)架構(gòu)和超參數(shù)的過程中沒有出現(xiàn)對驗證集的過擬合，也說明整體誤差對權(quán)重的變化比MAE更敏感.

圖5 使用蛋白質(zhì)全原子結(jié)構(gòu)訓練的模型的學習曲線Fig.5 The learning curve during the training processes with all-atomic protein structures

2.2 對SASA 的預測利用訓練得到的神經(jīng)網(wǎng)絡(luò)對蛋白質(zhì)SASA 進行預測.對單原子SASA 的預測結(jié)果如圖6 所示，圖6a 給出了單原子SASA預測值與相應解析值的比較（藍色散點）和線性擬合（紅線），圖6b 為預測值的分段平均誤差.可以發(fā)現(xiàn)，預測值與解析值的誤差大致在0.5，沒有隨著標簽值的增加而增加.

圖6 (a)單原子SASA 的預測值和解析值對比；(b)單原子SASA 預測值的平均絕對誤差與SASA 的關(guān)系Fig.6 (a) The comparison between the predicted and actual SASA for a single atom，(b) the MAE of predicted atomic SASA for the cases with various atomic SASA

選取2000 個蛋白質(zhì)進行總SASA 的預測，其包含的原子全部屬于測試集，結(jié)果如圖7 所示.圖7a 給出了預測值SASApre和相應解析值SASAact的比較（藍色散點）和線性擬合（紅線），其中，S0=1；圖7b 為預測值與相應解析值的相對誤差絕對值與體系原子數(shù)N之間的關(guān)系.結(jié)果表明，預測誤差的絕對值在100以內(nèi)的蛋白質(zhì)超過97%，測試集蛋白質(zhì)的總預測值比標簽平均高出約7.5，約為一個原子的平均暴露面積，說明網(wǎng)絡(luò)的預測存在約+0.08%的整體誤差.對測試集誤差取絕對值后除以各自標簽，得到平均相對誤差：且該相對誤差δ與原子數(shù)N不存在相關(guān)性（如圖7b 所示）.為了說明本方法的精確度，選取三種計算SASA 的傳統(tǒng)方法GetArea［32］，F(xiàn)reeSASA［33］以及御茶水女子大學提供的計算網(wǎng)頁［34］，分別對33 個蛋白質(zhì)進行測試.結(jié)果顯示，三種方法對全部測試蛋白的計算值與ARVO 精確值的相對誤差平均達到1.93%.因此，本方法與傳統(tǒng)典型方法比較，誤差縮小了近一個數(shù)量級.

圖7 (a)蛋白質(zhì)總SASA 的預測值SASApre 和解析值SASAact 的對比；（b)預測相對誤差δ 與原子數(shù)N 的關(guān)系Fig.7 (a) The comparison between the predicted and actual SASA for proteins，(b) the predicted relative error for the proteins with various numbers of atoms

2.3 對解折疊蛋白質(zhì)的預測以上的預測蛋白質(zhì)處于天然的折疊態(tài)，為了進一步驗證網(wǎng)絡(luò)對蛋白質(zhì)SASA 的預測能力，將其應用到訓練集中未曾出現(xiàn)的解折疊態(tài)蛋白質(zhì)中，對測試集中被加熱解折疊的三個蛋白質(zhì)進行了測試，得到的結(jié)果如表1 所示.比較網(wǎng)絡(luò)對這三個蛋白在折疊態(tài)和解折疊態(tài)的預測誤差，可以看出，和折疊態(tài)相比，解折疊態(tài)的蛋白質(zhì)SASA 都有升高，這符合折疊態(tài)是蛋白質(zhì)自由能極小狀態(tài)的觀點.但其預測的相對誤差都有顯著增大，從0.5%擴大到1%左右，說明解折疊態(tài)中的原子環(huán)境與訓練集中的折疊態(tài)存在較大差異，表明該網(wǎng)絡(luò)對未曾接觸的數(shù)據(jù)的可遷移性還有欠缺.

表1 折疊和解折疊態(tài)蛋白預測誤差比較Table 1 The comparison between predicted errors for folded and unfolded proteins

2.4 網(wǎng)絡(luò)精度及效率該網(wǎng)絡(luò)的預測精度與訓練集大小Nsample的關(guān)系如圖8 所示，其中圖8a 和圖8b 分別刻畫驗證集的平均絕對誤差MAE和整體誤差losssum，圖8c 為測試蛋白質(zhì)體系的平均相對誤差Dprotein，此時的訓練集容量Nsample為1 萬到3300 萬.可以發(fā)現(xiàn)，預測精度與訓練集容量的對數(shù)呈線性關(guān)系，訓練集容量每增大一個數(shù)量級，MAE大約下降0.29 ?2.

圖8 神經(jīng)網(wǎng)絡(luò)預測精度和訓練集容量Nsample的關(guān)系Fig.8 The relations between the predicted accuracy and the size of training set

該神經(jīng)網(wǎng)絡(luò)的預測用時與使用C 語言的傳統(tǒng)ARVO 工具相比，具有極大的優(yōu)勢.圖9 對比了用兩種方法計算不同原子數(shù)N的蛋白質(zhì)所需的時間，其中，t0=1 s，本文提出的解析方法加速了約80 倍.兩種方法的時間復雜度大致都為O(n)，但神經(jīng)網(wǎng)絡(luò)方法將復雜度的系數(shù)降低了近兩個數(shù)量級.以典型的蛋白質(zhì)原子數(shù)1000 為例，神經(jīng)網(wǎng)絡(luò)預測用時約0.05 s.為了克服解釋語言Python運算速度過慢的缺點，在建立本地坐標系時使用numba 中的JIT 解釋器［35］.

圖9 不同方法對不同大小蛋白質(zhì)SASA 的預測用時對比Fig.9 The time to calculate SASA for various sizes of proteins by different algorithms

2.5 拓展應用除了使用全原子的蛋白構(gòu)型進行SASA 預測外，此方法也可應用到只有α 碳原子的蛋白質(zhì)結(jié)構(gòu)中.α 碳原子的SASA 指其所屬殘基所有原子SASA 之和SASAres.與全原子一樣，同樣存在大量α 碳原子的暴露面積為0，即其所屬殘基被完全包裹在蛋白質(zhì)內(nèi)部.將殘基按親疏水性質(zhì)分類，親水殘基的平均SASAres達到70.9，而疏水殘基平均只有35.4.親水殘基中被完全埋藏在內(nèi)部的僅占2.5%，而疏水殘基中這一比例達到12.2%，如圖10 所示.圖10a 和圖10b 分別為親水和疏水殘基的情形，平均值分別為70.9和35.4.因此，蛋白質(zhì)傾向于將疏水殘基藏在內(nèi)部，這也是其折疊的主要驅(qū)動力.

圖10 親水殘基(a)和疏水殘基(b) SASAres的直方圖統(tǒng)計Fig.10 Histograms of SASA for the hydrophilic (a) and hydrophobic (b) residues

完整結(jié)構(gòu)中的4120 萬個原子包含520 萬個α 碳原子，選取18范圍內(nèi)的α 碳近鄰原子的信息輸入網(wǎng)絡(luò)，除了角度、距離和半徑信息外，每個近鄰的殘基種類也通過字典編號并輸入.經(jīng)過嵌入層展開后，每個種類的索引被映射到20 維的向量中，并與角度、距離、半徑信息連接，流入與全原子網(wǎng)絡(luò)架構(gòu)相同的全連接層.在前五個全連接層后分別插入一個Dropout 層，其在網(wǎng)絡(luò)中引入噪聲，有助于降低過擬合［36］.典型蛋白質(zhì)包含約100個殘基，批次大小從全原子的1000 降為100.損失函數(shù)和學習率與全原子網(wǎng)絡(luò)相同.520 萬樣本中，訓練集占360 萬，驗證集和測試集各80 萬.對2000 個蛋白質(zhì)進行測試，得到平均相對誤差為：

個別樣本的誤差δres超過10%.如圖11 所示，圖11a 為預測值和相應解析值的比較（藍色散點）和線性擬合（紅線），圖11b 為預測值與相應解析值的相對誤差絕對值與蛋白質(zhì)鏈長Nres之間的關(guān)系.考慮到原子數(shù)與SASA 本就呈明顯的線性關(guān)系，給網(wǎng)絡(luò)僅提供α碳結(jié)構(gòu)還是欠缺了較多的信息.

圖11 (a)基于α 碳結(jié)構(gòu)的蛋白質(zhì)總SASA 預測值與解析值對比；(b)預測相對誤差δres 與蛋白質(zhì)鏈長的關(guān)系Fig.11 (a) The comparison between the predicted and actual SASA of proteins，(b) the predicted relative error for the proteins with various sizes

3 結(jié)論

深度學習的發(fā)展為研究多體相互作用帶來了新范式，將神經(jīng)網(wǎng)絡(luò)方法應用到蛋白質(zhì)SASA 的計算是一種新嘗試.本文提取單個原子環(huán)境信息并轉(zhuǎn)換到滿足對稱性要求的局域坐標系中，將數(shù)據(jù)集輸入深度神經(jīng)網(wǎng)絡(luò)，通過最小化損失函數(shù)來優(yōu)化網(wǎng)絡(luò)權(quán)重參數(shù)，使其自行學習多體相互作用.將訓練穩(wěn)定的網(wǎng)絡(luò)預測的結(jié)果與使用解析工具ARVO 計算的SCOPe 數(shù)據(jù)集中蛋白質(zhì)所含原子的各自暴露面積進行比較，發(fā)現(xiàn)在對單原子SASA 預測時，MAE和整體誤差分別為0.29和2.3%，對蛋白質(zhì)SASA 預測的平均相對誤差為0.26%.還分別利用解析方法和神經(jīng)網(wǎng)絡(luò)方法對隨機選擇的蛋白質(zhì)SASA（單個蛋白質(zhì)包含的原子數(shù)為100～10000 個）進行預測并記錄所用的時間，發(fā)現(xiàn)和解析方法相比，在保證合理的準確度的前提下，本文提出的神經(jīng)網(wǎng)絡(luò)方法的預測速率提升近兩個數(shù)量級.還對蛋白質(zhì)的α 碳粗?；Y(jié)構(gòu)進行了SASA 預測，精度比全原子下降了約一個量級，但與傳統(tǒng)近似計算方法的精度相當，且仍能夠保持較高的計算效率.因此，在利用機器學習對蛋白質(zhì)全原子和粗?；Y(jié)構(gòu)的SASA 預測中，未來可以構(gòu)建更精巧的網(wǎng)絡(luò)拓撲，更高效地提取原子構(gòu)象的重疊關(guān)系.這也為動力學模擬提供了可能的支撐，結(jié)合更多真實的物理限制，使構(gòu)象沿勢能面的特定方向演化.