計算機模擬技術在SARS?CoV?2病毒分子機制研究中的應用

劉 凱，許志旺，何王秋，常 珊，孔 韌

（江蘇理工學院生物信息與醫(yī)藥工程研究所，常州 213001）

0 引言

新冠病毒暴發(fā)至今已有兩年之久，但此刻仍不能放松警惕，新冠病毒已進化出多種“令人擔憂的變異毒株”。截至2021年8月，Delta變體都是有很高傳播率的菌株；而Mu 變體則是對血清中和的抗性最高。世界衛(wèi)生組織（WHO）在2021年11月26日舉行的緊急會議中提到，最新出現(xiàn)的Omicron 變異株正在全球肆虐，Omicron變體的突變比Delta 多，甚至被認為是目前最危險的冠狀病毒株，這引起了新一輪的感染高峰。

在病毒傳播過程中，S 蛋白中氨基酸的變化可能會顯著改變病毒抗原性和中和抗體的效力。隨著SARS?CoV?2 在世界各地的傳播，S蛋白的突變不斷被報道。據(jù)GISAID 數(shù)據(jù)庫報告，SARS?CoV?2 刺突蛋白中有930 個自然發(fā)生的錯義突變，SARS?CoV?2刺突蛋白中從ASP614到GLY614（D614G）的關鍵突變使SARS?CoV?2比原始菌株更具傳染性。然而，SARS?CoV?2 的大多數(shù)疫苗、檢測試劑和抗體都基于武漢參考菌株的S蛋白。在其他冠狀病毒中，錯義突變被證明對MERS CoV 和SARS CoV 中的中和抗體具有耐藥性。就HIV 而言，已知錯義突變會影響包膜糖蛋白表達、病毒感染性、改變中和敏感性，并通過中和抗體產(chǎn)生耐藥性。SARS?CoV?2的S 蛋白是受體利用的關鍵決定因素，S 蛋白上的氨基酸突變甚至會增加病毒宿主的范圍，在不同物種間進行傳播。

SARS?CoV?2 病毒通過S 蛋白與人體受體ACE2 結合［1］，S 蛋白由S1 和S2 亞基組成。S1 亞基有三個結構域：N 端結構域（NTD），C 端結構域（CTD）和受體結合結構域（RBD）。而S 蛋白的RBD 區(qū)主要負責ACE2的識別，它是導致宿主感染的一個重要因素［2］。Omicron變異株最大的特點則是存在大量突變和免疫逃避［3?4］，特別是S蛋白RBD區(qū)域的突變更加導致了逃避免疫的發(fā)生。

分子動力學是一種以分子力學為基礎的計算機仿真模擬實驗方法，主要用于模擬并觀測分子在不同環(huán)境下的運動變化［5?6］。分子動力學模擬（MD）是分子模擬中最接近實驗條件的模擬方法，能夠從原子層面給出體系的微觀演變過程，直觀地展示實驗現(xiàn)象發(fā)生的機理與規(guī)律，促使研究向著更高效、更經(jīng)濟、更有預見性的方向發(fā)展，因此，分子動力學模擬在生物，藥學，化學，材料科學的研究中發(fā)揮著越來越重要的作用［7］。分子動力學模擬是目前為體系提供準確數(shù)據(jù)的最流行方法之一，隨著新冠病毒不停突變，越來越多的人開始將分子動力學模擬方法運用到自己的蛋白體系研究當中，這為研究新冠蛋白結構的穩(wěn)定性以及構象變化提供了幫助。在分子動力學模擬之前一般需要做一些初始化操作，以計算某體系的結合自由能為例。前期需要文件的處理，生成拓撲文件和坐標文件以及能量最小化，一般流程如圖1所示。

圖1 計算結合自由能的一般流程

接下來，本文從分子動力學模擬在刺突蛋白和ACE2的結合以及刺突蛋白和中和抗體的結合中的應用等多個方面進行闡述，然后討論目前刺突蛋白突變面臨的挑戰(zhàn)，分析并探討刺突蛋白突變對突變體和人體細胞結合的發(fā)展情況。

1 用分子動力學模擬研究刺突蛋白和ACE2的結合

1.1 刺突蛋白和ACE2復合物的RMSD和RMSF

體系經(jīng)過分子動力學模擬之后，一般會進行均方根偏差（RMSD）和均方根波動（RMSF）的計算［8?9］，簡單概括RMSD 表示的是分子結構隨著模擬時間變化的程度，而RMSF值表示的是分子中每個原子運動后與原始參考位置的波動程度。

RMSD代表了模擬過程中蛋白系統(tǒng)的構象穩(wěn)定性［6,10］，因為它反映了蛋白質(zhì)結構中原子相對于起始結構的平均位移，其計算公式如式（1）所示，N為原子總數(shù)，δi是指某一幀第i個原子相對于參考構像的位移偏移量。而對其所有原子進行平方和除N 平均再開方后，求得的為某一時刻這一幀的整體RMSD值。

RMSF表示模擬后蛋白系統(tǒng)的殘基位置與原始位置相比的波動程度［11?12］，其計算公式如式（2）所示，其中引入了計算時間T，這里是得到時間段T所有時刻的位移平移量的平方和，最終求得的是某個原子在時間T內(nèi)相對于初始結構的RMSF。

今年成都大學胡建平團隊研究了SARS?CoV?2野生型和突變體Delta，Mu 和Omicron 與人體ACE2 受體的結合情況，胡建平團隊使用殘基的基本結構取自PDB 數(shù)據(jù)庫，選取PDB ID 分別為7KJ2 和7V7S 當作野生型和Delta 變體的初始結構，然后在此基礎上使用Pymol 誘變工具和HadDock 服務器進行突變，通過SWISS?MODEL服務器進行序列補全，就構建出了①四個全長S 蛋白三聚體，包括WT、Delta、Mu 和Omicron變體；②四種具有ACE2 受體的三聚體復合物；③具有TMPRSS2 的四個S 蛋白三聚體對接復合體模型。

初始文件形成后開始進行分子動力學模擬，該團隊使用ff14SB 力場和Amber19 軟件包對12個S 蛋白三聚體系統(tǒng)進行了MD 模擬。溶質(zhì)被放置在邊界為15.0 ? 的八面體盒中，其中溶劑效應由TIP3P 水模型表征。除了常規(guī)MD 模擬的300 K 溫度設置外，還增加了310 K 和320 K 選項，以評估溫度是否影響S 蛋白構象。在MD 模擬之前，進行了最陡下降和溶質(zhì)無約束最小化兩步能量優(yōu)化。最后，在能量最小化后進行兩階段100 ns MD 模擬。第一階段是5 ns 溶質(zhì)約束動力學，力常數(shù)為10 kcal·mol?1??2。第二個是95 ns無約束生產(chǎn)模擬，其中采用SHAKE算法來防止涉及非氫原子的化學鍵的破壞。

通過其實驗發(fā)現(xiàn)WT、Mu 和Omicron 變體在30 ns 后趨于穩(wěn)定，RMSD 分別為4.61 ?、4.84 ?和5.15 ?，而Delta 在65 ns 后才達到平衡，均方根偏差為6.75 ?。這表明這三種變體的構象穩(wěn)定性不如WT，尤其是Delta［13］。為了分析不同溫度對S蛋白結構穩(wěn)定性的影響，在三種不同溫度（即300 K、310 K 和320 K）下進行了比較MD模擬，WT 在300 K/310 K/320 K 下的平均RMSD值依次為4.61 ?/5.47 ?/5.05 ?，Delta、Mu 和Omicron 的RMSD 值分別為6.75 ?/5.34 ?/5.24 ?、4.84 ?/4.59 ?/5.09 ? 和5.15 ?/4.8 ?/4.26 ?。從結果可以看出Mu 最為穩(wěn)定而Omicron 結構也比WT 要穩(wěn)定，但是Delta 變種的穩(wěn)定性要比WT差。同樣從RMSF的結果也可以看出幾種變體穩(wěn)定性從大到小依次是Mu、Omicron 和Delta，但是都要比WT更加穩(wěn)定。

1.2 刺突蛋白和ACE2復合物的結合自由能

結合自由能是查看體系結合情況的一個核心指標，平衡體系后使用腳本計算S蛋白?ACE2相互作用的能量項，該方法基于具有泊松?玻爾茲曼表面積的分子力學（MM?PBSA）方法［7,14?15］。MM?PBSA計算結合自由能的思想如圖2所示［16］。

圖2 MM?PBSA熱力學循環(huán)情況

ΔGSolv,Rec、ΔGSolv,Lig和ΔGSolv,Comp分別為氣相中的受體、配體和復合物進入液相中時溶劑化能的變化；而ΔGSolv,bind和ΔGGas,bind分別為液相中和氣相中的結合自由能。

理想情況下，如路徑1使用在溶劑中的受體和配體結合直接計算結合自由能，但是在這些溶劑化狀態(tài)的模擬中，大部分能量貢獻來自溶劑?溶劑相互作用，并且總能量的波動將比結合能大一個數(shù)量級，因此結合能的計算需要非常長的時間才能收斂。所以更有效的方法是把熱力學循環(huán)下的各部分分開計算。路徑2的簡化過程是先將受體和配體在氣相中結合，然后再把結合物放在溶劑環(huán)境中［17?19］。從圖2 可以看出，兩者的路徑始末狀態(tài)相同，則圖中的能量關系如下［20］：

式（3）用能量類型進行表達可得：

式（4）中，焓變（ΔH）可以分解為氣相能（ΔEMM）和溶劑化能（ΔGSol）［21］，忽略熵（?TΔS）對ΔGBind的貢獻，焓變公式表達為

式（5）中，氣相能（ΔEMM）可被分解為靜電項（ΔEele）、范德華項（ΔEvdW）和內(nèi)能項（ΔEint）。ΔEint內(nèi)能由鍵能、鍵角能和二面角能組成。ΔEMM分解表達如下式：

式（5）中的ΔGSol是極性溶劑化能（ΔGpb）和非極性溶劑化能（ΔGnp）之和：

2021年，Zhou等［22］計算了野生和突變體SARS?CoV?2 RBD?hACE2復合物的分子力結合自由能。N439K?hACE2的結合能要高于野生型RBD?hACE2，野生型和突變型之間的靜電能量有顯著差異。為了保證模擬的可靠性，Zhou等［22］將MD模擬時間增加到200 ns，200 ns MD模擬的結合能顯示出與100 ns 類似的結果，N439K?hACE2 的結合能（?1597.25 ±179.57 KJ/mol）也高于野生型的結合能（?959.29±130.36 KJ/mol）。同時，Zhou等［22］用MMPBSA計算了RBM（the receptor?binding motif）和hACE2之間的結合能，兩個復合物的結合能分別為?860.94±99.45 KJ/mol和?398.06±111.76 KJ/mol。比較RBD?hACE2和RBM?hACE2的結合自由能，發(fā)現(xiàn)能量的變化主要集中在RBM?hACE2 區(qū)域，說明了S蛋白與hACE2的主要結合點位［23］。

2 用分子動力學模擬研究刺突蛋白和中和抗體的結合

Kwarteng 等［24］使用GROMACS v2021 和CH?ARMM36 力場（模擬野生型和D614G S?蛋白系統(tǒng)）。將蛋白質(zhì)放置在一個立方體盒子中居中，置于盒子邊緣1 nm 處，使用TIP3P 水模型溶劑化。向系統(tǒng)中加入適當?shù)碾x子以靜電中和分子系統(tǒng)。使用最速下降算法對系統(tǒng)進行能量最小化。系統(tǒng)的蛋白質(zhì)和非蛋白質(zhì)組分獨立耦合到vrescale恒溫器和各向同性Berendsen算法進行壓力耦合［25］。通過弱耦合至1 bar 的參考壓力維持壓力，使用LINCS 算法限制蛋白質(zhì)內(nèi)的鍵長。NVT 和NPT 平衡總共進行了400 ps。對所有分子系統(tǒng)進行了三次各自100 ns的獨立模擬［26］。

根據(jù)GROMACS 指令求出野生型和D614G S?蛋白的平均RMSD 值分別為1.09±0.04 nm 和1.26±0.1 nm。野生型和D614G S 蛋白的結構相似性和構象穩(wěn)定性存在顯著差異，實驗探究D614G S蛋白和野生型的結構相似性，并通過計算每個S 蛋白結構域的RMSD 來研究構象穩(wěn)定性［27］。監(jiān)測突變對S?蛋白構象穩(wěn)定性的影響有助于研究病毒與中和抗體之間的相互作用。計算了野生型、D614G 和開放態(tài)S 蛋白的NTD、RBD 和S2 結構域的RMSD［28］。結果顯示野生型和D614G S 蛋白的RBD 比開放態(tài)RBD 具有更差的構象穩(wěn)定性［29］。與野生型和開放狀態(tài)相比，D614G S 蛋白的S2 結構域顯示出最高的RMSD值。這表明，S蛋白骨架原子的平均位移主要由S2 結構域的構象驅動［30］。突變對S 蛋白結構域構象的影響在S2 結構域中比NTD 和RBD 更顯著?？偟膩碚f，骨架RMSD 和結構域特異性RMSD 表明D614G S 蛋白采用了獨特的構象，但似乎比野生型構象更偏向于開放態(tài)構象，而比較野生型和D614G S 蛋白的RBD 區(qū)域會發(fā)現(xiàn)，突變會影響RBD的結構動力學。

3 分子動力學模擬在SARS?CoV?2中的其他應用

3.1 分子動力學模擬在病毒酶中的應用

由于該流行病仍在全球范圍內(nèi)蔓延，而負責病毒復制的酶又在病毒變異中起重要作用，因此采用模擬計算方法來確定這種酶的潛在抑制劑具有十分重要的意義。Parise 等［31］研究了有吸引力的抗病毒核苷酸類似物RNA 聚合酶（RdRp）鏈終結者抑制劑。研究以分子動力學模擬為基礎，探討了內(nèi)源性合成的三磷酸核苷（ddhCTP）與天然核堿基三磷酸環(huán)（CTP）在RdRp中結合的重要方面。

ddhCTP 種類是毒蛇素抗病毒蛋白的產(chǎn)物，是先天免疫反應的一部分。ddhCTP中核糖3'?OH的缺失可能對其抑制RdRp 的機制有重要影響。Parise 等［31］用實驗方法建立了一個嵌入RdRp 的RNA 鏈的硅學模型，從低溫電子顯微鏡結構開始，利用光譜法獲得了RNA 序列的信息。該模型在穩(wěn)定MD模擬計算后，獲得的結果為三磷酸核苷的結合提供了更深入的見解，但是RdRp活性位點的分子機制仍然是難以捉摸的。

3.2 分子動力學模擬在發(fā)現(xiàn)藥物靶點中的應用

最近，Yu等［32］將TLR9人體蛋白列為研究治療SARS?CoV?2 藥物靶點的目標蛋白，但此類治療的特異性和療效尚不清楚。團隊通過由相互作用化學物質(zhì)搜索工具（STITCH）、京都基因和基因組百科全書（KEGG）、分子對接和病毒?宿主藥物相互作用組映射組成的框架，對重組藥物進行了研究。以氯喹（CQ）和羥氯喹（HCQ）為探針，探索與SARS?CoV?2 相關的相互作用網(wǎng)絡。47個藥物靶點與SARS?CoV?2網(wǎng)絡重疊，并富含TLR 信號通路。分子對接分析和分子動力學模擬確定了TLR9 人體蛋白與CQ 和HCQ 的直接結合親和力［33］。此外，該團隊還建立了SARS?CoV?2?人類?藥物?蛋白質(zhì)相互作用圖，不僅為SARS?CoV?2 感染和病毒機制提供關鍵的見解，而且在發(fā)現(xiàn)新藥物靶點方面起重要作用［34］。

3.3 分子動力學模擬在疫苗設計中的應用

隨著2019 年冠狀病毒大流行的暴發(fā)，Mas?somi 等［35］專注于尋找治愈冠狀病毒疾病的解決方案上。在所有的疫苗接種策略中，納米顆粒疫苗已被證明可以使用于免疫系統(tǒng)，并在單次劑量中提供給免疫系統(tǒng)最佳的免疫狀態(tài)。鐵蛋白是一種用于疫苗生產(chǎn)的納米顆粒蛋白，目前已用于實驗研究。此外，糖基化在抗體和疫苗的設計中起著至關重要的作用，是開發(fā)有效疫苗的重要因素。在這項計算研究中，鐵蛋白納米顆粒和糖基化是疫苗設計的兩個獨特方面，首次用于對改進的納米顆粒疫苗進行建模［36］。在這方面，進行了分子建模和分子動力學模擬，以構SARS?CoV?2 RBD 區(qū)?鐵蛋白納米顆粒疫苗的三個原子模型，包括未糖基化、糖基化和在鐵蛋白?RBD 界面用O?聚糖改進的疫苗。結果表明，當聚糖添加到鐵蛋白?RBD 界面時，鐵蛋白?RBSD 復合物變得更穩(wěn)定，并且可以實現(xiàn)該納米顆粒的最佳性能。如果通過實驗驗證，這些發(fā)現(xiàn)可以改進納米顆粒的設計，以抵抗所有微生物感染。

RBD 鐵蛋白納米顆粒是最有效的疫苗之一，可在單次劑量內(nèi)提供對病毒的最佳免疫［37］。在這項工作中，分子建模和MD模擬將自組裝鐵蛋白納米顆粒支架與糖基化相結合，以改進目前可用的鐵蛋白RBD 疫苗。對SARS?CoV?2特異性鐵蛋白–RBD 納米顆粒疫苗的三種狀態(tài)進行了建模，包括未糖基化、糖基化和在鐵蛋白–RBS 界面用O?聚糖改進［38］。結果表明，糖基化通常保持納米顆粒的穩(wěn)定性。通過引入包括額外糖基化位點的修飾環(huán)，穩(wěn)定性顯著提高。這些發(fā)現(xiàn)可能對改進目前可用的鐵蛋白–RBD 納米顆粒疫苗以及未來針對各種冠狀病毒科的納米顆粒疫苗設計至關重要［39］。

4 面臨的挑戰(zhàn)

本文講述了分子動力學模擬在SARS?CoV?2中的應用情況，特別是SARS?CoV?2 S 蛋白突變前后和人體ACE2 在SARS?CoV?2 S 蛋白突變前后和中和抗體的結合情況，通過分子動力學模擬計算得出RMSD、RMSF和結合自由能這幾項指標來說明體系的結合能力［40］，從近些年的研究中可以發(fā)現(xiàn)，分子動力學模擬是分析SARS?CoV?2病毒情況很流行一種手段，但是SARS?CoV?2 病毒的不斷突變，數(shù)不勝數(shù)的突變體接踵而至，這為研究整個SARS?CoV?2 病毒的突變情況帶來了很大的挑戰(zhàn)。

5 結語

SARS?CoV?2 的突變導致疫情不斷嚴重，計算機模擬技術在S蛋白突變逃逸、基于酶結構的虛擬篩選、疫苗設計等方面均得到了廣泛的應用。其中MD模擬技術通過對野生型和突變型的S 蛋白的RMSD、RMSF、結合自由能分析，發(fā)現(xiàn)突變可導致S 蛋白與ACE2 蛋白結合能力增強，其與中和抗體的結合能力減弱，從原子水平上解釋了免疫逃逸的機制。計算機模擬技術與實驗研究的進一步結合必將為新冠病毒的研究帶來新的突破。