(Note: These tutorials are meant to provide illustrative examples of how to use the AMBER software suite to carry out simulations that can be run on a simple workstation in a reasonable period of time. They do not necessarily provide the optimal choice of parameters or methods for the particular application area.)
Copyright Ross Walker 2015

Analysis of water thermodynamics using Grid Inhomogeneous Solvation Theory of Factor Xa active site. SECTION 1

By Romelia Salomon, Crystal Nguyen, Steven Ramsey, Jonathan D. Gough, Ross Walker, & Tom Kurtzman

GIST Analysis Example
Analysis of water thermodynamics using Grid Inhomogeneous Solvation Theory of Factor Xa active site.

By Romelia Salomon, Crystal Nguyen, Steven Ramsey, Jonathan D. Gough, Ross Walker, & Tom Kurtzman

1) Generating and Relaxing the Initial Structure

GIST solvation thermodynamic mapping is performed on explicit water molecular dynamics simulations, typically containing a solute of interest that has been solvated using Leap. In this tutorial we will characterize the binding site water of coagulation factor Xa (fxa, 1FJS.pdb). In order to build our solvated simulation system we first need to prep a pdb structure file obtained from the protein databand (pdb, www.rcsb.org) so that leap can correctly read the structure, see the Leap tutorial and tutorial B0 for more information on this subject. For simplicity's sake in this tutorial we provide an initial pdb structure, provided below, that has been prepared for leap. This structure originated from pdb:1FJS.

fxa.pdb

Please note that since this structure has been preprocessed we have skipped some steps, such as removing unwanted crystallographic waters or ligand molecules, you should be cognizant of these structural elements when constructing your system of interest. We then use the program Leap, included in AmberTools, to generate the necessary parameter and coordinate files.

In this tutorial we will solvate the system in TIP3P water with no ions to balance charges since we will be restraining all heavy atoms during the production simulation, it is recommended to add ions when simulating an unrestrained system. The following leap commands can be used to load the provided pdb structure, solvate the system and output an AMBER parameter file (.prmtop) and restart coordinate file (.rst7):

$AMBERHOME/bin/tleap
> source leaprc.protein.ff14SB
> source leaprc.water.tip3p
> p = loadpdb fxa.pdb
> solvatebox p TIP3PBOX 10
> saveamberparm p fxa.prmtop fxa.rst7
> quit

Leap outputs the following files which will be used as the paramater-topology and the initial coordinates for minimizing and equilibrating the system: fxa.prmtop, fxa.rst7.

Using these initial structure files, we equilibrate the system before performing the production MD simulation. This process entails minimizing, heating, and relaxing the system before performing the production run. Here the equilibration is performed in 6 sequential stages:

Minimize only the water while restraining the protein (20000 cycles)

Minimize the water and the protein hydrogen atoms (20000 cycles)

Heat the system to 50K while restraining the protein heavy atoms (NVT 0 to 50K)

Heat the system from 50K to 300K in increments of 50K for 200ps each

Relax the system, restraining the protein heavy atoms (NPT, 300K, 10ns)

Relax the system, restraining all atoms in the protein (NVT, 300K, 5ns)

min_sys-1.in

min_sys-2.in

md-heat-50K.in

minimize water
 &cntrl
   imin=1, ntmin=1, nmropt=0, drms=0.1,
   maxcyc=20000, ncyc=1500,
   ntx=1,	// not a restart
   irest=0,
   ntpr=100, ntwr=100, iwrap=1,
   ntf=1,	// SHAKE is off
   ntb=1,	 // constant volume pbc 
   cut=9.0,
   ntr=1, restraint_wt=100, // Cartesian restraint
   restraintmask='!:WAT',
 /

minimize all
 &cntrl
   imin=1, ntmin=1, nmropt=0, drms=0.1,
   maxcyc=20000, ncyc=1500,
   ntx=1, irest=0,
   ntpr=100, ntwr=100, iwrap=1,
   ntf=1, ntb=1, cut=9.0,
   ntr=1,
   restraint_wt=100,
   restraintmask='!:WAT & !@H',
 /

heat NVT 200ps
 &cntrl
  imin = 0,
  irest = 0, ntx = 1, 
  ntb = 1, ntc = 2, ntf = 2,
  cut = 9.0,
  tempi = 0.0, temp0 = 50.0,
  ntt = 3, gamma_ln = 2.0,
  nstlim = 10000, dt = 0.002,
  iwrap = 1,
  ntpr = 1000, ntwx = 1000, ntwr = 10000,
  nscm = 1000, 
  ntr=1, restraint_wt=100,
  restraintmask='!:WAT & !@H',
  ig=-1,
  ioutfm=1,
 /

md-heat-300K.in

md-eq-npt-10ns.in

md-eq-nvt-5ns.in

heat NPT 200ps
 &cntrl
  imin = 0,
  irest = 1, ntx = 5,
  ntb = 2, ntc = 2, ntf = 2,
  cut = 9.0,
  temp0 = 300.0,
  ntt = 3, gamma_ln = 2.0,
  nstlim = 10000, dt = 0.002,
  ntp = 1, taup = 0.5,
  iwrap = 1,
  ntpr = 1000, ntwx = 1000, ntwr = 10000,
  nscm = 1000,
  ntr=1, restraint_wt=100,
  restraintmask='!:WAT & !@H',
  ig=-1,
  ioutfm=1,
 /

NTP 10ns
 &cntrl
  imin=0,
  ntpr   = 10000,    ntwx   = 10000,
  ntwr   = 10000,
  iwrap=1,
  ntx=5, irest=1,
  ntf    = 2,       ntb    = 2,
  ntc    = 2,       cut    = 10.0,
  nstlim = 5000000,
  dt     = 0.002,
  nscm   = 1000,
  temp0  = 300.0, tempi = 300.0,
  ntt    = 3, gamma_ln = 2.0,
  ntp    = 1, taup   = 0.5,
  ntr=1, restraint_wt=100,
  restraintmask="!:WAT  & !@H=",
  ig=-1,
  ioutfm=1,
 /

equil NVT 5ns
 &cntrl
  imin=0,
  ntpr   = 10000,    ntwx   = 10000,
  ntwr   = 10000,
  iwrap=1,
  ntx=5, irest=1,
  ntf    = 2,       ntb    = 1,
  ntc    = 2,       cut    = 9.0,
  nstlim = 2500000,
  dt     = 0.002,
  nscm   = 1000,
  temp0  = 300.0,
  ntt    = 3, gamma_ln = 2.0,
  ntr=1, restraint_wt=100,
  restraintmask="!:WAT",
  ig=-1,
  ioutfm=1,
 /

Note that our example simulation is performed with strong heavy atom restraints of 100 kcal/Å^2 to ensure the protein does not translate or rotate. GIST requires a simulation where the protein or host does not translate or rotate, here 100 kcal/Å^2 is used to ensure limited mobility but end-users are welcome to apply any restraint strength that will prevent protein translation and rotation.

The following is an example of the commands to run the equilibration, using pmemd (Note that these steps will take significant time, we suggest using pmemd.mpi or pmemd.cuda for all steps other than minimization, please refer to the AMBER manual for information on utilizing these versions):

$AMBERHOME/bin/pmemd -O -i min_sys-1.in -o min1.out -p fxa.prmtop -c fxa.rst7 -r min_fxa-1.rst7

$AMBERHOME/bin/pmemd -O -i min_sys-2.in -o min2.out -p fxa.prmtop -c min_fxa-1.rst7 -r min_fxa-2.rst7 -ref min_fxa-1.rst7

$AMBERHOME/bin/pmemd.cuda -O -i md-heat-50K.in -o heat1.out -p fxa.prmtop -c min_fxa-2.rst7 -r md_heat-50K.rst7 -x heat1.nc -ref min_fxa-2.rst7

$AMBERHOME/bin/pmemd.cuda -O -i md-heat-300K.in -o heat2.out -p fxa.prmtop -c md_heat-50K.rst7 -r md_heat-300K.rst7 -x heat2.nc -ref md_heat-50K.rst7

$AMBERHOME/bin/pmemd.cuda -O -i md-eq-npt-10ns.in -o equil-npt.out -p fxa.prmtop -c md_heat-300K.rst7 -r md_eq-npt.rst7 -x equil-npt.nc -ref md_heat-300K.rst7

$AMBERHOME/bin/pmemd.cuda -O -i md-eq-nvt-5ns.in -o equil-nvt.out -p fxa.prmtop -c md_eq-npt.rst7 -r md_eq-nvt.rst7 -x equil-nvt.nc -ref md_eq-npt.rst7

After running these simulations, we get the following files for the different stages (due to random seeds your result files may differ slightly):

Minimizing only the water gives the output file min1.out and the restart file min_fxa-1.rst7

Minimizing the water and protein Hydrogen atoms give the output file min2.out and the restart file min_fxa-2.rst7

Heating under NPT conditions from 0-50K gives the output file heat1.out the trajectory heat1.nc and the restart file md_heat-50K.rst7
Heating under NPT conditions from 250-300K gives the output heat2.out the trajectory heat2.nc and the restart file md_heat-300K.rst7

Relaxing the system under NPT conditions gives the output file equil-npt.out the trajectory equil-npt.nc and the restart file md_eq-npt.rst7

Relaxing the system under NVT conditions gives the output file equil-nvt.out the trajectory equil-nvt.nc and the restart file md_eq-nvt.rst7

If you run the calculations as stated above, you won't get exactly the same files, due to the use of different random seeds, but an equivalent initial configuration

CLICK HERE TO RETURN TO THE TITLE PAGE CLICK HERE TO GO TO SECTION 2