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

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

0 引言

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

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

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

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

圖1 計(jì)算結(jié)合自由能的一般流程

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

1 用分子動(dòng)力學(xué)模擬研究刺突蛋白和ACE2的結(jié)合

1.1 刺突蛋白和ACE2復(fù)合物的RMSD和RMSF

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

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

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

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

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

通過其實(shí)驗(yàn)發(fā)現(xiàn)WT、Mu 和Omicron 變體在30 ns 后趨于穩(wěn)定，RMSD 分別為4.61 ?、4.84 ?和5.15 ?，而Delta 在65 ns 后才達(dá)到平衡，均方根偏差為6.75 ?。這表明這三種變體的構(gòu)象穩(wěn)定性不如WT，尤其是Delta［13］。為了分析不同溫度對(duì)S蛋白結(jié)構(gòu)穩(wěn)定性的影響，在三種不同溫度（即300 K、310 K 和320 K）下進(jìn)行了比較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 ?。從結(jié)果可以看出Mu 最為穩(wěn)定而Omicron 結(jié)構(gòu)也比WT 要穩(wěn)定，但是Delta 變種的穩(wěn)定性要比WT差。同樣從RMSF的結(jié)果也可以看出幾種變體穩(wěn)定性從大到小依次是Mu、Omicron 和Delta，但是都要比WT更加穩(wěn)定。

1.2 刺突蛋白和ACE2復(fù)合物的結(jié)合自由能

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

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

ΔGSolv,Rec、ΔGSolv,Lig和ΔGSolv,Comp分別為氣相中的受體、配體和復(fù)合物進(jìn)入液相中時(shí)溶劑化能的變化；而ΔGSolv,bind和ΔGGas,bind分別為液相中和氣相中的結(jié)合自由能。

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

式（3）用能量類型進(jìn)行表達(dá)可得：

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

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

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

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

2 用分子動(dòng)力學(xué)模擬研究刺突蛋白和中和抗體的結(jié)合

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

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

3 分子動(dòng)力學(xué)模擬在SARS?CoV?2中的其他應(yīng)用

3.1 分子動(dòng)力學(xué)模擬在病毒酶中的應(yīng)用

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

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

3.2 分子動(dòng)力學(xué)模擬在發(fā)現(xiàn)藥物靶點(diǎn)中的應(yīng)用

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

3.3 分子動(dòng)力學(xué)模擬在疫苗設(shè)計(jì)中的應(yīng)用

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

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

4 面臨的挑戰(zhàn)

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

5 結(jié)語

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