利用Markov模型研究Aβ突變體的構象轉變過程

姚星宇, 劉瓔銳, 韓葳葳, 萬由衷

(1. 吉林大學 白求恩第三醫(yī)院, 長春 130033; 2. 吉林大學 生命科學學院, 長春 130012)

阿爾茲海默癥(AD)是一種進行性神經退行性疾病, 主要影響大腦的記憶、 思維和行為功能, 是老年人最常見的癡呆癥之一[1]. 阿爾茲海默癥的臨床表現包括記憶力減退, 尤其是新近記憶受損. 阿爾茲海默癥的發(fā)病機制尚不完全清楚, 但有幾個關鍵因素被認為是導致疾病的原因之一[2-3], 其中一個因素是β-淀粉樣蛋白(Aβ)的異常聚集[3]. 對阿爾茲海默癥的治療方法展開深入研究已引起人們廣泛關注[4]. Aβ是一種在大腦中產生的蛋白質片段, 正常情況下會被清除掉. 而在阿爾茲海默癥患者中, Aβ會形成斑塊, 干擾神經元的正常功能[5]. 此外, 在阿爾茲海默癥患者的大腦中, Tau蛋白會異常地聚集成纖維纏結, 干擾神經元之間的正常信號傳遞[6]. 炎癥反應、 氧化應激和神經元死亡也與阿爾茲海默癥的發(fā)病機制有關.

由于Aβ是由跨膜淀粉樣前體蛋白(APP)經β-分泌酶與γ-分泌酶共同作用下水解, 并由39～43個殘基結合而成. 因此存在一種淀粉質沉淀假說[7-8], 內容是Aβ聚集形成的寡聚物具有神經毒性, Aβ在大腦中的產生速率大于清除速率, 導致該神經毒性寡聚物聚集致病[9]. 因此, 對Aβ聚集過程進行抑制, 即可消除神經毒性物質的產生, 從而對阿爾茲海默癥起到一定治療作用[1,10].

研究表明, Aβ聚集的過程中, 有由α-螺旋或無規(guī)卷曲到β-折疊的構象轉變現象, 且該現象被證實為聚集過程的關鍵步驟[11]. 因此, 通過殘基突變抑制β-折疊的形成, 即可抑制Aβ聚集成為神經毒性物質[12-14]. 其中, 將4種甘氨酸(G25,G29,G33,G37)全部突變?yōu)楸彼? 甚至可以幾乎完全消除β-折疊結構的產生[15-16]. 但該構象轉變發(fā)生的十分迅速, 通過常規(guī)實驗方法無法進行跟蹤式檢測, 難以對其機理進行更進一步研究. 因此, 本文通過使用計算生物學的方法進行實驗, 結合Markov模型和分子動力學模擬, 對野生型、 A4型和D7N型Aβ的構象變化過程進行研究.

1 實 驗

1.1 分子動力學模擬

計算生物學中的分子動力學模擬是將物理、 數學和化學結合的一門學科技術, 能使物質的宏觀和微觀性質相結合從而進行研究, 其主要原理是通過牛頓力學關系計算, 對分子體系的運動過程進行模擬仿真[17]. 分子動力學模擬的一般步驟包括:

1) 定義系統(tǒng). 確定模擬的分子系統(tǒng), 包括分子的種類、 初始位置和速度等.

2) 確定勢能場. 選擇適當的勢能場描述分子之間的相互作用, 如分子力場.

3) 運動方程求解. 根據牛頓運動定律和勢能場, 求解分子的運動方程, 通常使用數值積分方法.

4) 模擬時間演化. 通過迭代求解運動方程, 模擬分子在一定時間范圍內的運動和相互作用.

5) 數據分析. 對模擬結果進行分析, 如計算物理量和構建分子軌跡等, 以獲得關于分子系統(tǒng)的信息.

1.2 Markov模型

Markov模型是一種用于描述隨機過程的數學模型, 它基于Markov性質, 即未來狀態(tài)僅依賴于當前狀態(tài), 與過去的狀態(tài)無關. 構建Markov模型的一般步驟包括: 定義狀態(tài)空間, 確定系統(tǒng)的所有可能狀態(tài); 確定狀態(tài)轉移概率, 計算狀態(tài)之間的轉移概率, 即從一個狀態(tài)轉移到另一個狀態(tài)的概率; 得到狀態(tài)轉移矩陣, 將狀態(tài)和轉移概率表示為一個矩陣, 其中每個元素代表從一個狀態(tài)到另一個狀態(tài)的概率; 進行狀態(tài)預測, 利用狀態(tài)轉移矩陣進行狀態(tài)預測, 可以通過給定的初始狀態(tài)和轉移概率計算未來狀態(tài)的概率分布.

1.3 實驗流程總述

圖1為本文實驗流程. 首先獲取野生型Aβ、 兩種突變體A4型Aβ和D7N型Aβ的軌跡, 用數組數據表示軌跡, 再用時滯獨立組分分析對數據進行降維處理, 并用k-均值聚類獲得準確的聚類結果, 然后確定時滯時間(lag time), 構建Markov模型, 最后將分子軌跡進行PCCA(perron cluster cluster analysis)聚類, 將上述獲得的狀態(tài)結果和轉移概率矩陣用于通量分析.

圖1 實驗流程Fig.1 Flow chart of experiment

2 Aβ及突變體構象轉變過程

2.1 3個體系的分子動力學模擬

為研究野生型Aβ(PDB id: 2M4J)[18]、 兩種突變體A4型Aβ(G25,G29,G33,G37)和D7N型Aβ的動態(tài)構象變化過程, 先進行分子動力學(MD)模擬. 對每個體系進行10次1 000 ns的顯式水模型模擬以探索其動態(tài)行為. 選擇GROMACS 2018.3[19]作為MD軟件包運行, 并采用Gromos54a7力場[20]描述分子間相互作用. 通過GROMACS自帶的x2top組件, 生成3種Aβ的拓撲結構文件, 并將其嵌入到周期性邊界模擬盒中, 盒子的尺寸為10 nm×10 nm×10 nm. 模擬溫度設定為36 ℃. 為確保體系的穩(wěn)定性, 先進行100 ps的預平衡MD模擬, 采用Berendsen恒溫器[21]在微正則(NVT)系綜下進行. 之后在等溫等壓(NPT)系綜下, 保持體系的壓強為100 kPa(1 bar), 溫度為36 ℃, 對3個體系進行模擬. 范德華作用截斷半徑設為1 nm. 在所有模擬過程中, 采用2 fs的時間步長, 并每10 ps保存一次原子坐標, 共得到30條軌跡用于后續(xù)的分析和研究.

2.2 Aβ體系的時滯獨立組分分析

應用pyemma軟件包構建Markov模型的全過程. 首先計算野生型Aβ、 兩種突變體A4型Aβ和D7N型Aβ的相互作用軌跡, 在平行條件下構建10份軌跡, 每份模擬時長為90 ns, 時間間隔為0.1 ns. 將得到的10份軌跡分別導入pyemma軟件包. 導入時選擇處殘基骨架的phi角和psi角信息為軌跡信息, 并借此對軌跡的每幀進行量化, 將分子動力學模擬軌跡的信息用數據組表示, 實現從抽象化轉為具體化. 10條軌跡被分別轉化為10個二維數組. 之后, 分別對10個二維數組進行時滯獨立組分分析降維, 降維結果至少保留95%的方差信息. 最終降維至82維度. 隨后將降維結果進行k-均值聚類, 聚類結果如圖2所示. 由圖2可見, 野生型Aβ和A4型Aβ的軌跡前后變化較小, 僅有一小部分與整體分離. 其原因可能是野生型Aβ和A4型Aβ的α-螺旋基本不發(fā)生變化. 最終野生型Aβ和A4型Aβ選取合適的k值為30, D7N型Aβ選取合適的k值為45.

圖2 k-均值聚類結果Fig.2 k-Means clustering results

2.3 Aβ體系時滯時間數值的確定

采用msm.its方法計算系統(tǒng)弛豫時間, 將系統(tǒng)弛豫時間標度與時滯時間進行繪圖, 結果如圖3所示. 由圖3可見, 對于野生型Aβ, 當時滯時間約為0.2 ns時, 10個軌跡系統(tǒng)的弛豫時間尺度呈收斂趨勢. 由于軌跡系統(tǒng)的間隔時間設定為0.1 ns, 因此弛豫時間為2 steps. 由于A4型Aβ和D7N型Aβ軌跡系統(tǒng)的間隔時間設定為0.1 ns, 10個軌跡系統(tǒng)的弛豫時間尺度呈收斂趨勢, 時滯時間約為0.5 ns, 因此弛豫時間為5 steps. 明確時滯時間后, 開始構建Markov模型, 并對轉移概率進行計算, 得到轉移概率矩陣.

圖3 野生型Aβ(A), A4型Aβ(B)和D7N型Aβ(C)時滯時間確定結果Fig.3 Results of determining delay time of wild type Aβ (A), A4 type Aβ (B) and D7N type Aβ (C)

2.4 構建Markov模型和結果分析

2.4.1 構建Markov模型和CK檢測

根據獲得的聚類結果和時滯時間構建Markov模型, 并計算轉移概率矩陣. 為進一步驗證所構建模型的Markov性質, 使用CK(Chapman-Kolmogorov test)進行檢測, CK檢測結果如圖4所示. 在該過程中, 使用pyemma中的msm.cktest方法, 并將參數mlags設為4, 即CK檢測算法中的倍數k, 同時將nstates設置為5, 表示后續(xù)PCCA聚類的宏觀狀態(tài)數. 由圖4可見, 各坐標系的Markov模型預測與實驗前估測都近于擬合, 表明時滯時間的選取及Markov模型構建均較合理, 說明構建的Markov模型具有良好的Markov性質, 可以用于開展后續(xù)研究.

圖4 Aβ的CK檢測結果Fig.4 CK-test results of Aβ

2.4.2 野生型Aβ的β-折疊構象轉變

使用PCCA聚類可對軌跡結果進行可視化分析, 并將分子軌跡系統(tǒng)用宏觀狀態(tài)進行描述. 采用msm.pcca方法對野生型Aβ的β-折疊構象轉變進行研究. 在PCCA聚類結果中, SA表示模擬軌跡中構象的起始狀態(tài), SB表示模擬軌跡中構象的最終收斂狀態(tài), S1,S2,S3表示按時間順序出現的中間狀態(tài). 各聚類狀態(tài)如圖5所示, 紫色標記區(qū)域為主要研究的α-螺旋區(qū)域. 野生型Aβ的通量分析結果列于表1.

表1 野生型Aβ通量分析結果

由表1可見, 在野生型Aβ的β-折疊構象轉變過程中, SA→S3→SB概率為33.49%, 是最具有代表性的路徑, 其次是SA→S3→S2→SB路徑, 占總通量的28.71%. SA→S1→S3→SB構象轉變路徑占比也較大, 占總通量的15.79%. 分析其殘基具體轉變過程, 得出各狀態(tài)代表構象幀的具體二級結構如下: SA→S3→S2→SB與SA→S1→S3→SB兩種構象轉變途徑大致的二級結構轉變與SA→S3→SB的結構轉變基本類似, 都出現一個β-折疊的過程, 且通過觀察各宏觀狀態(tài)的具體構象發(fā)現, 整個過程經歷S3的一個中間微狀態(tài), 且S3在殘基中出現了π-螺旋結構.

2.4.3 A4型Aβ的β-折疊構象轉變

A4型Aβ宏觀狀態(tài)構象如圖6所示. 由圖6可見, A4型Aβ的各PCCA聚類宏觀狀態(tài)的特征分布極接近, 4個宏觀狀態(tài)間的構象差距較小, 且紫色強調處都保持著α-螺旋結構.

A4型Aβ的通量分析結果列于表2. 由表2可見, 在A4型Aβ的β-折疊構象轉變過程中, SA→SB概率為43.96%, 是最具有代表性的路徑. 通過觀察宏觀狀態(tài)的構象和二級結構的計算結果可見, 宏觀狀態(tài)之間的差異非常小. 從S1到SB的二級結構基本沒有發(fā)生變化, SB與初始狀態(tài)SA存在一定差異, 但S1到SB仍保持螺旋結構, 表明螺旋只有一部分發(fā)生了解旋. 發(fā)生解旋的部分主要為該螺旋結構的N端, 而C端仍保持螺旋結構. 綜上, A4型Aβ的4個微狀態(tài)均未出現β-折疊, 說明A4型突變有顯著抑制β-折疊構象的特性.

表2 A4型Aβ的通量分析結果

2.4.4 D7N型Aβ的β-折疊構象轉變

D7N型Aβ宏觀狀態(tài)構象如圖7所示. D7N型Aβ的通量分析結果列于表3. 由表3可見, 在D7N型Aβ的β-折疊構象轉變過程中, SA→S1→SB概率為50.51%, 是最具有代表性的路徑. 觀察二級結構可知, S2和SB兩個微狀態(tài)均出現β-折疊. 由圖7可見, 微狀態(tài)SA代表的圓型面積很小, 說明符合微狀態(tài)SA的構象很少, 可能意味著初始微狀態(tài)SA發(fā)生變化的時間很短, 更進一步證明了D7N型突變對β-折疊轉變過程的促進作用.

表3 D7N型Aβ的通量分析結果

圖7 D7N型Aβ宏觀狀態(tài)構象Fig.7 Macroscopic state conformation of D7N type Aβ

綜上可見, 在野生型Aβ的β-折疊構象轉變過程中, SA→S3→SB路徑出現概率為33.49%, 概率最大, 較具有研究意義. 整個過程經歷S3一個中間微狀態(tài), 且微狀態(tài)S3在殘基中出現了π-螺旋結構, 此結構為野生型Aβ中發(fā)現的獨有結構.

在A4型Aβ的β-折疊構象轉變過程中, 4個微狀態(tài)構象極相似, 幾乎可認為沒有發(fā)生變化, 完全由α-螺旋與無規(guī)卷曲組成. 其中SA→SB路徑出現概率最大, 為43.96%, 說明該構象轉變較大概率沒有經過任何中間狀態(tài)的改變, 更證實了構象幾乎無變化. 因此, A4型突變可極大程度抑制β-折疊構象轉變的發(fā)生, 通過抑制β-折疊抑制Aβ聚集為細胞毒性物質, 為開發(fā)抗阿爾茲海默癥藥物提供一定理論依據.

在D7N型Aβ的β-折疊構象轉變過程中, SA→S1→SB概率為50.51%, 概率最大, 較具有研究意義. 其中, 最終構象SB的α-螺旋近乎完全解旋, 為β-折疊的形成提供有利條件. 相比于野生型Aβ只發(fā)生一個區(qū)域的β-折疊, D7N型突變體發(fā)生了兩個區(qū)域的β-折疊, 說明D7N型突變可促進β-折疊特性. 而且D7N型Aβ中SA微狀態(tài)的圓面積很小, 表明符合SA微狀態(tài)的構象數量較少, 說明初始構象轉變的速度非?？? 因此D7N突變可能會導致更快的β-折疊構象形成, 從而在阿爾茨海默癥的發(fā)病機制中發(fā)揮重要作用. 該發(fā)現為深入理解D7N突變對蛋白質結構和功能的影響提供了有價值的線索.