SPONGE进阶教程：MetaDynamics的简单用例

前序课程1
前序课程2

场景简述
蛋白与配体对接后，采用通常的分子动力学模拟，通常会采样得到势阱附近的大量结构（如下图左L-P），即蛋白-配体结合构象的平均系综。如果想要探究蛋白-配体解离过程，需要克服一个解离能垒，去到L-P*，甚至更远的L+P解离状态。这时候静候模拟体系自然地运动过去是不现实的（对于实验结合强度为μmol级的药物（pIC50~6.0），对应势阱深度约为10kcal/mol，其平衡态和稳定态的热力学概率比约为400：1，更何况解离路径上会存在蛋白位阻和溶剂排斥，实际成功游离的概率很低，统计热力学上不具有代表性）。

对于稀有事件的采样，有许多增强采样算法，本案例简单引入常用的MetaDynamics算法（以下简称为MetaD），演示在SPONGE软件中如何在较短观测时间内，模拟蛋白-配体解离过程。
MetaD是一种基于集合变量（Collective Variables,以下简称CVs）的增强采样方法，上图的横轴反应坐标即是一个广义的抽象的CV概念，CV的左侧L-P对应蛋白-配体结合态（状态A），CV的右侧L+P对应蛋白-配体解离态（状态B），CV应该满足连续可导。对于描述蛋白-配体结合-解离过程，CV可以有许多不同的定义，常见的包括配体质心与蛋白质心的距离，配体与蛋白结合位点原子的距离，配体的溶剂化面积等等。

本案例定义了一个非常粗糙的CV作为演示，定义配体质心MC到蛋白中心附近的174号GLY残基Cα(atomid=2643)的距离为d。

开启MetaD功能首先需要在mdin_metad.txt文件中添加“cv_in_file = cv.txt”指令，指向CV配置文件cv.txt。

    The first line is the name of the md tast

#本文本文件通过命令行中加入"-mdin 本文本文件名"读入，如果该文件名为"mdin.txt"且在当前命令行的工作路径中，可不在命令行中显式输入

#本文本文件内的所有参数指令均可通过加"-"前缀在命令行里输入，例如 -cutoff 8.0 -skin 2.0



#分子模拟最基本相关设置

mode = NPT #可选Rerun, Minimization, NVE, NVT, NPT



dt = 1e-3 #模拟步长，单位ps

step_limit = 100000 #模拟总步数，通常需要长一点的模拟，本教程用了偏小的值

write_information_interval = 100 #保存中间信息的间隔

write_restart_file_interval = 100000  #保存重启坐标和速度的间隔



target_temperature = 300.0 #模拟的目标平衡温度，单位为开

thermostat = middle_langevin #热浴方法

barostat =  andersen_barostat #压浴设置

target_pressure = 1.0 #模拟的目标平衡压强，单位为bar


restrain #此处重新开启restrain是为了固定蛋白的大致坐标和基本朝向，这样配体的绝对坐标分布可以与蛋白的坐标绘制在同一个图像上，而不需要进行平移旋转变换。
{
    atom_id = restrain_ca.txt 
    weight = 0.1
}

    #输入

        #文件目录

        default_in_file_prefix = 1ae5/1ae5



        #坐标和速度

        coordinate_in_file = Pmin_coordinate.txt #使用NPT得到的坐标作为输入坐标

        velocity_in_file = Pmin_velocity.txt
	
cv_in_file = cv.txt #定义CV文件


    #输出

        mdout = mdout.txt   #记录能量轨迹的文件

        mdinfo = mdinfo.txt #记录参数相关的文件

        box = mdbox.txt       #记录盒子轨迹的文件

        crd = mdcrd.dat        #记录坐标轨迹的文件

        rst = restart              #记录重启坐标和速度的文件



#其他参数

end_pause = 1 #程序运行结束后是否需要输入任意键退出

device = 0 #使用的GPU设备

dont_check_input = 1 #忽略多余参数

让我们看看cv.txt定义了些什么：

com_of_ligand #定义配体质心，atom列为配体原子的atom_id
{
vatom_type = center_of_mass
atom = 3352 3353 3354 3355 3356 3357 3358 3359 3360 3361 3362 3363 3364 3365 3366 3367 3368 3369 3370 3371 3372 3373 3374 3375 3376 3377 3378 3379 3380 3381 3382 3383 3384 3385 3386 3387
}
cx #取出质心的x方向坐标，cy,cz分别是y方向、z方向的分量，方便后续绘图示意
{
CV_type = position_x
atom = com_of_ligand
}
cy
{
CV_type = position_y
atom = com_of_ligand
}
cz
{
CV_type = position_z
atom = com_of_ligand
}
dis # 定义配体质心与蛋白内部174号GLY残基Cα(atomid=2643)的距离为dis
{
CV_type = distance
atom = com_of_ligand 2643
}
print #在mdout.txt及屏幕输出添加感兴趣的字段
{
CV = cx cy cz dis
}
meta #开启MetaD
{
Ndim = 1 #CV维度为1维
CV = dis #采样的CV为两点相对距离dis
CV_minimal = 0 #采样CV的下限
CV_maximum = 50 #采样CV的上限，体系边长约为50A，由于周期性边界条件，内部最远距离应该是体对角线的一半sqrt(3)/2*50=43.3A, 上限高于此值已经足够。
CV_period = 0 #周期性CV需要定义，本案例CV没有周期性
CV_grid = 100 #CV等分为100个格点进行统计
welltemp_factor = 50 #MetaD的biasfactor，单位为kcal/mol
height = 0.6 #MetaD的biasheight，单位为kcal/mol
CV_sigma = 1.0 #CV的高斯峰宽
}

接下来运行带MetaD的模拟：

awk '{if($1=="CX"){print NR-2}}' 1ae5/1ae5_atom_type_name.txt > restrain_ca.txt
#将蛋白的Cα原子id生成一个列表，用于restrain
sed -i '/2643/d' restrain_ca.txt
#因为2643号原子会受到MetaD的bias效果，不加restrain更干净些
./SPONGE -mdin mdin_metad.txt

模拟结束之后，对结果进行统计：

awk '{if(NR>1){print $19,$20,$21,$14}}' mdout.txt >xyz_bias.txt
#汇总坐标信息和概率分布到文件xyz_bias.txt中，进行fes统计。

fes统计可以使用插件cyfes,跳转博客链接，有详细的cyfes原理、安装和使用说明。