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.

The residues titratable in BPTI are displayed in picture below:

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

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

8. Multiple Runs for Standard Deviation

 

References:

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

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