Go Back

Section - 1 - Molecular Simulation Using VMD ( Membrane Protein )

Tutorial on Membrane Protein Simulation [Click to View]

Section - 2 - Molecular Simulation Using Gromacs

[A] Gromacs Installation and Checking

Check if GROMACS is installed.

How? Type: gmx on terminal.

Downloading Protein structure file

We must download the protein structure file we will be working with. For this tutorial, we will utilize T4 lysozyme L99A/M102Q (PDB code 3HTB). Go to the RCSB website and download the PDB text for the crystal structure.

[B]Clean the PDB file

Command to extract protein
grep '^ATOM' 3htb.pdb > protein.pdb

Command to extract ligand
grep '^HETATM' 3htb.pdb | grep 'JZ4' > ligand.pdb

If you want to make a complex pdb cat protein.pdb ligand.pdb > complex.pdb

[C]Prepare the Protein topology with pdb2gmx

Now that the crystal waters are gone and we have verified that all the necessary atoms are present, the PDB file should contain only protein atoms, and is ready to be input into the first GROMACS module, pdb2gmx. The purpose of pdb2gmx is to generate three files:
1. The topology for the molecule.
2. A position restraint file.
3. A post-processed structure file.

Execute pdb2gmx by issuing the following command:
gmx pdb2gmx -f protein.pdb -o protein_proc.gro -water spc

[D]Prepare the ligand topology using external tools

Go to the PRODRG website:
(a) http://davapc1.bioch.dundee.ac.uk/cgi-bin/prodrg
(b) Upload ligand.pdb file.

Chirality (yes/no): maintain chiral centers in the input molecule (yes), or reconstruct them (no). In our case, we have no chirality, so leave this option set to its default.

Charges (full/reduced): full charges are for condensed-phase systems (43A1 force field); reduced are for "vacuum" simulations (43B1 force field). Use full charges for nearly all applications. The accuracy of 43B1 is debatable, anyway.

EM (yes/no): perform energy minimization (yes), or leave coordinates alone (no). In our case, we want to construct our protein-ligand complex from existing coordinates, so we do not want PRODRG to minimize our molecule in vacuo. Choose "no" for this option.

Edit the topology and itp of ligand generated.

Including the parameters for the JZ4 ligand in the system topology is very easy. Just insert a line that says #include "drg.itp" into topol.top after the position restraint file is included. The inclusion of position restraints indicates the end of the "Protein" molecule type section. Concatenate protein and ligand topologies into complex, after add 15 to complex.gro header (2nd line).
cat protein_proc.gro jz4.gro > complex.gro Insert #include "drg.itp" into topol.top after the position restraint file is included. After line number 10647.

TOP

[F] Define the box

Let's define the box using editconf:

At this point, the workflow is just like any other MD simulation. We will define the unit cell and fill it with water.
gmx editconf -f 1AKI_processed.gro -o 1AKI_newbox.gro -c -d 1.0 -bt cubic

Now that we have defined a box, we can fill it with solvent (water). Solvation is accomplished using solvate:
gmx solvate -cp 1AKI_newbox.gro -cs spc216.gro -o 1AKI_solv.gro -p topol.top

[G] Adding charges to the system

Download em.mdp
Use grompp to assemble a .tpr file, using any .mdp file.
gmx grompp -f em.mdp -c solv.gro -p topol.top -o ions.tpr

Now we have an atomic-level description of our system in the binary file ions.tpr. We will pass this file to genion:
gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -nn 6

[H] Energy Minimization

Now that the system is assembled, create the binary input using grompp using this input parameter file: First download em_real.mdp

gmx grompp -f em_real.mdp -c solv_ions.gro -p topol.top -o em.tpr
Run the minimization

We are now ready to invoke mdrun to carry out the EM:
gmx mdrun -v -deffnm em
Observe the convergence of the system

[I] Equilibration

Equilibrating our protein-ligand complex will be much like equilibrating any other system containing a protein in water. There are a few special considerations, in this case:
Applying restraints to the ligand Treatment of temperature coupling groups

[J] Restraining the Ligand

PRODRG did not generate a file analogous to posre.itp for our ligand, but GROMACS provides the means to do so with the genrestr module. Simply run genrestr on the jz4.gro file that we obtained from PRODRG:

gmx genrestr -f jz4.gro -o posre_jz4.itp -fc 1000 1000 1000

Now, we need to include this information in our topology. We can do this in several ways, depending upon the conditions we wish to use. If we simply want to restrain the ligand whenever the protein is also restrained, add the following lines to your topology in the location indicated:

	; Include Position restraint file
	#ifdef POSRES
	#include "posre.itp"
	#endif
	
	; Include ligand topology
	#include "drg.itp"
	
	; Ligand position restraints
	#ifdef POSRES
	#include "posre_jz4.itp"
	#endif
	
	; Include water topology
	#include "gromos43a1.ff/spc.itp"

[K] Thermostats

Proper control of temperature coupling is a sensitive issue. Coupling every molecule type to its own thermostatting group is a bad idea.

Make index using commands

gmx make_ndx -f em.gro -o index.ndx
Merge the "Protein" and "JZ4" groups with the following, where ">" indicates the make_ndx prompt:

		> 1 | 13
		> q
	

We can now set tc_grps = Protein_JZ4 Water_and_ions to achieve our desired "Protein Non-Protein" effect.

[L] Equilibration, Part 2

Once the NVT simulation is complete, proceed to NPT with this .mdp file:
gmx grompp -f npt.mdp -c nvt.gro -t nvt.cpt -p topol.top -n index.ndx -o npt.tpr
gmx mdrun -deffnm npt Copy nvt.mdp file
Run command
gmx grompp -f nvt.mdp -c em.gro -p topol.top -n index.ndx -o nvt.tpr
gmx mdrun -v -deffnm nvt

Once the NVT simulation is complete, proceed to NPT with this .mdp file:

Go Back

TOP