Task 8: Molecular Dynamics Simulations

Task description

A detailed task description can be found here.

Intro and selection of two mutants

The description given here was applied for our wild type (1J8U). However, the steps done for our two mutants R408W and P281L are the same.

We selected R408 because it is the most abundant mutation of PAH which is associated with phenylketonuria. It has a total frequency of 6.67% (see [1]).

The second mutation, P281L, was selected because it is the closed mutation to the binding site (HIS 285) and we hope that we are able to see some interesting things here which give us an explanation why this mutation is disease causing.

Preperation

Before we can start to generate all necessary files for our MD simulation we have to prepare our structure first. We employed the following steps:

1. Extracting the crystal water below 15 Å

repairPDB 1J8U.pdb -ssw 15 > water_below_15.out

We received the following output:

HETATM 2559  O   HOH A1008       2.996   9.738  16.094  1.00 14.91           O  
HETATM 2561  O   HOH A1010      -3.667  12.788   8.538  1.00 14.80           O  
HETATM 2572  O   HOH A1021       3.557  27.911  10.913  1.00 14.69           O  
HETATM 2879  O   HOH A1328      -5.335  24.271  14.490  1.00 14.88           O  
TER

2. Extract only the protein

repairPDB 1J8U.pdb -nosol > 1J8U_nosol.pdb

3. Extract the amino acid sequence and turn it into lower case letter to give it as a input for SCWRL

repairPDB 1J8U.pdb -seq > 1J8U_seq.txt

vim 1J8U_seq.txt
:%s/.*/\L&/g

4. Use SCWRL to complete the side chains of all amino acids

/apps/scwrl4/Scwrl4 -i 1J8U_nosol.pdb -s 1J8U_seq.txt -o 1J8U_nosol_after_SCWRL.pdb | tee scwrl_1J8U_nosol_after_scwrl.out

5. Remove the hydrogen atoms from the SCWRL output

repairPDB 1J8U_nosol_after_SCWRL.pdb -noh > 1J8U_nosol_after_SCWRL_no_h.pdb

6. concatenate the protein and the crystal water into one file

We just added the output of step 1 to the end of the PDB file 1J8U_nosol_after_SCWRL_no_h.pdb and named this file 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water.pdb

Create important GROMACS files

1. We have to create a TOP (topology) and a PAR (parameter) file with the following command:

pdb2gmx -f 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water.pdb -o 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water.gro -p 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water.top -water tip3p -ff amber03 -vsite hydrogens | tee pdb2gmx.out

2. Next we create a box around the protein and fill it with water

Create the box:

editconf -f 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water.gro -o 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water_box.gro -bt dodecahedron -d 1.0 | tee editconf.out

Fill it with water:

genbox -cp 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water_box.gro -cs spc216.gro -p 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water.top -o 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water_water.gro

3. In the next step we have to add the ions. To do so we need to use the grompp command, which is always used to prepare a system for more complicated steps. This command needs an mdp file where we can define several additional parameters.

Create parameter file addions.mdp:

vim addions.mdp	

And fill it with the following content:

integrator = steep
emtol = 1.0
nsteps = 500
nstenergy = 1
energygrps = System
coulombtype = PME
rcoulomb = 0.9
rvdw	 = 0.9
rlist = 0.9
fourierspacing = 0.12
pme_order = 4
ewald_rtol = 1e-5

Now, we execute it with this command:

grompp -v -f addions.mdp  -c 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water_water.gro -p 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water.top -o 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water_water.tpr | tee grommp.out

4. Now we add the solvent to the system and neutralize the charge that is part of the system due to charged amino acids. To do so we use genion:

genion -s 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water_water.tpr -o 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water_solv.pdb -conc 0.1 -neutral -pname NA+ -nname CL- | tee genion.out

During the call we were asked for a number we selected 13.

5. Next we have to adjust the SOL value and add the number of NA+ and CL- in our 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water.top

We had 25 NA+ and 21 CL- ions generated from our last step. So we had to reduce the number in the SOL field to 10170 - 46 = 10124 the .top file had in the end the following values:

[ molecules ]
; Compound        #mols
Protein_chain_A     1
SOL                 4
SOL             10124
NA                 25
CL                 21

6. To make sure our initial crystal water does not overlap with the added water we use repairPDB again.

repairPDB 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water_solv.pdb -cleansol > 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water_solv2.pdb

7. In the last line there is a REMARK tag with the number of removed water molecules. This has to be changed in the TOP file if larger than 0.

REMARK RM 0

Value is 0 so we did not have to change anything.

8. We need to extract restraints from the structures. These restraints are useful to reduce the simulation time by disallowing very fast vibrations such as seen for hydrogen atoms. The command we use is called genrestr.

genrestr -f 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water_solv2.pdb -o 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water_solv2.itp | tee genrestr.out

In the command line menu we choose number 1 for protein.

Minimization solvent

1. Creating a .mdp file for grompp to prepare the system for the solvent minimization

We created an file with the name minimize_solvent.mdp by using vim:

vim minimize_solvent.mdp

The content of the file is as follows:

define = -DPOSRES
integrator = steep
emtol = 1.0
nsteps = 500
nstenergy = 1
energygrps = System
coulombtype = PME
rcoulomb	 = 0.9
rvdw = 0.9
rlist = 0.9
fourierspacing = 0.12
pme_order = 4
ewald_rtol = 1e-5
pbc = xyz

2. execution of grompp to prepare the system for the MD run

grompp -v -f minimize_solvent.mdp -c 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water_solv2.pdb -p 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water.top -o 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water_solv2_min.tpr | tee grompp_minimize_solvent.out

3. Executing mdrun to minimize the solvent

mdrun -v -deffnm 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water_solv2_min -c 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water_solv2_min.pdb | tee mdrun_minimize_solvent.out

Minimization system

In this part we minimize the solvent AND the side chains of the protein.

1. Creating the position restraint file.

genrestr -f 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water_solv2_min.pdb -o 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water_solv2_min.itp | tee genrestr_minimize_system.out

We selected option 4 which says that we only put restraints on the backbone of the protein.

2. Create an .mdp file for the grompp call

We use the same .mdp from the solvent minimization step

3. Calling grompp to prepare the tpr file for our mdrun step

grompp -v -f minimize_solvent.mdp -c 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water_solv2_min.pdb -p 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water.top -o 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water_solv2_min2.tpr | tee grompp_minimize_system.out

4. Calling mdrun to minimize the system

mdrun -v -deffnm 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water_solv2_min2 -c 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water_solv2_min2.pdb | tee mdrun_minimize_system.out 

Equilibration of system

Heating up the System with constant volume and particle number (NVT)

Here we will heat up our system to a certain temperature. We will keep our particle number and volume constant, the pressure is changing.

1. creating a .mdp configuration file for grompp

We created a mdp file with the name nvt.mdp :

vim nvt.mdp

The content is as follows:

define = -DPOSRES
integrator = md
dt = 0.005 
nsteps = 10000
nstxout = 0
nstvout = 0
nstfout = 0
nstlog = 1000
nstxtcout	 = 0
nstenergy = 5
energygrps = Protein Non-Protein
nstcalcenergy = 5
nstlist = 10
ns-type = Grid
pbc = xyz
rlist = 0.9
coulombtype = PME
rcoulomb = 0.9
rvdw = 0.9
fourierspacing = 0.12
pme_order = 4
ewald_rtol = 1e-5      
gen_vel = yes   
gen_temp = 200.0 
gen_seed = 9999  
constraints = all-bonds
tcoupl = V-rescale		
tc-grps = Protein  Non-Protein 
tau_t = 0.1	0.1 	
ref_t = 298	298 
pcoupl = no

2. Creating .tpr file for our MD run by employing grompp

grompp -v -f nvt.mdp -c 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water_solv2_min2.pdb -p 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water.top -o 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water_nvt.tpr | tee grompp_equi_system_nvt.out

3. MD run

mdrun -v -deffnm 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water_nvt | tee mdrun_equi_system_nvt.out

Heating up the System with constant pressure and particle number (NPT)

Here we will heat up our system to a certain temperature. We will keep our particle number and pressure constant, the volume is changing.

1. Creating an .mdp file for pressure coupling NPT

We create a .mdp file with the following command:

vim npt.mdp

And the content:

define = -DPOSRES
integrator = md
dt = 0.005 
nsteps = 10000
nstxout = 0
nstvout = 0
nstfout = 0
nstlog = 1000
nstxtcout = 0
nstenergy = 5
xtc_precision = 1000
xtc-grps = System
energygrps = Protein Non-Protein
nstcalcenergy = 5
nstlist = 5
ns-type = Grid
pbc = xyz
rlist = 0.9
coulombtype = PME
rcoulomb = 0.9
rvdw = 0.9
fourierspacing = 0.12
pme_order = 4
ewald_rtol = 1e-5
tcoupl = V-rescale 
tc-grps = Protein Non-Protein 
tau_t = 0.1 0.1 
ref_t = 298 298 
pcoupl = Berendsen
Pcoupltype = Isotropic
tau_p = 1.0
compressibility = 4.5e-5
ref_p = 1.0
gen_vel = no
constraints = all-bonds

2. Creating .tpr file for our MD run by employing grompp

grompp -v -f npt.mdp -c 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water_nvt.gro -p 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water.top -o 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water_npt.tpr | tee grompp_equi_system_npt.out

3. MD run

mdrun -v -deffnm 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water_npt | tee mdrun_equi_system_npt.out

Production run

In this part we prepare the files for our final production run on the LRZ cluster.

1. Creating a .mdp parameter file for grompp

We created a file with the following command:

vim md.mdp

And this content:

integrator = md
tinit = 0
dt = 0.005
nsteps = 2000000
nstxout = 50000
nstvout = 50000
nstfout = 0
nstlog = 1000
nstxtcout = 1000
nstenergy = 1000
energygrps = Protein Non-Protein 
nstcalcenergy = 5
nstlist = 5
ns-type = Grid
pbc = xyz
rlist = 0.9
coulombtype = PME
rcoulomb = 0.9
rvdw = 0.9
fourierspacing = 0.12
pme_order = 4
ewald_rtol = 1e-5
tcoupl = V-rescale 
tc-grps = Protein Non-Protein 
tau_t = 0.1 0.1 
ref_t = 298 298 
pcoupl = Berendsen
Pcoupltype = Isotropic
tau_p = 2.0
compressibility = 4.5e-5
ref_p = 1.0
gen_vel = no
constraints = all-bonds
constraint-algorithm = Lincs
unconstrained-start = yes
lincs-order = 4
lincs-iter = 1
lincs-warnangle = 30
comm_mode = linear

2. Creating a .tpr file for our MD run

grompp -v -f md.mdp -c 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water_npt.gro -p 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water.top -o 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water_md.tpr | tee grompp_prduction_run_md.out

LRZ

Locally

1. Copy files to LRZ

We copied the following files to the LRZ:

• md.mdp
• 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water.top
• 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water_npt.gro
• 1J8U_nosol_after_SCWRL_no_h_merged_crystal_water_md.tpr

By using the following command:

scp -r . di69duj@lx64ia3.lrz.de:/home/cluster/pr58ni/di69duj

On the LRZ

1. Connecting to the LRZ

ssh di69duj@lx64ia3.lrz.de

2. Creating a working dictionary

mkdir md_wt

3. Copy all files to this working directory

cp *.* md_wt

And finally change to that directory:

cd md_wt

4. Creating a runfile with the name md_run.cmd

Creating the file with vim:

vim md_run.cmd

With the following content:

#!/bin/bash
#$-o $HOME/md_wt/mdrun.WT_1.out -j y
#$-N MD_WT
#$-S /bin/bash
#$-M erik.pfeiffenberger@gmail.com 
#$-l h_rt=32:00:00
#$-l march=x86_64
#$-pe mpi_32 32
. /etc/profile
cd $HOME/md_wt
$HOME/md_wt/mpirun -np 32 mdrun_mpi -v -deffnm $HOME/md_wt/1J8U_nosol_after_SCWRL_no_h_merged_crystal_water_md

5. Submitting job

qsub md_run.cmd

6. Waiting...

While waiting, we recommend a cup of earl grey tea with one quarter milk and two tea spoons of sugar. Enjoy :)