Gō models

The MartiniGō model is a method of maintaining secondary and tertiary structure using native contacts of proteins to create a Gō-like model between beads. In contrast to an elastic network, the Gō model applies non-bonded interactions between pairs of beads within a protein based on residue overlap and restricted chemical structural unit criteria.

The latest version of Martinize2 (version ≥ 0.10.0) implements the newest version of the Gō model. In this version of the Gō model, interactions are mediated through the addition of extra virtual sites on top of backbone beads in the protein. Interactions are in the form of Lennard-Jones interactions, which are written as an extra file to be included in the protein topology.

The Gō model is described in the help:

Virtual site based GōMartini:
  -go [GO]              Use Martini Go model. Accepts either an input file from the server, or just provide the flag to
                        calculate as part of Martinize. (default: None)
  -go-eps GO_EPS        The strength of the Go model structural bias in kJ/mol. (default: 9.414)
  -go-low GO_LOW        Minimum distance (nm) below which contacts are removed. (default: 0.3)
  -go-up GO_UP          Maximum distance (nm) above which contacts are removed. (default: 1.1)
  -go-res-dist GO_RES_DIST
                        Minimum graph distance (similar sequence distance) below which contacts are removed. (default: 3)
  -go-write-file [GO_WRITE_FILE]
                        Write out contact map to file if calculating as part of Martinize2. (default: None)
  -go-backbone GO_BACKBONE
                        Name of protein backbone bead on which to place a virtual interaction go site (default: BB)
  -go-atomname GO_ATOMNAME
                        Name of the virtual interaction go site atom (default: CA)

To add a Gō model to your protein, the first step is to calculate the contact map of your protein. The contact map can be obtained in two ways. Firstly, by uploading your pdb structure to the web server, and downloading the associated contact_map.out file. Alternatively, the contact map can be calculated directly without the need for any external processes. While the implementations of the contact algorithm are identical, the Martinize2 implementation may be relatively slow for larger systems. Typically for proteins with fewer than 1000 residues, the calculation of the contact map as part of Martinize2 will add up to a minute of extra calculation. Note that while the implementations of the main algorithm are identical, there may be small differences in the resulting contact map files due to assumptions the server makes about the format of input pdb files, which the implementation in Martinize2 overcomes. If you want to check the contact map that Martinize2 has calculated, you can write it out using the -go-write-file argument. While the contact map files may have small differences, it is likely that they will still result in the same non-bonded file outputs, a result of how symmetrical contacts are further identified in the definition of the Gō model.

The Gō model is added to the input system using the -go argument of martinize2. If you have used a contact map from the server, give the path to the contact map file as the argument:

martinize2 -f protein.pdb -o topol.top -x cg_protein.pdb -ff martini3001 -dssp -go contact_map.out

Otherwise the contact map is calculated as part of Martinize2 by just specifying the -go argument:

martinize2 -f protein.pdb -o topol.top -x cg_protein.pdb -ff martini3001 -dssp -go

Without any further additions, this will:

  1. Generate virtual sites with the atomname CA directly on top of all the backbone beads BB in your protein. The atomname and the name of the backbone bead are controlled using the -go-atomname and -go-backbone flags respectively. Each CA atom has an underlying atomtype (see below), which has a default name inherited from the -name flag.

  2. Use the contact map to generate a set of non-bonded parameters between specific pairs of CA atoms in your molecule with strength 9.414 kJ/mol (changed through the -go-eps flag).

  3. Eliminate any contacts which are shorter than 0.3 nm and longer than 1.1 nm, or are closer than 3 residues in the molecular graph. These options are flexible through the -go-low and -go-up flags.

  4. If the contact map finds any atoms within contact range defined, but are also within 3 residues of each other, then the contacts are removed. This is defined through the -go-res-dist flag.

As a result, along with the standard output of martinize2 (i.e. itp files for your molecules, a generic .top file, a coarse grained structure file), you will get two extra files: go_atomtypes.itp and go_nbparams.itp. The atomtypes file defines the new virtual sites as atoms for your system, and the nbparams file defines specific non-bonded interactions between them.

For example, go_atomtypes.itp looks like any other [ atomtypes ] directive:

[ atomtypes ]
molecule_1 0.0 0 A 0.00000000 0.00000000
molecule_2 0.0 0 A 0.00000000 0.00000000
molecule_3 0.0 0 A 0.00000000 0.00000000
molecule_4 0.0 0 A 0.00000000 0.00000000
molecule_5 0.0 0 A 0.00000000 0.00000000
...

Similarly, go_nbparams.itp looks like any [ nonbond_params ] directive (obviously, the exact parameters here depend on your protein):

[ nonbond_params ]
molecule_17 molecule_0_13 1 0.59354169 9.41400000 ;go bond 0.666228018941817
molecule_18 molecule_0_14 1 0.53798937 9.41400000 ;go bond 0.6038726468003999
molecule_19 molecule_0_15 1 0.51270658 9.41400000 ;go bond 0.5754936778307316
molecule_22 molecule_0_15 1 0.73815666 9.41400000 ;go bond 0.8285528398039018
molecule_22 molecule_0_18 1 0.54218134 9.41400000 ;go bond 0.6085779754055839
molecule_23 molecule_0_19 1 0.53307395 9.41400000 ;go bond 0.5983552758317587
...

To activate your Gō model for use in Gromacs, the martini_v3.0.0.itp master itp needs the additional files included. The additional atomtypes defined in the go_atomtypes.itp file should be included at the end of the [ atomtypes ] directive as:

[ atomtypes ]
...
TX1er 36.0 0.000 A 0.0 0.0
W  72.0 0.000 A 0.0 0.0
SW 54.0 0.000 A 0.0 0.0
TW 36.0 0.000 A 0.0 0.0
U  24.0 0.000 A 0.0 0.0

#ifdef GO_VIRT
#include "go_atomtypes.itp"
#endif

[ nonbond_params ]
P6    P6  1 4.700000e-01    4.990000e+00
P6    P5  1 4.700000e-01    4.730000e+00
P6    P4  1 4.700000e-01    4.480000e+00
...

Similarly, the nonbonded parameters should be included at the end of the [ nonbond_params ] directive:

...
TX2er  SQ1n  1 3.660000e-01    3.528000e+00
TX2er  TQ1n  1 3.520000e-01    5.158000e+00
TX1er   Q1n  1 3.950000e-01    1.981000e+00
TX1er  SQ1n  1 3.780000e-01    3.098000e+00
TX1er  TQ1n  1 3.660000e-01    4.422000e+00

#ifdef GO_VIRT
#include "go_nbparams.itp"
#endif

Then in the .top file for your system, simply include #define GO_VIRT along with the other files to be included to active the Gō network in your model.

As a shortcut for writing the include statements above, you can simply include these files in your master martini_v3.0.0.itp file with the following commands in a bash shell:

sed -i "s/\[ nonbond_params \]/\#ifdef GO_VIRT\n\#include \"go_atomtypes.itp\"\n\#endif\n\n\[ nonbond_params \]/" martini_v3.0.0.itp

echo -e "\n#ifdef GO_VIRT \n#include \"go_nbparams.itp\"\n#endif" >> martini_v3.0.0.itp

The Gō model should then be usable in your simulations following the general protein tutorial. But careful! While the Gō model specifies nonbonded interactions, the interactions are only defined internally for each molecule. This means that if you have multiple copies of your Gō model protein in the system, the Gō bonds are still only specified internally for each copy of the molecule, not truly as intermolecular forces in the system as a whole. For more detail on this phenomenon, see the paper by Korshunova et al..

Multiple Gō models in the same system

If you have several different proteins that you have martinized with their own Gō models, then several extra steps are required to ensure that they can be simulated together.

Firstly, when the proteins are coarse grained with Martinize2, the -name flag must be used to ensure that all the virtual sites created for the purposes of the Gō model have unique names and atomtypes. For example:

martinize2 -f protein.pdb -x cg.pdb -o topol.top -go contact_map.out -name my_protein

will ensure that the atoms created to apply Gō sites to are names my_protein_{0..n} for a protein of n residues.

Having martinized the proteins in this way, all the go_atomtypes.itp and go_nbparams.itp files generated for each protein should be concatenated into a single file, which may then be included in the master force field file as described above.

Note on multiple Gō models

While the Gō model is expressed as a set of interactions between Gō sites on a protein, interactions are not generally extended over all copies of a protein. That is, if a simulation is set up with several copies of a multimeric protein where the monomers are held together by a Gō model, Gromacs will not permit the multimers to “fall apart” for monomers to find other monomers that were initially in other complexes. This limitation is a result of the nature of the Gromacs itp format, with the Gō interactions described within each [ moleculetype ] directive. This issue is discussed more extensively in a recent paper by Korshunova et al..

Visualising Gō networks

If you want to look at your Gō network in VMD to confirm that it’s been constructed in the way that you’re expecting, the MartiniGlass package can help write visualisable topologies to view.