Enhanced Sampling of Molecular Dynamics ... - ACS Publications

Subscriber access provided by UNIV OF DURHAM

Article

Enhanced Sampling of Molecular Dynamics Simulations of a Polyalanine Octapeptide: Effects of the Periodic Boundary Conditions on Peptide Conformation Kota Kasahara, Shun Sakuraba, and Ikuo Fukuda J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.7b10830 • Publication Date (Web): 13 Feb 2018 Downloaded from http://pubs.acs.org on February 22, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Enhanced Sampling of Molecular Dynamics Simulations of a Polyalanine Octapeptide: Effects of the Periodic Boundary Conditions on Peptide Conformation Kota Kasahara1,*, Shun Sakuraba2, and Ikuo Fukuda3,* 1. College of Life Sciences, Ritsumeikan University, 1-1-1 Noji-higashi, Kusatsu, Shiga 5258577, Japan 2. Graduate School of Frontier Sciences, The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8561, Japan 3. Institute for Protein Research, Osaka University, 3-2 Yamada-oka, Suita, Osaka 565-0871, Japan

* To whom correspondence should be addressed. Kota Kasahara, Tel: +81-77-566-1111 (4413); Email: [email protected]. Ikuo Fukuda, Tel: +81-6-6879-4311 (9286); Email: [email protected]

ACS Paragon Plus Environment

1

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 28

ABSTRACT:

We investigate the problem of artifacts caused by the periodic boundary conditions (PBC) used in molecular simulation studies. Despite the long history of simulations with PBCs, the existence of measurable artifacts originating from PBCs applied to inherently non-periodic physical systems remains controversial. Specifically, these artifacts appear as differences between simulations of the same system but with different simulation-cell sizes. Earlier studies have implied that, even in the simple case of a small model peptide in water, sampling inefficiency is a major obstacle to understanding these artifacts. In this study, we have resolved the sampling issue using the replica exchange molecular dynamics (REMD) enhanced-sampling method to explore PBC artifacts. Explicitly solvated zwitterionic polyalanine octapeptides with three different cubic-cells, having dimensions of L = 30, 40, and 50 Å, were investigated to elucidate the differences with 64 replica × 500 ns REMD simulations using the AMBER parm99SB force field. The differences among them were not large overall, and the results for the L = 30 and 40 Å simulations in the conformational free energy landscape were found to be very similar at room temperature. However, a small but statistically significant difference was seen for L = 50 Å. We observed that extended conformations were slightly overstabilized in the smaller systems. The origin of these artifacts is discussed by comparison to an electrostatic calculation method without PBCs.


2

Page 3 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Introduction Nowadays, molecular simulation is an important tool for understanding the physical, chemical, and biological features of macroscopic systems, especially with regard to the microscopic descriptions of the component particles and their interactions. In particular, electrostatic interactions are critical to maintaining physical structures, generating chemical properties, and performing biological functions. Unlike van der Waals interactions, electrostatic interactions are long-range, so the simple straight cutoff scheme does not work well. To handle this issue, periodic boundary conditions (PBCs) are typically used, where infinitely many copies of the target system confined in a unit (typically cubic) cell are arranged to fill a three-dimensional space. Here, the PBCs are not defined to be a simple toroidal boundary condition. Instead, the coulomb interactions between individual charges in the unit cell and all the charges in the all cells (unit and every copied cell) are taken into account, in principle, without an interaction cutoff. This can be achieved using the Ewald method1, where the long-range interactions are handled by calculations in Fourier space, assuming that the periodicity of the target molecular system is infinitely repeated. This method is a perfect match for inherently periodic target systems, such as a perfect crystal2, and it may also be applied to approximately periodic target systems, such as isotropic liquids3. This method has been successfully applied to a number of systems with such properties, and many physical properties have been revealed4,5. However, it is not trivial to determine whether this method applying PBCs is applicable to inherently non-periodic systems6,7. Specifically, imagine that we want to obtain conformational information about a single peptide in water (see Figure S1). If we use a PBC, we must consider the interactions between the original peptide (viz., the peptide in the unit


3


Page 4 of 28

cell) and other molecules, including not only the water molecules but also the (imaginary) peptides in the copied cells. Considering infinitely many such interactions between the original and imaginary peptides is not desirable for our purpose. In addition, although the interactions between the water molecules in the original cell and the peptides in the copied cells are not required, they are still accounted for in the PBC system. One may think that such unfavorable interactions from molecules in the neighboring copied cells would be weaker if the unit cell size were larger, but whether their strengths become greater or smaller after accounting for infinitely many cells is unknown. In summary, the finite-size effects arising from possible PBC artifacts must be clarified. In 2000, Weber et al.8 studied the conformation of an explicitly solvated alanine octapeptide (Ala8) using constant-temperature molecular dynamics (MD) to investigate the finite-size effect. They modeled a peptide with charged termini (zwitterionic) in three differently sized unit cells (the cell dimensions, L, were 20, 30, and 40 Å) subject to threedimensional (3D) PBCs. The simulations were performed with the GROMOS96 force field and SPC water molecules9. Based on the investigation of the root-mean-squared deviation (RMSD) and root-mean-squared fluctuation (RMSF) of the peptide in the 2 ns simulations at room temperature, with the help of continuum electrostatics analysis10–12, they concluded that the initial canonical α-helical structure is artificially stabilized with decreasing the cell size. The peptide in the large cell (L = 40 Å) quickly unfolded and became more compact. In contrast, Villarreal et al.13 observed no significant differences between the results of simulations using the three cell sizes (L = 20, 30, and 40 Å) involving the same zwitterionic peptide in explicit water in 20 ns MD simulations. The simulations were performed using the GROMOS87 force field and SPC/E water molecules14. For all cell sizes, the initial α-helical structure was lost over the 20 ns period, and the total duration of the


4

Page 5 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


termini contacts were comparable. They also conducted many short (0.5–1 ns) MD simulations and found a variety of trajectories that showed no significant bias in the unfolding velocity based on the cell size. Although Villarreal et al. discussed that the results by Weber et al. arose from insufficient sampling, attributed to the 2 ns simulation, they also admitted that a 20 ns simulation might not be sufficient to judge the existence of artifacts. They also suggested that investigations with larger system sizes and more accurate methods are necessary. Reif et al.15 also pointed out that a 20 ns simulation is still insufficient to achieve reversible folding– unfolding transitions in this system. In addition, while the work by Villarreal et al. addresses the kinetic bias by observing the difference in kinetics in unfolding process, it does not fully address the problem of a possible thermodynamic bias in the Ewald scheme with PBCs. Thus, Reif et al. concluded that the observations made by Villarreal et al. do not formally contradict the results of the previous continuum electrostatics analysis by Weber et al.8 For these reasons, it is still unclear whether the thermodynamic bias of the finite-size effect exists under the Ewald PBC, even for such a simple system. As suggested by these studies, one of the most promising ways to clarify the existence of artifacts is to conduct sufficiently long simulations using accurate methods and larger cells, followed by analysis of the thermodynamic properties from several viewpoints. Thus, we performed enhanced sampling for this system using 64 replica × 500 ns replica-exchange molecular dynamics (REMD) simulations. REMD is an established technique for efficiently exploring the conformational space of molecular systems16. While a canonical simulation at a constant temperature often becomes trapped into some local minima in the rugged energy landscape, a REMD simulation can escape from the local minima by exchanging temperatures (or coordinates) with other replica systems via Monte Carlo steps. Thus, the REMD method can


5


Page 6 of 28

prevent the artifacts arising from insufficient sampling. In addition, we employed more reliable simulation settings than those of previous reports, e.g., thermostatting the system via the Nosé-Hoover method and applying a longer cutoff distance for non-bonded interactions. Using these results, we analyzed the conformational ensemble of the peptide via the free energy (potential of mean force) based on the radius of gyration (Rg), the principal component analyses (PCA) of the interatomic distances, and the RMSD of the Cα atoms from that of an ideal helix structure, in addition to analyses of the secondary structure content. To further our analysis of simulation artifacts due to PBCs, we applied an alternative electrostatic calculation method without imposing PBCs. Specifically, to clarify whether the difference can be attributed to the PBCs, we applied the zero-multipole summation method (ZMM)17,18, which computes the electrostatic energy with a smooth cutoff function and, thus, does not require PBCs (Figure S1). The idea behind the ZMM is based on a key feature of the coulomb potential: it can have a positive or negative value, according to the charges of the two atoms, indicating a repulsive or attractive interaction. This is in contrast to the van der Waals and gravitational potentials which decay monotonically with increasing interatomic distance. This feature should strongly restrict the charge configuration (distribution). In particular, for a biological or condensed matter system, a charge configuration in which the electrostatic interactions are well canceled out is more stable. This simple observation of electrostatic neutrality, which is irrelevant to the PBC, is utilized in the ZMM and is achieved by using a mathematical formulation to define a pairwise function and the electrostatic energy. We applied the ZMM to measure any PBC artifacts. Note that Hünenberger et al.12 used the implicit continuum electrostatics model to detect or anticipate PBC artifacts. This is different from our approach via the ZMM, using explicit solvent models instead.


6

Page 7 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Other approaches to avoid PBC artifacts should also be addressed here. Caillol et al.19,20 have developed a method using three-dimensional hyper-sphere geometry, instead of the three-dimensional cubic PBC. Thus, neither the minimum image convention nor calculations in Fourier space are necessary. With a focus on biomolecular systems, Lin et al.21,22 combined the explicit and implicit solvent representations and have developed an image-based reaction field method, which treats a solute molecule surrounded with explicit solvent water molecules inside a spherical cavity and treats the solvent outside the cavity as a dielectric continuum generating reaction fields. Thus, it reduces the artificial periodicity and avoids the continuum approximation near the solute. More direct approaches, in the sense that the target system is surrounded with effective particles that play the role of a thermal bath, have also been developed23,24. To reveal the finite-size effects arising from PBC artifacts, we analyzed the results of enhanced-sampling REMD simulations of an explicitly solvated alanine octapeptide using the smooth particle mesh Ewald (PME) method, which is an accelerated version of the Ewald method using fast Fourier transforms (FFT) for the calculations in the Fourier space25. The influence of the size of the periodic boundary cell was evaluated by comparing the results from simulations of three system sizes, i.e., the cell dimensions were ca. 30, 40, and 50 Å. The simulation conditions are detailed in the Theoretical Methods section. As shown in the Results section, the current study with 192 µs of total simulation time revealed that the PBC artifacts in this system are not significant. However, some artifacts were observed, and these are discussed in detail in the Discussion section. Our concluding remarks are given in the Conclusion section.


7


Page 8 of 28

Methods Simulation Settings For all MD simulations, GROMACS26 ver. 4.6.3. was used. For the PME, a fourth-order interpolation function was applied with a maximum FFT grid spacing of 1.2 Å. The damping factor, α, was determined so that the parameter ewald_rtol (10-5 was used) in GROMACS is larger than the erfc(αrc)/rc, where erfc is the complementary error function, and rc is the cutoff radius of the real-space part evaluation. The AMBER parm99SB force field27 and the TIP3P water model28 were applied for the potential parameters. The cutoff radii for the Lennard-Jones (LJ) potential and the real-space part of the PME were both 12 Å. The bond lengths of hydrogen atoms were constrained by the LINCS algorithm29 and the numerical integration time step was 2.0 fs. The energies and atomic coordinates were recorded every 1.0 ps.

Simulation Systems The initial structure of the alanine octapeptide (Ala8) was generated by using MODELLER software30 without any template, and the termini were not capped but were charged (zwitterionic). The peptide chain was bathed in a water box, and the system was electrically neutral. Three kinds of cubic systems were studied: small, medium, and large, having cell dimensions (L) of 30, 40, and 50 Å, respectively. The individual systems were relaxed as follows. Energy minimizations with the steepest descent and conjugate gradient methods were first applied. Then, a simulation under NPT ensemble with the Berendsen barostat (1 atm) was performed for 2 ns. Position restraints for heavy atoms of the peptide and gradual heating from 10 to 300 K were applied for the first 1.0 ns of this NPT simulation. This relaxation process was performed with the PME. The structures of each system just after the relaxation are shown in Figure S2, and details of the systems are summarized in Table S1.


8

Page 9 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


REMD Simulations The REMD simulations were performed with three different sizes of the cubic cell, having dimensions of ca. 30, 40, and 50 Å (see Table S1), referred to as PME30, PME40, and PME50, respectively, from now on. The temperatures of the 64 parallel runs ranged from 290 to 557 K (Table S2). The temperatures were controlled by the Nosé-Hoover thermostat31,32 to produce reliable NVT canonical distributions at individual temperatures33. After 0.5 ns runs without any replica-exchange trials, REMD simulations were performed for 500 ns for each replica with the replica-exchange trials at every 500 steps (1 ps).

Remark on the Simulation Conditions The simulation conditions described here are not identical to those of Weber et al.8, where the GROMOS96 force field, SPC water molecules9, the Berendsen thermostat (NVT), and the P3M (particle particle-particle-mesh) electrostatic calculation method34 were used, and that of Villarreal et al.13, where the GROMOS87 force field, SPC/E water molecules14, Berendsen thermo/barostat (NPT), and the PME electrostatic calculation method were used. We have instead used the conditions mentioned above, which are considered to be standard or well tested. We did not assess the cubic cell with L = 20 Å, which was analyzed in the earlier studies because it is too small to handle the 12-Å cutoff in our simulation protocol. The 9 or 10-Å cutoff of the LJ interaction in these previous two studies seems insufficient to properly treat the nature of the van der Waals interaction in the PBC35. We believe that the conditions used here are sufficiently accurate to yield a standard and reliable molecular simulation at present. It should also be noted that the differences in simulation conditions from the previous studies do not affect the comparison to investigate the finite-size effects because the same conditions were used for the three target cell sizes in the comparison.


9


Page 10 of 28

ZMM Method Because the principle and derivation of the ZMM are detailed in Refs. [17,18], a quick review is provided in Supporting Information Text S1 for convenience. The accuracy and efficiency of the ZMM have been validated in applications to fundamental systems such as ionic systems17,36 and bulk water18,37, as well as heterogeneous bimolecular systems including DNA38 and membrane proteins39, and has also been used for motor protein40,41, DNA–protein complexes42,43, flexible molecular docking44–46, antibody modeling47, intrinsically disordered regions of protein48. For the implementation of the ZMM, we used our modified version of GROMACS (ver. 4.6.3.)26,49 (the modified version is freely available at https://github.com/shunsakuraba/gromacs). Currently, the zero-dipole condition with the damping parameter of α = 0 was used. The atom-based cutoff50 was used for the ZMM with the cutoff radius of 12 Å, which was the same cutoff radius as the LJ potential. The three different cell sizes were also used for the REMD simulations with the ZMM, instead of the PME, and the conditions L = ca. 30, 40, and 50 Å are referred to as ZMM30, ZMM40, and ZMM50, respectively (see Table S1).

Analysis Methods As a post-simulation analysis, the equilibrium ensembles at target temperatures were obtained using the multistate Bennett acceptance ratio method51 with the last 300 ns of the MD trajectories of seven temperatures centered at the target temperature. To characterize the structural properties of each snapshot, several quantities were analyzed. The radius of gyration (Rg) was calculated using the g_gyrate program implemented in the GROMACS suite. The RMSD of Cα atoms from the reference structure, which was an ideal helix structure, was analyzed and denoted as RMSDhelix. The PCA were performed for interatomic distances between all pairs of Cα atoms. The first and second principal component axes are denoted by PC1 and PC2, respectively. The secondary structure elements


10

Page 11 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and backbone dihedral angles were analyzed using the DSSP program52. The end-to-end distance, Dee, was defined as the distance between the Cα atoms of the N- and C-terminal residues.

Results Figure 1 shows the free energy landscapes at 300 K projected on the RMSDhelix–Rg space. All the simulations consistently generated a unimodal landscape, the basin of which was shallow and covered a wide range of Rg with a high RMSDhelix regime, and included a variety of denatured conformations (Figures 1A–C). The free energy landscape has a triangular geometry, which should mainly be due to elastic deformations to the helix structure and additional perturbations, as detailed in Supporting Information Text S2. One of the most stable structures formed an extended conformation (Figure 1D) compatible with polyproline II. A “shoulder” point of this basin was at a low Rg regime (Rg is ca. 4.7 Å and RMSDhelix is ca. 4 Å), and the structures in this region have closed termini, viz., a bent conformation with a β-turn stabilized by a salt-bridge between the charged termini (Figure 1E). The finite-size effect was not significant, and only a slight difference was observed around the shoulder. The helix conformations were not stable in these ensembles.


11


Page 12 of 28

Figure 1. The free energy landscapes in RMSDhelix–Rg space for A) PME30, B) PME40, and C) PME50. D) and E) indicate snapshots of conformations from the center of the basin and the shoulder point, respectively, which are indicated by dashed lines in panel A.

As a result of the PCA with respect to the interatomic distances, the contributions of PC1, PC2, and PC3 were 73.62%, 11.90%, and 4.222%, respectively, and, thus, we assumed that the energy landscape on the PC1–PC2 space captures the characteristics of the ensembles. PC1 resulted in a very high correlation with Rg, and the Pearson correlation coefficient (PCC) was 0.990. As shown in Figure 2, the landscapes on the PC1–PC2 space were very similar for all the three cell sizes and also show a unimodal distribution. The basin at PC1 ~10 and PC2 ~0 and the basin at PC1 ~-15 and PC2 ~2 correspond to the conformations D and E, respectively, shown in Figure 1.


12

Page 13 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2. The free energy landscapes on the PC1-PC2 plane for A) PME30, B) PME40, and C) PME50.

In all the conditions, only subtle populations of secondary structure elements were observed at room temperature (Figure 3). The α-helix contents were ca. 0.4–0.5%, and the 310-helix contents were ca. 3–4%. The turn and bend conformations were relatively highly populated, ca. 15–16% and 25–26%, respectively. There was a weak trend that larger cells showed a higher population of secondary structures.


13


Page 14 of 28

Figure 3. Secondary structure contents for A) α-helix, B) β-bridge, C) 310-helix, D) turn, and E) bend conformations, corresponding to H, B, G, T, and S in the definitions of the DSSP program, respectively. F) Contents without any secondary structure. The horizontal axis represents the amino acid position in Ala8. The blue, red, and green lines represent the results of the small (L = ca. 30 Å), medium (L = ca. 40 Å), and large (L = ca. 50 Å) cell sizes, respectively.

The free energy landscapes projected on the dihedral angles, viz., the φ-ψ plane (Figure 4), show that the most stable basin was located around the typical polyproline II parameters, φ = -75°, ψ = +145°, in all three conditions. Again, the three cell sizes generated similar conformational ensembles.


14

Page 15 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 4. The free energy landscape on the φ-ψ plane for A) PME30, B) PME40, and C) PME50. The horizontal and vertical axes denote the dihedral angles φ and ψ, respectively.

The free energy landscapes derived from the REMD simulations using the ZMM on the RMSDhelix-Rg plane and the PC1-PC2 plane are shown in Figures S3 and S4, respectively. The secondary structure contents for the ZMM is shown in Figure S5. The free energy landscape on the φ-ψ plane for the ZMM is shown in Figure S6. Significant deviations for the system size were not observed. Deviations between the PME and ZMM, especially in the free energy along Rg, are discussed in the next section. Note that the computational efficiency of the ZMM was greater than that of the PME (Table S3), mainly because the ZMM does not require calculations in the Fourier space. Both the results of the PME and ZMM simulations were based on the sufficient convergence of the sampling. In fact, the time course of Rg at room temperature shows that conformations in a wide range of Rg were steadily sampled between 200 and 500 ns (Figure S7). In addition, the averaged acceptance probabilities of the replica exchange trials were sufficiently high (at least 33.8%; Table S4) to efficiently traverse the temperature space. The average time required for a round trip in temperature space, which is the time required to


15


Page 16 of 28

travel from the lowest/highest temperature to the same temperature through the highest/lowest temperature, was 1.93–6.93 ns (Table S5), sufficiently short compared to the total simulation length. Discussion Stable Conformations Two reasons for the small populations of α-helix conformations were found. First, the helix stability in such a small model peptide is considered to be sensitive to the force field parameters, where we have applied the AMBER parm99SB. In fact, low helix contents were obtained for a polyalanine model in the AMBER parm99SB force field compared with other AMBER force fields, as reported in a benchmark using the reaction field method53. In contrast, an α-helix structure was presented as the most stable basin of an Ala10 model54 when using a hybrid AMBER force field, which is a weighted average of AMBER parm94 and parm99 parameters55, in a simulation with the multicanonical extendedensemble MD56 along with the cell-multipole expansion method57. Secondly, the conformational ensemble at equilibrium is influenced by the N- and C-terminal capping of the polypeptide, as well as the ion concentration of the solvent. Our simulation used the charged Ala residues as the termini without any capping groups and pure water as the solvent. Note that NMR and circular dichroism experiments have demonstrated that Ala7 capped by diaminobutyric acid and ornithine residues prefer the polyproline II conformation58. Although the capping condition was different, our simulation results agreed with this experimental outcome according to the free energy landscapes on the φ-ψ plane (Figure 4). This result indicates that our results possibly provide realistic conformations. Regardless of the agreement with experiment, the most important finding is that the three cell sizes generated similar conformational ensembles.


16

Page 17 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Finite-Size Effect Details A 1D free energy landscape with respect to Rg emphasizes the slight differences between the three conditions (Figures 5 and S8). The root-mean-squared error (RMSE; [PCC]) of the free energy along the Rg between PME50 and the other conditions are shown in Table 1. The artifacts arising from the simulation-length limitation can be considered to be negligible because the differences between the first and last 150 ns of the trajectories were smaller than the differences between the PME50 and individual methods, as indicated in Table S6 and Figure S9. Relatively large differences were observed at the “shoulder of the basin” or “the second basin” (i.e., Rg ~4.7 Å in Figure 5; see also Figure S8, which shows the error bars defined by two different statistical methods). Clearer differences were observed at a higher temperature (353 K; Figure 5B) than at room temperature; the RMSE [PCC] of the free energy along Rg with PME50 was larger than that in the room temperature, except PME40, as indicated in Table 1. The main difference, which appeared at the shoulder point in the low Rg regime, concerns the stability of compact conformations. Compact conformations tended to be stabilized in the large system with both electrostatic schemes (PME50 and ZMM50) and also in the medium size system with the ZMM (ZMM40). Assuming that a larger solvation box should result in fewer artifacts, the implied overstabilization of extended (polyproline II) conformations in the smaller boxes can be interpreted as the result of the finite-size effect. Regarding the system size convergence, the ZMM results imply the convergence at the size of L = 50 Å, since the difference of the results between L = 40 Å and L = 50 Å were very small, observed from Figures 5 and S3–6. In addition, Poisson-Boltzmann calculations59 to peptide conformations support the convergence, as detailed in Supporting Information Text S3. The overstabilization could be due to the attractive forces between the charged termini and the “copied” termini that are located in the nearest-neighbor image cell with opposite charge (dipole–dipole interactions). Rg is related to the dipole moment of the Ala8 molecule,


17


Page 18 of 28

as shown in Figure S10, where the correlation coefficient between them was 0.888. Figure 5 shows that the PME with medium cell-size (PME40) is similar to that with the small cell-size (PME30), while the ZMM with medium cell-size (ZMM 40) is similar to that with the large cell-size (ZMM50). In addition, Figure 5 shows that the ZMM with the small and medium cell-sizes (ZMM30 and ZMM40) at the room temperature produced similar landscapes to that of PME50 compared to the PME (PME30 and PME40). These tendencies were similar for a 1D free energy landscape with respect to the end-to-end distance, Dee (see Figure S11; the clearly seen second basin takes the structure represented in Figure 1E) and a 2D free energy landscape with respect to Dee and Rg (see Figure S12). These results suggest that the conformational properties derived via ZMM, which does not use PBCs, are less susceptible to the finite-size effect. In fact, the artificial dipole–dipole interactions can be off by the cutoff length (12 Å) of the ZMM when the cell dimension L is larger than about 35 Å because Dee is approximately 23 Å for the extended conformation. The temperature dependencies were also investigated, by taking advantage of generalized sampling methods such that the information for many temperatures can be obtained from one simulation. The free energy landscapes along Rg at 20 different temperatures are shown in Figure S13. As increasing the temperature, a population of compact conformations increases and the landscape becomes more flat. Figure S14 shows the difference of the free energy at the second basin from that of the PME50. Overall, while the differences generated by ZMM40 were smallest and less relevant to the temperature, the differences generated by PME30 were larger than PME40, ZMM30, and ZMM40 and increase as increasing the temperature. Higher temperature should provide diverse structures, as implied by the results of L = 50 Å, but relatively extended structures still take high populations in the smallest cell L = 30 Å.


18

Page 19 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5. The free energy landscape along Rg at A) 300 K and B) 353 K. The solid and dashed lines indicate the PME and ZMM, respectively. The red, green, and blue lines indicate the results of the small (L = 30 Å), medium (L = 40 Å), and large (L = 50 Å) cell sizes, respectively. The inset shows a close-up plot around the basin.

Table 1. Differences in the 1D Free Energy along Rg between PME50 and Others. temperature (K)

300

parameter

PCC

RMSE

353 353 PCC

RMSE

PME30 PME40 ZMM30 ZMM40 ZMM50

0.967 0.973 0.980 0.998 ,-0.998

0.0591 0.0535 0.0422 0.0183 0.0179

0.965 0.988 0.983 0.998 0.996

0.0609 0.0423 0.0443 0.0210 0.0238

*The landscapes were compared in the region with the value ≤ 1.0 kcal mol-1.

Comparison with Previous Studies Weber et al. found that, in the 1 ns MD simulations, the small solvation box (L = 20 Å) stabilizes the helix conformation of Ala8 and that the α-helical conformation was corrupted only in the large (L = 40 Å) system, where an extended structure was never observed although it may be stabilized by the dipole–dipole artifacts, especially for small systems8. In contrast, our results show subtle differences between the three cell sizes (L = ca. 30, 40, and 50 Å) for the PME, as well as the ZMM, and a very low population of the αhelical conformation, as well as the high population of the extended conformation, in all


19


Page 20 of 28

conditions. These discrepancies are mainly due to differences in the sampling simulation length and the differences in the force field. On the other hand, similarities between the results of Weber et al. and ours can also be found. First, the finite-size effect appeared, although it is exhibited as a small difference in the free energy at the shoulder point in our results. Secondly, the finite-size effect can be considered as an artificial dipole–dipole interaction, which is due to the conformations of relatively longer end-to-end distances. Thirdly, and intimately related, the relatively compact conformations are stable in larger cells. Our work also supports the previous hypothesis that the

compact

conformation

(previously

called

the

“circle

conformation”)

is

thermodynamically favored, which also, surprisingly, indicates that the continuum model can yield similar results based on smaller data. Our thermodynamic results also agree with the kinematical results based on 20 ns MD simulations by Villarreal et al.13. That is, no significant differences were observed between the results of simulations of different cell sizes, and the initial α-helical structure was not stable. Villarreal et al. also found that the structure folds to a hairpin structure in the L = 40 Å cell, which may also be compatible with our results. More diverse investigations could be carried out as follows. Because the dipole– dipole interaction artifact should be weakened by using capped polypeptide models, the evaluation of these models would be helpful. Reif et al.15 concluded that PBC artifacts are absent in the modeling of a β-heptapeptide in methanol solvent, although the screening of the electrostatic interactions should be weaker because of the lower dielectric environment, compared with the media with high dielectric constant, such as water, where the interactions are strongly damped. Our results do not directly extend to the behavior of the peptide in lower dielectric constant solvents, but a reinvestigation of this topic may lead us new insights. It has been shown60,61 that a solution with the lowest ionic strength (only the minimum


20

Page 21 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


number of ions were added to neutralize the system) minimizes the artifacts. While our current simulation system includes no ions, it would be interesting to investigate whether the use of a solvent with a high ion concentration generates larger differences between the cell sizes similar to the previous studies, or smaller differences due to the screening effects of the additional ions.

Conclusion Motivated by earlier, pivotal work8 and intimately relating prior12 and subsequent works13,15, the conformational ensembles of an explicitly solvated alanine octapeptide were compared among the three different sizes of periodic boundary cells. In all conditions, the simulated conformational ensembles support the absence of dramatic artifacts. All results were consistent in that a high population of the polyproline II conformation was observed. However, a slight difference was also observed concerning the stability of low Rg conformations (Figure 5). The PME30 and PME40 conditions stabilized high Rg extended conformations, which could be due to artificial dipole–dipole interactions between neighboring periodic cells. This study has been possible because of the significant increase in computer power and software development of the last decade. For example, Monticelli and Colombo discussed the sampling insufficiency on their constant temperature MD of an explicitly solvated 20-residue peptide for 500 ns simulations, which was the longest simulation in 200462. The present environment could, thus, help us to reconsider and reveal the artifacts arising from different electrostatic methods, not only the Ewald method. A complete review of the literature until 2009 discussing the artifacts of electrostatic methods is given in Ref. [15], concerning dipolar liquids, free solvated ions, solvated ion pairs, ionic systems, peptides,


21


Page 22 of 28

proteins, nucleic acids, and lipid membranes. For a more technical review of electrostatic methods, see Ref. [63], and a recent excellent review is given in Ref. [64]. Note that the problem addressed in the current study and our conclusions cannot be straightforwardly extended to other problems. However, this study is an important step forward, and, we hope that, just as earlier studies served as influential references (for example, Rocklin et al.65) for the investigation of the finite-size effects on the binding free energies of charged species, as well as for developing appropriate correction schemes to eliminate these effects, this study will help to reveal fundamental features of the electrostatic methods and to develop more efficient methodologies for molecular simulations.

ASSOCIATED CONTENT The following file is available free of charge. SupportingInformation.pdf includes Supporting Information Texts S1–S3, Figures S1–S14, Tables S1–S6.

ACKNOWLEDGMENTS KK, SS, and IF are funded by the Japan Society for Promotion of Sciences, Grant-in-Aid, Grant Numbers: JP16K18526, JP16K17778, and JP17K05143, respectively. This work was supported by the “Development of core technologies for innovative drug development based upon IT” from Japan Agency for Medical Research and development, AMED. The simulations were performed at the supercomputers of ACCMS, Kyoto University and Research Institute for Information Technology, Kyushu University. Computational resources were provided by the HPCI System Research Project (Project IDs: hp150065 and hp160165). The post-simulation analyses were performed with the


22

Page 23 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


supercomputer system provided by the National Institute of Genetics, Research Organization of Information and Systems, Japan. This work was performed in part under the Cooperative Research Program of Institute for Protein Research, Osaka University, CR-17-05. We

thank

Prof. Haruki Nakamura for fruitful discussions.

ABBREVIATIONS MD, molecular dynamics; PME, particle mesh Ewald; ZMM, zero-multipole summation; PCA, principal component analysis; FFT, fast Fourier transform; PCC, Pearson correlation coefficient; RMSD, root mean squared deviation; RMSF, root mean squared fluctuation; RMSE, root mean squared error; LJ, Lennard-Jones.

REFERENCES (1) (2) (3) (4)

(5)

(6) (7) (8)

(9)

Ewald, P. P. Die Berechnung Optischer Und Elektrostatischer Gitterpotentiale. Annalen der Physik 1921, 369 (3), 253–287. Tosi, M. P. Cohesion of Ionic Solids in the Born Model; Solid State Physics; Elsevier, 1964; Vol. 16, pp 1–120. Adams, D. J. On the Use of the Ewald Summation in Computer Simulation. J Chem Phys 1983, 78 (5), 2585–2590. Woodcock, L. V.; Singer, K. Thermodynamic and Structural Properties of Liquid Ionic Salts Obtained by Monte Carlo Computation. Part 1.—Potassium Chloride. J Chem Soc, Faraday Trans 1971, 67 (0), 12–30. Lewis, J. W. E.; Singer, K.; Woodcock, L. V. Thermodynamic and Structural Properties of Liquid Ionic Salts Obtained by Monte Carlo Computation. Part 2.— Eight Alkali Metal Halides. J Chem Soc, Faraday Trans 1975, 71 (0), 301–312. King, G.; Warshel, A. A Surface Constrained All‐Atom Solvent Model for Effective Simulations of Polar Solutions. J Chem Phys 1989, 91 (6), 3647–3661. Sagui, C.; Darden, T. A. Molecular Dynamics Simulations of Biomolecules: LongRange Electrostatic Effects. Annu Rev Biophys Biomol Struct 1999, 28 (1), 155–179. Weber, W.; Hünenberger, P. H.; McCammon, J. A. Molecular Dynamics Simulations of a Polyalanine Octapeptide Under Ewald Boundary Conditions: Influence of Artificial Periodicity on Peptide Conformation. J Phys Chem B 2000, 104 (15), 3668–3675. Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. Interaction Models for Water in Relation to Protein Hydration. In Intermolecular Forces; The Jerusalem Symposia on Quantum Chemistry and Biochemistry; Springer, Dordrecht, 1981; Vol. 14, pp 331–342.


23


(10) (11) (12)

(13)

(14) (15)

(16) (17) (18)

(19)

(20)

(21)

(22)

(23)

(24)

(25) (26)

(27)

Page 24 of 28

Harvey, S. C. Treatment of Electrostatic Effects in Macromolecular Modeling. Proteins 1989, 5 (1), 78–92. Warshel, A.; Papazyan, A. Electrostatic Effects in Macromolecules: Fundamental Concepts and Practical Modeling. Curr Opin Struct Biol 1998, 8 (2), 211–217. Hünenberger, P. H.; McCammon, J. A. Ewald Artifacts in Computer Simulations of Ionic Solvation and Ion–Ion Interaction: a Continuum Electrostatics Study. J Chem Phys 1999, 110 (4), 1856–1872. Villarreal, M. A.; Montich, G. G. On the Ewald Artifacts in Computer Simulations. the Test-Case of the Octaalanine Peptide with Charged Termini. J Biomol Struct Dyn 2005, 23 (2), 135–142. Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The Missing Term in Effective Pair Potentials. J Phys Chem 1987, 91 (24), 6269–6271. Reif, M. M.; Kräutler, V.; Kastenholz, M. A.; Daura, X.; Hünenberger, P. H. Molecular Dynamics Simulations of a Reversibly Folding Β-Heptapeptide in Methanol: Influence of the Treatment of Long-Range Electrostatic Interactions. J Phys Chem B 2009, 113 (10), 3112–3128. Sugita, Y.; Okamoto, Y. Replica-Exchange Molecular Dynamics Method for Protein Folding. Chem Phys Lett 1999, 314 (1), 141–151. Fukuda, I. Zero-Multipole Summation Method for Efficiently Estimating Electrostatic Interactions in Molecular System. J Chem Phys 2013, 139 (17), 174107. Fukuda, I.; Kamiya, N.; Nakamura, H. The Zero-Multipole Summation Method for Estimating Electrostatic Interactions in Molecular Dynamics: Analysis of the Accuracy and Application to Liquid Systems. J Chem Phys 2014, 140 (19), 194307. Caillol, J. M.; Levesque, D. Numerical Simulations of Homogeneous and Inhomogeneous Ionic Systems: an Efficient Alternative to the Ewald Method. J Chem Phys 1991, 94 (1), 597–607. Caillol, J. M. Numerical Simulations of Coulomb Systems: a Comparison Between Hyperspherical and Periodic Boundary Conditions. J Chem Phys 1999, 111 (14), 6528–6537. Lin, Y.; Baumketner, A.; Deng, S.; Xu, Z.; Jacobs, D.; Cai, W. An Image-Based Reaction Field Method for Electrostatic Interactions in Molecular Dynamics Simulations of Aqueous Solutions. J Chem Phys 2009, 131 (15), 154103. Lin, Y.; Baumketner, A.; Song, W.; Deng, S.; Jacobs, D.; Cai, W. Ionic Solvation Studied by Image-Charge Reaction Field Method. J Chem Phys 2011, 134 (4), 044105. Brünger, A.; Brooks, C. L.; Karplus, M. Stochastic Boundary Conditions for Molecular Dynamics Simulations of ST2 Water. Chem Phys Lett 1984, 105 (5), 495– 500. Li, X.; E, W. Variational Boundary Conditions for Molecular Dynamics Simulations of Crystalline Solids at Finite Temperature: Treatment of the Thermal Bath. Phys Rev B 2007, 76 (10), 104107. Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J Chem Phys 1995, 103 (19), 8577–8593. Pronk, S.; Páll, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D.; et al. GROMACS 4.5: a HighThroughput and Highly Parallel Open Source Molecular Simulation Toolkit. Bioinformatics 2013, 29 (7), 845–854. Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Comparison of Multiple Amber Force Fields and Development of Improved Protein Backbone Parameters. Proteins 2006, 65 (3), 712–725.


24

Page 25 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(28)

(29)

(30)

(31) (32) (33) (34) (35)

(36)

(37)

(38)

(39)

(40)

(41)

(42)

(43)

(44)

Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J Chem Phys 1983, 79 (2), 926–935. Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: a Linear Constraint Solver for Molecular Simulations. J Comput Chem 1997, 18 (12), 1463– 1472. Eswar, N.; Webb, B.; Marti-Renom, M. A.; Madhusudhan, M. S.; Eramian, D.; Shen, M.-Y.; Pieper, U.; Sali, A. Comparative Protein Structure Modeling Using Modeller. Curr Protoc Bioinformatics 2006, 5.6.1–5.6.30. Nosé, S. A Molecular Dynamics Method for Simulations in the Canonical Ensemble. Mol Phys 1984, 52 (2), 255–268. Hoover, W. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys Rev, A 1985, 31 (3), 1695–1697. Morishita, T. Fluctuation Formulas in Molecular-Dynamics Simulations with the Weak Coupling Heat Bath. J Chem Phys 2000, 113 (8), 2976–2982. Computer Simulation Using Particles; Institute of Physics Publishing: Bristol, 1988. Piana, S.; Lindorff-Larsen, K.; Dirks, R. M.; Salmon, J. K.; Dror, R. O.; Shaw, D. E. Evaluating the Effects of Cutoffs and Treatment of Long-Range Electrostatics in Protein Folding Simulations. PLoS ONE 2012, 7 (6), e39918. Fukuda, I.; Yonezawa, Y.; Nakamura, H. Molecular Dynamics Scheme for Precise Estimation of Electrostatic Interaction via Zero-Dipole Summation Principle. J Chem Phys 2011, 134 (16), 164107. Fukuda, I.; Kamiya, N.; Yonezawa, Y.; Nakamura, H. Simple and Accurate Scheme to Compute Electrostatic Interaction: Zero-Dipole Summation Technique for Molecular System and Application to Bulk Water. J Chem Phys 2012, 137 (5), 054314. Arakawa, T.; Kamiya, N.; Nakamura, H.; Fukuda, I. Molecular Dynamics Simulations of Double-Stranded DNA in an Explicit Solvent Model with the ZeroDipole Summation Method. PLoS ONE 2013, 8 (10), e76606. Kamiya, N.; Fukuda, I.; Nakamura, H. Application of Zero-Dipole Summation Method to Molecular Dynamics Simulations of a Membrane Protein System. Chem Phys Lett 2013, 568-569, 26–32. Nishikawa, Y.; Oyama, T.; Kamiya, N.; Kon, T.; Toyoshima, Y. Y.; Nakamura, H.; Kurisu, G. Structure of the Entire Stalk Region of the Dynein Motor Domain. J Mol Biol 2014, 426 (19), 3232–3245. Kamiya, N.; Mashimo, T.; Takano, Y.; Kon, T.; Kurisu, G.; Nakamura, H. Elastic Properties of Dynein Motor Domain Obtained From All-Atom Molecular Dynamics Simulations. Protein Eng Des Sel 2016, 29 (8), 317–325. Kasahara, K.; Fukuda, I.; Nakamura, H. A Novel Approach of Dynamic Cross Correlation Analysis on Molecular Dynamics Simulations and Its Application to Ets1 Dimer-DNA Complex. PLoS ONE 2014, 9 (11), e112419. Kasahara, K.; Shiina, M.; Fukuda, I.; Ogata, K.; Nakamura, H. Molecular Mechanisms of Cooperative Binding of Transcription Factors Runx1–CBFβ–Ets1 on the TCRα Gene Enhancer. PLoS ONE 2017, 12 (2), e0172654. Higo, J.; Dasgupta, B.; Mashimo, T.; Kasahara, K.; Fukunishi, Y.; Nakamura, H. Virtual-System-Coupled Adaptive Umbrella Sampling to Compute Free-Energy Landscape for Flexible Molecular Docking. J Comput Chem 2015, 36 (20), 1489– 1501.


25


(45)

(46)

(47) (48)

(49)

(50) (51) (52)

(53)

(54)

(55)

(56)

(57)

(58)

(59) (60)

(61)

Page 26 of 28

Dasgupta, B.; Nakamura, H.; Higo, J. Flexible Binding Simulation by a Novel and Improved Version of Virtual-System Coupled Adaptive Umbrella Sampling. Chem Phys Lett 2016, 662, 327–332. Bekker, G.-J.; Kamiya, N.; Araki, M.; Fukuda, I.; Okuno, Y.; Nakamura, H. Accurate Prediction of Complex Structure and Affinity for a Flexible Protein Receptor and Its Inhibitor. J Chem Theory Comput 2017, 13 (6), 2389–2399. Nishigami, H.; Kamiya, N.; Nakamura, H. Revisiting Antibody Modeling Assessment for CDR-H3 Loop. Protein Eng Des Sel 2016, 29 (11), 477–484. Iida, S.; Mashimo, T.; Kurosawa, T.; Hojo, H.; Muta, H.; Goto, Y.; Fukunishi, Y.; Nakamura, H.; Higo, J. Variation of Free‐Energy Landscape of the P53 C‐Terminal Domain Induced by Acetylation: Enhanced Conformational Sampling. J Comput Chem 2016, 37 (31), 2687–2700. Sakuraba, S.; Fukuda, I. Performance Evaluation of the Zero-Multipole Summation Method in Modern Molecular Dynamics Software. arXiv 2017, https://arxiv.org/abs/1704.07071. Fukuda, I.; Nakamura, H. Non-Ewald Methods: Theory and Applications to Molecular Systems. Biophys Rev 2012, 4 (3), 161–170. Shirts, M. R.; Chodera, J. D. Statistically Optimal Analysis of Samples From Multiple Equilibrium States. J Chem Phys 2008, 129 (12), 124105. Kabsch, W.; Sander, C. Dictionary of Protein Secondary Structure: Pattern Recognition of Hydrogen‐Bonded and Geometrical Features. Biopolymers 1983, 22 (12), 2577–2637. Thompson, E. J.; DePaul, A. J.; Patel, S. S.; Sorin, E. J. Evaluating Molecular Mechanical Potentials for Helical Peptides and Proteins. PLoS ONE 2010, 5 (4), e10056. Ikebe, J.; Umezawa, K.; Kamiya, N.; Sugihara, T.; Yonezawa, Y.; Takano, Y.; Nakamura, H.; Higo, J. Theory for Trivial Trajectory Parallelization of Multicanonical Molecular Dynamics and Application to a Polypeptide in Water. J Comput Chem 2010, 32 (7), 1286–1297. Kamiya, N.; Watanabe, Y. S.; Ono, S.; Higo, J. AMBER-Based Hybrid Force Field for Conformational Sampling of Polypeptides. Chem Phys Lett 2005, 401 (1-3), 312– 317. Nakajima, N.; Nakamura, H.; Kidera, A. Multicanonical Ensemble Generated by Molecular Dynamics Simulation for Enhanced Conformational Sampling of Peptides. J Phys Chem B 1997, 101 (5), 817–824. Ding, H.-Q.; Karasawa, N.; Goddard, W. A. Atomic Level Simulations on a Million Particles: the Cell Multipole Method for Coulomb and London Nonbond Interactions. J Chem Phys 1992, 97 (6), 4309–4315. Shi, Z.; Olson, C. A.; Rose, G. D.; Baldwin, R. L.; Kallenbach, N. R. Polyproline II Structure in a Sequence of Seven Alanine Residues. Proc. Natl. Acad. Sci. U.S.A. 2002, 99 (14), 9190–9195. Lu, Q.; Luo, R. A Poisson–Boltzmann Dynamics Method with Nonperiodic Boundary Condition. J Chem Phys 2003, 119 (21), 11035–11047. Bergdorf, M.; Peter, C.; Hünenberger, P. H. Influence of Cut-Off Truncation and Artificial Periodicity of Electrostatic Interactions in Molecular Simulations of Solvated Ions: a Continuum Electrostatics Study. J Chem Phys 2003, 119 (17), 9129–9144. Kastenholz, M. A.; Hünenberger, P. H. Influence of Artificial Periodicity and Ionic Strength in Molecular Dynamics Simulations of Charged Biomolecules Employing Lattice-Sum Methods. J Phys Chem B 2004, 108 (2), 774–788.


26

Page 27 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(62)

(63) (64) (65)

Monticelli, L.; Colombo, G. The Influence of Simulation Conditions in Molecular Dynamics Investigations of Model Β-Sheet Peptides. Theor Chem Acc 2004, 112 (3), 145–157. Koehl, P. Electrostatics Calculations: Latest Methodological Advances. Curr Opin Struct Biol 2006, 16 (2), 142–151. Cisneros, G. A.; Karttunen, M.; Ren, P.; Sagui, C. Classical Electrostatics for Biomolecular Simulations. Chem. Rev. 2014, 114 (1), 779–814. Rocklin, G. J.; Mobley, D. L.; Dill, K. A.; Hünenberger, P. H. Calculating the Binding Free Energies of Charged Species Based on Explicit-Solvent Simulations Employing Lattice-Sum Methods: an Accurate Correction Scheme for Electrostatic Finite-Size Effects. J Chem Phys 2013, 139 (18), 184103.


27


Page 28 of 28

TOC Graphic


28

Enhanced Sampling of Molecular Dynamics ... - ACS Publications

Recommend Documents