夏 翾，馬 帥，王 勤，李曉琴

(北京工業(yè)大學(xué)生命科學(xué)與生物工程學(xué)院，北京100124)

蛋白酶是一種天然的生物催化劑，具有較高的催化效率和底物特異性，被廣泛應(yīng)用于工業(yè)、食品業(yè)、醫(yī)藥業(yè)及人們的日常生活中［1-2］。蛋白酶發(fā)揮其高活性和高催化性是需要特定的條件作為保障的，如適當(dāng)?shù)腜H值、溫度、離子濃度等［3］。其中，溫度這一因素尤為重要。蛋白酶只有在一定的溫度范圍內(nèi)才能保持較高的活性，超出了這個(gè)溫度范圍，酶的活性就會(huì)降低甚至消失，嚴(yán)重影響它們的實(shí)際應(yīng)用。因此，了解蛋白酶的熱穩(wěn)定機(jī)制進(jìn)而提高蛋白酶的熱穩(wěn)定性［4］，對(duì)擴(kuò)大酶的應(yīng)用范圍與探究蛋白質(zhì)結(jié)構(gòu)和功能關(guān)系均具有非常重要的意義。

蛋白質(zhì)的熱穩(wěn)定性是由多種作用共同參與的，不同的蛋白質(zhì)可能有不同的耐熱機(jī)制，在這些作用中比較重要的有:氨基酸組成、氫鍵、鹽橋、疏水作用、包埋效應(yīng)、環(huán)和N、C端的穩(wěn)定性等［5］。目前，提高枯草桿菌蛋白酶熱穩(wěn)定性的方法主要有兩種，一種是改變酶的催化環(huán)境，如調(diào)節(jié)PH值、增加有機(jī)溶劑、底物保護(hù)等，這種方法不能從本質(zhì)上改變枯草桿菌蛋白酶的耐熱性，因此在工業(yè)應(yīng)用上會(huì)受到很多的限制;另一種是對(duì)蛋白酶進(jìn)行殘基突變，在蛋白質(zhì)的薄弱環(huán)節(jié)引入有利殘基，加強(qiáng)局部或整體的彼此相互作用，從而增加蛋白的熱穩(wěn)定性。

殘基突變是提高蛋白酶熱穩(wěn)定性的有效方法，如何評(píng)判殘基突變對(duì)蛋白酶熱穩(wěn)定性的影響?目前，對(duì)蛋白酶嗜熱性突變效果的判定主要有兩種方法，一種是生物實(shí)驗(yàn)的方法［6］，即將突變成功的蛋白通過蛋白質(zhì)合成技術(shù)制造出來，然后通過逐漸升溫來驗(yàn)證其突變效果的優(yōu)劣。該方法有較高準(zhǔn)確性與可靠性，但是會(huì)耗費(fèi)大量的時(shí)間和人力。另一種是生物信息學(xué)計(jì)算的方法［7-9］，即計(jì)算ΔΔGu值來評(píng)定嗜熱性突變的效果。ΔGu(蛋白質(zhì)解折疊自由能)是蛋白質(zhì)熱動(dòng)力學(xué)參數(shù)中最重要的參數(shù)，也是反應(yīng)蛋白質(zhì)熱穩(wěn)定性的一個(gè)通用指標(biāo)，因此多數(shù)方法都是利用生物信息學(xué)計(jì)算的方法來模擬突變前后ΔGu的變化，即ΔΔGu(式1)。當(dāng)ΔΔGu＞0時(shí)，表明該突變對(duì)蛋白質(zhì)的熱穩(wěn)定性有正效應(yīng)(即嗜熱性突變效果好)，反之有負(fù)效應(yīng)(即嗜熱性突變效果差)。

2005年，Korkegian等直接從蛋白質(zhì)結(jié)構(gòu)入手預(yù)測(cè)蛋白質(zhì)的ΔGu，該方法主要是利用一些經(jīng)驗(yàn)數(shù)據(jù)或方法來分析蛋白質(zhì)突變前后其結(jié)構(gòu)或能量變化，來推斷其 ΔGu的變化［10］。2006年，Saraboji等利用統(tǒng)計(jì)數(shù)據(jù)分析預(yù)測(cè)蛋白質(zhì)的ΔGu，即用一些統(tǒng)計(jì)學(xué)的手段對(duì)突變體數(shù)據(jù)進(jìn)行分析，得到一些統(tǒng)計(jì)學(xué)的規(guī)律，以此來判定突變前后的 ΔΔGu［11］。2008 年，Capriotti等提出機(jī)器學(xué)習(xí)方法，即利用蛋白質(zhì)的序列或結(jié)構(gòu)特征直接從數(shù)據(jù)入手，建立計(jì)算模型來預(yù)測(cè)ΔΔGu［12］。該方法較生物實(shí)驗(yàn)方法節(jié)約了時(shí)間和人力，但是其仍具有計(jì)算復(fù)雜、計(jì)算量大的缺點(diǎn)。

分子動(dòng)力學(xué)(MD)模擬是研究蛋白質(zhì)動(dòng)力學(xué)性質(zhì)進(jìn)而揭示其生物學(xué)功能的重要工具，利用MD模擬結(jié)果提取有效參數(shù)、建立蛋白酶嗜熱性改造突變效果評(píng)判方法是本文研究的重點(diǎn)。本文以研究資料頗為豐富的枯草桿菌蛋白酶BPN＇為研究對(duì)象，選取實(shí)驗(yàn)上已證明的能提高其熱穩(wěn)定性的突變體來進(jìn)行研究，運(yùn)用分子動(dòng)力學(xué)模擬的方法，通過對(duì)突變前后蛋白質(zhì)多種能量參數(shù)的對(duì)比分析，提取了多個(gè)有效參數(shù)，并進(jìn)一步建立了枯草桿菌蛋白酶嗜熱性改造單突變效果評(píng)判方法。

1 材料

1.1 枯草桿菌蛋白酶

枯草桿菌蛋白酶是芽孢桿菌屬細(xì)菌所分泌的胞外堿性蛋白酶，屬絲氨酸蛋白水解酶類［13-16］。它催化水解蛋白質(zhì)為氨基酸，在有機(jī)溶劑中也能催化多肽的合成。常見的枯草桿菌蛋白酶有 BPN＇、Carlsbery、E、Savinase、BL 等，它們的一級(jí)氨基酸序列非常類似，晶體結(jié)構(gòu)也基本相似［17］?？莶輻U菌蛋白酶BPN＇與E的區(qū)別是43個(gè)氨基酸殘基發(fā)生替換，Carlsberg與E在序列上的區(qū)別是89個(gè)氨基酸殘基替換再加上第56位殘基缺失，BL與BPN＇的區(qū)別是在275個(gè)氨基酸殘基上有103個(gè)殘基替換再加上4個(gè)殘基(160～163)缺失。

枯草桿菌蛋白酶具有重要的應(yīng)用價(jià)值，被廣泛應(yīng)用于洗滌劑、制革及絲綢等工業(yè);此蛋白的結(jié)構(gòu)與功能研究已有較長(zhǎng)的歷史，其蛋白結(jié)構(gòu)、生化特性、催化機(jī)理等均有了較清晰的闡釋。目前，在該酶的嗜熱性改造方面有豐富的突變資料可供參考;另外，解析度為1.8 ?的該蛋白酶的晶體結(jié)構(gòu)數(shù)據(jù)可以從PDB數(shù)據(jù)庫(kù)中獲取?；谝陨显?，我們將枯草桿菌蛋白酶作為研究對(duì)象。

以枯草桿菌蛋白酶BPN＇作為研究對(duì)象，它共有275個(gè)殘基，空間結(jié)構(gòu)如圖1所示。

圖1 枯草桿菌蛋白酶空間結(jié)構(gòu)Fig.1 Spatial structure of Subtilisin BPN'

圖中桶狀結(jié)構(gòu)代表α螺旋，固形箭頭代表β折疊，小球代表催化三聯(lián)體(Asp32、His64和Ser221)。圖中可見β折疊位于蛋白質(zhì)內(nèi)部形成疏水核心，α螺旋位于結(jié)構(gòu)表面，圍繞在β折疊周圍，整個(gè)蛋白質(zhì)呈一個(gè)球狀。

1.2 枯草桿菌蛋白酶突變位點(diǎn)的選取

選取6個(gè)實(shí)驗(yàn)上已證明的能提高蛋白質(zhì)熱穩(wěn)定性的突變點(diǎn)作為研究對(duì)象，它們分別為:S161C、G166R［18］、G169A［19］、S188P［20］、K213R［21］以 及P239R［22］。表1為野生型枯草桿菌蛋白酶 BPN＇及其6個(gè)突變體的Tm值及60℃下的半衰期時(shí)間，這兩個(gè)值越大說明結(jié)構(gòu)的熱穩(wěn)定性越強(qiáng)。由表中數(shù)據(jù)可以看出，六個(gè)突變體的熱穩(wěn)定性均強(qiáng)于野生型。表中Wild代表野生型枯草桿菌蛋白酶。

2 方法

2.1 枯草桿菌蛋白酶突變體的同源模建

從PDB數(shù)據(jù)庫(kù)中下載解析度為1.8 ?的枯草桿菌蛋白BPN＇晶體結(jié)構(gòu)(PDB:2SIC)，將此結(jié)構(gòu)的E鏈作為野生型樣本。利用SWISS-MODEL進(jìn)行同源模建，得到表1所列的6個(gè)突變體的結(jié)構(gòu)。

表1 枯草桿菌蛋白酶及其突變體的信息統(tǒng)計(jì)Table 1 Information statistics of subtilisin BPN'and its mutants

利用野生型及6個(gè)突變體的結(jié)構(gòu)信息，分別計(jì)算主鏈結(jié)構(gòu)的二面角，繪制圖2所示的拉氏構(gòu)象圖。圖A為野生型蛋白的拉氏構(gòu)象圖，圖B、圖C、圖D、圖 E、圖 F、圖 G 分別為 S161C、G166R、G169A、S188P、K213R、P239R突變體的拉氏構(gòu)象圖。圖中黑色區(qū)域表示構(gòu)象允許區(qū)域，在該區(qū)域內(nèi)的任何成對(duì)二面角所規(guī)定的構(gòu)象都是立體化學(xué)所允許的?；疑珔^(qū)域表示不完全允許區(qū)，表示在這個(gè)區(qū)域內(nèi)的二面角所規(guī)定的主鏈構(gòu)象雖然是立體化學(xué)可允許的，但是不夠穩(wěn)定，白色區(qū)域是不允許區(qū)域。在圖中可以看到，構(gòu)象角絕大部分都處于構(gòu)象允許區(qū)域內(nèi)，說明這些突變體的構(gòu)象均處于合理的狀態(tài)。

圖2 枯草桿菌蛋白酶及其突變體所對(duì)應(yīng)的拉氏構(gòu)象圖Fig.2 The ramachandran figure of subtilisin BPN'and its mutants

2.2 分子動(dòng)力學(xué)模擬實(shí)驗(yàn)的設(shè)計(jì)

使用NAMD軟件對(duì)野生型和6個(gè)突變體蛋白進(jìn)行分子動(dòng)力學(xué)模擬，選用CHARMM27全原子力場(chǎng)。預(yù)處理過程中，首先刪掉晶體中的氫原子，再使用VMD軟件自帶的psfgen程序包補(bǔ)足氫原子，測(cè)得該晶體結(jié)構(gòu)的凈電荷為零;然后給該結(jié)構(gòu)添加一個(gè)厚度為10 ?的TIP3P型水分子層，最終得到的水盒子的尺寸是 65 ? ×66 ? ×68 ?。

預(yù)平衡過程中，選用NPT系綜。在1個(gè)標(biāo)準(zhǔn)大氣壓下，使用SHAKE算法約束氫鍵鍵長(zhǎng)，截?cái)喟霃綖?2 ?，非鍵相互作用在14 ?處截?cái)?。?duì)1個(gè)野生型及6個(gè)突變體共進(jìn)行了以下4次分子動(dòng)力學(xué)模擬過程:

(1)溫度恒定為370 K，步長(zhǎng)為2 fs/step，運(yùn)行100 000步，分子動(dòng)力學(xué)模擬時(shí)間為0.2 ns，此過程命名為370 K_0.2 ns;

(2)溫度恒定為370 K，步長(zhǎng)為2 fs/step，運(yùn)行1 000 000步，分子動(dòng)力學(xué)模擬時(shí)間為2 ns，此過程命名為370 K_2 ns;

(3)溫度恒定為370 K，步長(zhǎng)為1 fs/step，運(yùn)行100 000步，分子動(dòng)力學(xué)模擬時(shí)間為0.1 ns，此過程命名為370 K_0.1 ns;

(4)初始溫度為300 K，步長(zhǎng)為2 fs/step，運(yùn)行100 000步，每2 000步升高1 K，終止溫度為350 K，分子動(dòng)力學(xué)模擬時(shí)間為0.2 ns，此過程命名為300 K_0.2 ns。

其中370K_0.2ns作為標(biāo)準(zhǔn)組，其他3組作為對(duì)照組。

2.3 基于模擬結(jié)果的參數(shù)計(jì)算

每個(gè)野生型及其突變體分子動(dòng)力學(xué)模擬完成后均輸出一個(gè)“.log”日志文件，包含了模擬過程中的壓力、溫度以及本實(shí)驗(yàn)中需要用到的能量等信息。配置文件中設(shè)置每1 000步(即2 ps)輸出一次能量信息，而我們需要計(jì)算的是在模擬的整個(gè)過程中或者某一特定時(shí)間段內(nèi)的平均能量，因此具體方法如下:利用VMD軟件的能量計(jì)算腳本 namdstats.tcl，結(jié)合輸出的.log日志文件，提取蛋白體系的以下平均能量信息:BOND(鍵能)、ANGLE(鍵角能)、DIHED(二面角能)、ELECT(靜電能)、VDW(范德華力能)、KINETIC(動(dòng)能)。

3 結(jié)果與分析

3.1 模擬結(jié)果

依照上述方法，我們共得到4組動(dòng)力學(xué)模擬后的能量結(jié)果(見表2～表5)。每組表格中上半部分為提取的能量值，下半部分為每個(gè)突變體與野生型能量的差值以及能量的變化趨勢(shì)。

首先，以溫度為1個(gè)標(biāo)準(zhǔn)大氣壓下溫度為370 K、時(shí)間為0.2 ns的模擬作為標(biāo)準(zhǔn)組。表2是該組的模擬能量結(jié)果，我們提取的是每個(gè)結(jié)構(gòu)從模擬第一步到最后一步的能量平均值。由表中數(shù)據(jù)可以看出，6個(gè)突變體與野生型酶相比，每種平均能量值均出現(xiàn)了相同的變化趨勢(shì)，即 BOND增大、ANGLE增大、DIHED減小、ELECT減小、VDW 增大、KINETIC增大。

表2 370 K_0.2 ns模擬平均能量(kcal/mol)Table 2 Average energy data simulated-370 K_0.2 ns(kcal/mol)

標(biāo)準(zhǔn)組的模擬溫度為370 K，模擬時(shí)間為0.2 ns，這樣高溫、短時(shí)間的模擬并沒有使野生型及突變體動(dòng)力學(xué)模擬達(dá)到平衡，為了檢驗(yàn)非平衡態(tài)的能量統(tǒng)計(jì)結(jié)果的可靠性，我們?cè)跇?biāo)準(zhǔn)組370 K_0.2 ns模擬的基礎(chǔ)上進(jìn)行了較長(zhǎng)時(shí)間的模擬，溫度不變，模擬時(shí)間增加為標(biāo)準(zhǔn)組的10倍，即2 ns。圖3為該模擬下野生型及其突變體的RMSD曲線圖。可以看到，在1.0 ns～1.5 ns的時(shí)間段，7個(gè)結(jié)構(gòu)的RMSD曲線均呈現(xiàn)平穩(wěn)狀態(tài)，說明在此時(shí)間段內(nèi)達(dá)到了平衡狀態(tài)。因此我們提取了這段時(shí)間內(nèi)的能量平均值，并計(jì)算突變體與野生型的能量差值(見表3)?？梢钥吹剑碇型蛔凅w的6種平均能量值有著與標(biāo)準(zhǔn)組相同的變化趨勢(shì)。

圖3 枯草桿菌蛋白及其突變體2ns分子動(dòng)力學(xué)模擬過程中的RMSD曲線Fig.3 The RMSD curves ofSubtilisin BPN'and its mutants in 2ns in the molecular dynamic simulation

表3 370 K_2 ns模擬平均能量差值(kcal/mol)Table 3 D-value of average energy data simulated-370 K_2 ns(kcal/mol)

從表2、表3的結(jié)果看，模擬時(shí)間的長(zhǎng)短以及體系是否平衡對(duì)能量變化趨勢(shì)并沒有影響，為進(jìn)一步檢驗(yàn)?zāi)M時(shí)間對(duì)結(jié)果的影響，我們?cè)跍囟炔蛔兊那闆r下，進(jìn)行了更短時(shí)間的模擬，模擬時(shí)間縮減為標(biāo)準(zhǔn)組的一半，即0.1 ns，其能量差值見表4。表4中所有能量的變化都較標(biāo)準(zhǔn)組呈現(xiàn)相同的趨勢(shì)。

表4 370 K_0.1 ns模擬平均能量差值 (kcal/mol)Table 4 D-value of average energy data simulated-370 K_0.1 ns(kcal/mol)

最后，我們檢驗(yàn)了溫度對(duì)模擬結(jié)果的影響。在標(biāo)準(zhǔn)組的基礎(chǔ)上我們做了一組溫度均勻變化的模擬300 K_0.2 ns，時(shí)間同標(biāo)準(zhǔn)組設(shè)置一樣，模擬過程中將溫度從300 K均勻升高到350 K。能量差值結(jié)果見表5。由表中數(shù)據(jù)可以看出，在溫度變化的模擬條件下突變體的平均能量還是出現(xiàn)了與標(biāo)準(zhǔn)組一樣的變化趨勢(shì)。

表5 300 K_0.2 ns模擬平均能量差值 (kcal/mol)Table 5 D-value of average energy data simulated-300 K_0.2 ns(kcal/mol)

3.2 結(jié)果分析

從上面4個(gè)表可以看出，突變體與野生型酶相比，均出現(xiàn)了BOND增大、ANGLE增大、DIHED減小、ELECT減小、VDW增大、KINETIC增大的趨勢(shì)。研究表明，分子柔性的降低會(huì)帶來嗜熱性的增強(qiáng)［23］，因此突變體的分子柔性均低于野生型，正是由于柔性降低、剛性增強(qiáng)，所以會(huì)使得BOND(鍵能)增大、ANGLE(鍵角能)增大以及 DIHED(二面角能)減小;而所有突變體的ELECT(靜電能)都比野生型小，這與以前的實(shí)驗(yàn)中所說引入帶電殘基增強(qiáng)蛋白質(zhì)熱穩(wěn)定性是相一致的，靜電能越小穩(wěn)定性越強(qiáng)［24］;而VDW(范德華能)增大說明分子的極性也增大了，表明蛋白質(zhì)的嗜熱性與分子極性也密切相關(guān);突變體的KINETIC(動(dòng)能)雖然比野生型都有所增大，但變化范圍很小，說明這個(gè)能量可能與熱穩(wěn)定性的關(guān)系并不大，但我們?nèi)詴簳r(shí)保留這一項(xiàng)，等待后續(xù)實(shí)驗(yàn)進(jìn)行進(jìn)一步驗(yàn)證。

另外可以看到表中有一些數(shù)據(jù)與其他數(shù)據(jù)相差較大，我們將這些數(shù)據(jù)稱為極端數(shù)據(jù)，并用特殊格式將這些極端數(shù)據(jù)標(biāo)記出來。除去這些極端數(shù)據(jù)，我們得到以下突變體與野生型相比的能量變化范圍(單位:kcal/mol):ΔBOND 為［ 2.66，8.57］，ΔANGLE 為［3.32，32.18］，ΔDIHED 為［-16.68，-2.66］，ΔELECT 為［-1 375.41，-286.86］，ΔVDW 為［81 362.83，87 027.74］，ΔKINETIC 為［0.26，2.89］。

為進(jìn)一步分析能量變化與蛋白嗜熱性的關(guān)系，我們利用SPSS分析軟件分別對(duì)以上4組模擬中每個(gè)突變體的6種能量差值與其對(duì)應(yīng)的Tm或半衰期進(jìn)行了斯皮爾曼(Spearman)相關(guān)性分析，99%置信度下的雙側(cè)檢驗(yàn)的計(jì)算結(jié)果見表6，其中表的上半部分是G166R、G169A、S188P的Tm值與各個(gè)能量差值的絕對(duì)值之間的相關(guān)系數(shù)，下表下半部分是S161C、K213R、P239R的Half life與各個(gè)能量差值之間的相關(guān)系數(shù)。相關(guān)系數(shù)越大表示相關(guān)性越顯著，正數(shù)表示正相關(guān)，負(fù)數(shù)表示負(fù)相關(guān)。

表6 能量差值與TmHalf life的相關(guān)系數(shù)Table 6 The correlation coefficient between the energy difference and TmHalf life

從表6可以看出:ΔBOND、ΔELECT與突變體的Tm相關(guān)性最強(qiáng)，3組模擬時(shí)間大于0.1 ns的相關(guān)系數(shù)均為1，同時(shí)與突變體的半衰期也分別有一組相關(guān)系數(shù)為1，總體平均值為0.75。說明 ΔBOND、ΔELECT對(duì)Tm及半衰期的變化最為敏感;ΔANGLE略低于ΔBOND、ΔELECT與突變體的Tm及半衰期的相關(guān)性，相關(guān)系數(shù)為1的共有 3組，平均為0.687 5;ΔVDW與突變體的Tm及半衰期的相關(guān)系數(shù)均為正值，有兩組模擬時(shí)間大于0.1 ns的相關(guān)系數(shù)為1，總體均值為0.5，相關(guān)性一般;而ΔDIHED和ΔKINETIC與突變體的Tm及半衰期的相關(guān)性不穩(wěn)定，說明ΔDIHED和ΔKINETIC對(duì)Tm及半衰期的變化趨勢(shì)不敏感。

3.3 判定方法

從以上分析可以看出，在4組不同條件的模擬過程中，突變體的6個(gè)平均能量值與野生型相比均出現(xiàn)了相同的變化趨勢(shì)，結(jié)合相關(guān)性的計(jì)算結(jié)果，ΔBOND、ΔELECT、ΔVDW、ΔANGLE、ΔDIHED、ΔKINETIC可以作為評(píng)判枯草桿菌蛋白酶嗜熱性改造單突變效果評(píng)判的6個(gè)有效參數(shù)，其中ΔBOND、ΔELECT、ΔANGLE、ΔVDW 為特征有效參數(shù)，并采用370 K_0.2 ns模擬結(jié)果來提取上述參數(shù)。

枯草桿菌蛋白酶嗜熱性改造單突變效果評(píng)判方法為:以標(biāo)準(zhǔn)組370 K_0.2 ns的模擬條件參數(shù)為準(zhǔn)，分別對(duì)需要判定嗜熱改造效果的突變體及野生型進(jìn)行分子動(dòng)力學(xué)模擬，利用模擬結(jié)果計(jì)算上述4個(gè)特征有效參數(shù)，如果ΔBOND、ΔVDW、ΔANGLE大于0，同時(shí)ΔELECT小于0，則可以判斷該枯草桿菌蛋白酶突變體嗜熱性提高，為有效突變，否則為無效突變。

4 檢驗(yàn)與討論

4.1 枯草桿菌蛋白酶其他突變點(diǎn)檢驗(yàn)

我們從Protherm數(shù)據(jù)庫(kù)中選定了2個(gè)實(shí)驗(yàn)上已經(jīng)獲取的枯草桿菌蛋白酶嗜熱突變體作為檢測(cè)樣本，它們分別是Q206C和N218S。野生型及2個(gè)突變體在ProTherm數(shù)據(jù)庫(kù)的記載見表7。表中2個(gè)突變體的Tm值都大于野生型，說明突變體的嗜熱性要強(qiáng)于野生型。另外我們還選取了一個(gè)穩(wěn)定性小于野生型的突變體P239G，實(shí)驗(yàn)測(cè)得此突變體在60℃時(shí)的半衰期比野生型減少了8.5 min。

表7 枯草桿菌蛋白酶及其嗜熱突變體在Protherm數(shù)據(jù)庫(kù)中的記錄Table7 Records of subtilisin BPN'and its mutants in protherm

根據(jù)3.3節(jié)提出的判定原則，我們對(duì)以上選取的兩個(gè)有效突變體、一個(gè)無效突變體以及野生型進(jìn)行分子動(dòng)力學(xué)模擬，野生型命名為Wild。模擬結(jié)果見表8所示。

表8 檢驗(yàn)樣本模擬的平均能量Table 8 Average energy data simulated-test samples

表8下半部分第一列變化趨勢(shì)為有效突變體Q206C及N218S的能量變化結(jié)果，第二列為無效突變體P239G的結(jié)果。從數(shù)據(jù)可以看出，Q206C和N218S這兩個(gè)突變體的特征有效參數(shù)的取值滿足有效突變條件，并且能量變化的數(shù)值也處于上文中所提到的范圍之間，可以判斷Q206C和N218S為有效突變;而P239G的能量變化趨勢(shì)正好與之相反，為無效突變。

4.2 討 論

本文以枯草桿菌蛋白酶BPN＇為研究對(duì)象，利用4組不同條件分子動(dòng)力學(xué)模擬結(jié)果，對(duì)其及其嗜熱突變體的能量變化趨勢(shì)進(jìn)行對(duì)比分析，計(jì)算并提取了6個(gè)有效參數(shù)，并結(jié)合相關(guān)性分析結(jié)果，確定了4個(gè)特征有效參數(shù)，建立了判定枯草桿菌蛋白酶BPN＇嗜熱性突變效果的理論方法。研究發(fā)現(xiàn)，不同的模擬溫度與時(shí)間對(duì)突變前后各個(gè)能量的變化趨勢(shì)并沒有影響。而且在較短的模擬時(shí)間條件下，即模擬沒有達(dá)到平衡的狀態(tài)下，有效參數(shù)及特征有效參數(shù)具有相同的變化趨勢(shì)。本文提出的對(duì)枯草桿菌蛋白酶嗜熱效果改造的判定方法非常簡(jiǎn)單有效，這將為以后篩選有效的突變體提供捷徑。

利用短時(shí)間(時(shí)間小于等于0.2 ns)分子動(dòng)力學(xué)模擬結(jié)果提取有效參數(shù)，在P239位點(diǎn)對(duì)應(yīng)的兩個(gè)突變體P239R和P239G上均出現(xiàn)了ΔVDW取值不同于正常ΔVDW取值范圍的離群值，具體結(jié)果見表2、4、5、及表 8。其中 P239R 為有效突變，對(duì)應(yīng)ΔVDW為正值;P239G為無效突變，對(duì)應(yīng)ΔVDW為負(fù)值，盡管ΔVDW的變化趨勢(shì)是對(duì)的，但P239R和P239G對(duì)應(yīng)的ΔVDW明顯超出ΔVDW的取值范圍。這可能與脯氨酸(Pro)特殊的性質(zhì)有關(guān)。由于脯氨酸是非極性氨基酸，而精氨酸(Arg)是極性帶正電的氨基酸，因而引入帶正電的Arg可以在增加蛋白的穩(wěn)定性的同時(shí)使得分子極性大大增強(qiáng)，導(dǎo)致P239R的范德華能較其他幾個(gè)能量呈現(xiàn)出更大的增長(zhǎng)，造成極端數(shù)據(jù)的出現(xiàn)。而甘氨酸(Gly)是非極性不帶電荷氨基酸，在P239G處為何也會(huì)出現(xiàn)極端數(shù)據(jù)，需要更多的實(shí)驗(yàn)來進(jìn)行統(tǒng)計(jì)與驗(yàn)證。

此外，本文選取的試驗(yàn)樣本僅為枯草桿菌蛋白酶，對(duì)其他蛋白酶是否具有普適性還需進(jìn)一步研究與探討。

References)

［1］ NASCIMENTO W C A，MARTINS M L L.Production and properties of an extracellular protease from thermophilic Bacillussp ［J］. Brazilian Journalof Microbiology，2004，35(1-2):91-96.

［2］ ARGOS P，ROSSMAN M.G，GRAU U.M，ZUBER H.Thermal stability and protein structure ［J］.Biochemistry，1979，18(25):5698-5703.

［3］ ARNOLD F H，GIVER L，GERSHENSON A，ZHAO H.Directed evolution ofmesophilic enzymesinto their thermophilic counterparts ［J］.Ann N Y Acad Sci，1999，870:400-403.

［4］ LEHMANN M，LOCH C，MIDDENDORF A，et al.The consensus conceptforthermostability engineering of proteins:further proof of concept［J］. Protein Engineering，2002，15(5):403-411.

［5］ SEELIGER D，GROOT B L.Proteinthermostability calculations using alchemical free energy simulations［J］.Biophysical Journal，2010，98(10):2309-2316.

［6］ 田健，王平.理性設(shè)計(jì)提高蛋白質(zhì)熱穩(wěn)定性的研究進(jìn)展［J］.生物技術(shù)進(jìn)展，2012，4(2):233-239.TIAN Jian，WANG Ping.Recent advances in the rational design to improve the proteinthermostability［J］.Current Biotechnology，2012，4(2):233-239.

［7］ GROMIHA M M，UEDAIRA H，AN J.ProTherm:Thermodynamic Database for Proteins and Mutants:developments in version 3.0［J］.Nucleic Acids Res.，2002，30(1):301-302.

［8］ GRIBENKO A V，PATEL M M，LIU J，et al.Rationalstabilization of enzymes by computational redesign of surface charge-charge interactions［J］.Proc.Natl.Acad.Sci.，2009，106(8):2601-2606.

［9］ JOO J C，PACK S P，KIM Y H，et al.Thermostabilization of bacillus circulans xylanase:Computational optimization of unstable residues based on thermal fluctuation analysis［J］.Journal of Biotechnology，2011，151(1):56-65.

［10］ KORKEGIAN A，BLACK M.E，BAKER D.Computational thermostabilization of an enzyme［J］.Science，2005，308(5723):857-860.

［11］ SARABOJI K，GROMIHA M.Average assignment method for predicting the stability of protein mutants［J］.Biopolymers，2006，82(1):80-92.

［12］ CAPRIOTTI E，F(xiàn)ARISELLI P.A neural-network-based method for predicting proteinstability changes upon single point mutations［J］.Bioinformatics，2004，20(Suppl 1):i63-68.

［13］畢汝昌，儲(chǔ)乃明.枯草桿菌蛋白酶與蛋白質(zhì)工程［J］.生物化學(xué)與生物物理進(jìn)展，1991，18(5):329-334.BI Ruchang，CHU Naiming.Subtilisin and protein engineering［J］.Progress in Biochemistry and Biophysics，1991，18(5):329-334.

［14］王凡強(qiáng)，馬美榮，王正祥，等.枯草桿菌蛋白酶基因工程的研究進(jìn)展［J］.生物工程進(jìn)展，2000，20(2):41-44.WANG Fanqiang，MA Meirong，WANG Zhengxiang，et al.The research progress of Subtilisin gene engineering［J］.Progress in Biotechnology，2000，20(2):41-44.

［15］ MCKENNEY J M，KOREN M J，KEREIAKES D J.Safety and efficacy of a monoclonal antibody to proprotein convertase subtilisin/kexin type 9 serine protease，SAR236553/REGN727，in patients with primary hypercholesterolemia receiving ongoing stable atorvastatin therapy［J］. Journalofthe American College of Cardiology，2012，59(25):2344-2353.

［16］ GUMBINER B，CHANDRASEKHAR U，JOH T.The effects of multiple dose administration of RN316(PF-04950615)，a humanized IgG2a monoclonal antibody bindingproprotein convertase subtilisin kexin Type 9，in hypercholesterolemic subjects［J］.Circulation，2012，126(suppl):A13524-A13524.

［17］ KIMBER I，BASKETTER D A.Categorisation of protein respiratory allergens:The case ofSubtilisin［J］.Regulatory Toxicology and Pharmacology，2014，68(3):488-492.

［18］ZHAO H M，ARNOLD F H.Directed evolution converts subtilisin E into a functional equivalent of thermitase［J］.Protein Engineering，1999，12(1):47-53.

［19］ PANTOLIANO M W，WHITLOW M，WOOD J F，et al.Large increases in general stability for subtilisin BPN＇through incrementalchangesin the free energy of unfolding［J］.Biochemistry，1989，28(18):7205-213.

［20］ S?TTLER A，KANKA S，MAURER K，et al.Thermostable variants of subtilisin selected by temperature-gradient gel electrophoresis［J］.ELECTROPHORESIS，1996，17(4):784-92.

［21］CUNNINGHAM B C，WELLS J A.Improvement in the alkaline stability of subtilisin using an efficient random mutagenesisand screening procedure［J］. Protein Engineering，1987，1(4):319-325.

［22］ TAKAGI H，MORINAGA Y，IKEMURA H，et al.The Role of Pro-239 in the Catalysis and Heat Stability of Subtilisin E［J］.Journal of Biochemistry，1989，105(6):953-956.

［23］ RADESTOCK S，GOHLKE H.Protein rigidity andthermophilic adaptation［J］.Proteins:Structure，F(xiàn)unction，and Bioinformatics，2011，79(4):1089-1108.

［24］CHAN C H，WILBANKS C C，MAKHATADZE G I，et al.Electrostatic contribution of surface charge residues to the stability of a thermophilic protein:Benchmarking experimental and predicted pKa values［J］.PloS One，2012，7(1):e30296.