Constant pH simulation of BPTI Implicit(amber)

This tutorial describes how to run constant pH simulation by Amber, we will be using BPTI as an example. This tutorial will be using implicit solvent model.

1. Structure Preparation

1.1 PDB Conversion

copy 6PTI.pdb to current directory,

pdb4amber -i 6PTI.pdb -o 6pti_fixed.pdb -d --constantph

 

==================================================
Summary of pdb4amber for file 6PTI.pdb
==================================================----------Chains
The following (original) chains have been found:
A---------- Alternate Locations (Original Residues!)
The following residues had alternate locations:
ARG_39
ASP_50
The alternate coordinates have been discarded.
Only the first occurrence for each atom was kept.---------- Non-Standard Residues
The following non-standard residue names in the original PDB file
are not recognized by Amber and have been written to the separate
file 6pti_fixed_nonprot.pdb
PO4---------- Water
73 water molecules have been removed
and stored in the file 6pti_fixed_water.pdb---------- Constant pH Simulation
ASP --> AS4: 2
GLU --> GL4: 2
HIS --> HIP: 0---------- Cysteines in Disulfide Bonds (Renumbered Residues!)
CYS_5 - CYS_55: S-S distance = 2.062611 Ang.
CYS_14 - CYS_38: S-S distance = 2.045531 Ang.
CYS_30 - CYS_51: S-S distance = 1.998788 Ang.
The above CYS have been renamed to CYX in the new PDB file.
Disulfide bond CONECT cards have been added to the new PDB file.---------- Missing Heavy Atoms (Renumbered Residues!)
GLY_57 misses 3 heavy atom(s)

 

This will rename titratable residues according to the table below.

#Amino Acid Titratable Residue Name Reference pKa
Aspartate AS4 4.0
Glutamate GL4 4.4
Histidine HIP 6.6
Cysteine CYS(CYX if DS) 8.5
Lysine LYS 10.4
Tyrosine TYR 9.6

The residues titratable in BPTI are displayed in picture below:

BPTI_titratable residues

all titratable residues includes(from file 6pti.cpin):


RESNAME='System: Unknown','Residue: AS4 3','Residue: GL4 7','Residue: TYR 10',
'Residue: LYS 15','Residue: TYR 21','Residue: TYR 23','Residue: LYS 26',
'Residue: TYR 35','Residue: LYS 41','Residue: LYS 46','Residue: GL4 49',
'Residue: AS4 50'

 

edit 6pti_fixed.pdb,remove all HETATM,CONECT

1.2 Prepare the Topology File and Input  Coordinate File

edit tleap.in

source leaprc.constph
# Load the PDBs
6pti = loadPDB 6pti_fixed.pdb# Make the disulfide bonds
bond 6pti.5.SG 6pti.55.SG
bond 6pti.14.SG 6pti.38.SG
bond 6pti.30.SG 6pti.51.SG# Save topology files
saveAmberParm 6pti 6pti.parm7 6pti.rst7
# Quit
quit

 

tleap -f tleap.in

This will add hydrogens and four missing heavy atoms  CA,C,O,OXT on residue 57 GLY, total 898 atoms.

2. PH Input File Preparation

Create input files for pH simulation using cpinutil.py, this step will let the simulation engine know which residues to titrate during simulation.

cpinutil.py -p 6pti.parm7 -o 6pti.cpin

Let's take a look at cpin file

STATEINF(0)%FIRST_ATOM=41, STATEINF(0)%FIRST_CHARGE=0,
23 STATEINF(0)%FIRST_STATE=0, STATEINF(0)%NUM_ATOMS=16, STATEINF(0)%NUM_STATES=5,

FIRST_ATOM of first titratable residue is 41(start from 1), "0" of FIRST_CHARGE is index of charge array which is -0.4157. FIRST_STATE=0, 0 is index of STATENE array, which is 0 for deprotonated state, energy for state 5 is also 0.NUM_ATOMS is the total number of atoms in that residue(from 41 to 56).NUM_STATES=5, the total states of this residue is 5 including 1 deprotonated state and 4 protonated states(syn- and anti- for both oxygens).

3. System Minimization

min.mdin

Minimization to prepare for constant pH MD
&cntrl
  imin = 1, ! Turn on minimization
  igb = 2, ! Use the GB model CpHMD was parametrized for
  saltcon = 0.1, ! Use the salt conc. CpHMD was parametrized for
  maxcyc = 1000, ! Total number of minimization cycles
  icnstph = 1, ! Turn on constant pH for implicit
  ntcnstph = 100000, ! Never attempt to change prot. states
  cut = 1000.0, ! No cutoff
  ntr = 1, ! Turn on positional restraints
  restraint_wt = 10, ! 10 kcal/mol/A**2 restraint force constant
  restraintmask = '@CA,C,O,N' ! Restraints on the backbone atoms only
/

mpirun -np 12 pmemd.MPI -O -i min.mdin -p 6pti.parm7 -c 6pti.rst7 -o 6pti.min.mdout -r 6pti.min.rst7 -ref 6pti.rst7 -cpin 6pti.cpin

mdout_analyzer.py 6pti.min.mdout

2016-05-10 16:14:41 的屏幕截图

4. Heating the System

heat.mdin

Implicit solvent constant pH initial heating mdin
 &cntrl
  imin=0, irest=0, ntx=1,
  ntpr=500, ntwx=500, nstlim=1000000,
  dt=0.002, ntt=3, tempi=100,
  temp0=300, tautp=2.0, ig=-1,
  ntp=0, ntc=2, ntf=2, cut=30,
  ntb=0, igb=2, tol=0.000001,
  nrespa=1, saltcon=0.1, icnstph=1,
  solvph=7.5, ntcnstph=100000000,
  gamma_ln=5.0, ntwr=500, ioutfm=1,
  nmropt=1, ntr=1, restraint_wt=2.0,
  restraintmask='@CA,C,O,N',
 /
 &wt
  TYPE='TEMP0', ISTEP1=1, ISTEP2=500000,
  VALUE1=10.0, VALUE2=300.0,
 /
 &wt TYPE='END' /

mpirun -np 12 pmemd.MPI -O -i heat.mdin -c 6pti.min.rst7 -p 6pti.parm7 -ref 6pti.min.rst7 -cpin 6pti.cpin -o 6pti.heat.mdout -r 6pti.heat.rst7 -x 6pti.heat.nc

5. Equilibrating the System

equil.mdin

Implicit solvent constant pH molecular dynamics
 &cntrl
  imin=0, irest=1, ntx=5,
  ntpr=1000, ntwx=1000, nstlim=1000000,
  dt=0.002, ntt=3, tempi=300,
  temp0=300, tautp=2.0, ig=-1,
  ntp=0, ntc=2, ntf=2, cut=30,
  ntb=0, igb=2, saltcon=0.1,
  nrespa=1, tol=0.000001, icnstph=1,
  solvph=4.5, ntcnstph=5, gamma_ln=5.0,
  ntwr=10000, ioutfm=1,
 /

mpirun -np 12 pmemd.MPI -O -i equil.mdin -p 6pti.parm7 -c 6pti.heat.rst7 -cpin 6pti.cpin -o 6pti.equil.mdout -cpout 6pti.cpout -r 6pti.equil.rst7 -x 6pti.equil.nc -cprestrt 6pti.equil.cpin

The cprestrt file should be used as the cpin file for the next production simulation.

6. Constant pH Production Run

md.mdin

Implicit solvent constant pH molecular dynamics
 &cntrl
  imin=0, irest=1, ntx=5,
  ntpr=1000, ntwx=1000, nstlim=1000000,
  dt=0.002, ntt=3, tempi=300,
  temp0=300, tautp=2.0, ig=-1,
  ntp=0, ntc=2, ntf=2, cut=30,
  ntb=0, igb=2, saltcon=0.1,
  nrespa=1, tol=0.000001, icnstph=1,
  solvph=4.5, ntcnstph=5, gamma_ln=5.0,
  ntwr=10000, ioutfm=1,
 /

mkdir 0(1,2,3,4,5,6,7,8,9,10,11,12,13,14)

edit pH_0.mdin (pH_1.mdin......)

Change solvph=<ph>

mpirun -np 12 pmemd.MPI -O -i pH_0.mdin -p 6pti.parm7 -c 6pti.equil.rst7 -cpin 6pti.equil.cpin -o 0/6pti.md1.mdout -cpout 0/6pti.md1.cpout -r 0/6pti.md1.rst7 -x 0/6pti.md1.nc -cprestrt 0/6pti.md1.cpin

7. Calculate pKas

We use cphstats  to parse cpout file and calculation pKa values.

 for ph in $(seq 0 14); do cphstats -i 6pti.cpin ./${ph}/6pti.md1.cpout -o pH${ph}_calcpka.dat --population pH${ph}_populations.dat; done

This will generate 2 files for each pH: pH<ph>_calcpka.dat  & pH<ph>_populations.dat

amber_pka.pl 0_14

titration_curves.dat

# pH ASP3 GLU7 TYR10 LYS15 TYR21 TYR23 LYS26 TYR35 LYS41 LYS46 GLU49 ASP50
0.00 0.007 0.026 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.004
1.00 0.064 0.098 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.008 0.035
2.00 0.696 0.866 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.071 0.275
3.00 0.936 0.993 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.256 0.704
4.00 0.996 1.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.813 0.941
5.00 1.000 1.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.976 0.991
6.00 1.000 1.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.996 1.000
7.00 1.000 1.000 0.002 0.001 0.000 0.000 0.000 0.000 0.014 0.002 0.999 1.000
8.00 1.000 1.000 0.029 0.011 0.000 0.000 0.007 0.000 0.072 0.016 1.000 1.000
9.00 1.000 1.000 0.142 0.070 0.005 0.000 0.056 0.006 0.341 0.121 1.000 1.000
10.00 1.000 1.000 0.550 0.493 0.036 0.054 0.367 0.021 0.869 0.557 1.000 1.000
11.00 1.000 1.000 0.808 0.882 0.260 0.088 0.819 0.272 0.881 0.920 1.000 1.000
12.00 1.000 1.000 0.985 0.985 0.675 0.557 0.976 0.614 0.986 0.987 1.000 1.000
13.00 1.000 1.000 0.998 0.998 0.959 0.913 0.997 0.493 1.000 0.999 1.000 1.000
14.00 1.000 1.000 1.000 1.000 0.995 0.995 1.000 0.925 1.000 1.000 1.000 1.000

 

gnuplot titration_curves.gpi

bpti_titration_amber_IM

8. Multiple Runs for Standard Deviation

pK1/2 of titrable residues in BPTI
#res name trial1 trail2 trial3 trial4 trial5 ave\pmstd
3 ASP 1.76 1.84 2.14 1.83 1.75 1.86\pm0.14
7 GLU 1.54 1.55 2.36 2.11 1.57 1.83\pm0.34
10 TYR 9.96 9.49 10.10 9.62 10.02 9.65\pm0.40
15 LYS 10.03 10.12 10.09 10.09 10.08 10.08\pm0.03
21 TYR 11.58 11.56 11.56 11.55 11.51 11.55\pm0.02
23 TYR 11.92 11.91 11.23 11.46 11.41 11.59\pm0.28
26 LYS 10.27 10.27 10.27 10.28 10.29 10.28\pm0.01
35 TYR 12.18 12.99 12.42 12.75 12.32 12.53\pm0.30
41 LYS 9.27 9.07 9.43 9.28 9.58 9.33\pm0.17
46 LYS 9.90 9.72 9.89 9.82 9.87 9.84\pm0.06
49 GLU 3.41 3.63 3.73 3.55 3.51 3.57\pm0.11
50 ASP 2.53 2.23 2.69 2.21 2.05 2.34\pm0.23

 

References:

http://ambermd.org/tutorials/advanced/tutorial18/

 http://archive.ambermd.org/201510/0064.html