一种磷酸甘油酸激酶与底物结合模式的筛选方法

摘要

一种磷酸甘油酸激酶与底物结合模式的筛选方法涉及一种计算筛选磷酸甘油酸激酶与底物结合模式的方法。本发明提供了一种磷酸甘油酸激酶与底物相互作用方式探究的方法。通过结合量子化学计算、分子模拟、分子动力学构象收集以及基于构象的结合自由能计算，实现准确筛选磷酸甘油酸激酶与底物的结合模式，在预测磷酸甘油酸激酶1(PGK1)与潜在的金盏花苷底物(CE)的结合模式上取得了良好的效果。

说明书

技术领域

本发明提供了一种磷酸甘油酸激酶与底物结合模式的筛选方法。通过结合量子化学计算、分子模拟、分子动力学构象收集以及基于构象的结合自由能计算，实现准确筛选磷酸甘油酸激酶与底物的结合模式。基于该结合模式，底物结合对磷酸甘油酸激酶的调控效果与实验预期一致，在预测磷酸甘油酸激酶1(PGK1)与潜在的金盏花苷底物(CE)的结合模式上取得了良好的效果。

背景技术

磷酸甘油酸激酶1(PGK1)(2.7.2.3)是糖酵解过程中第一个产生ATP的酶，它催化底物adenosine diphosphate(ADP)从1，3-二磷酸甘油酯获得磷酸基团从而生成adenosinetriphosphate(ATP)和3-磷酸甘油酯。PGK1的功能受到缺氧诱导因子1α(hypoxia-inducible factor-1α)的调控。除了在糖酵解过程中发挥重要作用，PGK1还与肿瘤的生物过程、血管再生、DNA复制和损伤修复紧密相关。研究PGK1与底物的底物相互作用，对探究底物，特别是有药用价值的底物对上述相关生物学过程的作用具有重要的意义。

酶与底物的相互作用主要包括以下几个方面：1，底物对酶活性位点的识别和结合；2底物与酶结合后产生的极性和非极性相互作用；对于变构调节的酶，酶与底物的相互作用还体现在酶结合底物后产生的构象变化。传统的分子模拟手段，极大的简化了底物的构型，使用隐式溶剂模型导致溶剂化作用的显著低估，使用有限制的柔性范围忽略蛋白和底物构型取向对结合模式产生的影响，使分子模拟手段仅用于提供蛋白与底物配体的可能结合模式，而不能较为准确的预测具体的结合模式。分子动力学模拟以及结合自由能的计算则能够很好的补充分子模拟手段的不足。由于分子动力学手段在模拟过程中尽可能多地进行蛋白和底物相互作用构象的收集，通过稳定平衡态结合自由能的计算来弥补分子模拟手段在预测效果上的不足，更加准确的筛选和鉴定出符合理论要求的蛋白-底物结合模式。另外，基于分子动力学的模拟过程能够提供更为充分的蛋白和底物相互作用的细节，包括统计的热力学平均、极性和非极性相互作用等。

发明内容

本发明针对单一分子模拟研究手段对PGK1与底物结合模式预测的不足，提出了一种结合动力学模拟的新方法，能够实现更加准确的酶-底物结合模式筛选。

本发明采用如下技术方案：利用量子化学计算手段优化被预测底物的几何构型，获得构型依赖的电子态信息和自建力场参数；将电子态信息用于分子模拟，经多次实现，收集可供分析的蛋白结合底物的结合模型；使用自建力场参数进行各结合模型的分子动力学平衡优化，收集各结合模型的稳态的平衡构型，并在此基础上计算各结合模型的结合自由能；根据结合自由能的取值，选出能量最低的结合模型结构。技术路线如附图1所示。

一种磷酸甘油酸激酶与底物结合模式的筛选方法，通过量子化学计算、分子模拟、分子动力学计算以及结合自由能计算实现结合模式的筛选，结合方式表述为底物结合酶的体系中体系原子的原子类型、原子坐标信息、原子拓扑结构信息和原子电荷信息。

量子化学计算：

1)被预测底物的结构信息准备：将通过数据库查询获得的底物原子种类和原子坐标信息通过量子化学软件gauss 09D01在B3LYP密度泛函和6-31+g*的基组条件下进行能量优化计算，得到底物分子的每个原子经优化过后的坐标信息以及当前坐标下原子的电荷信息；

2)酶PGK1的结构信息准备：通过数据库查询获得的包含PGK1酶氨基酸序列信息以及氨基酸中各原子的坐标信息的初始结构文件使用Amber 16计算软件中的tleap程序应用Amber力场参数ff14SB计算PGK1所包含的氨基酸中各原子的电荷信息；

分子模拟，底物以不同结合模式结合入PGK1的模拟体系建立：

1)分子模拟计算底物结合PGK1的方式：使用量子化学计算或权利要求2中步骤1)计算得到的底物分子的原子坐标信息、当前坐标下原子的电荷信息，步骤2)中通过数据库查询得到的PGK1氨基酸序列和原子坐标信息及经Amber中tleap程序计算得到的PGK1所包含的氨基酸中各原子的电荷信息及作为输入文件，应用分子模拟开源程序antodock vina计算底物在PGK1经文献确认的结合位点和结合区域内的结合方式，结合方式包含底物的原子类型、原子坐标；

2)被预测底物不同结合模式的原子电荷计算：使用步骤1)中计算得到的底物在PGK1中的结合模式(包含底物的原子类型和原子坐标)文件作为输入参数，通过量子化学软件gauss 09D01在B3LYP密度泛函和6-31+g*的基组条件下计算得到当前原子类型所对应坐标中的原子电荷信息；

3)创建被预测底物在分子动力学模拟中所需要的力场参数：使用步骤2)中计算得到的底物原子类型所对应坐标中的原子电荷信息，通过Amber 16中的antechamber程序拟合计算RESP信息，并以RESP信息为输入参数文件，通过Amber 16中的parmchk2程序计算用于分子动力学模拟所需要的底物在PGK1中的结合模式的原子类型信息、原子半径信息、原子与原子之间的键长参数、键角参数、二面角参数以及1-4长程相互作用参数。

4)底物以不同结合模式结合入PGK1的模拟体系建立：根据量子化学计算或权利要求2中步骤2)中通过数据库查询得到的PGK1氨基酸序列和原子坐标信息以及权利要求3步骤1)计算得到的底物在PGK1中的结合模式(包含底物的原子类型和原子坐标)文件为输入文件，通过shell脚本进行合并，形成底物以不同方式结合PGK1的结合模式的原子信息和原子坐标信息文件,通过Amber中的tleap程序，应用Amber力场参数ff14SB和步骤3)中计算得到的底物在PGK1中的结合模式的原子类型信息、原子半径信息、原子与原子之间的键长参数、键角参数、二面角参数以及1-4长程相互作用参数，

计算生成底物以不同结合模式结合入PGK1的模拟体系，建立用于分子动力学模拟计算的底物与PGK1不同方式的结合模式所需要的原子类型、原子坐标、拓扑结构以及电荷信息的参数文件。

分子动力学计算过程为：

用于模拟的体系包括：

A，对照组：模拟参数文件来源为权利要求2中步骤2)产生的、没有结合底物的PGK1酶体系

B，实验组:模拟参数文件来源于权利要求3中步骤4)中产生的底物以不同结合模式结合入PGK1的模拟体系；

使用分子模拟或权利要求3中得到的底物与PGK1不同方式的结合模式所需要的原子类型、原子坐标、拓扑结构以及电荷信息的参数文件作为输入文件，通过Amber中的Sander.MPI程序进行分子动力学计算模拟，收集底物与PGK1不同方式的结合模式在2纳秒模拟时间内的动力学构象信息集，构象信息表述为底物与酶结合后的体系中所有原子的类型信息、原子在计算设定的时间点(间隔0.04皮秒)的空间坐标信息和特定时间点的原子电荷信息。

结合自由能计算过程：

使用分子动力学计算或权利要求4中得到的底物与PGK1不同方式的结合模式在2纳秒模拟时间内的动力学构象信息集(包含底物与酶结合后的体系中所有原子的类型信息、原子在特定时间点的空间坐标信息和特定时间点的原子电荷信息)作为输入参数，通过Amber中的MMPBSA.py.MPI计算底物与PGK1不同方式的结合模式在特定模拟时间的模拟过程中的平均结合自由能。

模式的筛选：

使用权利要求5中计算得到的底物与PGK1不同方式的结合模式在特定模拟时间的模拟过程中的平均结合自由能为筛选数据来源，选择数据中值为负且绝对值最大的值所对应的底物结合PGK1的结合方式，结合方式表述为权利要求3中产生的相应的体系所有原子的原子类型、原子坐标信息、原子拓扑结构信息和原子电荷信息。

本发明的有益效果：以量子化学理论为基础建立关于被预测底物的电子态结构信息和自建立场，提升了分子模拟过程的准确性，使之给出的可能结合模型更加贴近实际；使用分子动力学模拟对各结合模型进行平衡态的优化更多的关注了酶类等生物大分子由于构象因素所产生的变化；基于结合自由能进行筛选判定是理论计算的常用指标，是建立在稳态平衡条件下的统计平均，更符合蛋白与底物结合的实际热力学过程，因而能够实现更为准确的预测筛选。

附图说明

图1一种磷酸甘油酸激酶与底物结合模式的筛选方法的技术路线示意图；

图2实施例1中被预测底物CE的初始结构示意图；

图3实施例1中所有的CE结合PGK1的模式(Alt_CE_M1～Alt_CE_M8)与对照体系(Apo)经2ns动力学模拟的均方根偏差(RMSD)情况；

图4实施例1中所有CE结合PGK1的模式(M1～M8)通过结合自由能计算得到的结合自由能取值情况(单位kcal/mol)；

图5实施例1中确定的CE结合PGK1的模式(Alt_M7)的可视化表述；

图6实施例2中CE结合PGK1的确定模式(Alt_M7)与对照体系(Apo)经20ns动力学模拟的均方根偏差(RMSD)情况；

图7实施例2中CE结合PGK1的确定模式(Alt)与对照体系(Apo)波动性对比结果；

具体实施方式

实施例用于进一步说明本发明，但本发明不限于实施例。

实施例1(PGK1与CE结合模式的筛选确定)：

1，被预测底物CE的结构信息准备

被预测底物CE的结构信息来源于开源数据库(https://pubchem.ncbi.nlm.nih.gov)，查询ID：176079，其二维结构如附图2所示。将该结构所包含的原子类型及其对应的坐标信息导入量子化学计算软件gauss 09D01(M.J.Frisch,etal.,Gaussian 09,revision D.1.2010:Gaussian,Inc.,Wallingford,CT.)，选择密度泛函B3LYP方法在6-31+g*的基组条件下(基组参数来自https://www.basissetexchange.org/，查询关键字:6-31+G*)进行计算量子化学的优化，优化得到能量收敛的计算结果，该计算结果包含被预测底物CE经过量子化学计算优化之后每个原子所具有的新的坐标信息以及当前坐标下每个原子的电荷信息；

2，酶PGK1的结构信息准备

PGK1的初始结构信息来源于蛋白质晶体结构数据库(http://www.rcsb.org/pdb/home/home.do)，查询ID：4O33。初始结构信息包含PGK1的氨基酸序列信息以及氨基酸中各原子的坐标信息。使用专用于蛋白质分子计算模拟的Amber力场参数ff14SB(Maier,J.A.,et al.,ff14SB:Improving the Accuracy of Protein Side Chain and BackboneParameters from ff99SB.Journal of Chemical Theory and Computation,2015.11(8):p.3696-3713.)通过分子动力学计算软件Amber 16(D.A.Case,et al.,AMBER 16.2017:University of California,San Francisco.)中的tleap程序及ambpdb程序计算PGK1所包含的氨基酸中各原子的电荷信息，并与PGK1的氨基酸序列信息及氨基酸中各原子的坐标信息一同用于后续的计算和模拟。使用tleap程序和ambpdb程序计算PGK1所包含的氨基酸中各原子的电荷信息所使用的命令如下：

tleap–s–f leaprc.ff14SB

loadpdb system＝loadpdb PGK1.pdb

saveamberparm system receptor.top receptor.crd

ambpdb–p receptor.top–c receptor.top–mol2>receptor.mol2

其中，PGK1.pdb为包含PGK1的氨基酸序列信息以及氨基酸中各原子的坐标信息的初始结构文件；receptor.mol2为包含PGK1氨基酸中各原子的电荷信息的文件。receptor.top为未结合配体的Apo空白对照体系原子类型信息、原子坐标信息和原子电荷信息的文件，receptor.crd为未结合配体的Apo空白对照体系所包含的原子类型信息及拓扑结构信息文件。

3，分子模拟计算CE结合PGK1的方式

文献报道(Tian,Y.,et al.,Calenduloside E Analogues ProtectingH9c2Cardiomyocytes Against H2O2-Induced Apoptosis:Design,Synthesis andBiological Evaluation.Frontiers in Pharmacology,2017.8:p.862)认为，CE结合PGK1之后会对PGK1催化底物发生磷酸化产生抑制作用，而PGK1的抑制性底物结合与PGK1上的388号赖氨酸(Lys388)有关。将PGK1蛋白的重心设置为坐标原点，Lys388的重心坐标为x＝7.885,y＝34.077,z＝22.892，以其为中心的，长宽高各为

4，被预测底物CE不同结合模式的原子电荷计算

步骤3获得的8组包含CE原子类型和原子坐标信息的计算结果CE_${i}.mol(i∈[1，8])导入量子化学计算软件gauss 09D01，使用与步骤1相同的条件进行计算量子化学计算，得到的计算结果是当前原子类型和原子坐标所对应的原子电荷信息，保存8个包含CE原子类型、原子坐标信息以及坐标信息对应的原子电荷信息的文件CE_${i}.log(i∈[1，8])。

5，创建被预测底物CE在分子动力学模拟中所需要的力场参数

被预测底物CE相对于氨基酸来说，在用于蛋白质模拟的分子力学力场参数中缺失模拟所需的必要信息，需要根据CE当前坐标信息下的电荷信息来计算拟合RESP信息，并通过RESP信息以及CE坐标信息来建立分子动力学模拟所需要的力场参数。将步骤4中保存的8组包含CE原子类型、原子坐标信息以及坐标信息对应的原子电荷信息的文件CE_${i}.mol(i∈[1，8])通过Amber 16中的antechamber程序拟合计算RESP信息，使用的命令为：

antechamber–i CE_${i}.log fi gout–o Alt_CE_${i}.prepin–fo prepi–cresp–nc 0–rn CEM–s 2–pf y

Alt_CE_${i}.prepin(i∈[1，8])为RESP拟合计算的输出结果，包含CE当前坐标下的RESP电荷参数。RESP计算的输出结果继续通过Amber 16中的parmchk2程序计算8组CE原子类型、原子坐标以及RESP信息用于分子动力学模拟计算所需要的原子类型信息、原子半径信息、原子与原子之间的键长参数、键角参数、二面角参数以及1-4长程相互作用参数，使用的命令为：

parmchk2–i Alt_CE_${i}.prepin–f prepi–o Alt_CE_${i}.frcmod–a Y

Alt_CE_${i}.frcmod(i∈[1，8])为parmchk2程序的计算输出结果，包含8组CE的原子类型信息、原子半径信息，原子与原子之间的键长参数、键角参数、二面角参数以及1-4长程相互作用参数信息。

6，CE以不同结合模式结合入PGK1的模拟体系建立

通过shell命令语句步骤3中计算得到的8组包含CE的原子类型和原子坐标信息的计算结果添加入步骤2中得到的PGK1的氨基酸序列信息以及氨基酸中各原子的坐标信息文件中，形成CE以不同方式结合PGK1的结合模式的原子信息和原子坐标信息文件，所使用的shell命令语句为：

其中，Alt_CE_M${i}.pdb(i∈[1，8])为8个CE以不同方式结合PGK1的结合模式包含原子信息和原子坐标信息的文件。将这些信息通过Amber 16的tleap程序计算，使用步骤5中产生的用于补充8个CE以不同方式结合PGK1的结合模式中关于CE的原子类型信息、原子半径信息，原子与原子之间的键长参数、键角参数、二面角参数以及1-4长程相互作用参数，生成用于分子动力学模拟计算的8个CE以不同方式结合PGK1的结合模式所需要的原子类型、原子坐标、拓扑结构以及电荷信息的参数文件，tleap程序的执行命令语句为：

source/apps/bio/amber16/dat/leap/cmd/oldff/leaprc.ff14SB

source leaprc.gaff

loadamberparams frcmod.ionsjc_tip3p

loadamberprep Alt_CE_${i}.prepin

loadamberparams Alt_CE_${i}.frcmod

ligand＝loadpdb ligand_Alt_CE_${i}.pdb

complex＝loadpdb deWAT_Alt_CE_${i}.pdb

system＝loadpdb Alt_CE_M${i}.pdb

set default PBRadii mbondi2

saveamberparm ligand ligand_Alt_CE_M${i}.top ligand_Alt_CE_M${i}.crd

saveamberparm complex deWAT_Alt_CE_M${i}.top deWAT_Alt_CE_M${i}.crd

solvateoct system TIP3PBOX 12.0

addions system Cl-0

saveamberparm system Alt_CE_M${i}.top Alt_CE_M${i}.crd

其中，Alt_CE_M${i}.crd(i∈[1，8])为8个CE以不同方式结合PGK1的结合模式所包含原子类型信息、原子坐标信息和原子电荷信息的文件，Alt_CE_M${i}.top(i∈[1，8])为8个CE以不同方式结合PGK1的结合模式所包含的原子类型信息及拓扑结构信息文件。

7，分子动力学模拟

用于模拟的体系包括：

C，对照组(Apo)：模拟参数文件来源为步骤2产生的、没有结合底物CE的PGK1酶体系

D，实验组(Alt_CE_M${i}(i∈[1，8])):模拟参数文件来源于步骤6中产生的CE以不同结合模式结合入PGK1的模拟体系

模拟计算过程：

以步骤2和步骤6中产生的包含体系原子类型信息、原子坐标信息和原子电荷信息的.top文件和所包含体系原子类型信息及拓扑结构信息的.crd文件为结构输入文件，通过Amber 16的一般动力学程序Sander.MPI并行程序进行下述计算：

1),能量最小化过程：

a),使用最陡下降和共轭梯度下降1：1组合的方式对溶剂环境和PGK1的氢原子进行真空条件下的能量最小化计算优化，计算所使用的约束条件文件min1.in包含的命令具体为：

Initial minimization of solvent and ions

&cntrl

imin＝1,

maxcyc＝5000,

ncyc＝2500,

ntb＝1,

ntr＝1,

cut＝10

/

Hold the protein fixed

100.0

RES 1 418

END

END

执行Sander.MPI并行程序计算的命令行为：

Sander.MPI–O–i min1.in–o min1_Alt_CE_M${i}.out–c Alt_CE_M${i}.crd\

-p Alt_CE_M${i}.top–r min1_Alt_CE_M${i}.rst–x Alt_CE_M${i}.crd

Sander.MPI–O–i min1.in–o min1_receptor.out–c receptor.crd\

-p receptor.top–r min1_receptor.rst–x receptor.crd

b),使用a)中产生的、新的包含体系原子类型信息及拓扑结构信息的.rst文件和原始的包含体系原子类型信息、原子坐标信息和原子电荷信息的.top文件继续进行对氢原子进行优化，计算所使用的约束计算的脚本文件min2.in所包含的命令信息：

Initial minimization of H system

&cntrl

imin＝1,

maxcyc＝5000,

ncyc＝2500,

ntb＝1,

ntr＝1,

cut＝10,

restraintmask＝':1-418&！@H＝',

restraint_wt＝5.0,

/

执行Sander.MPI并行程序计算的命令行为：

Sander.MPI–O–i min2.in–o min2_Alt_CE_M${i}.out–c min1_Alt_CE_M${i}.rst\

-p Alt_CE_M${i}.top–r min2_Alt_CE_M${i}.rst–x min1_Alt_CE_M${i}.rst

Sander.MPI–O–i min2.in–o min2_receptor.out–c min1_receptor.rst\

-p receptor.top–r min2_receptor.rst–x min1_receptor.rst

c),使用b)中产生的、新的包含体系原子类型信息及拓扑结构信息的.rst文件和原始的包含体系原子类型信息、原子坐标信息和原子电荷信息的.top文件继续进行对侧链进行优化，计算所使用的约束计算的脚本文件min3.in所包含的命令信息：

Initial minimization of main chain

&cntrl

imin＝1,

maxcyc＝5000,

ncyc＝2500,

ntb＝1,

ntr＝1,

cut＝10,

restraintmask＝':1-418@C,CA,N',

restraint_wt＝1.0,

/

执行Sander.MPI并行程序计算的命令行为：

Sander.MPI–O–i min3.in–o min3_Alt_CE_M${i}.out–c min2_Alt_CE_M${i}.rst\

-p Alt_CE_M${i}.top–r min3_Alt_CE_M${i}.rst–x min2_Alt_CE_M${i}.rst

Sander.MPI–O–i min3.in–o min3_receptor.out–c min2_receptor.rst\

-p receptor.top–r min3_receptor.rst–x min2_receptor.rst

d),使用c)中产生的、新的包含体系原子类型信息及拓扑结构信息的.rst文件和原始的包含体系原子类型信息、原子坐标信息和原子电荷信息的.top文件继续进行对整个体系的优化，计算所使用的约束计算的脚本文件min4.in所包含的命令信息：

System mimimazation

&cntrl

imin＝1,

maxcyc＝10000,

ncyc＝5000,

ntb＝1,

ntr＝0,

cut＝10,

/

执行Sander.MPI并行程序计算的命令行为：

Sander.MPI–O–i min4.in–o min4_Alt_CE_M${i}.out–c min3_Alt_CE_M${i}.rst\

-p Alt_CE_M${i}.top–r min4_Alt_CE_M${i}.rst

Sander.MPI–O–i min4.in–o min4_receptor.out–c min3_receptor.rst\

-p receptor.top–r min4_receptor.rst

2),体系的升温过程：由于最小化过程是在真空条件下进行的，而模拟的最终目的是在室温/生理条件下还原酶与底物的结合作用模式，因此需要升温的模拟过程a),使用1)能量最小化优化的最后一步d)中产生的、新的包含体系原子类型信息及拓扑结构信息的.rst文件和原始的包含体系原子类型信息、原子坐标信息和原子电荷信息的.top文件进行等容(NVT)缓慢升温优化，将模拟计算的温度升至300K，计算所使用的约束计算的脚本文件equil1.in所包含的命令信息：

Initial equibilization with protein fixed

&cntrl

imin＝0,

irest＝0,

ntx＝1,

ntb＝1,

cut＝10,

ntr＝1,

ntc＝2,

ntf＝2,

tempi＝0.0,

temp0＝300.0,

ntt＝3,

ig＝-1,

gamma_ln＝1.0,

nstlim＝100000,

dt＝0.002,

ntpr＝1000,

ntwx＝1000,

ntwr＝10000,

ioutfm＝1,

/

Keep protein fixed with weak restrains

10.0

RES 1 418

END

END

执行Sander.MPI并行程序计算的命令行为：

Sander.MPI–O–i equil1.in–o equ1_Alt_CE_M${i}.out–c min4_Alt_CE_M${i}.rst\

-p Alt_CE_M${i}.top–r equ1_Alt_CE_M${i}.rst–x equ1_Alt_CE_M${i}.mdcrd\

-r min4_Alt_CE_M${i}.rst

Sander.MPI–O–i equil1.in–o equ1_receptor.out–c min4_receptor.rst\

-p receptor.top–r equ1_receptor.rst–x equ1_receptor.mdcrd\

-r min4_receptor.rst

b),使用a)中产生的、新的包含体系原子类型信息及拓扑结构信息的.rst文件和原始的包含体系原子类型信息、原子坐标信息和原子电荷信息的.top文件进行定温(300K)优化，以保证温度升高的过程不会对体系产生破坏,计算所使用的约束计算的脚本文件equil2.in所包含的命令信息：

Equibilization with temperature reached its level under constantvolume

&cntrl

imin＝0,

irest＝1,

ntx＝5,

ntb＝1,

cut＝10,

ntr＝1,

ntc＝2,

ntf＝2,

tempi＝300.0,

temp0＝300.0,

ntt＝3,

ig＝-1,

gamma_ln＝1.0,

nstlim＝100000,

dt＝0.002,

ntpr＝1000,

ntwx＝1000,

ntwr＝10000,

ioutfm＝1,

restraintmask＝':1-418@C,CA,N',

restraint_wt＝0.1

/

执行Sander.MPI并行程序计算的命令行为：

Sander.MPI–O–i equil2.in–o equ2_Alt_CE_M${i}.out–c equ1_Alt_CE_M${i}.rst\

-p Alt_CE_M${i}.top–r equ2_Alt_CE_M${i}.rst–x equ2_Alt_CE_M${i}.mdcrd\

-r equ1_Alt_CE_M${i}.rst

Sander.MPI–O–i equil2.in–o equ2_receptor.out–c equ1_receptor.rst\

-p receptor.top–r equ2_receptor.rst–x equ2_receptor.mdcrd\

-r equ1_receptor.rst

3),稳定的平衡态模拟：

使用2)体系升温优化的最后一步b)中产生的、新的包含体系原子类型信息及拓扑结构信息的.rst文件和原始的包含体系原子类型信息、原子坐标信息和原子电荷信息的.top文件进行2纳秒(nanoseconds,ns)的300K条件下恒温恒压的一般动力学模拟过程，计算所使用的约束计算的脚本文件prod.in所包含的命令信息：

system equibilization restart from last one

&cntrl

imin＝0,

irest＝1,

ntx＝5,

ntb＝2,pres0＝1.0,ntp＝1,

taup＝2.0,

cut＝10,

ntr＝0,

ntc＝2,ntf＝2,

tempi＝300.0,

temp0＝300.0,

ntt＝3,

ig＝-1,

gamma_ln＝1.0,

nstlim＝1000000,

dt＝0.002,

ntpr＝400,

ntwx＝400,

ntwr＝10000,

ioutfm＝1,

/

执行Sander.MPI并行程序计算的命令行为：

Sander.MPI–O–i prod.in–o prod_Alt_CE_M${i}.out–c equ2_Alt_CE_M${i}.rst\

-p Alt_CE_M${i}.top–r prod_Alt_CE_M${i}.rst–x prod_Alt_CE_M${i}.mdcrd

Sander.MPI–O–i prod.in–o prod_receptor.out–c equ2_receptor.rst\

-p receptor.top–r prod_receptor..rst–x prod_receptor.mdcrd

计算得到的prod_Alt_CE_M${i}.mdcrd(i∈[1，8])为8个CE以不同方式结合PGK1的结合模式在2ns模拟过程中所有被采集时间点的动力学构象信息集，prod_receptor.mdcrd为Apo空白对照组在2ns模拟过程中所有被采集时间点的动力学构象信息集。构象信息包含体系的原子类型信息、原子在特定时间点的空间坐标信息和特定时间点的原子电荷信息。衡量2ns模拟过程是否达到稳定平衡的标准是计算被模拟体系中α碳原子在2ns模拟过程中位置偏移情况的统计平均值是否随着模拟时间的延长而趋于稳定(在一个较小的范围内变化)，通过Amber 16的cpptraj程序来计算，计算所使用的约束计算的脚本文件cpptraj.in所包含的命令信息:

parm../../Alt_CE_M${i}.top

trajin../../prod_Alt_CE_${i}.mdcrd

reference../../Alt_CE_M${i}.crd

rms reference mass out rms_Alt_CE_M${i}.rms:1-418@C,CA,N

parm../../receptor.top

trajin../../prod_receptor.mdcrd

reference../../receptor.crd

rms reference mass out rms_receptor.rms:1-417@C,CA,N

得到的结果rms_Alt_CE_M${i}.rms(i∈[1，8])为8个CE以不同方式结合PGK1的结合模式α碳原子在2ns模拟过程中位置偏移情况的均方根偏差(RMSD)随模拟时间的变化值，rms_receptor.rms为Apo空白对照体系α碳原子在2ns模拟过程中位置偏移情况的均方根偏差(RMSD)随模拟时间的变化值，数值单位为埃(Angstrom)，绘图显示，结果如图3.所有的CE结合模式模型经过2ns的分子动力学优化，均获得相对稳定的热力学状态。

8，不同结合模式的结合自由能的计算

将步骤7中获得的8个CE以不同方式结合PGK1的结合模式以在2ns模拟过程中所有被采集时间点的动力学构象信息集prod_Alt_CE_M${i}.mdcrd(i∈[1，8])、步骤6中获得的8个CE以不同方式结合PGK1的结合模式包含体系原子类型信息、原子坐标信息和原子电荷信息的Alt_CE_M${i}.top(i∈[1，8])文件作为输入参数文件，通过Amber16中的MMPBSA.py.MPI程序计算8个CE以不同方式结合PGK1的结合模式体系在2ns模拟过程中的平均结合自由能，计算所使用的约束计算的脚本文件mmpbsa.in所包含的命令信息：

执行MMPBSA.py.MPI并行程序计算的命令行为：

MMPBSA.py.MPI–O–i mmpbsa.in-o MMPBSA_M${i}.dat-sp Alt_CE_M${i}.top

-cp deWAT_Alt_CE_M${i}.top-rp receptor.top-lp ligand_Alt_CE_M${i}.top

-y prod_Alt_CE_M${i}.mdcrd

计算得到的MMPBSA_M${i}.dat(i∈[1，8])8个CE以不同方式结合PGK1的结合模式在2ns分子动力学模拟过程中平均结合自由能的计算结果。将结合模型用M1～M8表示，将平均结合自由能的结果列出，如图4所示，平均结合自由能的计算结果单位为kcal/mol。

9，CE结合PGK1的模式筛选

由步骤8中计算得到的8个CE以不同方式结合PGK1的结合模式在2ns分子动力学模拟过程中平均结合自由能的计算结果MMPBSA_M${i}.dat(i∈[1，8])显示如图4可知，M7结合模式在2ns的分子动力学模拟中所具有的平均结合自由能为负值，且绝对值最大，因此该结合模式被确定为CE结合PGK1的结合方式。该结合模型可以表述为步骤6产生的，Alt_CE_M7.top和Alt_CE_M7.crd所确定的体系原子的原子类型、原子坐标信息、原子拓扑结构信息和原子电荷信息。宏观描述为：位于PGK1上的357号天门冬氨酸(ASP 357)，360号缬氨酸(VAL 360)，388号赖氨酸(LYS 388)，364号丝氨酸(SER 364)，363号苏氨酸(THR 363)，385号苏氨酸(THR 385)，366号甘氨酸(GLY366)，361号赖氨酸(LYS 361)，356号甲硫氨酸(MET356)，362号丙氨酸(ALA 362)，359号缬氨酸(VAL 359)，389号缬氨酸(VAL 389)，358号谷氨酸(GLU 358)及368号异亮氨酸(ILE 368)参与形成PGK1与CE相互作用的网络，底物CE能够与其中的Lys388与Asp357形成氢键相互作用。CE在PGK1中结合的空间取向以及上述氨基酸在PGK1上的分布情况由图5所示，其中，Lys388与Asp357与CE的相互作用显示为图5左侧部分，其他重要氨基酸的分布显示为图5右侧部分。

实施例2(PGK1与CE结合模式的验证)：

1，已确定的结合模式体系长程动力学模拟

实施例1步骤7计算得到的，步骤9确定的第7号体系经过2ns动力学模拟计算新的包含体系原子类型信息及拓扑结构信息的prod_Alt_CE_M7.rst文件和原始的包含体系原子类型信息、原子坐标信息和原子电荷信息的Alt_CE_M7.top文件，步骤7计算得到的Apo空白对照体系经过2ns动力学模拟计算新的包含体系原子类型信息及拓扑结构信息的prod_receptor.rst文件以及步骤2得到的包含体系原子类型信息、原子坐标信息和原子电荷信息的receptor.top文件为输入文件，通过Amber 16的一般动力学程序Sander.MPI并行程序，继续进行20ns的300K条件下恒温恒压的一般动力学模拟过程，计算所使用的约束计算的脚本文件prod_20ns.in所包含的命令信息

system equibilization restart from last one

&cntrl

imin＝0,

irest＝1,

ntx＝5,

ntb＝2,pres0＝1.0,ntp＝1,

taup＝2.0,

cut＝10,

ntr＝0,

ntc＝2,ntf＝2,

tempi＝300.0,

temp0＝300.0,

ntt＝3,

ig＝-1,

gamma_ln＝1.0,

nstlim＝10000000,

dt＝0.002,

ntpr＝2000,

ntwx＝2000,

ntwr＝2000,

ioutfm＝1,

/

执行Sander.MPI并行程序计算的命令行为：

Sander.MPI–O–i prod.in–o prod_Alt_CE_M7_20ns.out–c prod_Alt_CE_M7.rst\

-p Alt_CE_M7.top–r prod_Alt_CE_M7_20ns.rst–x prod_Alt_CE_M7_20ns.mdcrd\

Sander.MPI–O–i prod.in–o prod_receptor_20ns.out–c prod_receptor.rst\

-p receptor.top–r prod_receptor_20ns.rst–x prod_receptor_20ns.mdcrd

计算得到的prod_Alt_CE_M7_20ns.mdcrd为第7号体系经20ns模拟所有被采集时间点的动力学构象信息集，prod_receptor_20ns.mdcrd为空白对照体系经20ns模拟所有被采集时间点的动力学构象信息集。构象信息包含体系的原子类型信息、原子在特定时间点的空间坐标信息和特定时间点的原子电荷信息。通过Amber 16的cpptraj程序来计算体系中α碳原子在20ns模拟过程中位置偏移情况的均方根偏差RMSD是否随着模拟时间的延长而趋于稳定，计算所使用的约束计算的脚本文件cpptraj.in所包含的命令信息:

parm../../Alt_CE_M7.top

trajin../../prod_Alt_CE_M7_20ns.mdcrd

reference../../Alt_CE_M7.crd

rms reference mass out rms_Alt_CE_M7_20ns.rms:1-418@C,CA,N

parm../../receptor.top

trajin../../prod_receptor_20ns.mdcrd

reference../../receptor.crd

rms reference mass out rms_receptor_20ns.rms:1-418@C,CA,N

得到的结果rms_Alt_CE_M7_20ns.rms为第7号体系经20ns模拟α碳原子位置偏移情况的均方根偏差(RMSD)随模拟时间的变化值,rms_receptor_20ns.rms为Apo空白对照体系经20ns模拟α碳原子位置偏移情况的均方根偏差(RMSD)随模拟时间的变化值，数值单位为埃(Angstrom)，，结果一同绘图对比，如图6，可知，第七号结合模型体系在整个20ns动力学过程表现出了比较小的波动性，体系保持了相对稳定的热力学状态。

2，CE结合PGK1导致的PGK1构型变化

以步骤1中计算得到的第7号体系经20ns模拟所有被采集时间点的动力学构象信息集prod_Alt_CE_M7_20ns.mdcrd以及Apo空白对照组经20ns模拟所有被采集时间点的动力学构象信息集prod_receptor_20ns.mdcrd作为输入文件，应用Amber软件中的cpptraj程序，计算该体系结合CE之后PGK1中各氨基酸在20ns模拟过程中的波动性数值与没有结合CE的对照组Apo体系的计算结果进行对比，数值单位为埃(Angstrom)，计算过程使用的约束计算脚本文件fluct.cpptraj包含的命令为：

parm../Alt_CE_M7.top

trajin../prod_Alt_CE_M7_20ns.mdcrd

average avg_Alt_M7.pdb:1-418

parm../Alt_CE_M7.top

trajin../prod_Alt_CE_M7_20ns.mdcrd

parm avg_Alt_M7.pdb

reference avg_Alt_M7.pdb avg_Alt_M7.pdb

rms reference

atomicfluct out Fluct_ref_M7.dat byres

parm../receptor.top

trajin../prod_receptor_20ns.mdcrd

average avg_receptor.pdb:1-417

parm../receptor.top

trajin../prod_receptor_20ns.mdcrd

parm avg_receptor.pdb

reference avg_receptor.pdb avg_receptor.pdb

rms reference

atomicfluct out Fluct_ref_receptor.dat byres

计算得到的结果Fluct_ref_M7.dat为第7号体系PGK1中各氨基酸在20ns模拟过程中的波动性计算结果，单位为埃(Angstrom)，Fluct_ref_receptor.dat为Apo体系的计算结果，对比绘图，结果如图7所示。由结果可知，PGK1在结合CE后，体系中氨基酸的自由度与未结合任何底物的Apo体系相比均有所减小，显示出波动性被抑制的状态，预示着酶的功能由于结构变化的抑制而不能正常发挥作用，产生宏观现象中酶活性受到抑制的效果，这一结果与实验测定的CE的结合抑制PGK1的作用相一致。

由步骤1和步骤2中的结果可知，应用本发明的筛选方法确定的CE结合PGK1的结合模式，与实验结果具有高度的一致性，达到了准确地计算筛选的目的，效果良好。