Amber学习第五天：模拟含有非标准残基的溶剂化蛋白
前面所及的都是标准的残基，所以我们不用为非标准的残基创建的单元和提供我们自己的参数。这章我们将介绍一种方法对于非标准的残基。以含有铜离子的质体蓝素蛋白为例。需要作的是：
1)铜离子由四个氨基酸帮定His37, Cys84, His87 和 Met92。我们将修饰这些氨基酸并提供新键类型的参数
2)应该保留晶体本身的水分子。但是痒的位子应该被指定所以用Xleap来增加丢失的质子。所以在运行模拟前要进行最小化这些位置。
3)不寻常的pdb文件含有显性的质子层。这种情况我们去除不能识别的质子并让leap来增加他们。
4)质子化的蛋白有－9电荷，我们需要增加9个Na+抗性离子来中和这个系统。
stage1 --pdb文件的编辑：
我们将用PDB文件：1PLC - 1PLC.pdb.
这个文件的前面写道：waters #187 and #183 form a disordered pair and so both should not be present.
（REMARK 4 1PLC 83
REMARK 4 HOH 187 AND HOH 183 FORM A DISORDERED PAIR AND BOTH ARE NOT 1PLC 84
REMARK 4 PRESENT SIMULTANEOUSLY. 1PLC 85）
所以我们除去水187并保留183，而且也除去在leap中不能用的数据(如最后的CONECT内容)。正常的情况下我们应该增加TER cards在没有结合到链上的残基，但在水溶剂的情况下不用这样做，因为这个被定义在WAT unit中。注：如果用不同的溶剂，它是必须的增加TER。
由于pdb件不区分参与键或其他东西的半胱氨酸残基，需要编辑半胱氨酸残基。对于有规律的质子化残基的名字为CYS;对于去质子化或绑定到金属离子的名字为CYM；对于涉及到二硫键或其它键的为CYX。由于在文件中的半胱氨酸84涉及到金属离子
所以需要改变CYS到CYM。The same is true for histidine residues which can be protonated in the 'delta' position (HID), the epsilon position (HIE) or at both (HIP). Fortunately in plastocyanin this is fairly easy since there are only two histidine residues (37 and 87) both of which are bonded to the copper via the 'delta' nitrogen. Thus they must both be protonated on the epsilon nitrogen. We will thus change both histidine residue names (37 and 87) from HIS to HIE. Leap will then add the correct number of protons in the correct positions.修饰后的文件为：1PLC_mod.pdb
接下来去除非标准的氢原子。它是可能的用质子化来纠正非标准的氢原子，从NMR转到PDB，实验显示在NMR结构的氢位置是不能相信的。因此最好的方法是去除所有的质子(一般在13列或14)并且让leap在标准的位置增加它们。产生的文件为：1PLC_mod2.pdb
除去质子：grep -v '.............H' 1PLC_mod.pdb > 1PLC_mod2.pdb. 模拟需要的氢原子和水分子将在leap中添加
接下来做的是怎样处理铜离子。简单的处理一下
编辑pdb文件把CU改为CUA；
原始的pdb文件含有一些氨基残基可选择用的参数：如LYS30
最后在铜离子和初始晶体水分子间加一个TER。(作用是不让铜离子成为蛋白链的一部分而使蛋白不稳定)。
stage2--创建非标准的CUA单元
有三种方法：
(1):最简单是利用Antechamber来创建，(2)在xleap中编辑CU残基。(3)创建一个新的库文件非标准残基，特别对于辅酶和NAD，我们用第三个。
我们可以现导入1PLC_mod_final.pdb文件来来检测残基的识别。
$ tleap
>source leaprc.ff99
>1PLC=loadpdb 1PLC_mod_final.pdb
会返回 Creating new UNIT for residue: CUA sequence: 100
Created a new atom named: CU within residue: .R<CUA 100>
total atoms in file: 894
Leap added 922 missing atoms according to residue templates:
922 H / lone pairs
The file contained 1 atoms not in residue templates
在模板中没有CUA残基。
创建CUA单元
最快的方法是从1PLC_mod_final.pdb 文件中剪切出来并保存为一个pdb文件cua.pdb。复制一句话
载入文件：
$ tleap
>source leaprc.ff99
> CUA = loadpdb cua.pdb
> edit CUA
接下来是指定这个原子的类型和电荷
需要计算铜原子的电荷和其有关系的残基的电荷，由于铜原子能使静电荷发生变化。在AMBER中应用Restrained Electrostatic Potential Method (RESP)来做。在这里就不计算了(教程在AMBER website)，假设铜原子为＋１价且不影响周围的残基。
在Manipulation上选择select button并点击分子，会发生颜色的变化。让后去Edit menu Edit->Edit Selected Atoms：在TYPE中填入CU,在CHARGE中填入1.0000. 然后选择Table->Save and Quit，关闭编辑窗口。在命令行中键入：
>decs MET (看这个新的单元头原子和尾原子，被用来联系蛋白链)
但是当键入decs CUA时，会发生：
UNIT name: CUA
Head atom: null
Tail atom: null
因为铜原子不是蛋白的一部分。如果创建了蛋白上氨基的新单元，你需要用这个命令显示那个原子是头原子和那个原子是尾原子。
保存完整的库文件：
> saveoff CUA cua.lib
stage3 --载入蛋白和创建库文件
在xleap中看缺失参数：
为了让xleap识别新的残基在导入1PLC_mod_final.pdb文件时，我们需要现载入新的库文件：
$AMBERHOME/exe/xleap -s -f $AMBERHOME/dat/leap/cmd/leaprc.ff99
>loadoff cua.lib
> 1PLC = loadpdb 1PLC_mod_final.pdb
接下来确保于铜原子作用的所有键被定义。我们需要增加在铜原子于半胱氨酸（84）硫原子的键，铜原子于甲硫氨酸(92)硫原子的键和两个组氨酸(37 & 87)的三角洲氮于铜原子的键。
> bond 1PLC.37.ND1 1PLC.100.CU
> bond 1PLC.87.ND1 1PLC.100.CU
> bond 1PLC.84.SG 1PLC.100.CU
> bond 1PLC.92.SD 1PLC.100.CU
注：DESC命令可以列出的残基和原子名称，因此可以识别结合的原子，如desc 1PLC.CUA.
接下来加水（a truncated octahedral box of TIP3P water.）
> solvateoct 1PLC TIP3PBOX 12
中性化我们的系统：
> addions 1PLC Na+ 9 （check 1PLC）
现在保存 prmtop and inpcrd 文件，发现xleappppp不能识别CU
> saveamberparm 1PLC 1PLC.prmtop 1PLC.inpcrd
是因为在标准立场中不存在CU，通过检查1PLC会发现由于溶剂的作用使大部分失去作用，丢失的参数于新类型CU有关。
> check 1PLC
所以我们需要增加这些参数到AMBER立场中，在关闭xleap前因该保存1PLC库以便重复以上的步骤。
> saveoff 1PLC 1PLC.lib
stage4 --创建Prmtop and Inpcrd 文件
创建立场文件：plc.frcmod
# modifications to force field for poplar plastocyanin
MASS
CU 65.36
BOND
NB-CU 70.000 2.05000 #kludge by JRS
CU-S 70.000 2.10000 #kludge by JRS
CU-SH 70.000 2.90000 #for pcy
CT-SH 222.000 1.81000 #met(aa)
ANGLE
CU-NB-CV 50.000 126.700 #JRS estimate
CU-NB-CR 50.000 126.700 #JRS estimate
CU-NB-CP 50.000 126.700 #JRS estimate
CU-NB-CC 50.000 126.700 #JRS estimate
CU-SH-CT 50.000 120.000 #JRS estimate
CU-S -CT 50.000 120.000 #JRS estimate
CU-S -C2 50.000 120.000 #JRS estimate
CU-S -C3 50.000 120.000 #JRS estimate
NB-CU-NB 10.000 110.000 #dac estimate
NB-CU-SH 10.000 110.000 #dac estimate
NB-CU-S 10.000 110.000 #dac estimate
SH-CU-S 10.000 110.000 #dac estimate
CU-SH-CT 50.000 120.000 #JRS estimate
CT-CT-SH 50.000 114.700 #met(aa)
HC-CT-SH 35.000 109.500
H1-CT-SH 35.000 109.500
CT-SH-CT 62.000 98.900 #MET(OL)
DIHE
X -NB-CU-X 1 0.000 180.000 3.000
X -CU-SH-X 1 0.000 180.000 3.000
X -CU-S -X 1 0.000 180.000 3.000
X -CT-SH-X 3 1.000 0.000 3.000
NONBON
CU 2.20 0.200
$AMBERHOME/exe/xleap -s -f $AMBERHOME/dat/leap/cmd/leaprc.ff99
> loadamberparams plc.frcmod
> loadoff 1PLC.lib
现在可以创建Prmtop and Inpcrd 文件;
> saveamberparm 1PLC 1PLC.prmtop 1PLC.inpcrd
用 ambpdb 创建pdb文件:
$AMBERHOME/exe/ambpdb -p 1PLC.prmtop < 1PLC.inpcrd > 1PLC.inpcrd.pdb
可以在VMD中观察这个蛋白了。
We could now use these files to run plastocyanin simulations. You can try this yourself if you want. Just remember that this is in explicit solvent and a peridiodic box so you use periodic boundary conditions.
You will also need to minimise the system initially to remove bad contacts. I would then heat it over 20ps from 0 to 300K with constant volume periodic boundaries before moving to a long equilibration at 300 K
with constant pressure
