scispace - formally typeset
Search or ask a question
Journal ArticleDOI

Software News and Update Reconstruction of Atomistic Details from Coarse-Grained Structures

TL;DR: 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 is presented.
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 CO 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. (C) 2010 Wiley Periodicals, Inc. J Comput Chem 31 : 1333-1343, 2010

Summary (4 min read)

Introduction

  • Molecular dynamics (MD) simulations have been successfully used to investigate biomolecular systems, such as proteins, DNA, or lipid membranes. [1] [2] [3].
  • The main benefit of these models is their computational efficiency due to the reduced number of interaction sites.
  • Thus, tools are desirable 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.
  • The latter examples demonstrate the suitability of the reconstruction algorithm to condensed-phase biomolecular systems.
  • A short conclusive section ends this article.

Methods

  • The authors goal is to generate a low-energy atomistic (AA) ensemble that underlies its corresponding CG system.
  • Initially, AA particles are positioned close to their reference CG beads.
  • Finally, the restraints are gradually removed to yield a relaxed atomistic system.
  • 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 temperature is gradually decreased from a high starting value to the desired target temperature, thus allowing the system to rapidly cross energy barriers and find a low-energy minimum at the end of the simulation.
  • The crucial parameters are the starting temperature and the cooling rate; too low initial temperatures or too rapid cooling will not yield properly equilibrated final structures.
  • This is possible because the system optimizes in the course of the simulation, and thus large forces occur less frequently.
  • 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.

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.
  • Thus, every AA water molecule that moves out of the cutoff radius is driven back toward the center of mass.
  • This procedure is iteratively applied to every water molecule.

Implementation

  • The authors implemented the reconstruction algorithm into the program mdrun (v3.3), which is the main MD engine of the program package GROMACS.
  • 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.
  • 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.
  • Such an unwanted minimum can be, for example, the cis rotamer of an amide bond in a polypeptide backbone.
  • The user can define a threshold for the dihedral angle force constant and thereby select certain dihedrals; each selected dihedral can then be restrained to a desirable value, for example, the one observed in the X-ray crystal structure of a protein.

Simulation Details

  • All simulations were carried out using the GROMACS simulation package (version 3.3.1) 5 and the implementation of the reconstruction algorithm described in this work.
  • In the CG simulations, version 2 of the MARTINI forcefield [8] [9] [10] was used together with a 20-fs integration time step.
  • The other simulation parameters were set to the standard values described in the original publications.
  • 31 The force field parameters for cholesterol were taken from Ref. 32.
  • 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.

NCQ Peptide

  • 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.
  • First, the authors 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 .
  • 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 .
  • 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).
  • 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.

Properties of Reconstructed Ensembles

  • Next, the authors more closely characterized the properties of the reconstructed atomistic ensembles of the NCQ peptide.
  • The coarse amino acid-to-1 mapping yields an ensemble in which three dihedral states are populated (magenta curve): the reference state at DihA = −60 , the next Journal of Computational Chemistry DOI 10.1002/jcc Journal of Computational Chemistry DOI 10.1002/jcc minimum at +60 , and a third state at 180 .
  • The black solid line in Figure 3C shows that all three dihedral states are populated.
  • As described earlier, the authors have optimized the reconstruction procedure for the MARTINI mapping and found that it also works well for finer mapping schemes.
  • To summarize, the reconstructed atomistic ensembles correctly represent the possible dihedral states and energy minima.

DPPC and Cholesterol

  • The authors next goal was to investigate whether and how the parameters of the reconstruction procedure that were optimized for the NCQ tripeptide can be transferred to other molecules as well.
  • The reference structures for the definition of the restraints during the reconstruction simulations were taken from atomistic equilibrium MD simulations of the single molecules in vacuum; the CG representations were constructed via MARTINI mapping.
  • By this means, a proper ensemble of reconstructed cholesterol molecules was generated even within shorter annealing times and from lower initial temperatures (data not shown).
  • As shown in Figure 7 , the dipole moment distributions of the two systems are different: the correlation between individual dipoles of water molecules that are grouped together by the restraining potential yields a dipole distribution (dashed line) that is slightly shifted with respect to the one of free SPC (solid line).
  • Note that the simulations were carried out in the NVT ensemble, i.e., at constant box volume.

WALP Peptide in Solvated DPPC Bilayer

  • In the previous sections, the authors have tested and optimized a set of parameters for the restrained annealing simulations, which turned out to be well suited for a number of different biomolecules in the gas phase and for bulk water.
  • Finally, their aim is to demonstrate that their reconstruction procedure can also be applied to a real-life problem, i.e., a condensed-phase biomolecular system with several components.
  • Figure 8 shows a part of the simulation system before (B) and after (C) the reconstruction.
  • The same value was obtained from a free 10-ns atomistic simulation of the same system.
  • Another successful reconstruction was performed for a number of different lipid membranes in the presence of the antimicrobial peptide magainin, as described elsewhere.

Summary and Conclusions

  • The authors describe the implementation of an algorithm to reconstruct atomistic details from CG structures into GROMACS, a free and highly efficient MD program package.
  • Second, a SA MD procedure is used, during which the atomistic and CG systems are coupled via restraints.
  • The authors systematically tested and optimized the reconstruction procedure for single molecules in the gas phase, such as the NCQ tripeptide, the DPPC lipid, and cholesterol.
  • These reconstructions yielded proper low-energy atomistic ensembles underlying the CG representations.
  • Additionally, because the annealing is merely used to optimize low-energy structures, and not to calculate dynamic properties of the system, the mass of the hydrogen atoms can be increased, which would allow for a larger integration time step.

Did you find this useful? Give us your feedback

Figures (10)

Content maybe subject to copyright    Report

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

Citations
More filters
Journal ArticleDOI
TL;DR: The focus of this work is on the coarse-grained MARTINI force field, for which mapping definitions are written to allow conversion to and from the higher-resolution force fields GROMOS, CHARMM, and AMBER, and to andFrom a simplified three-bead lipid model.
Abstract: The conversion of coarse-grained to atomistic models is an important step in obtaining insight about atomistic scale processes from coarse-grained simulations. For this process, called backmapping or reverse transformation, several tools are available, but these commonly require libraries of molecule fragments or they are linked to a specific software package. In addition, the methods are usually restricted to specific molecules and to a specific force field. Here, we present an alternative method, consisting of geometric projection and subsequent force-field based relaxation. This method is designed to be simple and flexible, and offers a generic solution for resolution transformation. For simple systems, the conversion only requires a list of particle correspondences on the two levels of resolution. For special cases, such as nondefault protonation states of amino acids and virtual sites, a target particle list can be specified. The mapping uses simple building blocks, which list the particles on the di...

540 citations

Journal ArticleDOI
TL;DR: An overview of some of the more popular CG models used in biomolecular applications to date, focusing on models that retain chemical specificity, are provided.
Abstract: Computational modeling of biological systems is challenging because of the multitude of spatial and temporal scales involved. Replacing atomistic detail with lower resolution, coarse grained (CG), beads has opened the way to simulate large-scale biomolecular processes on time scales inaccessible to all-atom models. We provide an overview of some of the more popular CG models used in biomolecular applications to date, focusing on models that retain chemical specificity. A few state-of-the-art examples of protein folding, membrane protein gating and self-assembly, DNA hybridization, and modeling of carbohydrate fibers are used to illustrate the power and diversity of current CG modeling.

464 citations

Journal ArticleDOI
TL;DR: The state of the art in the field of realistic membrane simulations is reviewed and the current limitations and challenges ahead are discussed.
Abstract: Cell membranes contain a large variety of lipid types and are crowded with proteins, endowing them with the plasticity needed to fulfill their key roles in cell functioning. The compositional complexity of cellular membranes gives rise to a heterogeneous lateral organization, which is still poorly understood. Computational models, in particular molecular dynamics simulations and related techniques, have provided important insight into the organizational principles of cell membranes over the past decades. Now, we are witnessing a transition from simulations of simpler membrane models to multicomponent systems, culminating in realistic models of an increasing variety of cell types and organelles. Here, we review the state of the art in the field of realistic membrane simulations and discuss the current limitations and challenges ahead.

427 citations

Journal ArticleDOI
09 Apr 2012-ACS Nano
TL;DR: This study provides, for the first time, the self-assembly mechanism and the molecular organization of FF dipeptide nanostructures.
Abstract: Nanostructures, particularly those from peptide self-assemblies, have attracted great attention lately due to their potential applications in nanotemplating and nanotechnology. Recent experimental studies reported that diphenylalanine-based peptides can self-assemble into highly ordered nanostructures such as nanovesicles and nanotubes. However, the molecular mechanism of the self-organization of such well-defined nanoarchitectures remains elusive. In this study, we investigate the assembly pathway of 600 diphenylalanine (FF) peptides at different peptide concentrations by performing extensive coarse-grained molecular dynamics (MD) simulations. Based on forty 0.6-1.8 μs trajectories at 310 K starting from random configurations, we find that FF dipeptides not only spontaneously assemble into spherical vesicles and nanotubes, consistent with previous experiments, but also form new ordered nanoarchitectures, namely, planar bilayers and a rich variety of other shapes of vesicle-like structures including toroid, ellipsoid, discoid, and pot-shaped vesicles. The assembly pathways are concentration-dependent. At low peptide concentrations, the self-assembly involves the fusion of small vesicles and bilayers, whereas at high concentrations, it occurs through the formation of a bilayer first, followed by the bending and closure of the bilayer. Energetic analysis suggests that the formation of different nanostructures is a result of the delicate balance between peptide-peptide and peptide-water interactions. Our all-atom MD simulation shows that FF nanostructures are stabilized by a combination of T-shaped aromatic stacking, interpeptide head-to-tail hydrogen-bonding, and peptide-water hydrogen-bonding interactions. This study provides, for the first time to our knowledge, the self-assembly mechanism and the molecular organization of FF dipeptide nanostructures.

250 citations

Journal ArticleDOI
TL;DR: This paper presents a probabilistic analysis of the response of the immune system to E.G.D.’s, a type of “spatially aggregating” disease, which is known as E.g. “cell reprograming”.
Abstract: G. Andreś Cisneros,† Mikko Karttunen,‡ Pengyu Ren, and Celeste Sagui* †Department of Chemistry, Wayne State University, Detroit, Michigan 48202, United States ‡Department of Chemistry and Waterloo Institute for Nanotechnology, University of Waterloo, 200 University Avenue West, Waterloo, Ontario, Canada N2L 3G1 Department of Biomedical Engineering, BME 5.202M, The University of Texas at Austin, 1 University Station, C0800, Austin, Texas 78712-1062, United States Department of Physics, North Carolina State University, Raleigh, North Carolina 27695, United States

238 citations

References
More filters
Journal ArticleDOI
13 May 1983-Science
TL;DR: There is a deep and useful connection between statistical mechanics and multivariate or combinatorial optimization (finding the minimum of a given function depending on many parameters), and a detailed analogy with annealing in solids provides a framework for optimization of very large and complex systems.
Abstract: There is a deep and useful connection between statistical mechanics (the behavior of systems with many degrees of freedom in thermal equilibrium at a finite temperature) and multivariate or combinatorial optimization (finding the minimum of a given function depending on many parameters). A detailed analogy with annealing in solids provides a framework for optimization of the properties of very large and complex systems. This connection to statistical mechanics exposes new information and provides an unfamiliar perspective on traditional optimization problems and methods.

41,772 citations

Journal ArticleDOI
TL;DR: The dynamical steady-state probability density is found in an extended phase space with variables x, p/sub x/, V, epsilon-dot, and zeta, where the x are reduced distances and the two variables epsilus-dot andZeta act as thermodynamic friction coefficients.
Abstract: Nos\'e has modified Newtonian dynamics so as to reproduce both the canonical and the isothermal-isobaric probability densities in the phase space of an N-body system. He did this by scaling time (with s) and distance (with ${V}^{1/D}$ in D dimensions) through Lagrangian equations of motion. The dynamical equations describe the evolution of these two scaling variables and their two conjugate momenta ${p}_{s}$ and ${p}_{v}$. Here we develop a slightly different set of equations, free of time scaling. We find the dynamical steady-state probability density in an extended phase space with variables x, ${p}_{x}$, V, \ensuremath{\epsilon}\ifmmode \dot{}\else \.{}\fi{}, and \ensuremath{\zeta}, where the x are reduced distances and the two variables \ensuremath{\epsilon}\ifmmode \dot{}\else \.{}\fi{} and \ensuremath{\zeta} act as thermodynamic friction coefficients. We find that these friction coefficients have Gaussian distributions. From the distributions the extent of small-system non-Newtonian behavior can be estimated. We illustrate the dynamical equations by considering their application to the simplest possible case, a one-dimensional classical harmonic oscillator.

17,939 citations

Journal ArticleDOI
TL;DR: A set of simple and physically motivated criteria for secondary structure, programmed as a pattern‐recognition process of hydrogen‐bonded and geometrical features extracted from x‐ray coordinates is developed.
Abstract: For a successful analysis of the relation between amino acid sequence and protein structure, an unambiguous and physically meaningful definition of secondary structure is essential. We have developed a set of simple and physically motivated criteria for secondary structure, programmed as a pattern-recognition process of hydrogen-bonded and geometrical features extracted from x-ray coordinates. Cooperative secondary structure is recognized as repeats of the elementary hydrogen-bonding patterns “turn” and “bridge.” Repeating turns are “helices,” repeating bridges are “ladders,” connected ladders are “sheets.” Geometric structure is defined in terms of the concepts torsion and curvature of differential geometry. Local chain “chirality” is the torsional handedness of four consecutive Cα positions and is positive for right-handed helices and negative for ideal twisted β-sheets. Curved pieces are defined as “bends.” Solvent “exposure” is given as the number of water molecules in possible contact with a residue. The end result is a compilation of the primary structure, including SS bonds, secondary structure, and solvent exposure of 62 different globular proteins. The presentation is in linear form: strip graphs for an overall view and strip tables for the details of each of 10.925 residues. The dictionary is also available in computer-readable form for protein structure prediction work.

14,077 citations

Journal ArticleDOI
TL;DR: A new implementation of the molecular simulation toolkit GROMACS is presented which now both achieves extremely high performance on single processors from algorithmic optimizations and hand-coded routines and simultaneously scales very well on parallel machines.
Abstract: Molecular simulation is an extremely useful, but computationally very expensive tool for studies of chemical and biomolecular systems Here, we present a new implementation of our molecular simulation toolkit GROMACS which now both achieves extremely high performance on single processors from algorithmic optimizations and hand-coded routines and simultaneously scales very well on parallel machines The code encompasses a minimal-communication domain decomposition algorithm, full dynamic load balancing, a state-of-the-art parallel constraint solver, and efficient virtual site algorithms that allow removal of hydrogen atom degrees of freedom to enable integration time steps up to 5 fs for atomistic simulations also in parallel To improve the scaling properties of the common particle mesh Ewald electrostatics algorithms, we have in addition used a Multiple-Program, Multiple-Data approach, with separate node domains responsible for direct and reciprocal space interactions Not only does this combination of a

14,032 citations

Journal ArticleDOI
TL;DR: The software suite GROMACS (Groningen MAchine for Chemical Simulation) that was developed at the University of Groningen, The Netherlands, in the early 1990s is described, which is a very fast program for molecular dynamics simulation.
Abstract: This article describes the software suite GROMACS (Groningen MAchine for Chemical Simulation) that was developed at the University of Groningen, The Netherlands, in the early 1990s. The software, written in ANSI C, originates from a parallel hardware project, and is well suited for parallelization on processor clusters. By careful optimization of neighbor searching and of inner loop performance, GROMACS is a very fast program for molecular dynamics simulation. It does not have a force field of its own, but is compatible with GROMOS, OPLS, AMBER, and ENCAD force fields. In addition, it can handle polarizable shell models and flexible constraints. The program is versatile, as force routines can be added by the user, tabulated functions can be specified, and analyses can be easily customized. Nonequilibrium dynamics and free energy determinations are incorporated. Interfaces with popular quantum-chemical packages (MOPAC, GAMES-UK, GAUSSIAN) are provided to perform mixed MM/QM simulations. The package includes about 100 utility and analysis programs. GROMACS is in the public domain and distributed (with source code and documentation) under the GNU General Public License. It is maintained by a group of developers from the Universities of Groningen, Uppsala, and Stockholm, and the Max Planck Institute for Polymer Research in Mainz. Its Web site is http://www.gromacs.org.

13,116 citations

Frequently Asked Questions (16)
Q1. What are the contributions in "Reconstruction of atomistic details from coarse-grained structures" ?

The authors 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. 

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. 

because the annealing is merely used to optimize low-energy structures, and not to calculate dynamic properties of the system, the mass of the hydrogen atoms can be increased, which would allow for a larger integration time step. 

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. 

In about 8% of the cases, a structure was obtained in which one of the two stereocenters is inverted, whereas a structure with both stereocenters inverted occurred at a probability of only 2%. 

Because of the random initial placement of the atoms close to their reference CG beads (see Methods), there is a chance that, e.g., the methyl group bound to C10 is initially on the “wrong” side of the plane spanned by the other three bound groups. 

The distributions obtained with the 2- to-1 (dashed red curve) and 1-to-1 mapping (solid red curve) have even lower average energies and are very sharp. 

a system composed of the WALP20 transmembrane peptide embedded in a solvated DPPC bilayer was reconstructed from its CG representation. 

To generate a starting structure for the subsequent reconstruction simulations, a 40-ns CG simulation (MARTINI force field) of a system containing a WALP20 peptide embedded in a membrane bilayer composed of 112 DPPC lipids, solvated by 1186 water beads, was carried out at T = 323 K in the NPT ensemble (p = 1 bar). 

the authors investigated how the properties of the ensemble of reconstructed atomistic structures depend on two crucial parameters of the SA: the total annealing time, ttot , and the initial temperature, Tinit . 

a system composed of a WALP20 transmembrane peptide embedded in a solvated DPPC bilayer was transformed from a CG to an atomistic representation. 

A three-step approach is used to optimize low-energy ensembles: First, the atomistic particles are positioned close to their reference CG beads. 

By this means, the reconstruction method allows to check and validate the results and predictions obtained with CG models against atomistic models, thereby combining the efficiency of the former with the accuracy of the latter. 

Although the helical structure of the peptide is already fullyformed in about 20 ps, the potential energy of the system (dashed line in Fig. 9) indicates that also in this case, simulation times longer than about 60 ps are required to obtain a low-energy structure. 

The different mapping schemesapplied were an amino acid-to-1 mapping, i.e., each CG bead represents a complete amino acid; the MARTINI mapping, in which on an average four heavy atoms are mapped onto one CG bead; and a 2-to-1 mapping, i.e., each CG bead represents two heavy atoms. 

The second set of reconstruction simulations (Set2) started from a snapshot taken from an equilibrium MARTINI CG simulation at 300 K.