Our example is a beta-hairpin and we will base our simulations on the OV contact map (you may want to obtain more information from our manuscript here):

Preparing the necessary files:

- 1gb1.pdb: Original protein G file with 56 residues and 60 different models

- beta-hairpin.pdb: residues 41−56 of the protein G.

- martini_v2.2.itp: MARTINI force field version 2.2

To create the contact map use the go_martinize.py script in the following way (Note that we need at least version 2.7 of python):

./go_martinize.py -f beta-hairpin.pdb -o system-vaccum.top -x 1GB1-CG.pdb -dssp  ./dssp -p backbone -ff elnedyn22 -go -goepsilon 9.414

Note: the “-ff elnedyn22” is included in the command line as contacts are between Cα atoms and the parameterization of elnedyn for the MARTINI force-field is used.

On the parameters:

-o create a typical GROMACS topology file in vaccum

-x create the coarse grained structure in this case for beta-hairpin.pdb according to MARTINI representation. BB atom coincides with the position of CA atoms.

-dssp when using the -dssp option you'll need the dssp binary (here as ./dssp), which determines the secondary structure classification of the protein backbone from the structure

-ff read the Martini force field in this case version 2.2, different ff has different mapping or parameters. The GoMartini runs for v.2 martini_v2.2.itp)

-go -goepsilon: imposes Go contacts (-go) and defines the strength (-goepsilon) of the Go contacts (epsilon_ij of the LJ potential).  The units are in kJ/mol, so 9.414 corresponds to epsilon_ij=λ*epsilon with epsilon=6.276 kJ/mol (a typical value for the hydrogen bond energy in proteins). We have shown that this value gives similar results as in the case of ELNEDIN or all-atom simulations. The value is not universal and we advice before starting a large scale simulation to test values for λ in the range [1,1.5] in the case of MARTINI for RMSD/RMSF or any other structural property. However, we should note that overall the method is not very sensitive to λ. It all depends on how much precise you would like to be depending on the problem you would like to address.

Note: go_martinize.py script here works only for a single protein chain. If you would like to use multi-chain system, simply name all chains with the same letter and then remove by hand from the topology (itp file) the bonds and angles between different chains that should not be there!

Once you execute the script it will create a file Protein_A.itp which contains the list of native contacts. Search for [ pairs ] in the file and you will find the contacts list based on the OV contact map. To include rCSU contacts, please, use the server here.

You may use the version of Gromacs that suits you the best. This tutorial follows Gromacs version 4.6.5.

All files needed to run in GROMACS are in the RUN folder.

1) Energy minimization in vaccum for CG structure [set in system-vaccum.top, first line the corresponding entry for martini ff, as #include "martini_v2.2.itp" ]

- editconf -f 1GB1-CG.pdb -d 1.0 -o 1GB1-CG.gro : convert CG pdb to GRO format in nanometer with box size

- grompp -f minimization-vaccum.mdp -p system-vaccum.top -c 1GB1-CG.gro -o minimization-vaccum.tpr

- mdrun -deffnm minimization-vaccum -v

2) Solvate the minimized CG structure. If you want, you can change your simulation box size Lx,Ly and Lz by modifying

the last line in the .gro file. Then the solvation will take in consideration the new size.

- genbox -cp minimization-vaccum.gro -cs water-box-CG_303K-1bar.gro -vdwd 0.21 -o system-solvated.gro 

- grompp -f minimization.mdp -c system-solvated.gro -p system.top -o minimization.tpr

- mdrun -deffnm minimization -v

3) NPT equilibration with Parrinello-Rahman barostat

- grompp -f equilibration.mdp -c minimization.gro -p system.top -o equilibration.tpr

- mdrun -deffnm equilibration -v

4) Production RUN :), for 10 ns in this case.

- grompp -f dynamic.mdp -c equilibration.gro -p system.top -o dynamic.tpr

- mdrun -deffnm dynamic -v