Software News and Update Reconstruction of Atomistic Details from Coarse-Grained Structures
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
Citations
540 citations
464 citations
427 citations
250 citations
238 citations
References
41,772 citations
17,939 citations
14,077 citations
14,032 citations
13,116 citations
Related Papers (5)
Frequently Asked Questions (16)
Q2. Why are the AA particles randomly placed around the corresponding CG beads?
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.
Q3. Why is the annealing used to optimize low-energy structures?
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.
Q4. How many nanoseconds can be used to simulate a molecule?
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.
Q5. What is the probability of a structure with both stereocenters inverted?
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%.
Q6. Why is the methyl group bound to C10 on the wrong side of the plane?
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.
Q7. What are the distributions obtained with the 2- to-1 mapping?
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.
Q8. What was the first step in the reconstruction of a CG system?
a system composed of the WALP20 transmembrane peptide embedded in a solvated DPPC bilayer was reconstructed from its CG representation.
Q9. How was the reconstruction of a WALP20 peptide performed?
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).
Q10. What are the properties of the ensemble of reconstructed atomistic structures?
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 .
Q11. What is the atomistic representation of a WALP20 peptide?
a system composed of a WALP20 transmembrane peptide embedded in a solvated DPPC bilayer was transformed from a CG to an atomistic representation.
Q12. What is the main method used to optimize low-energy ensembles?
A three-step approach is used to optimize low-energy ensembles: First, the atomistic particles are positioned close to their reference CG beads.
Q13. What is the purpose of the reconstruction method?
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.
Q14. How long is the simulation required to obtain a low-energy structure?
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.
Q15. What mapping scheme was used to map the dihedral angle of asparagine?
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.
Q16. What was the first set of reconstruction simulations?
The second set of reconstruction simulations (Set2) started from a snapshot taken from an equilibrium MARTINI CG simulation at 300 K.