General Purpose Coarse-Grained Force Field

Other tutorials

Backward

The method

Mapping

Backmapping aquaporin 1

Polarizable water

CHARMM-GUI Martini Maker

Backward

Backmapping or reverse coarse-graining or fine-graining a coarse-grained (CG) structure requires a correspondence between the two models; i.e., for atomistic and CG: which atoms make up which bead. Actually, an atom can in principle contribute to several beads. A backmapping 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 used in this tutorial is backward^[1], developed by Tsjerk Wassenaar. Backward 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.

First, we go over the scripts used in the backward^[1] method. Second, we describe the structure and setup of the needed CG to fine-grained mapping files, taking the mapping file for valine as an example. Last, we demonstrate the method by backmapping a CG membrane embedded aquaporin 1, as described in^[2].

More extensive discussion and examples, including tutorial material can be found in the paper by Wassenaar et al.^[1], and Supporting Material to that paper.

The method

The backward program is available here (download the zip-file for the version of Mapping files you need and possibly download another version of the initram.sh script), and consists of three scripts and a number of CG to fine-grained mapping definition files. The scripts are backward.py, initram.sh and in the Mapping directory __init__.py and a number of .map files. The .map files describe the CG to fine-grained mapping and a file needs to be provided for each molecule (or its building blocks) in your system. The __init__.py script interprets the .map files. The backward.py script performs the actual backmapping and initram.sh is a bash wrapper that runs a series of minimizations and molecular dynamics steps, using the fine-grained force field to push the initial backmapped structure to one that satisfies the fine-grained force field.

Mapping

A requirement for the procedure to work is that the subdirectory Mapping contains definitions for how the atomic positions are generated from the CG positions. The subdirectory Mapping contains a file for each residue and/or molecule that can be backmapped, named for the atomistic target force field, e.g. val.oplsaa.map for a valine residue targeted to OPLS-AA. Updates of the Martini and/or other force fields may results in changes in the mapping files, so beware you may need to make changes and/or find the appropriate mapping files. You can download mapping files for each of the lipids in the lipidome library, for example the mapping file for cholesterol. The structure of a .map file is explained below for the valine residue. The file val.oplsaa.map reads:

[ molecule ]
VAL
[ martini ]
BB SC1
[ mapping ]
oplsaa
[ atoms ]

    1   N     BB
    2   H     BB

    3   CA    BB

    4   HA    BB

    5   CB    SC1 BB
    6   HB    SC1 BB
    8   CG1   SC1
    9   HG11  SC1
   10   HG12  SC1
   11   HG13  SC1
   12   CG2   SC1
   13   HG21  SC1
   14   HG22  SC1
   15   HG23  SC1
   16   C     BB
   17   O     BB
[ chiral ]
  CB     CA    N    C
  HB     CA    N    C
[ chiral ]
  HA   CA  N CB  C  ;L-Val
 ;HA   CA  N C  CB  ;D-Val
[ out ]
  CG2  CB CG1 CA
  HG21 CB CG1 CA
  HG22 CB CG1 CA
  HG23 CB CG1 CA

Directives analogous to gromacs topologies contain specifications that build the atomistic structure from the CG positions. The [ molecule ] directive contains the name of the residue or molecule. The [ martini ] directive contains the names of the CG beads in the Martini model: valine has two beads called BB and SC1. The [ mapping ] directive contains the name of the object model. Note that this directive may contain multiple object models. If you do not care for the naming convention of different force fields, you can use the same mapping file for the CHARMM36 and OPLS-AA/L force fields, because these are both all-atom models which in the GROMACS implementation also use the same order of the atoms (if not the same names). Thus, the mapping files for the methylated terminal ends explicitly state that they can be used for mapping to both OPLS-AA and CHARMM36 force fields.

The [ atoms ] directive contains the index numbers and names of the atoms in the object model and their relation to the CG beads. Note that a single atom may be in a relation with more than one CG bead. The back-mapping procedure starts by putting each atom that is related to a single bead on the position of that bead. If an atom is related to more than one bead, it will be placed on the weighted average position of the beads listed. It is allowed to list the same bead multiple times; thus the line

4  OE1  BB  BB  BB  SC1  SC1

places the fourth atom (with name OE1) of the residue on the line connecting the BB and SC1 beads at 2/5 of the distance between the beads, starting at the BB bead. This mechanism is a simple aid to position atoms already at fairly reasonable starting positions. Using the -kick flag displaces all atoms randomly after their initial placement. Note that the script applies a random kick to atoms that are initially put at exactly the same place, e.g. because they are defined by the position of a single bead. Thus, no two atoms will be on top of each other. Switching on an atomistic force field usually results in an error if two interacting atoms are at exactly the same place.

The backward procedure implements a few other more sophisticated mechanisms to place atoms and some are used in the valine residue. Implementation can be found in the file Mapping/__init__.py. The [ chiral ] directive generates stereochemically specifically arranged groups of atoms. As is seen for valine, two types of input can be provided. In the first instance of the [ chiral ] directive, four atoms are listed on a line. The first atom is the atom to be placed (named A in Figure 1e, chiral*) on the basis of the positions of the other three. The Figure shows how vectors are defined from the positions of the other three particles to construct the position of the first atom. In the second instance of the [ chiral ] directive, five atoms are listed on a line. This corresponds to the construction shown in Figure 1d, chiral. The [ out ] directive may be used to spread out several equivalent atoms (as on a CH3-group) away from the center of the bead in a reasonable manner, as shown in Figure 1c, out. Again, be aware that atoms initially placed on the same spot, are randomly displaced; therefore, using the same reference atoms as in the example still leads to different positions for the generated atoms.

Figure 1 | Figure from^[1]. Prescriptions for placing atom A on the basis of the positions of other atoms (upper case letters B, C, D, E). The other atoms define vectors (lower case letters a, b, c, d, e, p, q) that are used to calculate the position of A. A bar is used over the vectors to denote their normalization. The x denotes taking the outer product of two vectors.

Backmapping aquaporin 1

Here we backmapp a CG membrane embedded aquaporin 1 into CHARMM36 atomistic coordinates, as described in^[2]. The files for this part of the tutorial are available in aquaporin_backmap.tar.gz. Missing residues were added to aquaporin 1 and it was converted to Martini CG coordinates, solvated in a CG POPC bilayer with ions and polarizable water. Then simulated for 100 ns at the CG level with position restrains on the protein, see^[2] for method details and CG_posre.gro for final coordinates. Note, without position restraints on the protein the CG protein might (depending on the protein in question and the CG force field used) evolve to far away from a possible fine-grained structure, rendering backmapping impossible.

We are going to use the initram-v5.sh script, which calls backward.py and then runs a series of minimization and equilibrium simulations to get the final backmapped structure. To run the script we need the following:

* The CG structure to backmap, provided in CG_posre.gro

* A complete fine-grained force field corresponding to all the CG molecules in CG_posre.gro. Here we use CHARMM36, see all .itp files provided and topol.top, which contains the molecules in the same order they are present in CG_posre.gro and with the same names. Note, water and ions can be skipped in the .top files as they are automatically detected by backward.py.

* A .map file in the Mapping directory for all residues and molecules to be backmapped (water and ions can also be skipped here as their definitions are included in backward.py). The .map files are collected in a directory called Mapping. NOTE that the coarse-grained structure used here uses an older Martini model for oleoyl chains, represented by 5 beads, rather than the current version of Martini that uses 4-bead representation of oleoyl chains. Therefore, you will need the 'older' backward.zip. file and the 'newer' script initram-v5.py.

The initram.sh script uses the gromacs package so a proper version needs to be sourced, i.e. a version of initram.sh that is compatible with the gromacs version you are using.

Run the script using the following command:

$ ./initram-v5.sh -f CG_posre.gro -o aa_charmm.gro -to charmm36 -p topol.top

After successful backmapping the aa_charmm.gro file will contain POPC membrane embedded aquaporin 1 as well as the solvent at fully atomistic resolution according to the CHARMM36 force field. Figure 2 illustrates the backmapping of a few POPC molecules. When running the backmapping script please keep in mind that initram.sh generates a significant number of temporary files so backmapping in a separate directory can be a good idea and that backward.py used random kicks to initially displace the atoms so rerunning the same command can give different results (and even though some runs might results in an error others may not).

Figure 2 | Backmapped all-atom representation of a Martini POPC membrane. Figure from^[2].

Polarizable water

In this part of the tutorial we will convert water in an existing Martini-system to polarizable water. In the polarizable Martini paper^[3] the model is described as follows:

"The polarizable CG water consists of three particles instead of one in the standard Martini force eld. The central particle W is neutral and interacts with other particles in the system by means of the Lennard-Jones interactions, just like the standard water particle. The additional particles WP and WM are bound to the central particle and carry a positive and negative charge of +0.46 and -0.46, respectively. They interact with other particles via a Coulomb function only, and lack any LJ interactions. The bonds W-WP and W-WM are constrained to a distance of 0.14 nm. The interactions between WP and WM particles inside the same CG water bead are excluded, thus these particles are "transparent" toward each other. As a result the charged particles can rotate around the W particle. A harmonic angle potential with equilibrium angle phi = 0 rad and force constant K_h = 4.2 kJ mol^-1 rad^-2 is furthermore added to control the rotation of WP and WM particles and thus to adjust the distribution of the dipole momentum."

If you want to solvate a new system with polarizable water, follow the steps in the lipid or protein Martini tutorials. Instead of a normal waterbox and .mdp/.itp files use the waterbox containing polarizable water and polarizable Martini .mdp and .itp files (available from the Martini website: cgmartini.nl). Note that minimizing a polarizable Martini system requires some tweaking, as described below. We encourage you to get all downloads yourself, but if you insist on taking a short-cut, you can download an archive with just the basic files here and with worked example output here.

If you have an existing system with normal Martini water and want to change to polarizable water you may use the python script triple-w.py. For this example we convert dppc_bilayer.gro/dppc_bilayer.top available here.

1. Create a new gro file:

$ python triple-w.py dppc_bilayer.gro

The python-script triple-w.py adds positive and negative sites at a small distance to every central water bead in a .gro file.

2. Adapt the .top file. Make sure the polarizable version of the particle definition file  martini_v2.2refP.itp is included if you want to use the refined polarizable water model, otherwise include the older version. An example of an .mdp file can be downloaded here. Make sure your .mdp is compatible with the polarizable water version; this one use the reaction-field electrostatics, but the refined polarizable model is particularly suitable for long-range electrostatics (coulomb_type = PME). The .itp files for the lipids and possibly other molecules do not have to be changed.

3. If polarizable water is used in combination with proteins or peptides, all AC1 and AC2 beads have to be replaced by normal C1 and C2 beads. AC1 and AC2 are obsolete in polarizable Martini.

4. Minimize the system. For polarizable water to minimize without problems, it is SOMETIMES necessarry to change the constraints to stiff bonds. Using the ifdef-statement in the .itp file, this can be set using a define value in the .mdp file (have a look at the bottom of martini_v2.P.itp to see how it works).

5. Generate the input files for the minimization run:

$ gmx grompp -f martini_v2P.x_new-rf-em.mdp -c dppc_bilayer_PW.gro -p dppc_bilayer.top -o em -maxwarn 1

6. Proceed with equilibration (e.g. using martini_v2P.x_new-rf-eq.mdp) and production (e.g. using martini_v2P.x_new-rf-md.mdp) runs. Note, the mdp files linked to here have not been tested extensively (August 2017). They appear to give a stable simulation of the DPPC bilayer in polarizable water, but the scope of these parameters has not been validated systematically in the systems described in the original paper.   

CHARMM-GUI Martini Maker

The CHARMM-GUI input generator now supports the Martini force field^[4]. The Martini Maker supports building solution, micelle, bilayer, and vesicle systems, and systems with randomly distributed lipids. The Martini Maker has all the common Martini lipids as well as supporting variations of the Martini model (polar and nonpolar Martini, Dry Martini, and ElNeDyn). The CHARMM-GUI has an easy to use step-by-step user interface, please explore the Martini Maker input generator at http://www.charmm-gui.org/.

Following is a very short tutorial on how to make a Dry Martini vesicle using Martini Maker.

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

[2] Pluhackova, K., Wassenaar, T. A., and Böckmann, R. A. (2013) Molecular Dynamics Simulations of Membrane Proteins, in Methods in Molecular Biology (Rapaport, D., and Herrmann, J. M., Eds.), pp 85–101. Humana Press.

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

[4] Qi, Y., Ingólfsson, H.I., Cheng, X., Lee, J., Marrink, S.J., and Im, W. (2015) CHARMM-GUI Martini Maker for coarse-grained simulations with the Martini force field. JCTC, doi:10.1021/acs.jctc.5b00513.