Introduction

The correct treatment of long-range electrostatic interactions, which is a sticking point for biomolecular simulations, has been the focus of modelling methods for at least the past decade. The discussion on the advantages and disadvantages of various approaches has continued (1, 2, 3, 4) since 1995, when the first molecular dynamics simulations with the particle mesh Ewald technique (PME) appeared. Although the use of PME is essential for nucleic acids (5), for proteins and peptides, which are not such highly charged molecules as polynucleotides, the choice of technique may be subject to question.

Electrostatic interactions are treated in molecular mechanics by pairwise calculation using the Coulomb law. However, with the increasing size of the systems involved, and especially by the inclusion of explicit water molecules, the calculation of all possible interactions, if we still use this method, would take more than 90% of computational time and would make simulations unfeasible. One possibility to overcome this problem is to truncate the forces at a certain distance called the cutoff. It is well known that this method introduces significant nonphysical effects especially at and also beyond the cutoff distance (6). Several methods of avoiding these artifacts (e.g., shifting and switching functions) have been developed (7, 8).

Substantial progress has occurred with the introduction of the particle mesh Ewald technique (9, 10), which is able to run stable simulations even on highly charged nucleic acid systems. This method, which was originally developed for the analysis of periodic crystalline systems (11), calculates electrostatic potential beyond chosen cutoff distance by the solution of a Poisson-Boltzmann equation for a solute in a continuum solvent bath. The electrostatic potential within cutoff distance is calculated as Coulomb forces.

In the systems under periodic boundary conditions (PBC), a combination of these algorithms provides us with a good approximation of environmental influence on the solute. Pairwise interactions of atoms separated by a distance that is less than cutoff are treated by exact calculation and the forces beyond the cutoff radius are approximated by the PME method. However, what seems unavoidable for nucleic acids may not be necessary for proteins and peptides. Using PME usually results in stable trajectories (1, 2), but the question remains if the stability of the trajectory is equivalent to its correctness (4) (despite the fact that the term ?corectness of the trajectory? itself could be rather speculative). It has been shown that in the case of proteins PME may bring unnatural bindings and overstabilization of the system (12). We note that the solution of the above question is not only important for computer simulations themselves, but also for correct refinement and interpretation of experimental data, especially that obtained by NMR spectroscopy and X-ray crystallography.

In the present study, we wish to compare both PME and cutoff methods. We have selected a short peptide, endothelin-1, for which structural data from both X-ray crystallography and NMR spectroscopy is available. The peptide hormone endothelin-1 (13, 14) belongs to the family of vasoconstrictors. It is composed of 21 residues (see Fig. 1). Even if it is partially constrained by two disulfide bridges, it is still flexible enough to exhibit significant conformational variability.

Figure 1: Scheme of amino acid sequence and disulfide bridging pattern for endothelin-1. Charged residues are denoted by + or -.

Method

Abbreviations

Four trajectories were run with the protocol described below, two starting with the NMR structure (13), two with the X-ray structure (14). Two methods of treating electrostatics were used, PME and cutoff. At the end of our study we decided to confirm our conclusions also by two trajectories with cutoff value of 15 Å. A list of trajectories and abbreviations is given in Table I.



Starting Geometries

Two experimentally obtained structures were downloaded from the PDB database (15), the NMR structure (1edp) (13) and the X-ray structure (1edn) (14) obtained with the resolution of 2.18 Å. Two of the 21 residues are considered as positively charged in a water environment (Cys1, Lys9), four of them are negatively charged (Asp8, Glu10, Asp18, Trp21), resulting in the overall charge of -2 for the whole molecule. Two disulfide bridges (Cys1-Cys15 and Cys3-Cys11) give the molecule the shape of a bicyclic structure with the tail corresponding to the C-end (further called the C-tail). Several secondary structure elements have been found in both experimental structures. Following the Kabsch and Sanders approach (16), there is an ?α-helix? in the 9-15 region of the NMR structure and a ?3₁₀-helix? in the 8-10 region of the X-ray structure. Several additional β- or γ-turns were found in both structures. The difference between the two structures is rather large as described in terms of the secondary structure elements. Although the X-ray study determined all residues of the peptide (in contrast to the NMR structure, in which four residues at the C-end were missing), it did not show the helical arrangement in the 9-15 region. It can therefore be stated that both experimental structures represent different parts of the conformational space and that they can serve as good starting points for conformational studies.

Initial Model Building: As the experimental structures were not fully determined, they had to be completed before the simulations were started. The X-ray structure required the addition of hydrogen atoms, which was done by the protonate module of the AMBER software package. The missing residues in the NMR structure (residues 18-21) were added by the xleap module of AMBER. This added part was then minimized while the rest of the molecule was heavily restrained. The above structures were used as input files for molecular dynamics simulations in the explicit solvent.

Molecular Dynamics: Molecular dynamics simulations were performed with AMBER-6.0 (17) software. The all-atom Cornell et al. force field (18), parm94, was used. The system was neutralized by adding two Na+ ions, which were positioned using the xleap module of AMBER. The solute was then embedded in a TIP3P water box (19), which extended approximately 10 Å in each direction from the solute. The number of water molecules was 3038 for the NMR and 2409 for the X-ray structure. All simulations were run with the SHAKE algorithm (20) (with a tolerance of 0.0005 Å) to constrain covalent bonds involving hydrogens, with periodic boundary conditions, a 2-fs time step, and a temperature of 300 K (Berendsen temperature coupling algorithm [with time constant of 0.2 ps (21)].

As specified above, the main attention was paid to long-range electrostatic interactions representation. Two simulations for each starting structure were performed at first -- with a cutoff distance of 10 Å and with a PME representation with a direct sum cutoff of 10 Å. After analysis of these four trajectories we decided to run an additional trajectory with cutoff distance of 15 Å on each starting structure. The common part of the equilibration procedure started with a 1000 step minimization with the solute constrained (solvent and ions allowed to move). The first step was followed by molecular dynamics (MD) at 300 K for 25 ps with the solute position fixed. Equilibration continued with 1000 steps of minimization, followed by a 3 ps MD, with harmonic restraints of 25 kcal/(mol·Å²) on the peptide in order to relax water molecules around the solute. This equilibration was followed by 5 rounds of 1000 step minimizations where the solute restraints were reduced by 5 kcal/(mol·Å:²) in each round. Finally, a 20 ps unrestrained MD simulation was performed, which differed in its treatment of long-range electrostatic interactions (PME or cutoff). During that phase, the systems were heated from 100 K to 300 K in the initial 2 ps and equilibrated for 18 ps. Then, production runs were initiated with equivalent simulation protocols without any restraints.

Standard analysis tools were used to process the trajectories, namely the carnal and ptraj modules of the AMBER package. The DSSP program (16) was used to determine secondary structure elements.

Free Energy Calculations: The total energy of a solvated molecule (22), as calculated by the MM-PBSA approach, can be written as E_tot = E_vac + ΔG_solv. Here E_vac represents a molecule?s energy in vacuum (gas-phase or E_gas) and ΔG_solv is the free energy of transferring the molecule from the vacuum into the solvent, i.e., solvation free energy. This means that E_tot does not include solute entropy. Solute entropic contributions were estimated using a subset of these structures (at 200 ps intervals) based on a harmonic approximation to the normal modes, and standard (quantum) formulas at 300 K. Only a subset of the configurations was minimized due to the relatively high computational cost, as described in our previous work (23, 24). Although this minimization causes the structures to distort somewhat, the calculated entropy is not tremendously sensitive to small structural perturbations.

Free energy was calculated using a modified MM-PBSA script distributed with AMBER 6.0. The Cornell et al. (parm94) charge set, PARSE van der Waals radii (25) and dielectric constants of 80 for external dielectric, and 1 for internal dielectric (as is consistent with the utilization of a nonpolarizable force field) (26) were used. For free energy of solvation applies: ΔG_solv = ΔG_PB + ΔG_nonpolar. The electrostatic component of the solvation ΔG_PB was estimated with the Poisson-Boltzmann electrostatic continuum method using the program Delphi II (27). The hydrophobic component of the solvation free energy ΔG_nonpolar was calculated from a surface area dependent term using Sanner?s MSMS program for solvent accessible surface area (SASA) calculation (28). The total energy was calculated over all the trajectories with the frequency of 10 snapshots and averaged over parts of interest.

Secondary Structure Assignment: Secondary structure elements were determined using the Kabsch and Sander?s (16) DSSP program.

Results and Discussion

As we stated before (23, 29), the development of MD simulation trajectories depends on starting geometry. This was confirmed also by the present study where the simulations started from two very different experimental conformations.

The RMSD Development

RMS deviations were calculated for heavy atoms as well as for backbone C, N, CA atoms only, and for the whole molecule as well as for C-tail only (residues 16-21). Averaged RMSD values are introduced in Table II, and RMSD development of the backbone atoms is pictured in Figure 2. It can be seen that three trajectories (CX, PX, PN) exhibit the same trend.




A

B
Figure 2: RMS deviations from starting structures evaluated for backbone C, N, and CA atoms during the trajectories ? part A: i) CX, ii) PX, iii) CN, iv) PN; part B: i) C15X, ii) C15N. Calculated for entire molecules (residues 1-21, gray) or for the C-tail only (residues 16-21, black). Averaged values are shown in Table II.

The average RMSD values for the whole molecule are very similar to the RMSD values for the C-tail residues (it applies for heavy atoms as well as for backbone atoms). For example, for the CX trajectory, the heavy atom average RMSD values for the entire molecule and for the C-tail are 5.25 Å and 5.31 Å, respectively, while for the backbone atoms, average RMSD values are 3.57 Å (entire molecule) and 3.00 Å (C-tail only). The situation is pictured in Figure 2. The overlapping curves imply that the C-tail is the most flexible part of the molecule in CX, PX, and PN trajectories. However, the situation is different in the CN trajectory. Here, quite a large difference between the RMS deviation of the whole molecule, and that of the C-tail residues only, is seen. It shows that not only the C-tail but also the remaining part of the molecule is very flexible in this trajectory. The CN trajectory appears to be exceptional among those remaining not only for its highest RMSD values, but also for their largest fluctuations, as can be seen from the B-factors displayed in Figure 3. It is interesting that the tail region is rather rigid in C15X trajectory, so its overall flexibility is given mainly by bicyclic part (as seen from low C-tail RMSD). It can be seen, for example, that this trajectory exhibits the lowest B-factors values for the tail residues -- see Figure 3. However, C15N exhibits similar tendencies in B-factors as CN.

Figure 3: B-factors evaluated over the entire trajectories. The B-factors were calculated as described in the literature (37) using the equation B_i = (8π²/3)R_i², where B_i and R_i are the B-factor and RMSD of atom i, respectively. The RMSD values were calculated by the ptraj module of AMBER. Compared to the remaining trajectories, the CN trajectory exhibits exceptionally high B-factors in the bicyclic part of the molecule (residues 1-15). This part of the molecule significantly influences the flexibility of the whole system. The remaining trajectories (CX, PX, PN) display a high increase of B-factors in the C-tail region when compared to the bicyclic part which means that the bicyclic region is a rather rigid part of the molecule in these simulations. This applies also for C15N trajectory for which the B-factors exhibit the same behavior as those for the CN trajectory. In contrary the C15X tail appears to be more rigid in comparison with the rest of the molecule. It corresponds with RMS deviations for this simulation, as mentioned before. Both PME trajectories have low B-factors in comparison to both cutoff trajectories.

Free Energy Calculations

With the help of the MM-PBSA method (25, 30, 31, 32) we have calculated the free energy for all MD simulations as described in the Method section. The averaged results are summarized in Table III and shown in Figure 4.

Table III shows that among the trajectories ran for 10 Å cutoff only the CN trajectory exhibits a decreasing trend in total energy E_tot when compared over the first 0.5 ns and last 2 ns. The remaining trajectories lead to structures with either higher E_tot (CX, PX) or approximately unchanged E_tot (PN). This is a consequence of the fact that E_tot is composed of E_vac and E_solv contributions that have a compensating effect, though their scale is different. It is generally valid for all four trajectories that E_vac decreases during MD (see Fig. 4). However, the solvation energy E_solv is able to revert this trend. As the contribution of E_solv is larger than that of E_vac, it may result in the increase of total energy E_tot. This phenomenon is best seen from the comparison of the CN and PN trajectory energies. While E_vac is always above zero in CN, it is well balanced here by E_solv to give low E_tot. In the case of PN, E_vac is low, but the structure is disfavored by E_solv, giving high E_tot. Both X-ray started trajectories exhibit an increasing trend in E_tot with decreasing E_vac and the opposite trend in E_solv. The reason why E_vac decrease is overridden by E_solv increase for three out of four trajectories is discussed in detail in the next section. The other dependency can be observed for different electrostatics treatment. While both trajectories with PME treatment of electrostatics give low E_vac energy, both cutoff trajectories show above-zero values in E_vac. This is balanced by E_solv, which is lower for cutoff trajectories than for PME ones. When the cutoff value was extended to 15 Å, E_tot decreased in both trajectories. While E_tot decreased only slightly in C15X trajectory, there was a rapid decrease by almost 100 kcal/mol in C15N run. This latter simulation gave also the lowest average E_tot value in the whole study. In order to see what happens if we include solute entropy, we calculated its contribution via harmonic vibration analysis (as described in the Method section) -- see Table III A. The average values from the last 2 ns of all the trajectories are almost identical, which means that the entropy does not change the relative order of the free energy values.

A


B

Figure 4: Free energy development as calculated by MM-PBSA approach for all trajectories. The total energy is marked by a thick solid line and it is composed of E_vac (upper, dashed line) and E_solv (lower, dotted line).

Interactions of Charged Ends

In order to understand the above mentioned trends we focused our attention on analysis of the E_vac energy, especially on its substantial decrease within the last 3 ns of the CX and PX trajectories (see Fig. 4). The explanation comes from E_vac decomposition into single terms in relation to the structural changes. Both ends of the peptide, the N- and C-termini, bear charges at their amino acid sidechains. This is caused by partial charges used by the force field and summarized over all of the residues. Cys1 (N-terminus) possesses a charge of +1 and Trp21 (C-terminus) has a charge of -1. Another negative charge is located at Asp18. As the C-terminal of the peptide is flexible it can move in such a way that it can reach close proximity to the N-end. This results in strong electrostatic interaction between opposite charges at the peptide ends, and in H-bond establishment between Cys1 and Trp21 or between Cys1 and Asp18. This kind of stabilization appears in three of the first four MD trajectories we ran. Hydrogen atoms at the N atom of Cys1 make H-bond contacts with O or OXT atoms of Trp21 (CX trajectory), with OD1 or O atoms of Asp18 (PX trajectory), and with OD1 or OD2 atoms of Asp18 (PN trajectory). No H-bond of this kind was observed in the CN trajectory. The above H-bond contacts were established either in the later part of the trajectory (CX, PX trajectories) or even in a very early stage, probably during the equilibration period (PN trajectory). We have observed that the creation of the H-bond contact strongly correlates with the change in E_vac energy, especially with its electrostatic non-bonded component (see Fig. 5).

A


B


C


D


E


F


Figure 5: Time evolution of the electrostatics component of E_vac in correlation with the distance between the C- and N-ends. Point representation shows the layout of the snapshots with respect to Energy and Distance common axes. The CX trajectory with distance Cys1@N-Trp21@C (A), PX trajectory with distance Cys1@N-Asp18@CA (B), CN trajectory with distance Cys1@N-Trp21@C (C), PN trajectory with distance Cys1@N-Asp18@CG (D), C15X trajectory with distance Cys1@H1-Asp18@OD2 (E), C15N trajectory with the distance Cys1@H1-Trp21@OXT (F).

In the first case where no interaction of the charged termini appeared (CN trajectory), E_vac conserves its above-zero value (compare with Fig. 4). However, as the E_solv balances E_vac by the opposite trend, but with higher deviations in the opposite direction, it results in the lowest total E_tot energy for the CN trajectory (as seen from Etot values in Table III). On the other hand, the contact between the charged ends lasts the entire PN trajectory causing very low E_vac but higher E_solv. All those tendencies together give the highest E_tot average value. Correlation graphs in Figure 5 show how spread out or concentrated the values are for the distance between the charged ends and for the energy (influenced mainly by the non-bonded electrostatic term) calculated for single snapshots along the trajectories. This is especially visible in the correlation part of Figure 5 (right bottom corners) where the corresponding points create clusters that are represented by charge stabilized structures. It is remarkable that the clusters are much denser for the PME than for the cutoff method. It implies that in both PME trajectories the movement of the ends is more restricted to a short distance, which results in lower electrostatic energy.

An important observation is that even if the charges at the terminal residues 1 and 21 were behind the cutoff distance in the experimental structures (i.e., at the beginning of the simulations, more than 13 Å in both cases), the distance between the charged residues 1 and 18 was shorter than the cutoff (around 7 Å) in both cases. It means that the system had in all cases a chance to place the charged ends closer to each other and keep them so for the entire trajectory. However, it only happened for the PME trajectories even if the free energy was higher than in the cutoff trajectories. A possible explanation is that the structure with the charged ends close to each other creates a minimum of free energy that is deep enough to be retained inside it within the time scale for which the MD was run. On the other hand, the cutoff trajectories were more driven by the solvation of the charged ends, which lead to free energy decrease for the CN trajectory. The above observation confirms the fact that PME restricts the structure more than the cutoff method.

In the trajectories with extended cutoff (C15X and C15N), stable interactions of charged ends also did not appear. The system in C15N trajectory had a chance to establish these H-bonds in the early part of the simulation (300-600ps, see Fig. 5F), but the chance did not lead to closing the structure. Such a close distance between charged ends did not appear in C15X simulation.

Secondary Structure Elements

We have analyzed the development of secondary structure elements during the trajectories, in an attempt to find a relationship with free energy values.

The main elements of secondary structure were identified by Kabsch and Sander?s method of ?turns?. H-bonded patterns found in the endothelin structure are formed in several secondary structures: bend, turn, β-bridge, and three types of helices -- 3₁₀-helix (3-turn), α-helix (4-turn), and 5-helix (5-turn). No β-sheet was found in any of the structures within all the trajectories. The overall picture is shown in Figure 6.

Both NMR trajectories, CN and PN, started from the structure with a well formed regular α-helix at the residues 9-15. Although the averaged structures do not exhibit this element (see Fig. 6 B and D), trajectory development shows that the α-helix sometimes appears between the residues 10-14. This is more apparent from the CN trajectory, where the structure keeps close to the regular α-helix between the residues 10-13. During the PN trajectory, the α-helix is broken in the early stage of the simulation, but the structure remains close to it during the entire simulation. The α-helix appears in the X-ray trajectories as well, although it was not present in the experimental structure. The most stable regular α-helix, 11-14, appears in the PX trajectory for about 3.5ns without any interuption. After that period the helix shortens to 11-13 and reverts to the 3₁₀-form. The regular α-helix appears from time to time also during the CX trajectory.

The more tightened 3₁₀-helix (which requires only three consecutive residues bound with H-bonds), is present in the experimental X-ray structure, and it appears in some modifications in all the trajectories. This initial arrangement is almost retained in both CX and PX trajectories, even if the 3₁₀-helix shifts by one residue (PX trajectory) or it becomes part of the larger 5-helix (CX trajectory). The C-end added to the NMR structure folds also into the 3₁₀-helix in CN and PN simulations. Three 3₁₀-helices are found in the PX average structure. They may be considered as parts of one helix at the 7-17 residues with two single-residue interruptions at residues 10 and 14. No secondary structure element is seen between the residues 1-6, probably because they are fixed by two disulphidic bridges.

Based on the above, it is difficult to decide whether cutoff or PME conserves secondary structures better. Even if the secondary structures have different flexibility and locations for each of the two methods, it is confirmed that both PME and cutoff methods are able to describe them in a proper manner. After inspection of the relationships between the history of secondary structure elements and free energy development during all the trajectories we have found that there is no apparent correlation. So we may state that the formation/breaking of secondary structure elements is not a leading force in free energy development and it does not depend on the electrostatic interaction energy calculation method used.

Conclusions

Four 10 ns and two 2 ns molecular dynamics trajectories were run on endothelin-1 in order to investigate the influence of two different treatments of long-range electrostatic interactions, cutoff and PME, on geometrical changes and free energy development during the trajectories. The simulations started from two different experimental structures that were obtained by X-ray crystallography and NMR spectroscopy. We have found that the charge-charge interaction between opposite charged peptide ends is the leading force governing development of the trajectories. The trajectory course is also strongly influenced by the starting structures. The peptide ends, C- and N- termini, were rather distant in both starting structures. Both trajectories that started from the X-ray structure needed approximately the same time to bring the charged ends into close proximity and stabilize them by H-bonds. In the case of the NMR starting structure, the PME trajectory lead to this kind of stabilization immediately. In contrast, the cutoff method never created stable H-bonds that would connect the opposite ends of the peptide and that would last a longer time, even though the peptide ends met at the appropriate distance several times during the simulation. This was confirmed also by simulations with larger cutoff. This is in agreement with the new NMR structure (33) that appeared when this paper was in revision.

The NMR experiments suggested (13) that the C-tail of the peptide is rather flexible. It undergoes an interconversion between at least two different conformers with lifetimes in the microsecond domain. The new NMR structure shows that the C-terminus is bent apart of N-terminus (open structure). Three of the six MD trajectories we ran converged to the structure with stable interaction between charged peptide ends (closed structure). Only the simulations with cutoff treatment of electrostatics lead to open structures. We may therefore conclude, seen from this point if view, that the most reliable simulations were those with cutoff electrostatics treatment. This is confirmed by the fact that also the X-ray structure does not exhibit such an interaction of terminals that would lead to a closed structure.

A frequently asked question concerning MD simulations is over the stability of trajectories. This is usually expressed by RMS deviation. It is known that the PME method produces more stable trajectories than the cutoff method in the sense of amplitude of motion (1, 2, 34). However, this does not directly imply that such trajectories are more reliable (4). It is also well known that using PME is essential for charged systems like nucleic acids (35). Since PME can lead to protein ?overstabilization? during MD simulations, it remains an open question whether it must be used for peptides and/or proteins. This ?overstabilization? could cause an artificial locking of geometry within a local minimum which would remain there for an unnaturally long time.

The data presented here does not suggest an unambiguous conclusion about which of the two methods of long-range electrostatics treatment should be used in MD simulations. However, we can certainly state that using cutoff may not always be catastrophic (as mentioned in ref. 2, 36). On the contrary, it may in some cases provide better results for peptides and proteins than PME does. This leads to a recommendation that the choice of electrostatics treatment should be made carefully and not automatically by choosing the PME method due to it being the most widely used.

Acknowledgement

The study was supported by the Ministry of Education, Youth and Sports of the Czech Republic, grant no. LN00A016. The financial support is gratefully acknowledged. We would like to thank the Brno Supercomputer Center for providing us with access to computer facilities, and Roger Turland for language corrections.

References and Footnotes

1. T. Fox and P. A. Kollman. Proteins 25, 315-334 (1996).
2. O. N. de Souza and R. L. Ornstein. J. Biomol. Struct. Dyn. 16, 1205-+ (1999).
3. M. Krol. J. Mol. Model. 9, 316-324 (2003).
4. R. Walser, P. H. Hunenberger, and W. F. van Gunsteren. Proteins 43, 509-519 (2001).
5. R. Svobodova-Varekova, J. Koca, and C.-G. Zhan. Int. J. Mol. Sci. 5, 154-173 (2004).
6. C. Sagui and T. A. Darden. Annu. Rev. Bioph. Biom. 28, 155-179 (1999).
7. P. J. Steinbach and B. R. Brooks. J. Comput. Chem. 15, 667-683 (1994).
8. R. J. Loncharich and B. R. Brooks. Proteins 6, 32-45 (1989).
9. T. Darden, D. York, and L. Pedersen. J. Chem. Phys. 98, 10089-10092 (1993).
10. U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. G. Pedersen. J. Chem. Phys. 103, 8577-8593 (1995).
11. P. Ewald. Ann. Phys. 64, 253-287 (1921).
12. P. H. Hunenberger and J. A. McCammon. J. Chem. Phys. 110, 1856-1872 (1999).
13. N. H. Andersen, C. P. Chen, T. M. Marschner, S. R. Krystek, and D. A. Bassolino. Biochemistry 31, 1280-1295 (1992).
14. R. W. Janes, D. H. Peapus, and B. A. Wallace. Nat. Struct. Biol. 1, 311-319 (1994).
15. H. M. Berman, J. Westbrook, Z. Feng, G. Gilliland, T. N. Bhat, H. Weissig, I. N. Shindyalov, and P. E. Bourne. Nucleic Acids Res. 28, 235-242 (2000).
16. W. Kabsch and C. Sander. Biopolymers 22, 2577-2637 (1983).
17. D. A. Case, D. A. Pearlman, J. W. Caldwell, T. E. Cheatham III, W. S. Ross, C. L. Simmerling, T. A. Darden, K. M. Merz, R. V. Stanton, A. L. Cheng, J. J. Vincent, M. Crowley, V. Tsui, R. J. Radmer, Y. Duan, J. Pitera, I. Massova, G. L. Seibel, U. C. Singh, P. K. Weiner, and P. A. Kollman. AMBER Version 6, University of California, San Francisco (1999).
18. W. D. Cornell, P. Cieplak, C. I. Bayly, I. R. Gould, K. M. Merz, D. M. Ferguson, D. C. Spellmeyer, T. Fox, J. W. Caldwell, and P. A. Kollman. J. Am. Chem. Soc. 117, 5179-5197 (1995).
19. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein. J. Chem. Phys. 79, 926-935 (1983).
20. J. P. Ryckaert, G. Ciccoti, and H. J. C. Berendsen. J. Comput. Phys. 23, 327-341 (1977).
21. H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. Dinola, and J. R. Haak. J. Chem. Phys. 81, 3684-3690 (1984).
22. A. Onufriev. AMBER Version 8, University of California, San Francisco (2004).
23. E. Fadrna, N. Spackova, R. Stefl, J. Koca, T. E. Cheatham III, and J. Sponer. Biophys. J. 87, 227-242 (2004).
24. R. Stefl, T. E. Cheatham, N. Spackova, E. Fadrna, I. Berger, J. Koca, and J. Sponer. Biophys. J. 85, 1787-1804 (2003).
25. D. Sitkoff, K. A. Sharp, and B. Honig. J. Phys. Chem. 98, 1978-1988 (1994).
26. N. Spackova, T. E. Cheatham, F. Ryjacek, F. Lankas, L. van Meervelt, P. Hobza, and J. Sponer. J. Am. Chem. Soc. 125, 1759-1769 (2003).
27. K. A. Sharp and B. Honig. Annu. Rev. Bioph. Bio. 19, 301-332 (1990).
28. M. F. Sanner, A. J. Olson, and J. C. Spehner. Biopolymers 38, 305-320 (1996).
29. E. Fadrna and J. Koca. J. Biomol. Struct. Dyn. 20, 715-731 (2003).
30. B. Honig and A. Nicholls. Science 268, 1144-1149 (1995).
31. J. Srinivasan, T. E. Cheatham, P. Cieplak, P. A. Kollman, and D. A. Case. J. Am. Chem. Soc. 120, 9401-9409 (1998).
32. B. Jayaram, D. Sprous, and D. L. Beveridge. J. Phys. Chem. B 102, 9571-9576 (1998).
33. H. Takashima, N. Mimura, T. Ohkubo, T. Yoshida, H. Tamaoki, and Y. Kobayashi. J. Am. Chem. Soc. 126, 4504-4505 (2004).
34. L. Monticelli and G. Colombo. Theor. Chem. Acc. 112, 145-157 (2004).
35. T. Darden, L. Perera, L. P. Li, and L. Pedersen. Struct. Fold. Des. 7, R55-R60 (1999).
36. J. Norberg and L. Nilsson. Biophys. J. 79, 1537-1553 (2000).
37. H. Resat and M. Mezei. Biophys. J. 71, 1179-1190 (1996).