ETH Library
Definition and testing of the
GROMOS force-field versions
54A7 and 54B7
Journal Article
Author(s):
Schmid, Nathan; Eichenberger, Andreas P.; Choutko, Alexandra; Riniker, Sereina; Winger, Moritz; Mark, Alan E.; van Gunsteren,
Wilfred F.
Publication date:
2011-07
Permanent link:
https://doi.org/10.3929/ethz-b-000038976
Rights / license:
In Copyright - Non-Commercial Use Permitted
Originally published in:
European Biophysics Journal 40(7), https://doi.org/10.1007/s00249-011-0700-9
This page was generated automatically upon download from the ETH Zurich Research Collection.
For more information, please consult the Terms of use.
ORIGINAL PAPER
Definition and testing of the GROMOS force-field versions 54A7
and 54B7
Nathan Schmid
•
Andreas P. Eichenberger
•
Alexandra Choutko
•
Sereina Riniker
•
Moritz Winger
•
Alan E. Mark
•
Wilfred F. van Gunsteren
Received: 18 January 2011 / Revised: 24 February 2011 / Accepted: 17 March 2011 / Published online: 30 April 2011
Ó European Biophysical Societies’ Association 2011
Abstract New parameter sets of the GROMOS biomo-
lecular force field, 54A7 and 54B7, are introduced. These
parameter sets summarise some previously published force
field modifications: The 53A6 helical propensities are
corrected through new u/w torsional angle terms and a
modification of the N–H, C=O repulsion, a new atom type
for a charged -CH
3
in the choline moiety is added, the
Na
?
and Cl
-
ions are modified to reproduce the free energy
of hydration, and additional improper torsional angle types
for free energy calculations involving a chirality change are
introduced. The new helical propensity modification is
tested using the benchmark proteins hen egg-white lyso-
zyme, fox1 RNA binding domain, chorismate mutase and
the GCN4-p1 peptide. The stability of the proteins is
improved in comparison with the 53A6 force field, and
good agreement with a range of primary experimental data
is obtained.
Keywords GROMOS 54A7 Force field
Secondary structure
Abbreviations
CM Chorismate mutase
FOX Fox1 RNA binding domain
GCN GCN4-p1 peptide
HEWL Hen egg-white lysozyme
PDB Protein Data Bank
RMSD Root-mean-square deviation
SPC Simple point charge
Introduction
Biomolecular simulation involves four major challenges:
(1) the force field must faithfully represent the atomic and
molecular interactions, (2) the conformational space must
be sampled in a manner which is both fast and efficient, (3)
a Boltzmann configurational ensemble must be generated in
order to reproduce thermodynamic quantities and (4)
appropriate experimental data must be available against
which the simulations can be validated. The quality of the
force field is perhaps the most important of these issues and
is the aspect addressed in this work. The interaction
between the atoms in a system must be described with
sufficient accuracy to reproduce the properties and mecha-
nisms underlying the processes of interest. In the first
generation of classical biomolecular force fields such as
AMBER (Weiner and Kollman 1981; Pearlman et al. 1995;
Cornell et al. 1995), CHARMM (Brooks et al. 1983;
MacKerell et al. 1995, 1998), OPLS-AA (Jorgensen and
Tirado-Rives 1988; Jorgensen et al. 1996) and GROMOS
(van Gunsteren et al. 1996) the parameters were chosen so
as to reproduce either spectroscopic or crystallographic
structural data. Subsequently, increasing computer power
Electronic supplementary material The online version of this
article (doi:10.1007/s00249-011-0700-9) contains supplementary
material, which is available to authorized users.
N. Schmid A. P. Eichenberger A. Choutko S. Riniker
W. F. van Gunsteren (&)
Laboratory of Physical Chemistry, Swiss Federal Institute
of Technology ETH, 8093 Zu
¨
rich, Switzerland
e-mail: wfvgn@igc.phys.chem.ethz.ch
M. Winger A. E. Mark
School of Chemistry and Molecular Biosciences,
The University of Queensland, Brisbane, Australia
123
Eur Biophys J (2011) 40:843–856
DOI 10.1007/s00249-011-0700-9
made the simulation of liquids possible and condensed-
phase thermodynamic data such as densities, energies and
free energies were included in the parameterisation, leading
to second-generation force fields. For example, the GRO-
MOS 45A4 (Daura et al. 1998; Schuler et al. 2001)
parameter set was parameterised against the thermo-
dynamic properties of aliphatic chains. In the subsequent
generation of the GROMOS force field the polar amino acid
side-chains and peptide backbone moiety were repara-
metrised. This resulted in the GROMOS 53A6 parameter
set (Oostenbrink et al. 2004, 2005). However, the GRO-
MOS 53A6 force field, in which the hydration properties of
amino acid analogues are in good agreement with experi-
ment, did not improve the stability of the dominant fold for
all peptides (Oostenbrink et al. 2005). Rather, short
a-helices were found to be less stable than expected. This
suggested that the dihedral-angle parameters of the back-
bone transferred from the earlier version of the force field
were no longer appropriate. Subsequently, a variety of
torsional-angle potential energy functions were investigated
by different workers. Cao et al. (2009) proposed a correc-
tion where the u=w torsional angle terms were repara-
metrised and a torsional cross term depending on the sum of
the u- and w-angles was added. In an alternative approach
leading to the 54A7 force field, the torsional angle terms
were reparametrised based on fitting to a large set of high-
resolution crystal structures (Xu et al. in preparation) and
the N–O non-bonded interactions between the peptide
nitrogen and oxygen atoms was adjusted to be less repul-
sive. In the present work we test these alternative approa-
ches using four different test systems: hen egg-white
lysozyme (HEWL, 129 residues), the fox1 RNA binding
domain (FOX, 88 residues and 7 RNA bases), chorismate
mutase (CM, 165 residues) and the GCN4-p1 peptide
(GCN, 16 residues). Because the results for the two men-
tioned modifications yielded similar results, the simpler one
was adopted in the new 54A7 parameter set. Recently,
improved parameters for the simulation of lipids were
reported (Poger et al. 2010). Incorporation in the GROMOS
53A6 force field required the definition of an additional
atom type, which led to the 54 atom types of the 54A7 set.
Finally, new Na
?
and Cl
-
parameters from Reif and
Hu
¨
nenberger (2010) were added and an additional improper
dihedral angle type was defined in order to facilitate free
energy calculations involving a change in chirality. The
definition of the changes from the 53A6 parameter set to the
54A7 parameter set is given in this work.
Definition of the GROMOS 54A7 force field
The GROMOS force field 54A7 is a modification of the
GROMOS 53A6 force field, with four modifications:
1. The torsional-angle energy term for the polypeptide
u- and w-dihedral angles is modified in conjunction
with a change of the combination prescription of the
C
12
van der Waals parameters for the atom type
pair O(IAC = 1)–N(IAC = 6):
(a) In the selection table for the repulsive van der
Waals C
12
1/2
(I,I) parameters, Table 8 of Oosten-
brink et al. (2004), the type for the O(IAC = 1)–
N(IAC = 6) pair (i.e. line 1, column 6) is changed
from ‘‘2’’ to ‘‘1’’. This means that a C
12
1/2
(O,O)
value of 1.000 9 10
-3
[kJ mol
-1
nm
12
]
1/2
for the
O-atom (IAC = 1) is selected for the interaction
with an N-atom (IAC = 6) compared with the
C
12
1/2
(O,O) value of 1.130 9 10
-3
[kJ mol
-1
nm
12
]
1/2
in 53A6.
(b) Four different dihedral torsional angle types
are added to Table 5 of Oostenbrink et al.
(2004):
Type code K
un
[kJ mol
-1
]
cos(d
n
) m
n
Example K
un
[kcal mol
-1
]
42 3.50 -1 2 –CHn–C– 0.84
43 2.80 ?1 3 –CHn–N– 0.67
44 0.70 -1 6 –CHn–N– 0.17
45 0.40 ?1 6 –CHn–C– 0.10
In the molecular topology building blocks for the a-pep-
tides and b-peptides the dihedral-angle type 39 (53A6) in
the backbone C–N–CA–C dihedral angle (a-peptide) or the
backbone C–N–CB–CA dihedral angle (b-peptide) is
changed to type 44 (54A7) and the same dihedral angle
with type 43 (54A7) is added. In addition, the dihedral-
angle type 40 (53A6) for the backbone N–CA–C–N dihe-
dral angle (a-peptide) or the backbone CB–CA–C–N
dihedral angle (b-peptide) is to be changed to type 45
(54A7), and the same dihedral angle with type 42 (54A7) is
added.
These changes increase the hydrogen-bonding interaction
between the N–H and the C=O groups in the polypeptide
backbone and bring the u- and w-angle distributions for a
number of proteins more in line with the preferences
observed in PDB protein structures.
2. A new van der Waals non-bonded atom type for a
charged -CH3 group (IAC = 54) is introduced in
Table 6 of Oostenbrink et al. (2004) in order to
increase the repulsion between the positively charged
-CH3 groups of the choline moiety and the negatively
charged -OM oxygen atoms of the phosphate moiety
844 Eur Biophys J (2011) 40:843–856
123
in dipalmitoylphosphatidylcholine (DPPC)-type lipids,
see Poger et al. (2010).
(a) In the normal van der Waals parameters table,
Table 7 of Oostenbrink et al. (2004), an addi-
tional atom type is added:
IAC Atom type [C6(I,I)]
1/2
[kJ mol
-1
nm
6
]
1/2
[C12(I,I)]
1/2
10
-3
[kJ mol
-1
nm
12
]
1/2
123
54 CH3p 0.09805 5.162 – –
(b) In the third-neighbour van der Waals parameters
table, Table 9 of Oostenbrink et al. (2004), an
additional atom type is added:
IAC Atom type [C6(I,I)]
1/2
[kJ mol
-1
nm
6
]
1/2
[C12(I,I)]
1/2
10
-3
[kJ mol
-1
nm
12
]
1/2
123
54 CH3p 0.08278 2.456 – –
(c) In the selection table for the repulsive van der
Waals C
12
1/2
(I,I) parameters, Table 8 of Oosten-
brink et al. (2004), the matrix is enlarged by one
column and one row to accommodate the new
atom type. For all pairs, type 1 is selected with
the exception of the OM(IAC = 2)-CH3p(IAC =
54) pair (line 2, column 54) where type 3 is
selected.
These changes increase the area per lipid for DPPC bilayers,
see Poger et al. (2010), by increasing the repulsion of their
phosphate oxygen which leads to the Lennard–Jones
parameters r = 0.3877 nm and = 0.3433 kJ mol
-1
for the choline CH
3
(CH3p)-phosphate oxygen (OM)
pair.
3. The van der Waals non-bonded interaction parameters
for the Na
?
and Cl
-
ions are taken from the set L
proposed by Reif and Hu
¨
nenberger (2010).
(a) In the normal van der Waals parameters table,
Table 7 of Oostenbrink et al. (2004), the param-
eters for Na
?
(IAC = 37) and Cl
-
(IAC = 38)
are set to:
IAC Atom type [C6(I,I)]
1/2
[kJ mol
-1
nm
6
]
1/2
[C12(I,I)]
1/2
10
-3
[kJ mol
-1
nm
12
]
1/2
123
37 NA
?
0.0088792 0.2700 0.2700 –
38 CL
-
0.11318 7.776 7.776 7.776
The third-neighbour CS6 and CS12 parameters were also
changed accordingly.
These changes bring the solvation properties of Na
?
and
Cl
-
more in line with experiment.
4. To facilitate the calculation of differences in free
energy involving changes in chirality, two additional
improper dihedral-angle types are defined, i.e. one
with -35° and one with 180° as energy minimum.
(a) The two new improper (harmonic) dihedral-angle
types are added to Table 4 of Oostenbrink et al.
(2004):
Type
code
K
nn
[kJ mol
-1
deg
-2
]
n
0n
[deg] Example
of usage
K
nn
[kcal
mol
-1
rad
-2
]
4 0.0510 180.0 Planar
groups
40
5 0.102 -35.26439 Tetrahedral
centres
80
Corresponding changes in the 53B6 force field for
in vacuo simulations yield the 54B7 force field for in vacuo
simulations.
Materials and methods
All simulations were performed using a modified version of
the GROMOS biomolecular simulation software (Christen
et al. 2005) in conjunction with the parameter set of the
GROMOS force field indicated: 45A4 (Schuler et al. 2001),
53A6 (Oostenbrink et al. 2004), Liu (Cao et al. 2009)or
54A7. Note that the CM system was not simulated using the
Liu modification. Also, the GCN system was only simulated
using the 53A6 and 54A7 force fields, meaning that a longer
timescale could be examined. The dihedral-angle potential
energy function in the GROMOS force field for the backbone
u- and w-angles is defined as
Eur Biophys J (2011) 40:843–856 845
123
V
GROMOS
ðu; wÞ¼
X
N
u
i¼1
K
i
1 þ cos m
i
u d
i
ðÞ½
þ
X
N
w
i¼1
K
i
1 þ cos m
i
w d
i
ðÞ½; ð1Þ
where N
u
and N
w
are the number of terms for one dihedral
angle, see Table 1. In the approach of Liu and coworkers, a
cross term that depends on the sum of the u- and w-angle,
V
cross
ðu; wÞ¼
X
N
u;w
i¼1
K
i
1 þcos m
i
u þ wðÞd
i
ðÞ½ð2Þ
is added to the potential energy function (Eq. 1), resulting
in the complete potential energy function for the backbone
u; w dihedral angles. The parameters of the different force
fields are summarised in Table 1. The 45A4 and 53A6
force fields use the same description of the torsional
potential energy term. In the 54A7 force field, these terms
are adjusted and the repulsive (C
12
) term of the Lennard–
Jones potential energy term is changed from type 2 to
type 1, which means that it is less repulsive.
The initial coordinates of the protein and RNA mole-
cules were taken from structures (Fig. 1) deposited in the
Protein Data Bank (PDB). The entry codes were 1AKI for
HEWL (Artymiuk et al. 1982), 2ERR for FOX (Auweter
et al. 2006), 2FP2 for CM (Okvist et al. 2006) and 2OVN
for GCN (Steinmetz et al. 2007). In the case of 2ERR and
2OVN, the first structure of the NMR bundle was taken. In
the case of 2FP2, only the first subunit of the dimeric
protein was taken. All hydrogens were (re)generated by the
GROMOS?? (Christen et al. 2005) program gch. Each
system was first energy minimised in vacuo, then the
protein plus RNA molecules in the case of FOX were
solvated in cubic boxes filled with simple point charge
(SPC) water (Berendsen et al. 1981) molecules. Periodic
Table 1 Force-field parameters
of the dihedral-angle term for
peptide backbone u- and
w-torsional angles in the
GROMOS force field
Oostenbrink et al. (2004) and
from Cao et al. (2009)
Term 45A4 and 53A6 (Oostenbrink et al. 2004) 54A7 Liu (Cao et al. 2009)
K/kJ mol
-1
m d K/kJ mol
-1
m d K/kJ mol
-1
m d
u 1.0 6 180° 2.8 3 0° 1.2 1 0°
0.7 6 180° 0.5 1 -120°
0.8 2 180°
1.0 3 0°
w 1.0 6 0° 3.5 2 180° 0.8 1 150°
0.4 6 0° 0.5 3 -90°
0.5 4 90°
u þ w 0.5 1 90°
1.1 2 0°
0.5 3 -135°
0.5 4 120°
Non-bonded C
12
(N,O) type 2 C
12
(N,O) type 1 C
12
(N,O) type 2
Fig. 1 Ribbon pictures of the four proteins investigated: a HEWL,
a-helices (residues 5–14, 25–34, 89–100 and 109–114), 3
10
-helices
(residues 80–83 and 120–123) and b-sheet (residues 43–45, 51–53
and 58–59). b FOX, a-helices (residues 21–32 and 58–69) and b-sheet
(residues 10–14, 38–40, 50–54, 72–73, 76–77 and 79–82). The
nucleic acid is shown in ball-and-stick representation. c CM, a-helices
(residues 4–17, 18–28, 35–52, 56–84, 85–89, 95–116, 117–121,
124–141 and 144–155). d GCN, a-helix (residues 1–16). Colour code:
a-helix (black), 3
10
-helix (green) and b-sheet (red)
846 Eur Biophys J (2011) 40:843–856
123
