Protein in water
Keeping in line with the overall Martini philosophy, the coarse-grained protein model groups 4 heavy atoms together in one CG bead. Each residue has one backbone bead and zero or more side-chain beads depending on the amino acid type. The secondary structure of the protein influences both the selected bead types and bond/angle/dihedral parameters of each residue as explained below. It is noteworthy that, even though the protein is free to change its tertiary arrangement, local secondary structure is predefined and thus imposed throughout a simulation. Conformational changes that involve changes in the secondary structure are therefore beyond the scope of Martini CG proteins.
Setting up a CG protein simulation consists basically of two steps:
1. converting an atomistic protein structure into a CG model;
2. generating a suitable Martini topology.
Both steps are done using the publicly available 'martinize.py' script, of which the latest version can be downloaded here.
The aim of this module is to set up a CG simulation of ubiquitin in a water box. In this part of the tutorial, some basic knowledge of gromacs commands is assumed and not all commands will be given explicitly. Please refer to the previous tutorials and/or the gromacs manual. After getting the atomistic structure of ubiquitin (1UBQ), you'll need to convert it into a CG structure and to prepare a Martini topology for it. Once this is done the CG structure can be minimized, solvated and simulated. The steps you need to take are roughly the following:
1. Download 1UBQ.pdb from the Protein Data Bank.
gzip -d 1UBQ.pdb.gz
2. The pdb-structure can be used directly as input for the martinize.py script, to generate both a structure and a topology file. Have a look at the help function (i.e. run martinize.py -h) for the available options. Hint: During equilibration it might be useful to have (backbone) position restraints. The final command might look a bit like this:
./martinize.py -f 1UBQ.pdb -o system.top -x cg_1UBQ.pdb
-dssp /wherever/dsspcmbi -p backbone -ff martini22
When using the -dssp option you'll need the dssp binary, which determines the secondary structure classification of the protein backbone from the PDB-structure. It can be downloaded from a CMBI website:
As an alternative, you may prepare a file with the required secondary structure yourself and feed it to the script:
./martinize.py -f 1UBQ.pdb -o system.top -x cg_1UBQ.pdb
-ss your-sec-struct-file -p backbone -ff martini22
NOTE: The martinize.py script might not work with older versions of python! Also python 3.x migth not work. We know it does work with versions: [2.6.x, 2.7.x] and does not work with [2.4.x]. If you have tested it with any version inbetween, we would like to hear from you.
3. If everything went well, the script generated three files: a coarse grain structure (.gro/.pdb), a master topology file (.top), and a protein topology file (.itp). In order to run a simulation you need two more: the Martini topology file (martini v2.1.itp) and a run parameter le (.mdp). You can get examples from the Martini website or from the protein-tutorial package. Don't forget to adapt the settings where needed!
4. Do a short (+/-10 steps!) minimization in vacuum. (Before you can generate the run-input files with grompp, you have to generate a box, for example using editconf).
5. Solvate the system with genbox (an equilibrated water box can be downloaded here) and minimize. Make sure the box-size is large enough (i.e. there is enough water around the molecule) and remember to use a larger van der Waals distance when solvating to avoid clashes, e.g.:
genbox -cp confout.gro -cs water-1bar-303K.gro -vdwd 0.21 -o solvated.gro
6. Do a short energy minimization and position restrained simulation. Since the martinize.py script already generated position restraints (thanks to the -p option), all you have to do is specify
define = -DPOSRES
in your .mdp file and add the appropriate number of water (W) beads to your .top file.
7. Start production run (without position restraints!).
8. . . .
9. PROFIT! What sort of analysis can be done on this molecule? Start by having a look at the protein with vmd (section Visualizing Martini systems using vmd).
A mapped fine-grained simulation of ubiquitin for comparison is available in the archive mapubq.tar.gz
Protein and elastic network
The aim of this module is to see how application of elastic networks can be combined with the Martini model to conserve tertiary and quartenary structures more faithfully without sacricing realistic dynamics of a protein. We offer two alternative routes. Please be advised that this is an active field of research and that there is as of yet no "gold standard". The first option is to generate a simple elastic network on the basis of a standard Martini topology. The second options is to use an ElNeDyn network. This second option constitutes quite some change to the Martini forcefield and thus is a different forcefield! The advantage is that the behavior of the method has been well described*. Both approaches can be set up using the martinize.py script and will be shortly described below.
First you should simulate a normal Martini CG protein without an elastic network and then see what changes when you use a CG topology with an elastic network.
1. Download HIV-1 Protease (1A8G.pdb).
2. Repeat steps 2-7 from the previous exercise. (You'll need to simulate the protein for hundreds of nanoseconds to see major changes in the structure, a sample long simulation is provided in the archive 1A8G_trajectories.tar.gz).
3. Visualize the simulation; look especially at the binding pocket of the protein: does it stay closed, open up, or what? What happens to the overall protein structure?
Martini + elastic network
The first option to help preserve higher-order structure of proteins is to add to the standard Martini topology extra harmonic bonds between non-bonded beads based on a distance cut-off. Note that in standard Martini, long harmonic bonds are already used to impose the secondary structure of extended elements (sheets) of the protein. The martinize.py script will generate harmonic bonds between backbone beads if the options -elastic is set. It is possible to tune the elastic bonds (e.g.: make the force constant distance dependent, change upper and lower distance cut-off, etc.) in order to make the protein behave properly. The only way to f8nd the proper parameters is to try different options and compare the behavior of your protein to an atomistic simulation or experimental data (NMR, etc.). Here we will use basic parameters in order to show
1. Use the martinize.py script to generate the coarse grain structure and topology as above. For the elastic network options use:
-elastic -ef 500 -el 0.5 -eu 0.9 -ea 0 -ep 0
This turns on the elastic network (-elastic), sets the elastic bond force constant to 500 kJ mol-1 nm-2 (-ef 500), the lower and upper elastic bond cut-o to 0.5 and 0.9 nm, respectively, (-el 0.5 -eu 0.9) and makes the bond strengths independent of the bond length (elastic bond decay factor and decay power, -ea 0 -ep 0, respectively; these are default). The elastic network is defined in the .itp file by an "ifdef" statement (have a look!). The "#define RUBBER BANDS" in the .top or .itp file switches it on. Note that martinize.py does not generate elastic bonds between i->i+1 and i->i+2 backbone bonds, as those are already connected by bonds and angles.
2. Proceed as before and start a production run. If you are using a Gromacs version older then 4.5.x you have to explicitly set the "define = -DRUBBER_BANDS" in the mdp-file. In 4.5.x and further this is automatically set in the top file.
The second option to use elastic networks in combination with Martini puts more emphasis on remaining close to the overall structure as deposited in the PDB than standard Martini does. The main difference from the standard way (used in the previous exercise) is the use of a global elastic network between the backbone beads to conserve the conformation instead of relying on the
angle/dihedral potentials and/or local elastic bonds to do the job. The position of the backbone beads is also slightly dierent. In standard Martinithe center-of-mass of the peptide plane is used as the location of the backbone bead, but in the elastic network implementation the positions of the C-atoms are used for this.
The martinize.py script automatically sets these options and sets the correct parameters for the elastic network. As the elastic bondstrength and the upper cut-off have to be tuned in an ElNeDyn network, these options can be set manually (-ef -eu).
1. Use the martinize.py script to generate the coarse grain structure and topology as above. Set the following options (-ef 500 and -eu 0.9 are the default options, so can be left away):
-elnedyn -ef 500 -eu 0.9
This selects the ElNeDyn topology options, with a force constant of 500 kJ mol-1nm-2 (-ef 500) and an upper bond length cut-o of 0.9 nm (-eu 0.9).
2. Proceed as before and start a production run.
Now you've got three simulations of the same protein with different elastic networks. If you do not want to wait, some pre-run trajectories can be found in: 1A8G_trajectories.tar.gz One of them might fit your needs in terms of structural and dynamic behavior. If not, there are an almost innite number of ways to further tweak the elastic network!
An easy way to compare the slightly dierent behaviors of the proteins in the previous three cases is to follow deviation/fluctuation of the backbone during simulation (and compare it to an all-atom simulation if possible). vmd provides a set of friendly tools to compute RMSD/RMSF, but it needs some tricks to be adapted to CG systems ("protein", "water", etc. keywords are not understood by vmd on CG structures). In this specific case, to visualize/select only the backbone you need to create a representation containing only the corresponding beads: "name BB". See section on Coarse Grain vizualization for more detailed howtos and advice.
*X. Periole, M. Cavalli, S.J. Marrink, M.A. Ceruso, Combining an Elastic Network With a Coarse-Grained Molecular Force Field: Structure, Dynamics, and Intermolecular Recognition, J. Chem. Theory Comput., 5, 2531-2543 (2009)