University of Groningen
Software News and Update Reconstruction of Atomistic Details from Coarse-Grained
Structures
Rzepiela, Andrzej J.; Schafer, Lars V.; Goga, Nicolae; Risselada, H. Jelger; De Vries, Alex H.;
Marrink, Siewert J.
Published in:
Journal of Computational Chemistry
DOI:
10.1002/jcc.21415
IMPORTANT NOTE: You are advised to consult the publisher's version (publisher's PDF) if you wish to cite from
it. Please check the document version below.
Document Version
Publisher's PDF, also known as Version of record
Publication date:
2010
Link to publication in University of Groningen/UMCG research database
Citation for published version (APA):
Rzepiela, A. J., Schafer, L. V., Goga, N., Risselada, H. J., De Vries, A. H., & Marrink, S. J. (2010). Software
News and Update Reconstruction of Atomistic Details from Coarse-Grained Structures.
Journal of
Computational Chemistry
,
31
(6), 1333-1343. https://doi.org/10.1002/jcc.21415
Copyright
Other than for strictly personal use, it is not permitted to download or to forward/distribute the text or part of it without the consent of the
author(s) and/or copyright holder(s), unless the work is under an open content license (like Creative Commons).
The publication may also be distributed here under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license.
More information can be found on the University of Groningen website: https://www.rug.nl/library/open-access/self-archiving-pure/taverne-
amendment.
Take-down policy
If you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediately
and investigate your claim.
Downloaded from the University of Groningen/UMCG research database (Pure): http://www.rug.nl/research/portal. For technical reasons the
number of authors shown on this cover page is limited to 10 maximum.
Software News and Update
Reconstruction of Atomistic Details from
Coarse-Grained Structures
ANDRZEJ J. RZEPIELA, LARS V. SCHÄFER, NICOLAE GOGA, H. JELGER RISSELADA,
ALEX H. DE VRIES, SIEWERT J. MARRINK
Groningen Biomolecular Sciences and Biotechnology Institute & Zernike Institute for Advanced
Materials, University of Groningen, Nijenborgh 4, 9747 AG Groningen, The Netherlands
Received 18 May 2009; Revised 12 August 2009; Accepted 14 August 2009
DOI 10.1002/jcc.21415
Published online 19 January 2010 in Wiley InterScience (www.interscience.wiley.com).
Abstract: We present an algorithm to reconstruct atomistic structures from their corresponding coarse-grained (CG)
representations and its implementation into the freely available molecular dynamics (MD) program package GROMACS.
The central part of the algorithm is a simulated annealing MD simulation in which the CG and atomistic structures are
coupled via restraints. A number of examples demonstrate the application of the reconstruction procedure to obtain low-
energy atomistic structural ensembles from their CG counterparts. We reconstructed individual molecules in vacuo (NCQ
tripeptide, dipalmitoylphosphatidylcholine, and cholesterol), bulk water, and a WALP transmembrane peptide embedded
in a solvated lipid bilayer. The first examples serve to optimize the parameters for the reconstruction procedure, whereas
the latter examples illustrate the applicability to condensed-phase biomolecular systems.
© 2010 Wiley Periodicals, Inc. J Comput Chem 31: 1333–1343, 2010
Key words: back-mapping; multiscale simulation; coarse-grained; molecular dynamics
Introduction
Molecular dynamics (MD) simulations have been successfully used
to investigate biomolecular systems, such as proteins, DNA, or
lipid membranes.
1–3
The huge computational effort involved in
conventional atomistic MD simulations currently limits accessible
simulation times to hundreds of nanoseconds and length scales to
tens of nanometers. However, most biomolecular processes occur at
slower time scales and, at the same time, often involve larger length
scales. To study such processes, coarse-grained (CG) models, in
which several atomistic particles are grouped together into effective
beads, have recently gained more and more popularity in the field.
4
The main benefit of these models is their computational efficiency
due to the reduced number of interaction sites. In this way, many
CG methods allow probing the structural dynamics of large systems
on time scales up to milliseconds and length scales up to hundreds
of nanometers. This large gain in efficiency, however, comes at the
cost of a reduced accuracy compared to atomistic (or all-atom, AA)
models due to the inherent simplifications. Thus, tools are desir-
able that allow switching between the different levels of resolution,
i.e., to back-map the atomistic structural ensemble that underlies
a CG representation, thereby combining the efficiency of CG with
the accuracy of AA models. This would allow CG models to live
up to their full potential by applying them to explore large regions
of phase space, followed by zooming in on the atomistic details of
interesting configurations.
In this contribution, we present our recent implementation of
an algorithm to reconstruct AA structures from their corresponding
CG representations into the fast and freely available MD program
package GROMACS.
5, 6
The basic idea of the reconstruction algo-
rithm is to carry out a simulated annealing (SA)
7
MD simulation
of an atomistic system that is coupled to its corresponding CG
system via restraints. During the SA, the system is cooled down
from a high initial temperature to a desired target temperature, thus
allowing the system to cross energy barriers and optimize under the
restraints. Finally, the coupling is gradually removed to ensure a
smooth relaxation of the final AA structure.
By means of a number of example applications presented in
this work, we systematically and extensively tested the reconstruc-
tion procedure for different molecule types. Proper reconstructed
ensembles were obtained using optimized parameters. Furthermore,
although the applications presented in this contribution are based
on the MARTINI CG force field
8–10
and the GROMOS96 AA
force field,
11, 12
the algorithm is general and can be applied in a
Correspondence to: S.J. Marrink; e-mail: s.j.marrink@rug.nl
Contract/grant sponsor: Netherlands Organisation for Scientific Research
(NWO), Veni grant; Contract/grant number: 700.57.404
Contract/grant sponsor: National Supercomputing Facilities (NCF); Con-
tract/grant numbers: NRG.2008.04, SH.002.08
© 2010 Wiley Periodicals, Inc.
1334 Rzepiela et al.
•
Vol. 31, No. 6
•
Journal of Computational Chemistry
straightforward manner to a wide range of force fields and molecule
types, because it does not rely on libraries of predefined fragments.
A number of techniques to reconstruct atomistic details from CG
structures have been proposed in the literature. Some of these meth-
ods
13–15
make use of a multiscale Hamiltonian exchange approach,
in which the system is at the same time represented at both levels of
resolution (and at several intermediate, i.e., mixed levels), and con-
tinuous attempts are made along an MD trajectory to switch between
neighboring representations. The adaptive resolution schemes use
a spatial compartmentalization of the simulation system and allow
an on-the-fly interchange between atomistic and CG representa-
tions.
16–20
Other methods are more similar to the approach presented
here in that they aim at reconstructing atomistic structures from
single CG structures.
21–27
Most of these methods differ from the
approach presented in this work by using libraries of predefined
fragments to construct the initial atomistic structure, followed by a
minimization and equilibration protocol. Additionally, these meth-
ods were often optimized only for the specific application at hand
and may therefore not be generally applicable and transferable to
other molecular systems and force fields.
The rest of this article is organized as follows. We set out to
describe the methodological basis of the algorithm and detail the
implementation. Additionally, we present details of the user inter-
face. Subsequently, a number of examples will serve to illustrate
the application of the algorithm, optimization of SA parameters,
and generation of unbiased distributions. As a first example, the
reconstruction algorithm was applied to obtain atomistic structural
ensembles of different biomolecules in vacuum: the NCQ tripep-
tide, the dipalmitoylphosphatidylcholine (DPPC) lipid molecule,
and cholesterol. Then, a box of atomistic water was reconstructed
from CG water. Finally, a system composed of the WALP20
transmembrane peptide embedded in a solvated DPPC bilayer
was reconstructed from its CG representation. The latter exam-
ples demonstrate the suitability of the reconstruction algorithm to
condensed-phase biomolecular systems. A short conclusive section
ends this article. In the appendix, we give a list of GROMACS
commands used to carry out the reconstruction simulations.
Methods
Our goal is to generate a low-energy atomistic (AA) ensemble that
underlies its corresponding CG system. This is achieved by means
of a three-step reconstruction algorithm. Initially, AA particles are
positioned close to their reference CG beads. Then, a SA proce-
dure is used, during which the AA system is coupled to the CG
system via harmonic restraints. Finally, the restraints are gradually
removed to yield a relaxed atomistic system. As a prerequisite, the
reconstruction algorithm requires the definition of a mapping of the
AA structure to the CG structure, i.e., a prescription of which AA
particles are represented by which CG beads. The mapping can be,
e.g., defined via the center of mass of AA particles (as done in this
work, see later) or via the positions of the C
α
atoms in an amino
acid chain.
Initial Placement of Atomistic Particles
To generate a starting configuration, the AA particles are randomly
positioned within a sphere around their corresponding CG beads.
In the MARTINI CG force field, on average four heavy atoms
are mapped onto one CG bead. For the applications based on the
MARTINI force field presented in this work, the radius of the
sphere was chosen as r
CG
= 0.3 nm, which roughly corresponds
to the typical van der Waals radius of a MARTINI bead. Such
an initial placement of AA particles close to their expected final
positions significantly speeds up the reconstruction procedure for
condensed-phase systems.
Restrained Simulated Annealing
The central part of the reconstruction algorithm is a SA MD protocol.
During the annealing, a restraining potential keeps the center of mass
of groups of AA particles close to their reference CG beads. The
complete potential U
tot
is described by eq. (1)
U
tot
= U
AA
+ U
restr
, (1)
where U
AA
represents the atomistic force field and U
restr
is a
harmonic potential,
U
restr
=
n
i=1
k
2
r
CG
i
− r
AA,com
i
2
.(2)
In eq. (2), r
CG
i
is the position of the reference CG bead i, r
AA,com
i
the
center of mass of the AA particles that are mapped to this bead, n
is the total number of CG beads, and k a restraining force constant.
The r
AA,com
i
are updated at every simulation step.
Using this restraining potential, a SA MD procedure is used
to generate low-energy structures. The temperature is gradually
decreased from a high starting value to the desired target tempera-
ture, thus allowing the system to rapidly cross energy barriers and
find a low-energy minimum at the end of the simulation. The cru-
cial parameters are the starting temperature and the cooling rate; too
low initial temperatures or too rapid cooling will not yield properly
equilibrated final structures.
Because the AA particles are initially placed randomly around
the corresponding CG beads, very large forces will occur at the
beginning of the simulation and cause numerical instabilities. To
ensure a stable simulation, these forces are reset to a specified
threshold; this threshold is then linearly increased during the anneal-
ing simulation to increase the sampling. This is possible because
the system optimizes in the course of the simulation, and thus
large forces occur less frequently. The initial value for the force
threshold, F
cap,0
= 15, 000 kJ mol
−1
nm
−1
, and its increase rate,
F
cap
= F
cap,0
+ At, with A = 100 kJ mol
−1
nm
−1
ps
−1
, were chosen
high enough to ensure that the relevant energy barriers are crossed,
and at the same time low enough to avoid instabilities.
Release of the Restraints
At the end of the SA, the resulting atomistic structure is still coupled
to the CG structure (cf. eq. (2)). However, in general, the (mapped)
AA and CG minimum-energy structures will deviate because of the
inherent differences between the force fields. Therefore, to release
the strain on the AA structure, U
restr
is smoothly removed at the end
Journal of Computational Chemistry DOI 10.1002/jcc
Reconstruction of Atomistic Details 1335
of the reconstruction procedure to yield the final low-energy AA
configuration.
The restrained SA simulation is carried out in the NVT ensemble,
i.e., without pressure coupling. Simulating at constant volume pre-
vents an expansion of the simulation box due to high initial forces.
Any differences in the densities of the CG and atomistic systems
can be taken into account at the end of the reconstruction by switch-
ing on pressure coupling after the release of the restraints. For CG
systems simulated with the MARTINI force field, this difference is
usually small.
Reconstruction of Solvent
A special situation arises for the reconstruction of solvents in which
several small molecules are mapped to one CG bead, such as, e.g.,
water in the MARTINI model, where one CG water bead effectively
represents four atomistic water molecules. To avoid water molecules
diffusing away from their associated CG bead, a potential of the form
U
restr,W
j
=
0 for r
ij
≤ r
CGW
k
W
2
(r
ij
− r
CGW
)
2
for r
ij
> r
CGW
(3)
is used in addition to the restraining potential described in eq. (2).
In eq. (3), r
ij
is the distance between the oxygen atom of AA water
j and the center of mass of the four AA waters that belong to CG
bead i, r
CGW
the cutoff radius, and k
W
the restraining force constant.
Thus, every AA water molecule that moves out of the cutoff radius is
driven back toward the center of mass. To ensure that the additional
external force does not lead to a net center of mass movement of
the four AA water molecules, a counter force is distributed among
the remaining three water molecules such that
4
j=1
F
j
= 0. This
procedure is iteratively applied to every water molecule.
Implementation
We implemented the reconstruction algorithm into the program
mdrun (v3.3), which is the main MD engine of the program package
GROMACS.
5
The conversion between the AA and CG structures is
done with the program g_fg2cg, which reads in both the AA and the
CG topologies. Here, in the AA topology file, an additional mapping
section is included that describes the correspondence between the
two levels of resolution. With this information, g_fg2cg can either
generate a CG structure from a given atomistic structure (in this case,
the former is uniquely defined by the latter through the mapping)
or generate an initial atomistic structure for the SA reconstruction
simulation. The input parameters (cf., Table 1) for the restrained SA
simulation are added to the MD-parameter (mdp) file and are read
in by the GROMACS preprocessor (grompp). During the restrained
SA simulation, the forces due to the external restraining potentials
(Eqs. (2) and (3)) are computed at each integration step and are
added to the forces derived from the original atomistic potential.
To simplify the mapping of large complex molecules, such
as proteins or polycarbohydrates, we have modified the program
pdb2gmx, which uses database files to generate a topology from a
structure file. In these database files, now also the CG/AA mapping
Table 1. Mdp-Parameters for Reconstruction Simulation.
Parameter mdp-option Recommended value
Initial capping force F
cap,0
cap_force 15,000 kJ mol
−1
nm
−1
Capping increase rate A cap_a 100 kJ mol
−1
nm
−1
ps
−1
Restraining force constant k fc_restr 12,000 kJ mol
−1
nm
−2
Radius of CG water r
CGW
r_CGW 0.21 nm
Water restraining force fc_restrW 400 kJ mol
−1
nm
−2
constant k
W
Nr of steps to release rel_steps 5000
restraints
Annealing method annealing single
Annealing time annealing_time 60 ps
Initial annealing temperature annealing_temp 1300 K
is defined for each building block, such as the individual amino
acids.
Dihedral angle states that are separated by high energy barri-
ers can lead to problems during the reconstruction, because the
system might be trapped in the unwanted minimum. Such an
unwanted minimum can be, for example, the cis rotamer of an
amide bond in a polypeptide backbone. A possible solution is to
add dihedral angle restraints to the topology. This can be done with
the program g_dihfix, which processes the structure and topology
files and selects those dihedral angles with a high energy bar-
rier. The user can define a threshold for the dihedral angle force
constant and thereby select certain dihedrals; each selected dihe-
dral can then be restrained to a desirable value, for example, the
one observed in the X-ray crystal structure of a protein. The code
described in this work can be downloaded from our web page under
http://md.chem.rug.nl/∼marrink/coarsegrain.html.
Simulation Details
All simulations were carried out using the GROMACS simulation
package (version 3.3.1)
5
and the implementation of the reconstruc-
tion algorithm described in this work. In the CG simulations, version
2 of the MARTINI forcefield
8–10
was used together with a 20-fs inte-
gration time step. The other simulation parameters were set to the
standard values described in the original publications.
9, 10
The atomistic simulations were carried out with a 2-fs integra-
tion time step, and the temperature was controlled by coupling to
aNo
´
se–Hoover thermostat (τ
T
= 0.1 ps).
30
Because of the ran-
dom initial placement of the atomistic particles (see Methods), no
constraints were applied in the reconstruction simulations, except
for SPC water. For the simulations of the NCQ tripeptide, the 53a6
parameter set of the GROMOS united atom force field
12
was used.
For DPPC, we used the parameters published by Berger et al.
31
The
force field parameters for cholesterol were taken from Ref. 32. In
the vacuum simulations, no cutoffs for the nonbonded interactions
were applied. The WALP20 peptide
33
was represented by the 43a2
parameter set of the GROMOS force field
11
and solvated with SPC
water.
34
The system was simulated within periodic boundary condi-
tions. Nonbonded interactions were calculated using a triple-range
cutoff scheme: interactions within 0.9 nm were calculated at every
time step from a pair list, which was updated every 20 fs. At these
Journal of Computational Chemistry DOI 10.1002/jcc
1336 Rzepiela et al.
•
Vol. 31, No. 6
•
Journal of Computational Chemistry
Figure 1. NCQ peptide. The coarse-grained structure is represented by
gray spheres; an ensemble of four atomistic structures is shown as col-
ored sticks. The N- and C-termini are capped with −NH
2
and −COOH
groups, respectively.
time steps, interactions between 0.9 and 1.5 nm were also calculated
and kept constant between updates. A reaction-field contribution
35
was added to the electrostatic interactions beyond this long-range
cutoff, with
r
= 62. In the simulations of bulk SPC water, a twin-
range cutoff scheme was used, with a single cutoff at 0.9 nm and a
pair list updated every 20 fs. Here, a long-range dispersion correction
was applied in addition to the reaction field.
Example Applications
NCQ Peptide
As a first application of our algorithm, we reconstructed atomistic
structures of an NCQ tripeptide in vacuum from CG structures. NCQ
is shown in Figure 1 and is composed of the amino acids asparagine,
cysteine, and glutamine. The application to NCQ has a twofold aim:
(i) to find optimal parameters for the annealing procedure and (ii)
to verify that proper ensembles are generated.
Parametrization of Annealing Procedure
Two sets of reconstruction simulations were carried out. The first
set (Set1) was initialized from a structure taken from an equilibrium
atomistic simulation at 300 K. To generate a reference CG struc-
ture for the definition of the restraints (see eq. (2)), the atomistic
structure was converted to its CG representation using the MAR-
TINI mapping. The second set of reconstruction simulations (Set2)
started from a snapshot taken from an equilibrium MARTINI CG
simulation at 300 K. For each combination of parameters, 1000
independent annealing simulations were carried out to generate an
ensemble of reconstructed structures.
First, we investigated how the properties of the ensemble of
reconstructed atomistic structures depend on two crucial parameters
of the SA: the total annealing time, t
tot
, and the initial temperature,
T
init
. In Figure 2A, the average potential energy of the ensemble
of final structures of Set1 is plotted as a function of the length
of the annealing simulation. As expected, E
pot
decreases with t
tot
and reaches a minimum at about −400.1 kJ mol
−1
after 60 ps
(standard deviation = 20.8 kJ mol
−1
, std error = 0.7 kJ mol
−1
).
For comparison, the average E
pot
obtained from a 400-ns equi-
librium atomistic MD simulation is −392.8 kJ mol
−1
(standard
deviation = 19.4 kJ mol
−1
, standard error = 0.8 kJ mol
−1
). For
short simulation times, the average energies are higher and the distri-
butions are broader due to a number of high-energy structures within
the ensemble. For longer simulation times, there are no such high-
energy structures. The corresponding average energies are lower
and the distributions are narrower, suggesting that the structures
within the ensemble are similar to each other. Indeed, as shown in
the inset, the rmsd with respect to the reference structure (i.e., the
structure taken from the equilibrium atomistic simulation, see ear-
lier) reaches its minimum value of 0.05 nm for simulation times
longer than 60 ps. To obtain such a converged ensemble, the initial
temperature T
init
needs to be chosen high enough to overcome the
relevant energy barriers during the simulation. Figure 2B shows that
low-energy ensembles are obtained for T
init
≥ 1300 K, which was
thus chosen for the simulations shown in Figure 2A.
Figure 2C shows how the potential energy of the reconstructed
ensemble depends on the force that couples the atomistic to the CG
system. For Set1, the potential energy of the reconstructed ensemble
depends only weakly on the restraining force constant (dashed line).
Very weak force constants do not limit the sampling to only the
lowest energy regions of the potential energy landscape, thus leading
to a slightly higher E
pot
. For Set2, i.e., the reference structure taken
from the equilibrium CG simulation, a lower force constant yields
ensembles with lower energies when compared with those obtained
at higher force constants (solid line). This observed strain energy
is due to the different regions of phase space sampled at the CG
and AA levels of resolution. The larger this mismatch, the higher
the strain. Thus, an atomistic system that is only weakly coupled to
the CG system is allowed to sample lower energy configurations,
which might be further away from the CG structure to which it is
restrained.
On the one hand, the aim is to generate an atomistic ensemble
that is close to its reference CG structure; on the other hand, strained
structures should of course be avoided. This can be achieved in a
two-step approach: first, an ensemble is generated using a rather
high restraining force constant during the SA. Then, in a second step,
the strain energy is released. Figure 2D shows the time evolution
of the potential energy (solid line) and the rmsd with respect to
the reference structure taken from the equilibrium CG simulation
(dashed line). First, a 60-ps SA was carried out, followed by 20-ps
equilibration at the final temperature of 300 K under the restraints.
Then, after 80 ps, the restraining potential was removed within a
time period of 10 ps. During this time, the system releases its strain
energy in the fast degrees of freedom by relaxing to a local minimum
(drop in E
pot
), which is structurally further away from the reference
structure (rise in rmsd). To release the additional strain energy in
slow degrees of freedom or in degrees of freedom that involve high
energy barriers requires more extensive sampling.
Journal of Computational Chemistry DOI 10.1002/jcc
