Hybrid model using virtual sites

The material offered in this page is LEGACY material. This means that it either pertains to older iterations of the Martini force field and/or that it is no longer actively supported. However, this does not mean that there aren’t cool things to be learned here! Have a look! If you are looking for the current iterations of the Martini force field lectures or hands-on tutorials, click here to go back!

Aim

Demonstrate how to combine different Classical Mechanical level models within GROMACS and run a HYBRID simulation. In general, this will involve providing user-defined tabulated potentials to GROMACS. A second aim is to demonstrate a back-mapping procedure to build (partially) atomistic structures from (fully/partially) coarse-grained models.

Background

Many properties of molecular systems are local, i.e. primarily determined by the molecule and its nearest neighbors. Simulations of small systems therefore often already provide a reasonable to good representation of macroscopic systems. Nevertheless, small system simulations may contain artifacts due to the boundary conditions placed on them. Also, the minimal size of the simulated systems may still prevent sufficient sampling of phase space. The idea of HYBRID models is that a detailed description of relatively distant molecules is not required and that enhanced sampling of a region of interest may be achieved by embedding it in a surroundings that represents the environmental influences properly but at a computationally less demanding level. This approach has long been applied to the combination of quantum chemistry and molecular mechanics models (QM/MM), but is much less well developed for combining molecular models at different levels of resolution. It should be kept in mind that developments in this area are still ongoing and that there is as of yet no standard or best practice for this type of simulations.

In this tutorial, we shall focus on the combination of an atomistic (AA, or fine-grained, FG) and a coarse-grained (CG) model. The coarse-grained model is one employing a particle — as opposed to a continuum — description of the surroundings. The standard machinery present in GROMACS allows a quite generic implementation of HYBRID particle-based models through using VIRTUAL sites for the particles at the coarser level in addition to the normal atoms at finer level, and possibly adding user-defined interactions between the particles in accordance with the appropriate expression for the total energy.

The correspondence between the atoms and the virtual sites is called the MAPPING of the atomistic to the CG model. A CG model by itself may use such mapping schemes to parameterize (parts of) the force field. In a purely hierarchical scheme, such as force matching, atomistic simulations are used to calculate the forces on the mapped centers, which are then averaged over many snapshots to define the interactions between the CG particles. In an empirical approach, such as the Martini model, the bonded interactions between CG particles are often refined based on atomistic simulations; the distributions of bond lengths, bond angles and torsional degrees of freedom at the CG level are compared to those from the mapped atomistic simulations and bonded parameters are adapted to get an overall reasonable agreement. Tutorials involving fine-tuning can be found here and in the lipid and polymer tutorials. The non-bonded interactions in the Martini model are nevertheless generic and based on parameterization to experimental data.

Hybrid OPLS-AA/L—Martini model

The material presented in this tutorial leans on the work originating in the Molecular Dynamics Group at the University of Groningen, in particular that by Tsjerk Wassenaar. The two main publications primarily concern the combination of GROMOS united atom and Martini force fields. Here, we show that one can combine OPLS-AA/L and Martini in the same fashion.

GROMACS implements a considerable number of standard interaction potentials, both bonded and non-bonded. It also enables the users to define their own interaction potentials. In addition, the code implements the generation and use of interaction centers (called VIRTUAL sites) whose positions depend in some geometrically well-defined way on the positions of two or more other particles. A hybrid model combines molecular models at different levels of resolution. The different models may use different types of interaction potentials, and may therefore not be compatible with the same non-bonded (and bonded) functional forms, necessitating a more complex set-up of hybrid model simulations than simulations at a single level of resolution. Here, a set-up will be demonstrated for a combination of the OPLS-AA/L peptide model with the Martini coarse-grained model. The peptide will interact internally at the atomistic resolution, while it interacts with solvent at the coarse-grained (CG) level. The solvent interacts with itself only at CG level. Sections 1-3 take you through setting up and running such a hybrid system for a single AA peptide in Martini water.

In Section 4 a tool is introduced to generate fully atomistic solvent configurations from the hybrid simulation that can serve as starting points for production runs at all-atom (AA) level. Such procedures are known as back-mapping techniques.

Section 5 uses the back-mapping tool to generate coordinates for a hybrid simulation from a fully coarse-grained snapshot. Here, a collection of peptides is back-mapped from the Martini model to the OPLS-AA/L model, while the water and coarse-grained ions are maintained at the Martini level.

All files are provided in the tar-ball hybrid-tutorial.tgz, which expands to the main directory hybrid-tutorial/OPLS-MARTINI. Paths will be given with respect to this directory. A directory with all the results is provided in the tar-ball hybrid-tutorial-worked.tgz, under hybrid-tutorial/WORKED.

1. The Models

A. Atomistic model: OPLS-AA/L

As atomistic model we will use the OPLS-AA/L force field.[1] Here, we will have a peptide interacting internally through this force field. The standard atomistic model can be built using the GROMACS tool pdb2gmx. Here, we will use a small tripeptide, TYR-VAL-TYR, with methylated terminal ends. A Protein Data Bank (PDB) file can be built very simply using e.g. the builder of pymol, or by finding an existing protein in the PDB with this sequence of residues and cutting those out. We wish to use neutral methylated end-groups, normally known in the PDB as ACE for the C-terminal end and NMA for the N-terminal end. To make this work with the standard library file for OPLS-AA/L in GROMACS, the residue name of the N-terminal end must be changed manually from NMA to NAC.

Hands-on

Go to the directory OPLSAA. The file YVY.pdb (Tyr-Val-Tyr), built using pymol, is available for you, and it incorporates the naming of the N-terminal end (NAC instead of NME, done by manual editing) to comply with the GROMACS implementation of the OPLS-AA/L force field. Build the atomistic topology file:

> gmx pdb2gmx -f YVY.pdb -ter -ignh

This command generates an OPLS-AA peptide topology, topol.top, after ENTERING the correct numbers (which may depend on the exact GROMACS version) for the options OPLS-AA/L force field, TIP3P for the water model, and None for both terminal ends. Check that the file topol.top produced makes sense.

B. Coarse-grained solvent model: Martini

As a solvent, the Martini model for water and hydrated ions will be used.[2] This model will also deal with the interactions between the peptide and the solvent and hydrated ions (see Section 2 for further details). Standard Martini files can be downloaded (but that is not necessary here because they are provided for you) from the Martini website. They cannot be used as such, however, in the hybrid set-up, and modifications will need to be made. The standard file for the Martini model which serves as the basis for the modifications can be found in the directory SOLVENT/martini_v2.2.itp.

2. The Combination of the Models

Models at different levels can be combined in a number of schemes. We describe the simplest approach, viz. we choose to administer all interactions as being at either AA or CG level. The particles making up the peptide will interact at both atomistic and coarse-grained levels. Within the peptide, the interactions should be treated according to the atomistic force field. The atoms of the peptide see the other atoms of the peptide as they would in a normal atomistic simulation. Interactions of the peptide with the solvent are to be CG interactions, i.e. treated entirely according to the Martini model. This clear-cut separation is possible for the system defined here because there are no Coulomb interactions between the AA and CG subsystems. Standard Martini water beads are single, neutral particles. (The proper treatment of electrostatic interactions across the two levels is by no means a finished area of research. Aspects are discussed in Section 5.)

The present set-up achieves the separation between AA and CG interactions by adding VIRTUAL sites to the atomistic topology, and making sure that these sites have no interaction with the atoms, nor with each other, but only with the solvent. The atomistic topology generated for the peptide must be extended with the virtual sites to build a hybrid topology for the peptide.

A. HYBRID TOPOLOGY FOR THE PEPTIDE

The hybrid topology for the peptide makes use of VIRTUAL SITES. These are described in the GROMACS manual in Section 4.7. In the present implementation, the virtual sites are defined as the centers of mass of a number of atoms. The virtual sites themselves are also particles belonging to the molecule. Thus, two entries need to be ADDED to the atomistic topology. The first extension is either an extension of the [ atoms ] directive containing the atomistic particles or a second [ atoms ] directive placed somewhere below the first one. In both cases, the atom numbering must be consecutive to the atoms already defined. The second extension is a [ virtual_sitesn ] directive, which defines how the added virtual sites are related to the normal atoms of the molecule. The additional virtual sites and the mapping can be added by hand, but a script is also available to do this for you. The conf.gro file of the peptide generated by pdb2gmx can be used to generate a Martini topology using an adapted martinize.py script (available in HYBRIDTOPOLOGY — the original script is downloadable from the Martini web-site) after slight modification of the terminal methylated groups. The standard Martini model does not implement the neutral methylated terminal ends. To distinguish them properly in the martinize.py script, a name change of the methyl carbon is required.

Hands-on

Go to the directory HYBRIDTOPOLOGY. Copy the standard pdb2gmx output file conf.gro produced in Section 1 to conf-martini.gro and edit the CH3 of the ACE residue to be named CTJ and the CH3 of the NAC residue to be named CTB. This change is for the benefit of making the adapted martinize-ter.py script recognize these neutral methylated termini.

> cp ../OPLSAA/conf.gro ./conf-martini.gro

> [vi/gedit] conf-martini.gro

Edit the conf-martini.gro to look like this:

Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
   70
    1ACE    CTJ    1  -1.244   0.089  -0.561
    1ACE   HH31    2  -1.272   0.187  -0.524
    1ACE   HH32    3  -1.199   0.101  -0.659
    1ACE   HH33    4  -1.333   0.029  -0.571
    1ACE     C     5  -1.146   0.024  -0.467
    1ACE     O     6  -1.176   0.001  -0.349
...
    5NAC     N    55  -0.107  -0.003   0.013
    5NAC     H    56  -0.153   0.028   0.100
    5NAC   CTB    57   0.034   0.012  -0.021
    5NAC   HH31   58   0.062  -0.063  -0.094
    5NAC   HH32   59   0.052   0.110  -0.062
    5NAC   HH33   60   0.094  -0.001   0.068
   1.42700   0.79900   0.81020

Now, execute the martinize-ter.py script (this is the adapted version):

> python martinize-ter.py -f conf-martini.gro -ff martini22 -multi all -nt -o martini.top -x hybrid.pdb

This generates an extension to the molecule definition of the tri-peptide, defining 12 new atom entries that are defined as virtual sites. The flag -ff selects the Martini force field version, the flag -multi all produces a topology extension for hybrid AA-CG simulations (omitting this flag, you get a topology for pure Martini CG simulations), the flag -nt ensures neutral termini, and the -x flag produces a starting structure including the virtual sites. The topology extension is written to the standard file Protein_A.itp, which should look like the text shown below:

; MARTINI (martini22) Multiscale virtual sites topology section for "Protein_A"

; Sequence:
; JYVYB
; Secondary Structure:
; CCCCC
[ atoms ]
   71   vN0     1   ACE   vBB   71 0.0000 ; C
   72   vP5     2   TYR   vBB   72 0.0000 ; C
   73  vSC4     2   TYR  vSC1   73 0.0000 ; C
   74  vSC4     2   TYR  vSC2   74 0.0000 ; C
   75  vSP1     2   TYR  vSC3   75 0.0000 ; C
   76   vP5     3   VAL   vBB   76 0.0000 ; C
   77  vAC2     3   VAL  vSC1   77 0.0000 ; C
   78   vP5     4   TYR   vBB   78 0.0000 ; C
   79  vSC4     4   TYR  vSC1   79 0.0000 ; C
   80  vSC4     4   TYR  vSC2   80 0.0000 ; C
   81  vSP1     4   TYR  vSC3   81 0.0000 ; C
   82   vN0     5   NAC   vBB   82 0.0000 ; C
;
; Coarse grained to atomistic mapping
;
#define mapping virtual_sitesn
[ mapping ]
   71     2    1    2    3    4    5    6
   72     2    7    8    9   10   26   27
   73     2   11   14   15   16
   74     2   17   18   21   22
   75     2   19   20   23   24   25
   76     2   28   29   30   31   42   43
   77     2   32   33   34   35   36   37   38   39   40   41
   78     2   44   45   46   47   63   64
   79     2   48   51   52   53
   80     2   54   55   58   59
   81     2   56   57   60   61   62
   82     2   65   66   67   68   69   70

Each virtual site is defined as the center of mass of a number of atoms of the tri-peptide. The virtual sites are given first, under the [ atoms ] directive, and are numbered starting from 71; the all-atom model tri-peptide has 70 atoms. The atom types given to the virtual sites are those of the Martini model for proteins and peptides, with a prefix v (e.g., the standard back-bone bead for a coil residue in Martini is P5; therefore the bead-type for the back-bone beads 72, 76, and 78 are vP5). The mapping is defined in the directive [ mapping ], which is made to be interpreted by grompp as a directive [ virtual_sitesn ]. The first number is the index number of the virtual site, the second number (2 for all sites) specifies the type of virtual site (center of mass), and the other numbers are the indices of the atoms that define the virtual site. You may cross-check the mapping with the atom definition given in the atomistic topology.

To weave the two models together, copy the topol.top obtained from the atomistic model to hybrid.top.

> cp ../OPLSAA/topol.top ./hybrid.top

> [vi/gedit] hybrid.top

Use an editor to add the Martini addition in Protein_A.itp at the appropriate place, which is just after the definition of the peptide ends. Excerpts of the file hybrid.top are shown below:

;     Command line:
;     gmx pdb2gmx -f YVY.pdb -ter -ignh
;     Force field was read from the standard GROMACS share directory.
;
; Include forcefield parameters
#include "oplsaa.ff/forcefield.itp"
[ moleculetype ]
; Name           nrexcl
Protein            3
[ atoms ]
;   nr       type resnr residue atom   cgnr     charge       mass typeB   chargeB     massB
; residue 20 ACE rtp ACE q 0.0
     1   opls_135   1   ACE    CH3    1      -0.18    12.011   ; qtot -0.18
     2   opls_140   1   ACE   HH31    1       0.06     1.008   ; qtot -0.12
...
    69   opls_140   5   NAC   HH32   25       0.06     1.008   ; qtot -0.06
    70   opls_140   5   NAC   HH33   25       0.06     1.008   ; qtot 0
[ bonds ]
...
[ pairs ]
...
[ angles ]
...
   69   67   70     1
[ dihedrals ]
...
...
   63   67   65   66     1   improper_Z_N_X_Y
; MARTINI (martini22) Multiscale virtual sites topology section for "Protein_A"
; Sequence:
; JVVVB
; Secondary Structure:
; CCCCC
[ atoms ]
   71   vN0     1   ACE   vBB   71 0.0000 ; C
...
   82   vN0     5   NAC   vBB   82 0.0000 ; C
;
; Coarse grained to atomistic mapping
;
#define mapping virtual_sitesn
[ mapping ]
   71     2     1   2    3    4    5    6
...
   82     2   65   66   67   68   69   70
; Include Position restraint file
...
[ system ]
; Name
Protein
[ molecules ]
; Compound       #mols
Protein             1

B. HYBRID FORCE FIELD FILES

We now have a topology for the hybrid peptide model. In general, further modifications are required before we can run a simulation.

  1. First of all, the virtual sites are assigned atom-types that belong to the CG model that are not known in the atomistic force field. They must be added to the atom-type list. Of course, the atom-type definitions must not overlap between the two models. For the combination of OPLS-AA/L and Martini this is no problem, but in general, you may need to rename the atom-types of either force field, whichever is more convenient.

  2. Also, interactions must be defined between the sets of the two different force fields: either they must be set to zero, or defined explicitly if required. In GROMACS, standard files for interactions are provided for a number of force fields in libraries, provided in directories, see the GROMACS manual, Section 5.8. The OPLS-AA/L force field definitions can be found in the directory YOUR-PATH-TO-GROMACS/share/gromacs/top/oplsaa.ff. We will need to adapt the file ffnonbonded.itp in this directory. Since we may not have administrative privileges, we will copy the entire directory to a place that is under our control, effect the required changes there and make sure those files are used in our simulations.

Hands-on Go to the directory FORCEFIELDS, and copy the standard library directory (or use the one already copied for you there):

> cp -R YOUR-PATH-TO-GROMACS/share/gromacs/top/oplsaa.ff .

> cd oplsaa.ff

Valid atom-types for any force field are defined in the file ffnonbonded.itp under the directive [ atomtypes ]. In the GROMACS implementation of the OPLS-AA/L force field there is a long list of atom types. For each atom type, the Lennard-Jones sigma and epsilon values are given in the 7th and 8th column of each entry. An excerpt of the file ffnonbonded.itp is given below to illustrate the point.

[ atomtypes ]
; full atom descriptions are available in ffoplsaa.atp
; name bond_type   mass   charge   ptype     sigma     epsilon
opls_001   C   6     12.01100     0.500   A   3.75000e-01 4.39320e-01 ; SIG
opls_002   O   8     15.99940    -0.500   A   2.96000e-01 8.78640e-01 ; SIG
...

The nonbonded interaction parameters can be calculated from the values entered using standard combination rules. GROMACS implements and geometric and arithmetic means, specified in GROMACS in the [ defaults ] directive, which can be found in the forcefield.itp file. In addition, the interaction parameters for pairs of types can be set explicitly, overriding the standard values. This is done in the [ nonbond_params ] directive in the file ffnonbonded.itp. Since OPLS-AA/L uses standard combination rules, the [ nonbond_params ] directive is missing from the file ffnonbonded.itp. You may skip the OPTIONAL part below if you like and continue at the end of the OPTIONAL part.

OPTIONAL (intermediate step)

Just adding the Martini types used for the virtual sites as valid OPLS-AA atom types (make sure you define them as type V, not A) with standard LJ parameters 0 and 0, allows you to do a simulation in vacuum with the virtual sites present, but not interacting with anything. You may attempt to add the virtual site types yourself by hand (you will need vN0, vP5, vSC4, vSP1, and vAC2), or copy the pre-prepared force field files in which this edit has been done for you.

Hands-on

Go to the VACUUM directory and either copy the given files:

cp -R ../FORCEFIELDS/hyb-vac.ff oplsaa.ff

or copy the standard files and edit the file ffnonbonded.itp to include the atom types for the virtual sites as described in the text above:

cp -R gromacs_path/share/gromacs/top/oplsaa.ff .
[vi/gedit] oplsaa.ff/ffnonbonded.itp

Next, copy your hybrid starting structure and edited topology there, and you should be able to run a vacuum simulation using the provided set-up file sd.mdp (using straight cut-off Coulomb interactions instead of PME and without pressure correction):

cp ../HYBRIDTOPOLOGY/hybrid.pdb .
cp ../HYBRIDTOPOLOGY/hybrid.top .
gmx grompp -p hybrid.top -c hybrid.pdb -f sd.mdp -maxwarn 1
gmx mdrun -v
vmd -e viz-sdrun.vmd

When visualizing with vmd, using the .vmd file, you will see the virtual sites, the positions of the Martini beads as green spheres. Due to the standard settings of vmd drawing bonds on the basis of distances, bonds will be drawn between the virtual sites and nearby atoms.

END OPTIONAL

The MARTINI force field has interaction levels that do not follow from a multiplication rule; it therefore defines a full matrix of interactions. The interaction matrix between Martini beads must be added to the nonbonded definition file in the [ nonbond_params ] directive of the file ffnonbonded.itp. Interactions must be defined between the sets of the two different force fields: either they must be set to zero, or defined explicitly if required. This can be done by hand by appropriately mixing the two files and extending the table with pairs of interactions. Finally, in general one may not want all molecules of a particular kind to interact at a particular level, but only a limited sub-set. For example, we may want to study the active site of a protein atomistically, but embedded in a CG protein and CG water. The virtual sites representing the CG active site must not interact with each other (the atomistic interactions take care of that), but should interact with all surrounding molecules in the appropriate CG manner. One could generate a list of exclusions for the virtual sites while using their normal CG atom types. A more general approach is to double the type of interactions, here CG, where a distinction is made between the types that do not interact with each other, here with prepended v, e.g. vP5 for a P5-type particle, and the normal CG types. Now v-particles can interact with all CG particles in the normal CG manner, whereas they do not interact at all amongst themselves if the LJ parameters between all v-particles are set to zero. This approach is followed here, and you may try to edit an appropriate file yourself, or use a minimal nonbonded file defined for you in FORCEFIELDS/hybrid.ff/ffnonbonded.itp.

You may skip the OPTIONAL part below if you like and use the pre-prepared force field files in which this edit has been done for you (FORCEFIELDS/hybrid.ff/ffnonbonded.itp) and continue at the end of the OPTIONAL part.

OPTIONAL (defining the set of CG and virtual particle interactions)

Hands-on

Go to the FORCEFIELDS directory and make a copy of the OPLS/AA force field files:

cp -R oplsaa.ff myhybrid.ff

Open the file myhybrid.ff/ffnonbonded.itp and add the Martini types used for the coarse-grained system and for the virtual sites as valid OPLS-AA atom types (make sure you define the ones for the virtual sites as type V, not A) with standard LJ parameters 0 and 0 (you will need P4, BP4, vN0, vP5, vSC4, vSP1, and vAC2). If you have gone through the previous OPTIONAL HANDS-ON, you will already have defined the virtual site atom types and you can use that again.

[gedit/vi] myhybrid.ff

Next, while still in the editor, add a directive [ nonbond_params ] and specify explicitly the non-bonded interaction matrix for interactions that should be non-zero. If a pair is not specified in the [ nonbond_params ] directive, standard combination rules will be used; thus by defining a particle type with LJ parameters 0 and 0 this type will always have no LJ interaction with any other particle type if the pair interaction is not explicitly defined in the [ nonbond_params ] directive. An example of a [ nonbond_params ] directive entry is shown below.

[ nonbond_params ]
; type1 type2  flag  sigma     epsilon
P4   P4   1        0.47e-00 0.50e+01 ; supra atrractive Martini interaction

When you are done, you may compare your file myhybrid.ff/ffnonbonded.itp with the file hybrid.ff/ffnonbonded.itp.

END OPTIONAL

NOTE that the standard Martini force field file available from the Martini website (and provided in the SOLVENT directory) uses C6 and C12 parameters instead of sigma and epsilon, as appropriate for the OPLS-AA/L force field. Thus, the Martini parameters must now also be specified in terms of sigma and epsilon (which is the way they are defined in Martini papers).

3. Setting up and running the hybrid system

We will now work toward a hybrid simulation of a single atomistic tripeptide solvated in Martini water.

A. Surrounding the solute with solvent

First, a starting structure must be prepared in which the peptide is solvated. GROMACS has a standard tool to solvate systems, gmx solvate. One can either use an equilibrated solvent box to fill an existing system containing a solute, or insert individual molecules one by one. After one or more molecules are inserted, a check for close contacts is made and new molecules that are too close to existing molecules are removed.

Hands-on

Go to the HYBRID directory. Take the hybrid solute co-ordinates, center them in the box and define a box-size of 3.5x3.5x3.5 nm, and fill this with Martini water beads. A co-ordinate file containing 400 equilibrated Martini waters (water.gro) is available from the Martini web-site, and is also provided for you in the SOLVENT directory.

> cp ../SOLVENT/water.gro .

> gmx editconf -f ../HYBRIDTOPOLOGY/hybrid.pdb -center 1 1 1 -box 3.5 3.5 3.5

> gmx solvate -cp out.gro -cs water.gro -o solvated.gro -radius 0.2

> cp ../HYBRIDTOPOLOGY/hybrid.top .

Edit the topology file (hybrid.top) to include the definition of the Martini water bead (define a moleculetype or use a #include; the file SOLVENT/MartiniWater.itp contains it) and specify the number of water beads just added to your system (right at the bottom in the [ molecules ] directive; the Martini water bead is called W). The number of water beads added should be clear from the output of the gmx solvate command.

> [vi/gedit] hybrid.top

Next, either copy the directory with pre-edited OPLS-AA/L force field files:

> cp -R ../FORCEFIELDS/hybrid.ff oplsaa.ff

or copy the directory with the standard OPLS-AA/L force field files and edit them to include the Martini-bead definitions (see the OPTIONAL hands-on section just prior to this section):

> cp -R gromacs_path/share/gromacs/top/oplsaa.ff .

> [vi/gedit] oplsaa.ff/ffnonbonded.itp

OPTIONAL (intermediate stage)

The system thus defined will run; you can test this (in an energy minimization run) by running:

```bash gmx grompp -p hybrid.top -c solvated.gro -f em.mdp -o em -maxwarn 1

gmx mdrun -v -deffnm em

vmd -e viz-em.vmd

END OPTIONAL

B. Generating proper tabulated potentials

Just by making all CG atom types known to the OPLS-AA/L force field and adding the LJ-parameters, the simulation will run, but it is not a proper hybrid simulation because the two force fields do not use the same cut-offs for the nonbonded interactions, nor — in the case of OPLS-AA/L and Martini — the same functional form for the Coulomb and LJ interactions. If you would run an MD simulation with the present set-up, the density of the system is incorrect because with standard LJ potentials you would not get the proper density of Martini water.

The next ingredient is thus making sure that the interactions are treated properly according to the appropriate force fields. This may or may not be entirely possible within GROMACS. A powerful tool that enables multiscale simulations are the so-called “tables” that GROMACS implements. Internally, potentials and forces are tabulated within GROMACS and the forces are determined by an interpolation scheme between the tabulated values. This is hidden from the user, but users may define their own interaction tables, by which they can input any pairwise nonbonded interaction potential. The next section briefly takes you through the setting up of such tables. They are described in the GROMACS manual, Section 6.9 and a more extensive tutorial is available from the GROMACS website.

Nonbonded interaction tables can be defined for the interactions between any number of groups. Here, we need the following set, which are provided for you in the FORCEFIELDS directory:

  1. interaction tables between the atomistic particles, which will be called AA:
  • the OPLS-AA/L force field with nonbonded cut-off of 1.0 nm, using a standard LJ potential contained in the file tables_AA_AA.xvg,
  1. interaction tables between coarse-grained particles, which will be called CG:
  • the MARTINI force field with nonbonded cut-off of 1.2 nm, using a modified LJ potential, contained in the file tables_CG_CG.xvg, and
  1. interaction tables between the two types (AA and CG) of particles:
  • for the example here, the potentials can be set to zero for all distances, but a file with these values must be present. We use the naming convention of GROMACS and let this interaction be the default, which means these interactions must be contained in the file table.xvg.

Each tables file contains 7 columns, the first column contains the distances between the pair of particles (in nanometer), and the next six columns implement three generic functions in three pairs of potential and force. Within GROMACS, the magnitude of the interaction or force is determined by the interpolated value taken from the tables file, multiplied by the appropriate Coulomb or LJ prefactors for the pair at hand. The structure of the tables files is thus:

distance   f(r)       f’(r)     g(r)       g’(r)       h(r)         h’(r)

and with heavily truncated numbers could start like this:

0.0000E+00 0.0000E+00 0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00 0.0000E+00
0.2000E-02 0.4999E+03 0.2499E+06 -0.1562E+17 -0.4687E+20 -0.2441E+33 0.1464E+37
....

For the first interaction (modeled on the Coulomb interaction), the prefactor for the functions f(r) and f’(r) is calculated by multiplication of the two charges and dividing by the dielectric constant (epsilon_r in the mdp-file). For the second (g) and third (h) functions, (modeled on the LJ interaction), the prefactors are calculated from the standard combination rules for the C6 and C12 parameters, respectively (which may themselves be calculated from the input sigma and epsilon values), or taken directly from the [ nonbonded_params ] list in the ffnonbonded.itp file. The distance spacing of the table files can be chosen at will but needs to be uniform. The tables should extend to at least the longest cut-off used in the interactions between the different types of particles, here 1.2 nm, but it is useful to make them extend beyond this — by default, tables are extended 1 nm beyond the cut-off (mdp option table-extension).

GROMACS provides a number of tables files for you, for example for the Lennard-Jones 6-9 potential, table6-9.xvg. They are provided in the YOUR-PATH-TO-GROMACS/share/gromacs/top directory, with an extension up to 3.0 nm. You will generally need to generate these files yourself, writing a script, a program, etc. Force-matching or related hierarchical coarse-grained force fields[3] generally use tabulated potentials derived from atomistic simulations. For the purpose of this tutorial, the tables files are provided for you in the FORCEFIELDS directory and you can simply copy them over. Technical information about the functional forms of the potentials used here is given in the appendix to a paper by Baron et al.[4]. NOTE that you also need a separate tablep.xvg file for 1-4 (pair) interactions.

OPTIONAL Hands-on

Go to the FORCEFIELDS directory. The tabulated potentials can be viewed using xmgrace, for example:

xmgrace -nxy table_CG_CG.xvg -p xmgr.par

Return to the HYBRID directory.

END OPTIONAL

A basic mdp-file for energy minimization in the hybrid scheme using user-defined potentials is shown below (and provided for you as min.mdp):

; RUN CONTROL PARAMETERS
integrator               = steep
nsteps                   = 100
; OUTPUT CONTROL OPTIONS
; Output frequency for coords (x), velocities (v) and forces (f)
nstxout                 = 1
nstvout                 = 1
nstenergy               = 1
; Selection of energy groups
energygrps              = AA CG
; NEIGHBORSEARCHING PARAMETERS
; Periodic boundary conditions: xyz, no, xy
pbc                     = xyz
; nblist cut-off      
rlist                   = 1.0
; OPTIONS FOR ELECTROSTATICS AND VDW
; Method for doing electrostatics
coulombtype             = user
rcoulomb                = 1.2
; Relative dielectric constant for the medium and the reaction field
epsilon_r               = 1
; Method for doing Van der Waals
vdwtype                 = user
rvdw                    = 1.2
; Seperate tables between energy group pairs
energygrp_table         = AA AA CG CG

Compared to standard settings, points to note in the hybrid set-up are:

  1. definition of the energy groups AA and CG; these are required because their interactions are treated differently, as specified in the energygroup table (see (iii) below):
energygrps = AA CG

defines two energy groups, which must be present in the index file (run.ndx) supplied by the user (see below);

  1. specification of user-defined potentials for both coulomb and vanderwaals interactions:
coulombtype = user

vdwtype = user

tells grompp to look for table files;

  1. designation of which user-defined table is used to handle which interaction:
energygrp_table = AA AA CG CG

tells grompp to look for table_AA_AA.xvg (first pair of arguments, AA AA) and table_CG_CG.xvg (second pair of arguments, CG CG) that define the interactions between the energy groups AA with AA and CG with CG, respectively. For all other interactions the data in the file table.xvg will be used.

The separate treatment of interactions is taken care of by defining energy groups as explained above. The energy groups must be known to the program and you must provide an index file to that end, because the energy groups needed are not standard groups known to GROMACS. You can use the conformation with added water beads to build the index file. The sequence shown below does that for you (lines preceded by > are commands within make_ndx; the numbering used in this sequence presumes that the standard groups that are defined number up to 13, where 13 are the W beads, and the group Protein is group number 1):

Hands-on (make sure you are in the HYBRID directory)

> gmx make_ndx -f solvated.gro

(the next series of commands is typed within make_ndx)

> a v*

> 13 | 14

> 1 &! 14

> name 15 CG

> name 16 AA

> q

(quits make_ndx and puts you back on the Linux command line)

> mv index.ndx run.ndx

The next sequence copies the tables files, does a minimization and short MD run in the OPLS-AA/L-Martini hybrid set-up.

> cp ../FORCEFIELDS/table* .

> gmx grompp -p hybrid.top -c solvated.gro -f min.mdp -n run.ndx -maxwarn 3

> gmx mdrun -v

> vmd -e viz-min.vmd

At room temperature, Martini water may freeze, which may cause problems. Water may be prevented from freezing by treating around 10% of the water beads as so-called “anti-freeze water”. You may do this quite simply by editing the topology file, hybrid.top. Add a line after the water beads, specifying WF nnn, where nnn is about 10% of the number of W beads. Reduce the number of W beads accordingly. You do not need to change the coordinates file. You must make sure the particle type BP4 is specified in the force field file ffnonbonded.itp.

> [vi/gedit] hybrid.top

> [vi/gedit] oplsaa.ff/ffnonbonded.itp

> gmx grompp -p hybrid.top -c confout.gro -f md.mdp -n run.ndx -maxwarn 1

> gmx mdrun -v

> gmx trjconv -center -pbc mol -o viz.xtc < convin

> vmd -e viz-run.vmd

Hopefully by now, the system is stable enough to run for a longer period of time and you can start production. OPLS-AA/L uses PME by default, and a possible set-up is given in pme.mdp. In my own experience, using constraints leads to problems, and the time-step is fairly limited. Still, there is a substantial reduction in the number of solvent particles, and the overall efficiency in sampling the peptide degrees of freedom is higher than with, say, the TIP3P model as solvent. As a reminder, these types of simulations are at an early stage of development and it is still very much open to debate what methodology is the best one to yield optimal reliability and efficiency. Feel free to try different integrators, settings for temperature and pressure coupling and changes to the topology. Treating the H-atoms as virtual sites and/or make them heavy may enhance the numerical stability of the system.

To my knowledge, the combination of OPLS-AA/L and Martini models has not been reported in the literature. Our group in Groningen has published on combining GROMOS united-atom models with Martini models.[5, 6] It is clear from that work, that the combination of different resolutions is not without artifacts, and especially the treatment of electrostatics needs careful attention. An attempt to incorporate this is discussed in Section 5.

4. Fine-graining a hybrid snapshot

Hybrid models promise the best of two worlds: efficient sampling of complex systems and high-level detail in the region of interest. It is nevertheless likely that in practice, some trade-off between accuracy and computational speed will be made causing artifacts near the boundary(ies) of the two (or more) levels. The boundary artifacts may reflect on the region of interest and a method to put in more detail based on the less detailed model can help to study such artifacts. Alternatively, such a method should provide a good starting structure for a simulation at the more detailed level. Approaches that sample long-term processes at coarse-grained level and investigate the details of the interesting intermediate structures after fine-graining have been described already in the literature.

Here, we will take a snapshot from the hybrid simulation of an OPLS-AA/L peptide in Martini water discussed in the previous section and put in TIP3P waters based on the positions of the CG water beads. The method also allows (partial) fine-graining purely CG snapshots, which we will do for a system of 100 tripeptides in Martini water. Several back-mapping schemes have been published. The present method, developed by Tsjerk Wassenaar, is in my view a flexible and user-friendly one.[7]

Back-mapping or reverse coarse-graining or fine-graining a (partially) coarse-grained structure requires a correspondence between the two models; specifically for AA and CG: which atoms make up which bead. Actually, an atom can in principle contribute to several beads. The hybrid set-up described in this tutorial defines virtual sites (the positions of the CG beads) as the centers of mass of a number of atoms. In this mapping scheme, each atom appears only once, but there is no technical reason why it could not appear in the definition of several virtual sites. A back-mapping protocol needs to know at least which atoms contribute to which bead. Existing schemes then use rigid building blocks anchored on the CG bead, or place the atoms randomly near the bead in an initial guess and the structure is relaxed based on the atomistic force field, usually by switching it on gradually. The method presented here allows for an intelligent, yet flexible initial placement of the atoms based on the positions of several beads, thereby automatically generating a very reasonable orientation of the groups of atoms with respect to each other. When back-mapping the hybrid model used in this tutorial we do not need this, because only water beads are fine-grained. When partially back-mapping only the peptides in the CG model to generate a hybrid model, we do. The more general procedure is described in a separate tutorial on the Martini website.

Back-mapping Martini water to TIP3P (or SPC) water

A Martini water bead models four atomistic water molecules. The backward.py script used here to fine-grain the hybrid snapshot generates four water molecules around a water bead. By default, any particle named W in the structure file is interpreted as a water bead and the script itself generates four 3-atom water molecules (which may then be used as SPC or TIP3P waters). For this exercise, a fully fine-grained starting structure can be prepared with a little editing.

Hands-on Go to the directory BACKMAPPING.

  1. make two copies of a hybrid snapshot
> cp ../HYBRID/confout.gro W.gro

> cp ../HYBRID/confout.gro peptide.gro
  1. remove all non-water (i.e. those that are not W or not WF) beads from the file W.gro, and edit the second line of the file to correctly state the number of atoms/beads in the file (number of waters added in section 3A minus 82). Make sure that if there is anti-freeze water, the bead is renamed from WF to W. You can do this within the editor, or after the edit as shown here. Alternatively, implement the backmapping of anti-freeze water in backward, where it is hard-coded.
> [vi/gedit] W.gro

> sed -i 's/WF/W /g' W.gro
  1. apply the back-mapping script backward.py to the edited file
> python backward.py -f W.gro -sol -kick 0.2 -o fgW.gro

The flag -sol tells the script to convert W beads into (standard 3-particle water) solvent atoms (four molecules for one bead), and the flag -kick 0.2 applies random displacements to the initially placed atoms, reducing the strong orientational correlations of the initial placement.

  1. remove the header (first two lines) and all virtual sites and W and WF beads and the last line from the file peptide.gro; this file now only contains the real atoms of the peptide
> [vi/gedit] peptide.gro
  1. insert the peptide atoms in the file containing the atomistic water, fgW.gro; you should add them near the top of the .gro file, just after the line stating the number of atoms in the file and before the solvent molecules. Also edit the number of atoms in the file (add 70, this is the number of atoms in the peptide).
> cp fgW.gro fg.gro

> [vi/gedit] fg.gro
  1. copy the atomistic topology file and edit it to add a SOL entry in the [ molecules ] directive. The number of SOL molecules is 4 times the number of W beads you had in the hybrid set-up.
> cp ../OPLSAA/topol.top .

> [vi/gedit] topol.top
  1. minimize the back-mapped structure and run a short relaxation/equilibration
> gmx grompp -p topol.top -c fg.gro -f min.mdp -o min

> gmx mdrun -v -deffnm min

> gmx grompp -p topol.top -c min.gro -f md.mdp -o md -maxwarn 1

> gmx mdrun -v -deffnm md

> gmx trjconv -center -pbc mol -s md.tpr -f md.xtc -o viz.xtc < convin

> vmd -e viz-run.vmd

Now you are all set to run production at fully fine-grained level.

5. An assembly of peptides including charged side chains and ions

In this part of the tutorial, we will address a more complex system in which electrostatic interaction between the FG and CG subsystems play a role. The aim of this part of the tutorial is to show how to set up hybrid systems. The best way of dealing with the interactions between CG and FG subsystems is by no means clear yet. In fact, the literature from our own group shows that proper treatment in the hybrid model is problematic. Therefore, this section should be regarded as an incentive for further research in this area.

In the context of high throughput screening simulations, all possible tripeptides were tested for their self assembly behavior using the Martini coarse-grained force field. Experimental testing subsequently confirmed that the type of self-assembly behavior of the peptides was predicted satisfactorily by the simulations. It may be interesting to have more detailed insight in the self assembly pricess and in the organization of the assemblies; this may be achieved by using one or more snapshots from the CG model as a basis for running either all-atom or hybrid simulations. Here, we will run a hybrid simulation of an assembly of 100 TYR-GLU-TYR (YEY) peptides in water that is the result of spontaneous aggregation with the Martini model. The spontaneous aggregation set-up is treated in the high throughput peptide screening tutorial.

Hands-on Go to the directory HYBRID-YEY-AGG. This will be the parent directory for hybrid simulations of 100 TYR-GLU-TYR (one-letter code abbreviation YEY) peptides in Martini water. The starting structure is taken from another Martini tutorial and is the result from a spontaneous self-assembly run of these peptides at the Martini level. The hybrid set-up aims at bringing in atomistic details and relaxation at the atomistic level of the aggregate formed thus far. A number of steps done for the TYR-VAL-TYR peptide need to be repeated and/or modified here. You may go through them all explicitly or use the results of some steps without doing them yourself.

A. Preparing the hybrid topology

You may follow the same procedure here for the YEY-peptide as for the YVY peptide earlier in this tutorial, or just take the topology that is generated in that way for you.

OPTIONAL Hands-on: preparing hybrid topology

Go to the directory HYBRID-YEY-AGG/HYBRIDTOPOLOGY. Use a structure editor, e.g. pymol to generate a pdb structure of a protected Ace-YEY-NMA peptide. Save as a PDB-file under the name YEY.pdb, edit it to remove the line “TER TYR”, change the residue name NME to NAC (see Hands-on in Section 1A) and from it, generate an OPLS-AA/L model for the peptide:

pymol

[vi/gedit] YEY.pdb

pdb2gmx -f YEY.pdb -ter -ignh

Edit the conf.gro file to comply with the naming of terminal methyl groups in the adapted martinize-ter.py script (see Section 2A), and use the script to generate the definitions of the virtual sites:

[vi/gedit] conf.gro

python martinize-ter.py -f conf.gro -ff martini22 -multi all -nt -o martini.top -x hybrid.pdb

Copy the AA topology and add the information for the virtual sites in Protein_A.itp to the AA topology.

cp topol.top myhybrid.top

[vi/gedit] myhybrid.top

END OPTIONAL

B. Partial back-mapping of the fully CG snapshot

We need to generate atomistic coordinates from the CG coordinates for the peptides only. This may be done using the backward.py script. However, since we are using virtual sites, we need to also make sure there are virtual site coordinates present. This is achieved here by introducing special mapping files for the first and last residues. The procedure is described in more detail in an appendix, which you may want to refer to later on. Also, we will need to use the polarizable Martini water model. In the polarizable Martini model, water consists of three particles, two of which bear opposite charges which are free to move on a circle around the central uncharged bead.[8] The LJ interaction resides on the central bead and the charges within the water bead do not have an electrostatic interaction, but there is an angle potential between the fixed-length bonds from the central particle to each of the charges. The model mimics the orientational polarization of a cluster of four water molecules. Each water bead from the fully CG model (that used single-bead Martini water) needs to be converted to three beads making up the polarizable water model. The script triple_w.py does this for you (see also here).

Hands-on Go to the directory HYBRID-YEY-AGG/PARTIALBACKMAP.

  1. convert the W beads in the fully CG snapshot provided to polarizable Martini water (each bead generates three particles). The result is written to the file TYR-GLU-TYR_eq_PW.gro:
> python triple_w.py TYR-GLU-TYR_eq.gro
  1. edit the CG snapshot with polarizable water beads to change residue names 1TYR to 1TYN and 3TYR to 3TYC. You can do this conveniently using two sed commands:
> cp TYR-GLU-TYR_eq_PW.gro tobackmap.gro

> sed -i 's/1TYR/1TYN/g' tobackmap.gro

> sed -i 's/3TYR/3TYC/g' tobackmap.gro
  1. apply the back-mapping script backward.py to the edited file
> python backward.py -f tobackmap.gro -from martini -to oplsaa -nt -o YEY-backmapped.gro

The file YEY-backmapped.gro should now contain atomistic coordinates for the 100 YEY peptides and coordinates for the Na-ions. The polarizable water beads present have not been backmapped to atomistic water and also they have not been kept or copied into the backmapped result. We must add them to the backmapped file to arrive at a system containing both the atomistic peptides and the coarse-grained water and ions. To this end, extract the polarizable water beads from the file TYR-GLU-TYR_eq_PW.gro.

> awk '{if (substr($2,1,1) == "W") print $0 }' < TYR-GLU-TYR_eq_PW.gro > W.gro
  1. edit the file YEY-backmapped.gro. Add the polarizable water beads below the Na-ions, and make sure the number of particles in the file is correctly stated on the second line of the file.
> [vi/gedit] YEY-backmapped.gro

For comparison, the file PartialBackmap.gro is provided for you.

C. Setting up the hybrid simulation

The treatment of the interactions is different from the example of the single peptide in water described in section 3B. In this case, electrostatic interactions between atomistic and coarse-grained system play a role. They are dealt with by using a different set of table-files. This is explained in more detail in an appendix, which you should consult once the simulation is up and running.

Hands-on Go to the directory HYBRID.

  1. get the partially back-mapped CG structure
> cp ../PARTIALBACKMAP/YEY-backmapped.gro system.gro
  1. copy and edit the file hybrid.top (you can use your own if you prepared myhybrid.top in the OPTIONAL part earlier on) to reflect the number and order of molecules (100 Protein molecules, 100 NA+ beads and 2,936 polarizable water molecules (PW). Put in #include statements to the Martini polarizable water (../../SOLVENT/PolMartiniWater.itp) and ion topologies (../../SOLVENT/MartiniSodium.itp).
> cp ../HYBRIDTOPOLOGY/hybrid.top .

> [vi/gedit] hybrid.top
  1. make sure that all the CG and virtual site types and their interactions are defined properly in the force field files and that the user-defined potentials (tables) are present
> cp -R ../../FORCEFIELDS/hybrid-charged.ff oplsaa.ff

> cp ../../FORCEFIELDS/table* .

> cp ../../FORCEFIELDS/gro-table_AA_AA.xvg table_AA_AA.xvg

> cp ../../FORCEFIELDS/gro-table_AA_CG.xvg table_AA_CG.xvg

> cp ../../FORCEFIELDS/mar-table_CG_CG.xvg table_CG_CG.xvg

> cp ../../FORCEFIELDS/mar-table_CG_VS.xvg table_CG_VS.xvg
  1. minimize the hybrid system using an improper hybrid forcefield (have a look at the settings for the coulomb and vdw interactions). This will hopefully easily heal the sometimes rather haphazardly placed atoms; you may try to improve on this by writing a better mapping file, see the backmapping tutorial. (Note that the virtual sites are put on the centers of mass before anything else is done by GROMACS itself, therefore their placement in the starting structure is not important as long as they are in the file!) Here, we use the soft-core potentials via the free-energy code to deal with atoms that are placed very close to each other and cause problems when all interactions are immediately switched on. Using the free-energy code slows down the calculations considerably, but we usually need only a few thousand steps to get a decent structure.
> gmx grompp -p hybrid.top -c system.gro -f em-sc.mdp -o em-sc.tpr -maxwarn 2

> gmx mdrun -deffnm em-sc -v

> gmx grompp -p hybrid.top -c em-sc.gro -f sd-sc.mdp -o sd-sc.tpr -maxwarn 2

> gmx mdrun -deffnm sd-sc -v
  1. make the index file including the definitions of AA, CG, and virtual site (VS) subsystems for the chosen treatment of the interactions (see Appendix D). Here, we assume that the standard groups run up to and including number 14 (PW is the final group, and the ions are group 13)
> gmx make_ndx -f sd-sc.gro

(the next series of commands is typed within make_ndx)

> a v*

> name 15 VS

> 1 & !15

> name 16 AA

> 0 & !15 & !16

> name 17 CG

> q

(you have now quit the make_ndx tool and your index file is saved as index.ndx)

> mv index.ndx run.ndx
  1. minimize again, now using the proper hybrid force field (have a look at the settings for the coulomb and vdw interactions and what table files are expected). If you are lucky, the following sequence will work. If not, you need to look at the details and try to find out what may cause problems. This may require quite a bit of grit, patience, and determination…
> gmx grompp -p hybrid.top -n run.ndx -c sd-sc.gro -f min.mdp -maxwarn 1

> gmx mdrun -v -o min.gro

> gmx grompp -p hybrid.top -n run.ndx -c min.gro -f sd1.mdp -maxwarn 1

> gmx mdrun -v -o sd1.gro

> gmx grompp -p hybrid.top -n run.ndx -c sd1.gro -f sd2.mdp -maxwarn 1

> gmx mdrun -v -o sd2.gro

> gmx grompp -p hybrid.top -n run.ndx -c sd2.gro -f md1.mdp -maxwarn 1

> gmx mdrun -v -o md1.gro

Note that it is inconvenient here to use the -deffnm option on mdrun because the table files should then have names corresponding to the name of the -deffnm option. If you want to keep the trajectories under specific names, move traj.xtc after each step to your name of choice.

Appendix: Special mapping files

In sequential multiscaling, resolution transformations are made between different levels of description: an atomistic snapshot is coarse-grained to run a CG model for computational efficiency, and when an interesting event or structure occurs during the CG simulation, the atomstic details are put back in to set up a more expensive simulation studying the details. A powerful tool to toggle between resolutions is provided by the backward procedure, implemented in python script backward.py. The package provides a set of mapping files, defining the correspondence for amino acids and selected other molecules (lipids mostly) between the Martini model and common united atom (gromos) and all atom (charmm, amber) models.

In the hybrid models described in this tutorial, peptides are described as partially interacting through atomistic force field interactions and partly through Martini-type interactions. This is achieved by adding virtual sites to the topology of the peptide. The virtual sites interact only with the coarse-grained environment. If we have a fully CG snapshot, we need to put in the atomistic details of the peptide molecule(s), INCLUDING the virtual sites. The standard mapping files do not provide this information. Here, we describe how to add the virtual sites in a simple manner. A second aspect in this tutorial is the fact that we use neutral (methylated) terminal residues ACE and NAC. There are no mapping files for these residues either. Both aspects can be dealt with simultaneously.

In the chosen hybrid set-up, all virtual sites are added after all atoms of the peptide. If we make a mapping file for the terminal resdue (here a Tyrosine), we can simply define as many extra atoms as we need virtual sites, and add them at the end of the normal (standard) atoms. In the mapping assignment, we can assign them all to the back-bone bead. They will all be placed near this bead in the backmapped structure but this should be no problem at all. The virtual sites are put on the centers of mass by GROMACS itself before anything else is done, e.g. energy minimization or the first dynamics step. All that is required is then to copy the mapping file for an existing amino acid and give a unique residue name to the terminal residue. Here, we chose TYC (for Tyrosine, C-terminal). In this new mapping file, we can also easily add atoms to account for the neutral terminal we created, here the C=O(Me) group. In this case, the Martini model did not use an extra bead for the extra atoms, but assumes these part of the final backbone bead (BB). Thus, before the extra atoms for the virtual sites, modify and add atoms for the methylated terminus.

Hands-on

Go to the directory HYBRID-YEY-AGG/PARTIALBACKMAP/Mapping. Copy the standard file for a tyrosine residue from martini to CHARMM to a new file that is to implement the special C-terminus. If you do not wish to override the file provided, choose a different name. Note that CHARMM and OPLS/AA use the same names but not the same numbering (sequence) of the atoms in amino acid residues.

> cp tyr.charmm36.map tyc.oplsaa.map

> [vi/gedit] tyc.oplsaa.map

In the editor, change the residue name from TYR to TYC. Change the entry in the [ mapping] directive to oplsaa. Change the order of the atoms according to the definition in the HYBRID-YEY-AGG/HYBRIDTOPOLOGY/hybrid.top file. At the end of the atom list, add extra atoms for the terminal methyl group (CH3, HH31, HH32, and HH33) first, then add atom definitions for all virtual sites of the entire peptide (think of suitable names yourself, or take them from the hybrid.top file). Assign all the new atoms and virtual sites to the BB bead in the mapping. You should have a total of 37 atoms for this residue.

The first residue, in our case also a tyrosine, is also extended with a methylated terminal end. This means, we now have to place atoms at the start of the amino acid. This then requires a renumbering of the other atoms. Copy the tyrosine residue mapping file once more, make edits to change the name of the residue to TYN, change the force field name, and change the order of the atoms in the aromatic ring and add the required atoms at the beginning of the [ atoms ] directive (CH3, HH31, HH32, HH33, and C1). The remaining (original) atoms must now be renumbered and you should end up with 27 atoms.

> cp tyr.charmm36.map tyn.oplsaa.map

> [vi/gedit] tyn.oplsaa.map

Finally, you will need a copy of the glutamate residue mapping file. Copy the one from charmm and change the name of the force field. The charmm mapping file has the acid group protonated. Here we want a charged side group. Delete the H-atom on the acid group and renumber the remaining atoms.

> cp glu.charmm36.map glu.oplsaa.map

> [vi/gedit] glu.oplsaa.map

If you have downloaded the backward package, you will also need to adapt the backward.py script. The downloaded version does not accept mapping to OPLS/AA, and it does not recognize the methylated termini. Compare the two scripts to see how these changes were implemented.

Appendix: Setting up hybrid electrostatics

The treatment of electrostatics in hybrid simulations is unlikely to be straightforward, especially if two force fields are combined that are not related to each other in a hierarchical scheme. The two types of force fields may not use the same functional form for the electrostatics. In atomistic models electrostatic potentials are usually based on unscreened Coulomb interactions, either cut off at some distance or treated by full Ewald-type electrostatics in which all charges interact with each other. The standard choice for OPLS-AA/L is Particle Mesh Ewald (PME), a numerical variant of Ewald electrostatics. Coarse-grained models may have no electrostatic term at all, but if they do, the electrostatic potential is usually of some screened form to account for the absence of orientational polarization by dipolar species. In the standard Martini model, the beads modeling water as a solvent do not bear any charge. Formal charges are given to the two beads making up the Zwitterionic lipid head groups, and to the bead(s) representing charged amino-acid side chains. In atomistic models the charge-charge interactions are screened by the explicit water reorientation which is absent in the CG model. Therefore, the Martini model uses a dielectric constant of 15 to reduce the interaction between the charged species. Also, the Martini model does not use Ewald-type electrostatics; instead the potential smoothly goes to zero at the cut-off distance. It should be noted that the standard Martini model for lipids does run stably under PME and that lipid bilayer properties are similar to those obtained with the standard screened and smoothed Coulomb potential. Thus, combining models may be possible if the electrostatics treatments of the two models can be brought to the same level in a practical sense, possibly with some adjustments of parameters. A published study of combining the GROMOS 53A6 force field with Martini[6] shows that it is not trivial to achieve a working model but at least a working starting model can be constructed.

This tutorial contains a possible way in which to combine the OPLS-AA/L force field with the polarizable water Martini model and Martini ions using a mixed model for electrostatics. This appendix explains in more detail how it is set up. In the polarizable Martini model, water consists of three particles, two of which bear opposite charges which are free to move on a circle around the central uncharged bead.[8] The LJ interaction resides on the central bead and the charges within the water bead do not have an electrostatic interaction, but there is an angle potential between the fixed-length bonds from the central particle to each of the charges. The model mimics the orientational polarization of a cluster of four water molecules. The model uses the same form of the Coulomb potential (shifted. i.e. with a smoothed cut-off at 1.2 nm, and no electrostatic interaction beyond that point), but with a relative dielectric constant of 2.5 instead of 15. In contrast, the interaction between partial charges of the atomistic OPLS-AA/L model use a relative dielectric constant of 1 because screening is explicit by TIP3P water molecules. In addition PME electrostatics is used in the atomistic model.

In the implementation used in this tutorial, PME electrostatics is not used (except in the model without charged peptides). Instead, the electrostatic interactions between the AA particles are treated with the Coulombic reaction field potential that is the standard in the GROMOS force fields. Here, we chose a cut-off of that potential of 1.2 nm (the same as for the Martini model; 1.4 nm is the standard cut-off for the GROMOS force fields), and a value of the dielectric constant of the surroundings, epsilon_rf in Gromacs settings files, of 78. There could be several possible other choices made for these settings. The point here is to show how you can set up these different choices. It is by no means claimed that the set-up given here is optimal! As far as I am aware, no satisfactory set-up has been published yet, so finding the proper way to treat this type of hybrid model is a nice challenge for everyone.

All interactions are implemented using tabulated user-specified potentials. You will need to do your book-keeping carefully. The system is divided into three groups:

  1. AA: atomistic particles, these are all atoms of the 100 peptide molecules that would also be there in a purely atomistic simulation

  2. VS: virtual sites, these are the virtual sites on the 100 peptide molecules that were added to take care of the interaction with the coarse-grained solvent and ions

  3. CG: the coarse-grained beads, the polarizable water beads and the ions

You first need to make sure all possible atom, CG bead and virtual site particle types are known. The file FORCEFIELDS/hybrid-charged.ff/ffnonbonded.itp has the minimum set implemented for YEY in polarizable water with Na-ions only (D, POL, Qd, vQa, vN0, vP5, vSC4, vSP1, and vAC2).

Compared to the simpler implementation of Section 3, we now distinguish between virtual sites and coarse-grained beads; they were a single group in the simpler implementation. Due to having three interaction groups we have more table files. Here is an overview of which interactions potentials act between the groups:

  1. AA-AA: table_AA_AA.xvg when building and using the topol.tpr, gro-table_AA_AA.xvg in the FORCEFIELDS directory. The table contains the normal cut-off Lennard-Jones potentials that are used in the OPLS force field. The cut-off is 1.0 nm. The Coulomb potential is the GROMOS reaction field potential (due to Tironi et al), with epsilon_r=1 and epsilon_rf=78, and cut-off rc=1.2 nm. This potential smoothly goes to zero at 1.2 nm and approximately mimics electrostatic screening due to a homogenous surrounding dielectric medium. Note that this electrostatics treatment is different from standard OPLS, where lattice-sum electrostatics is used (usually PME).

  2. AA-CG: table_AA_CG.xvg when building and using the topol.tpr, gro-table_AA_CG.xvg in the FORCEFIELDS directory. The table contains the normal cut-off Lennard-Jones potentials that are used in the OPLS force field. The cut-off is 1.2 nm. It should be noted, however, that in practice, the vdw interaction between atoms and beads is repulsive only and very short-ranged. The purpose of the short-range repulsion is to prevent oppositely charged AA and CG particles attracting each other to infinitely short (zero) distance. This short-range repulsive potential is achieved by a special settings of the non-bonded interaction parameters epsilon and sigma in the ffnonbonded.itp file. It has negative sigma specified. In this manner, the C6 parameter is zero and there is a non-zero C12 parameter, see Gromacs manual section 5.3.2. An example of an entry in the [ nonbond_params ] directive is shown below.

[ nonbond_params ]
; type1 type2  flag  sigma     epsilon
opls_135   Qd   1     -1.00e-01 0.25e+05 ; C6 = 0, C12 = 10E-07

The Coulombic part of the potential is the GROMOS reaction field potential (due to Tironi et al), i.e. the same functional form as used for the AA-AA interactions, with epsilon_rf=78, and cut-off rc=1.2 nm, but with er=1.45 instead of 1. Wassenaar et al.[6] showed that this was a reasonable choice for a Martini-GROMOS hybrid model in the sense that charged pairs of amino acids show a PMF as a function of their distance that is similar to the all-atom PMF. It is one of the areas, however, in which progress needs to be made to get to a better model.

  1. AA-VS: table.xvg when building and using the topol.tpr, table.xvg in the FORCEFIELDS directory. The standard because no this pair is not listed in the energygrp_table keyword in the settings (.mdp) file. All potentials here are zero everywhere. Computation efficiency is gained, however, if in the settings file the keyword energygrp-excl specifies this pair as excluded.

  2. CG-CG: table_CG_CG.xvg when building and using the topol.tpr, mar-table_CG_CG.xvg in the FORCEFIELDS directory. The table contains the shifted coulomb and shifted Lennard-Jones potentials of the Martini model with standard settings, i.e. cut-off 1.2 nm for both and a switching at 0.9 nm for the LJ potential, and epsilon_r=2.5 because we use the polarizable water model.

  3. CG-VS: table_CG-VS when building and using the topol.tpr, gro-table_CG_VS.xvg in the FORCEFIELDS directory. The table contains the shifted Lennard-Jones potentials of the Martini model with standard settings, i.e. cut-off 1.2 nm and a switching at 0.9 nm for the LJ potential. The Coulomb potential is everywhere zero because electrostatic interactions between AA and CG particles are calculated. Note that an alternative choice would be to treat the electrostatic interactions between peptides and surroundings at the CG level. In that case, the AA-CG electrostatic potential shoud be all zero (in other words, exclude all AA-CG interactions) and the CG-VS electrostatic potential should be the one for the Martini P version, as for CG-CG.

  4. VS-VS: table.xvg when building and using the topol.tpr, table.xvg in the FORCEFIELDS directory. The standard because no this pair is not listed in the energygrp_table keyword in the settings (.mdp) file. All potentials here are zero everywhere. Computation efficiency is gained, however, if in the settings file the keyword energygrp-excl specifies this pair as excluded.

NOTE that the different dielectric constants (1 for AA-AA, 1.45 for AA-CG, 2.5 for CG-CG) are implemented in the tabulated potentials. Therefore, the value of the epsilon_r setting is 1. There is no other way to differentiate amongst the screening of the electrostatic interactions.

The tabulated potentials must be prepared by you. In the FORCEFIELDS directory, you can find the FORTRAN programs used to prepare the tables files of this tutorial. You can compile them if you need to (use gfortran, but the binaries might work anyway), or modify them to implement other potentials. The input to the programs are also given. I trust that in future, multi-purpose, well-documented python tools will become available.

The tabulated potentials can be viewed using xmgrace, for example, when you are in the FORCEFIELDS directory:

> xmgrace -nxy mar-table_CG_CG.xvg -p xmgr.par

Tools and scripts used in this tutorial

References

[1] G.A. Kaminsky, R.A. Friesner, J. Tirado-Rives and W.L. Jorgensen, Evaluation and Reparametrization of the OPLS-AA Force Field for Proteins via Comparison with Accurate Quantum Chemical Calculations on Peptides, J. Phys. Chem. B 105, 6474-6487 (2001).

[2] S.J. Marrink, H.J. Risselada, S. Yefimov, D.P. Tieleman and A.H. de Vries, The MARTINI force field: coarse grained model for biomolecular simulations, J. Phys. Chem. B 111, 7812-7824 (2007).

[3] G.S. Ayton, W.G. Noid and G.A. Voth, Multiscale modeling of biomolecular systems: in serial and in parallel, Curr. Opin. Struct. Biol. 17, 192-198 (2007).

[4] R. Baron, D. Trzesniak, A.H. de Vries, A. Elsener, S.J. Marrink and W.F. van Gunsteren, Comparison of thermodynamic properties of coarse-grained and atomic-level simulation models, ChemPhysChem 8, 452-461 (2007).

[5] A.J. Rzepiela, M. Louhivuori, C. Peter and S.J. Marrink, Hybrid simulations: combining atomistic and coarse-grained force fields using virtual sites, Phys. Chem. Chem. Phys. 13, 10437-10448 (2011).

[6] T.A. Wassenaar, H.I. Ingolfsson, M. Prieß, S.J. Marrink and L.V. Schafer, Mixing MARTINI: Electrostatic Coupling in Hybrid Atomistic−Coarse-Grained Biomolecular Simulations, J. Phys. Chem. B 117, 3516-3530 (2013).

[7] T.A. Wassenaar, K. Pluhackova, R.A. Bockmann, S.J. Marrink and D.P. Tieleman, Going Backward: A Flexible Geometric Approach to Reverse Transformation from Coarse Grained to Atomistic Models, J. Chem. Theory Comput. 10, 676-690 (2014).

[8] S.O. Yesylevskyy, L.V. Schafer, D. Sengupta and S.J. Marrink, Polarizable Water Model for the Coarse-Grained MARTINI Force Field, PLOS Comput. Biol. 6, e1000810 (2010).