J Mol Model (2014) 20:2212 DOI 10.1007/s00894-014-2212-x
ORIGINAL PAPER
Molecular dynamics investigation of Helicobacter pylori chemotactic protein CheY1 and two mutants Ahmet Yildirim & Mustafa Tekpinar & Tsjerk A. Wassenaar
Received: 8 October 2013 / Accepted: 16 March 2014 # Springer-Verlag Berlin Heidelberg 2014
Abstract CheY is a chemotactic response regulator protein modulating the rotation direction of bacterial flagellar motors. It plays an important role in the colonization and infection of Helicobacter pylori (H. pylori), which is a common pathogen. Recently, the structure of CheY1 of H. pylori (HpCheY1) was solved, showing similarities and differences with CheY from E. coli. Here, we report 200 ns atomistic molecular dynamics (MD) simulations of HpCheY1 and two mutants. The results suggest that the surface of HpCheY1 has regions with increased affinity for Mg2+. In addition, wildtype HpCheY1 (WT HpCheY1) shows characteristic dynamics in helix 4, which is involved in FliM binding. This dynamics is altered in the D53A mutant and completely suppressed in the T84A mutant. The results are discussed in relation to the binding and function of HpCheY1. Keywords Chemotaxis protein . D53A mutant . H. pylori . Molecular dynamics simulation . Principal component analysis . T84A mutant
Introduction In 1982 Marshall and Warren identified H. pylori as a bacterium that grows in human gastric epithelial tissue and mucus A. Yildirim Department of Physics, Faculty of Science and Art, Siirt University, Siirt 56100, Turkey M. Tekpinar Department of Physics, Faculty of Science and Art, Yüzüncü Yıl University, Van 65100, Turkey T. A. Wassenaar (*) Department of Biology, University of Erlangen-Nürnberg, Staudtstr. 5, 91058 Erlangen, Germany e-mail: [email protected]
[1]. Its presence was found to be linked to gastritis, many forms of ulcers, as well as to gastric cancer [2–4]. Currently, it is known that H. pylori is present in approximately 50 % of the world population and that additional factors are involved in the development of H. pylori related diseases [5]. The bacterium is able to colonize and infect tissues through directed movement in response to a gradient of a chemical stimulus. This kind of motion, called chemotaxis, is affected through flagellae, which rotate by means of a flagellar motor [6]. A key role in the control of this motor in H. Pylori, as well as in other bacteria, is the chemotactic effector protein CheY. Upon activation through phosphorylation, this protein binds to the flagellar motor domain [7], causing the direction of rotation to invert, which changes the motion from 'running' to 'tumbling' [8]. In the case of H. pylori, the chemotaxis causes the bacterium to avoid the gastric acid environment and colonize the mucus layer, where the pH is higher [9]. Unlike E. coli, H. pylori has two CheY response regulators, here referred to as HpCheY1 and HpCheY2. HpCheY1 seems the main partner of FliM, while HpCheY2 acts as a phosphate sink to modulate the half-life of phosphorylated HpCheY1 [7, 10]. Recently, the crystal structure of activated HpCheY1 was solved. Beryllium fluoride (BeF3) was used as substitute for phosphate, because the lifetime of phosphate activated CheY1 is too short to allow structure elucidation. In addition to the structure of the wild-type protein (WT), crystal structures were solved for two mutations, namely D53A and T84A. The former of these cannot be activated by phosphorylation, while the latter displays increased activity [11]. Unfortunately, solving an apo-structure of HpCheY1 failed due to the sulfate present in the medium, which apparently also binds in the active site, resulting in an active conformation. From the available structures of E. coli CheY, which includes a FliM-bound structure [12], it is expected that the differences between the activated and the apo forms are small.
2212, Page 2 of 8
For such cases, MD simulations can be used to complement experimental methods and assess the conformation and dynamics of the unbound protein to contribute to a comprehensive characterization. In this work, we use atomistic scale MD simulations to characterize the conformation and dynamics of WT HpCheY1 and the two mutants D53A and T84A. D53 is the active residue for phosphorylation, and T84 is located near the binding site and hydrogen bonds to the phosphate. The results highlight regions on the surface of the protein with increased affinity for Mg2+ and show differences in the conformations as well as in the dynamics of wildtype and mutants. To compare the conformations and dynamics, an alternative approach of principal component analysis (PCA) was taken, using a newly developed, interactive method in PyMOL (manuscript in preparation). The results are discussed in relation to the mechanism of action of HpCheY1, with reflection on the role of dynamics.
Methods Simulations All simulations were performed with GROMACS version 4.5.4 [13, 14], in conjunction with the GROMOS96 53a6 united atom force field [15]. The integration time step used was 2 fs. Hydrogen bonds were treated as constraints using the LINCS algorithm [16], except for solvent molecules, which were constrained using the SETTLE algorithm [17]. Water was modeled using the simple point charge model [18]. Van der Waals interactions were calculated using a 1.4 nm cut-off. Electrostatic interactions were calculated using the particlemesh Ewald (PME) method [19, 20], with a grid spacing of 0.12 nm and a real space cut-off of 0.9 nm. All calculations were performed on a workstation with two Intel Xeon E52687 W processors and 64GB RAM. Starting structures were taken from Protein Data Bank with accession codes 3GWG (WT), 3H1F (D53A), and 3H1G (T84A) [11]. Bound water and Mg2+ ions present in the crystal structures were retained during the setup of the simulations. Structures were first energy minimized using a steepest descent algorithm, after which they were solvated in a periodic rhombic dodecahedron unit cell with a minimal distance between periodic images of 2 nm. The solvent consisted of a 0.15 M NaCl solution, corrected to compensate for the net charge of the protein. The resulting system for WT HpCheY1 contained 18,187 atoms in total, including 5645 water molecules, 17 sodium and 19 chloride ions. For the D53A mutant the solvated system measured 17,969 atoms in total, including 5574 water molecules, 17 sodium, and 18 chloride ions. The T84A mutant yielded a system consisting of 17,288 atoms in
J Mol Model (2014) 20:2212
total, including 5347 water molecules, 16 sodium and 18 chloride ions. The solvated systems were subjected to another round of steepest descent energy minimization, followed by 100 ps MD simulation in which position restraints were applied to the heavy atoms in the system. During this run, the system was weakly coupled to an external heat bath with a temperature of 300 K and a coupling time of 0.1 ps using a V-rescale thermostat [21]. This stage was followed by 100 ps position restrained MD in which the pressure was kept constant at 1 bar by coupling to a Berendsen barostat [22] with 1.0 ps and a compressibility of 4.5×10−5 bar. After this initial equilibration, productions runs were performed for a duration of 200 ns, with frames written to disk every 40 ps. Analysis Analysis was performed using a combination of the tools available from GROMACS, PyMOL [23] and auxiliary tools written for specific purposes. Principal component analysis (PCA) was performed using a new PyMOL module (manuscript in preparation). Post-hoc analysis of PCA results was performed in R [24]. The root mean square deviation To follow the relaxation over time, the root-mean-squared deviation (RMSD) was determined for each system. The RMSD was determined against the starting structure, as well as against the average structure from the last 100 ns and 20 ns of the simulations. The reason for this approach is that the RMSD with respect to the starting structure does not usually allow assessment of relaxation, as any plateau in RMSD corresponds to a volume of conformational space, with a size proportional to the RMSD. In other words, there may be drift in conformation, which is not reflected by a change in RMSD. Calculating the RMSD against the average structure does allow assessment of relaxation, because the average structure is an estimator of the mean structure, which is an ensemble property, and because the direct neighborhood of the average structure is associated with lower RMSD values and consequently smaller associated volume elements, resulting in a more direct relation between RMSD and conformational state. The RMSD of a given structure with respect to a reference structure rref, is calculated as RMSDðt Þ ¼
2 1=2 1 XN ref m r ð t Þ−r ; i i i¼1 i M
ð1Þ
where M=∑mi and ri(t) is the position of atom i at time t after least square fitting the structure to the reference structure.
J Mol Model (2014) 20:2212
Page 3 of 8, 2212
Distributions of Mg2+ For visualizing Mg2+ distributions around HpCheY1, first the molecular shaped unit cell representation was determined, which means that every particle is positioned in the periodic system where it is closest to a given reference group, in this case the protein [25]. Subsequently, all frames were fit to the WT starting structure, and the positions of the ions were projected around the reference structure.
Principal component analysis PCA in MD simulations starts with the construction of the positional covariance matrix. If x denotes a single structure consisting of n atoms represented as a single column vector of length 3n, then combining the structures sampled over time gives a data matrix X of dimensions 3n × m, with m being the number of frames. From this frame the covariance matrix can be calculated as S¼
1 −1 X I−ð j0 jÞ jj0 X0 : m−1
ð2Þ
Here I is the m × m identity matrix and j is a vector of ones with length m. The matrix (I−(j′j)−1jj′) is the idempotent centering matrix, and an apostrophe is used to indicate transposition. The covariance matrix S is decomposed into a matrix of eigenvectors P and a diagonal matrix of eigenvalues D to obtain the principal components: S ¼ PDP0 :
ð3Þ
Each principal component corresponds to a linear correlation between atomic displacements and can be regarded as a collective motion, although the real underlying motion may be non-linear and will then be dispersed over a number of components. We introduce here an alternative, hierarchical approach to use PCA more effectively to study conformation and dynamics in relation to mutations. The first stage consists of determining the average structures. These are used to determine the plane through conformational space in which the distinct states are most dispersed. The relaxation of the systems toward the respective end states then follows from the projections of the trajectories onto this plane. The second stage consists of PCA applied to each of the systems separately. This approach allows constructing a comprehensive framework to assess similarity and differences in structure and dynamics of different mutants.
PCA was applied to the last 100 ns of the simulations, as well as to the last 20 ns, allowing assessment of the robustness of the results in relation to the length of the trajectory.
Results and discussion The root mean square deviation To assess the stability of each structure, the RMSD was calculated for backbone atoms and for heavy atoms against the starting structure and against the mean structures obtained from the last 100 ns and the last 20 ns of the simulations as a function of time. These graphs are shown in Fig. 1. The RMSD against the starting structure, for both backbone (Fig. 1a) and all non-hydrogen atoms (Fig. 1b), increases over a time span of 15 ns to values of approximately 0.25 nm (backbone) or 0.35 nm (heavy atoms). The RMSD against the average obtained from the last 100 ns shows a steady decrease over 75 ns with little difference between backbone (Fig. 1c) and all heavy atoms (Fig. 1d). The T84A mutant reaches a stable state with little fluctuation, while the D53A mutant displays more flexibility. The wild type is the most dynamic, displaying a higher RMSD and showing larger fluctuations. A similar view is obtained from the RMSD against the average from the last 20 ns, which shows little difference with the other graphs for D53A and T84A, but displays significant fluctuations for wild type, temporarily reaching RMSD values against the average close to the value observed for the starting structure. Binding and distribution of magnesium The crystal structures of WT HpCheY1 and T84A each contain two magnesium ions. The D53A mutant contains only a single magnesium ion, as the binding site at D53 is removed. In Fig. 2, the positions of the magnesium ions over the whole length of the simulation are shown for all three systems to investigate the stability and distribution of Magnesium in and around HpCheY1. In the case of WT and T84A one ion remained bound for the whole length of the simulation, shown in red in Fig. 2. The other ion, as well as the single ion in the system containing the D53A mutant, diffused over the surface of the protein. This ion is shown in blue in Fig. 2. In the WT, there appear to be four sites on the protein surface with increased density, which are labeled B-E in Fig. 2. Removal of the D53 side chain causes a shift in density from site C to B, but also shifts the density from site E to an alternative site F near helix 5. In addition, a small cluster of density is observed above site B, which is labeled A, and which is located close to the binding site in the protein. Apparently, the loss of the internal magnesium ion causes a shift in the distributions in
2212, Page 4 of 8
J Mol Model (2014) 20:2212
Fig. 1 Root mean-square deviation (RMSD) over time. RMSD against starting structure for backbone atoms (a) and all heavy atoms (b). RMSD against average structures from the last 100 ns of each simulation for backbone atoms (c) and all heavy atoms (d). RMSD against average
structures from the last 20 ns of each simulation for backbone atoms (e) and all heavy atoms (f). Colors are black for wild type, magenta for D53A, and green for T84A
the direction of the binding site, suggesting a mechanism for catching and incorporating magnesium.
The three mean structures span a plane in conformational space. Two orthonormal vectors in this plane are easily obtained by performing PCA on those mean structures. Subsequently, the trajectories can be projected onto these vectors, which then allows assessment of the similarity or otherwise of the relaxation pathways and the ensembles. These projections are shown for the planes through the means from the last 100 ns and the last 20 ns in Fig. 4, where the starting structure and mean structures are indicated with large dots. From this Fig., it can be seen that the three structures relax in different directions and sample distinct areas in conformational space. Interestingly, the two mutants relax in opposite directions from the starting structure, while the WT relaxes in a direction almost orthogonal to the axis D53A-T84A. This is true for both views. Previously, it has been suggested that investigation of differences in structure and dynamics could be achieved by performing PCA on a combined trajectory [26]. However, that is only true if the dynamics orthogonal to the relaxation has much smaller variance than the dynamics in the plane of the
Conformations and dynamics To investigate the effects of the mutations on the conformation of HpCheY1, the mean structures from the last 100 ns and the last 20 ns of the simulations were determined, the latter of which are shown in Fig. 3. These structures show that there are small, but significant changes in the conformations, in particular in helix 4. In the T84A mutant, shown in green, this helix becomes more ordered compared to WT, while in the D53A mutant it becomes more disordered. In addition, the latter shows a shift in helices 1 and 5. Interestingly, the conformation of these helices in the D53A mutant is more similar to the crystal structures, which means that the helices are shifted in WT and T84A during the relaxation. The increased ordering of T84A is caused by E90, which is internalized to compensate for the loss of the threonine hydroxyl group. As a consequence, the loop K88-A89-E90-A91 becomes part of helix 4.
J Mol Model (2014) 20:2212
Fig. 2 Distributions of magnesium over time in and around HpCheY1 wildtype and mutants. WT HpCheY1 (a), D53A mutant (b), and T84A mutant (c). A total of six sites on the surface with increased Mg2+ density are observed, which are labeled A–F. The internal Mg2+ ion in WT (a) and T84A (c) is colored red and does not leave the binding site. In the WT, the main site on the surface is D. Removal of the internal binding site by mutation of D53 to alanine causes the density to shift toward B and A, as well as to a site F on the other side. Mutation of T84 to alanine retains the internal binding site, but also shifts the surface binding to site B
Fig. 3 Stereo image of HpCheY1 mean structures after equilibration. The image shows the mean structures from the last 20 ns of a 200 ns simulation of WT HpCheY1 (gray), D53A mutant (purple), and T84A mutant (green). The helices are numbered from N-terminus to C-
Page 5 of 8, 2212
relaxation, and the two are uncorrelated. In most cases, one of these conditions is not satisfied and the modes obtained mix relaxation and steady-state dynamics, making it difficult to comment on either. The procedure used here solves this issue and provides a comprehensive framework for investigating differences between states in terms of both structures and dynamics. From Fig. 4 it is clear that each system ends up in a distinct state. The overall picture does not change much when the last 100 ns is used for the analysis (Fig. 4a) or only the last 20 ns (Fig. 1b). In particular, the profiles for D53A and T84A hardly change, in agreement with the view obtained from the RMSD. For WT there is a difference in the projections of the trajectory around the mean. When using the mean structure from the last 100 ns, the projections are more dispersed and several clusters can be identified, which is not the case when using the average structure from the last 20 ns. This suggests that the last 20 ns alone are not sufficiently representative for the WT equilibrium ensemble, which is in accordance with the view obtained from the RMSD curves. For a first characterization of the equilibrium dynamics, ellipses are drawn around each mean structure, summarizing the in-plane dynamics of the last 100 (yellow) or 20 (orange) ns of each system. Again, D53A and T84A show significant overlap of the ellipses, suggesting that the collective motions in both regimes are similar, whereas the ellipses for WT show less overlap, indicating different mean projections and different main directions of motion. These ellipses also show that WT and D53A have differently oriented main radii, reflecting differences in the collective motions. In the T84A mutant, there seems to be little overall motion after relaxation, suggesting that the removal of the threonine side chains rigidifies the structure. Because the collective motions appear different in nature for each of the mutants, we decided to apply PCA to each system separately, of which the results are shown in Fig. 5. PCA was applied to the last 100 and last 20 ns from each
terminus. Residue E90 is highlighted with spheres. The image shows an extended helix 4 in T84A due to inward positioning of E90. In addition, helices 1 and 5 are shifted in D53A
2212, Page 6 of 8
J Mol Model (2014) 20:2212
Fig. 4 Projections of trajectories on the plane through mean structures. Projections on the plane determined from the mean structures from the last 100 ns (a) and on the plane determined from the last 20 ns (b). In both panels the starting structure is marked with s. The three 200 ns simulations are shown in gray/ black (wildtype), violet/purple (D53A mutant), and green/forest (T84A mutant). The darkest, thickest line marks the last 20 ns of each simulation, whereas the intermediate line marks the trajectory from 100 to 180 ns. The respective mean structures are indicated with a yellow dot, while orange ellipses are drawn to mark the principal directions of the motions in the plane in the last 100 ns, and yellow ellipses mark the principal directions in the plane in the last 20 ns
trajectory, corresponding to the thicker, darker parts of the trajectories drawn in Fig. 4. The traces of the covariance matrices from the last 100 ns and the last 20 ns are given in Table 1. The results from the last Fig. 5 Collective dynamics in WT HpCheY1 and mutants D53A and T84A. Visualization of the first principal component from simulations of HpCheY1. The loadings are drawn as rainbowcolored bicones, centered on the respective mean structures. The residues corresponding to D53 (magenta) and T84 (green) are highlighted and the side chain of E90 is shown as spheres. The images show the first principal component determined from the last 100 ns simulation for WT (a), D53A (b), and T84A (c), as well as the first component determined from the last 20 ns of simulation for WT (d), D53A (e), and T84A (f)
100 ns show that the wild type is the most dynamic, although the effect is more severe for T84A. If the PCA is performed on the last 20 ns only, the total variance for wild type is found to be lower than for D53A, but the values for D53A and T84A
J Mol Model (2014) 20:2212 Table 1 Total variance and first three eigenvalues
Page 7 of 8, 2212
Time
Structures
Total variance 2
Last 100 ns
Last 20 ns
WT D53A T84A WT D53A T84A
are unchanged. This actually shows that the part of the simulation, even when it is commonly considered to be in an equilibrium state, still has a large influence on the results from analysis of conformations and dynamics. The mean structure (Fig. 3) suggests that the decreased motility of T84A is caused by the extension and stabilization of helix 4, as discussed above. The first principal component obtained from T84A also only captures 8 % of the total variance, which suggests that the dynamics is mostly governed by thermal fluctuations, rather than by collective dynamics. For both WT and D53A, approximately 25–30 % of the variance is captured by the first component, suggesting a certain degree of collectivity. The collectivity can be assessed from Fig. 5. The top panel shows that in WT there is collectivity in the motions of helix 4 and the preceding loop. In particular, there seems to be a certain degree of rotational motion of the helix, related to the dynamics of the loop. This is observed in the analysis of the last 100 ns as well as of the last 20 ns, showing that, despite the differences, there is qualitative agreement between the regimes. In the D53A mutant, the rotational motions are largely absent. Interestingly, the region that shows collective dynamics in wild type, helix 4, forms part of the FliM interacting region. This suggests that the altered dynamics profiles of the D53A and T84A mutants could explain part of the differences observed in vivo, aside from effects on the phosphate binding. This applies in particular to T84A, which displays higher activity. In turn, this raises the suggestion that modulation of the dynamics by binding of phosphate plays a role in the binding of HpCheY1 and FliM, or even more broadly in binding of CheY to FliM in bacteria. This possibility will be further explored in a future study.
Conclusions We performed MD simulations to characterize the conformation and dynamics of apo WT HpCheY1, D53A mutant, and T84A mutant. The results suggest that the apo structures relax to distinct states, mainly characterized by changes in the FliM binding region helix 4. The T84A mutant is characterized by lengthening of the helix and stabilization due to internalization
Eigenvalue 1 2
Eigenvalue 2 2
Eigenvalue 3
nm
nm (%)
nm (%)
nm2 (%)
29.85 19.06 13.73 17.23 19.76 13.18
10.33 (35 %) 2.98 (16 %) 0.83 (6 %) 4.28 (25 %) 4.58 (23 %) 1.11 (8 %)
3.34 (11 %) 1.39 (7 %) 0.74 (5 %) 0.97 (6 %) 2.28 (12 %) 0.97 (7 %)
1.32 (4 1.23 (6 0.64 (5 0.90 (5 1.04 (5 0.70 (5
%) %) %) %) %) %)
of E90 to compensate the loss of the T84 hydroxyl group in the active site. The different states are also characterized by differences in the dynamics, with the wild type showing a certain degree of collectivity of backbone and side chain motions in helix 4 and the preceding loop. These motions are suppressed in the D53A mutant and diminished in the T84A mutant. Future studies should show whether such changes in dynamics are also invoked by phosphate activation and what role this dynamics plays in FliM binding and CheY mediated flagellar motor control.
References 1. Marshall BJ, Warren JR (1984) Unidentified curved bacilli in the stomach of patients with gastritis and peptic ulceration. Lancet 1: 1311–1315 2. Labigne A, de Reuse H (1996) Determinants of Helicobacter pylori pathogenicity. Infect Agents Dis 5:191–202 3. McColl KEL (1997) Helicobacter pylori: clinical aspects. J Infect 34: 7–13 4. Riegg SJ, Dunn BE, Blaser MJ (1995) Microbiology and pathogenesis of Helicobacter pylori. In: Blaser MJ et al (eds) Infections of the gastrointestinal tract. Raven, New York, pp 535–550 5. Covacci A, Telford JL, Del Giudice G, Parsonnet J, Rappuoli R (1999) Helicobacter pylori virulence and genetic geography. Science 284:1328–1333 6. Scharf BE, Fahrner KA, Turner L, Berg HC (1998) Control of direction of flagellar rotation in bacterial chemotaxis. Proc Natl Acad Sci 95:201–206 7. Foynes S, Dorrell N, Ward SJ, Stabler RA, McColm AA, Rycroft AN, Wren BW (2000) Helicobacter pylori possesses two CheY response regulators and a histidine kinase sensor, CheA, which are essential for chemotaxis and colonization of the gastric mucosa. Infect Immun 68:2016–2023 8. Berg HC, Brown DA (1972) Chemotaxis in E. coli analysed by threedimensional tracking. Nature 239:500–504 9. Schreiber S, Konradt M, Groll C, Scheid P, Hanauer G, Werling HO, Josenhans C, Suerbaum S (2004) The spatial orientation of Helicobacter pylori in the gastric mucus. Proc Natl Acad Sci 101: 5024–5029 10. Jimenez-Pearson MA, Delany I, Scarlato V, Beier D (2005) Phosphate flow in the chemotactic response system of Helicobacter pylori. Microbiology 151:3299–3311 11. Lam KH, Ling TK, Au SW (2010) Crystal structure of activated CheY1 from Helicobacter pylori. J Bacteriol 192:2324–2334
2212, Page 8 of 8 12. Lee SY, Cho HS, Pelton JG, Yan D, Henderson RK, King DS, Huang LS, Kustu S, Berry EA, Wemmer DE (2001) Crystal structure of an activated response regulator bound to its target. Nat Struct Mol Biol 81:52–56 13. Berendsen HJC, van der Spoel D, van Drunen R (1995) Gromacs: a message-passing parallel molecular dynamics implementation. Comput Phys Commun 91:43–46 14. Lindahl E, Hess B, van der Spoel D (2001) GROMACS 3.0: a package for molecular simulation and trajectory analysis. J Mol Model 7:306–317 15. Oostenbrink C, Villa A, Mark AE, van Gunsteren WF (2004) A biomolecular force field based on the free enthalpy of hydration and solvation: the GROMOS force-field parameter sets 53A5 and 53A6. J Comput Chem 25:1656–1676 16. Hess B, Bekker H, Berendsen HJC, Fraaije JGEM (1997) LINCS: a linear constraint solver for molecular simulations. J Comput Chem 18:1463–1472 17. Miyamoto S, Kollman PA (1992) Settle—an analytical version of the shake and rattle algorithm for rigid water models. J Comput Chem 13:952–962 18. Berendsen HJC, Postma JPM, Gunsteren WFV, Hermans J (1981) Interaction models for water in relation to protein hydration. In: Pullman B (ed) Intermolecular forces, vol 1. Reidel, Dordrecht, pp 331–342
J Mol Model (2014) 20:2212 19. Darden T, York D, Pedersen L (1993) Particle mesh Ewald: an N-log (N) method for Ewald sums in large systems. J Chem Phys 98: 10089–10092 20. Essmanm U, Perera L, Berkowitz ML, Darden T, Lee H, Pedersen LG (1995) A smooth particle mesh Ewald method. J Chem Phys 103: 8577–8593 21. Bussi G, Donadio D, Parrinello M (2007) Canonical sampling through velocity-rescaling. J Chem Phys 126:14101–14107 22. Berendsen HJC, Postma JPM, Gunsteren WFV, DiNola A, Haak JR (1984) Molecular dynamics with coupling to an external bath. J Chem Phys 81:3684–3690 23. DeLano WL, Bromberg S (2004) PyMOL user’s guide. DeLano Scientific LLC, San Francisco 24. R Development Core Team (2008) R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna 25. Wassenaar TA (2006) The molecular dynamics of sense and sensibility in processing and analysis of data. Faculty of Mathematics and Natural Sciences, University of Groningen, Groningen 26. Van Aalten DMF, Amadei A, Linsse ABM, Eijsink VGH, Vriend G, Berendsen HJC (1995) The essential dynamics of thermolysin: confirmation of the hinge‐bending motion and comparison of simulations in vacuum and water. Proteins 22:45–54
ORIGINAL PAPER
Molecular dynamics investigation of Helicobacter pylori chemotactic protein CheY1 and two mutants Ahmet Yildirim & Mustafa Tekpinar & Tsjerk A. Wassenaar
Received: 8 October 2013 / Accepted: 16 March 2014 # Springer-Verlag Berlin Heidelberg 2014
Abstract CheY is a chemotactic response regulator protein modulating the rotation direction of bacterial flagellar motors. It plays an important role in the colonization and infection of Helicobacter pylori (H. pylori), which is a common pathogen. Recently, the structure of CheY1 of H. pylori (HpCheY1) was solved, showing similarities and differences with CheY from E. coli. Here, we report 200 ns atomistic molecular dynamics (MD) simulations of HpCheY1 and two mutants. The results suggest that the surface of HpCheY1 has regions with increased affinity for Mg2+. In addition, wildtype HpCheY1 (WT HpCheY1) shows characteristic dynamics in helix 4, which is involved in FliM binding. This dynamics is altered in the D53A mutant and completely suppressed in the T84A mutant. The results are discussed in relation to the binding and function of HpCheY1. Keywords Chemotaxis protein . D53A mutant . H. pylori . Molecular dynamics simulation . Principal component analysis . T84A mutant
Introduction In 1982 Marshall and Warren identified H. pylori as a bacterium that grows in human gastric epithelial tissue and mucus A. Yildirim Department of Physics, Faculty of Science and Art, Siirt University, Siirt 56100, Turkey M. Tekpinar Department of Physics, Faculty of Science and Art, Yüzüncü Yıl University, Van 65100, Turkey T. A. Wassenaar (*) Department of Biology, University of Erlangen-Nürnberg, Staudtstr. 5, 91058 Erlangen, Germany e-mail: [email protected]
[1]. Its presence was found to be linked to gastritis, many forms of ulcers, as well as to gastric cancer [2–4]. Currently, it is known that H. pylori is present in approximately 50 % of the world population and that additional factors are involved in the development of H. pylori related diseases [5]. The bacterium is able to colonize and infect tissues through directed movement in response to a gradient of a chemical stimulus. This kind of motion, called chemotaxis, is affected through flagellae, which rotate by means of a flagellar motor [6]. A key role in the control of this motor in H. Pylori, as well as in other bacteria, is the chemotactic effector protein CheY. Upon activation through phosphorylation, this protein binds to the flagellar motor domain [7], causing the direction of rotation to invert, which changes the motion from 'running' to 'tumbling' [8]. In the case of H. pylori, the chemotaxis causes the bacterium to avoid the gastric acid environment and colonize the mucus layer, where the pH is higher [9]. Unlike E. coli, H. pylori has two CheY response regulators, here referred to as HpCheY1 and HpCheY2. HpCheY1 seems the main partner of FliM, while HpCheY2 acts as a phosphate sink to modulate the half-life of phosphorylated HpCheY1 [7, 10]. Recently, the crystal structure of activated HpCheY1 was solved. Beryllium fluoride (BeF3) was used as substitute for phosphate, because the lifetime of phosphate activated CheY1 is too short to allow structure elucidation. In addition to the structure of the wild-type protein (WT), crystal structures were solved for two mutations, namely D53A and T84A. The former of these cannot be activated by phosphorylation, while the latter displays increased activity [11]. Unfortunately, solving an apo-structure of HpCheY1 failed due to the sulfate present in the medium, which apparently also binds in the active site, resulting in an active conformation. From the available structures of E. coli CheY, which includes a FliM-bound structure [12], it is expected that the differences between the activated and the apo forms are small.
2212, Page 2 of 8
For such cases, MD simulations can be used to complement experimental methods and assess the conformation and dynamics of the unbound protein to contribute to a comprehensive characterization. In this work, we use atomistic scale MD simulations to characterize the conformation and dynamics of WT HpCheY1 and the two mutants D53A and T84A. D53 is the active residue for phosphorylation, and T84 is located near the binding site and hydrogen bonds to the phosphate. The results highlight regions on the surface of the protein with increased affinity for Mg2+ and show differences in the conformations as well as in the dynamics of wildtype and mutants. To compare the conformations and dynamics, an alternative approach of principal component analysis (PCA) was taken, using a newly developed, interactive method in PyMOL (manuscript in preparation). The results are discussed in relation to the mechanism of action of HpCheY1, with reflection on the role of dynamics.
Methods Simulations All simulations were performed with GROMACS version 4.5.4 [13, 14], in conjunction with the GROMOS96 53a6 united atom force field [15]. The integration time step used was 2 fs. Hydrogen bonds were treated as constraints using the LINCS algorithm [16], except for solvent molecules, which were constrained using the SETTLE algorithm [17]. Water was modeled using the simple point charge model [18]. Van der Waals interactions were calculated using a 1.4 nm cut-off. Electrostatic interactions were calculated using the particlemesh Ewald (PME) method [19, 20], with a grid spacing of 0.12 nm and a real space cut-off of 0.9 nm. All calculations were performed on a workstation with two Intel Xeon E52687 W processors and 64GB RAM. Starting structures were taken from Protein Data Bank with accession codes 3GWG (WT), 3H1F (D53A), and 3H1G (T84A) [11]. Bound water and Mg2+ ions present in the crystal structures were retained during the setup of the simulations. Structures were first energy minimized using a steepest descent algorithm, after which they were solvated in a periodic rhombic dodecahedron unit cell with a minimal distance between periodic images of 2 nm. The solvent consisted of a 0.15 M NaCl solution, corrected to compensate for the net charge of the protein. The resulting system for WT HpCheY1 contained 18,187 atoms in total, including 5645 water molecules, 17 sodium and 19 chloride ions. For the D53A mutant the solvated system measured 17,969 atoms in total, including 5574 water molecules, 17 sodium, and 18 chloride ions. The T84A mutant yielded a system consisting of 17,288 atoms in
J Mol Model (2014) 20:2212
total, including 5347 water molecules, 16 sodium and 18 chloride ions. The solvated systems were subjected to another round of steepest descent energy minimization, followed by 100 ps MD simulation in which position restraints were applied to the heavy atoms in the system. During this run, the system was weakly coupled to an external heat bath with a temperature of 300 K and a coupling time of 0.1 ps using a V-rescale thermostat [21]. This stage was followed by 100 ps position restrained MD in which the pressure was kept constant at 1 bar by coupling to a Berendsen barostat [22] with 1.0 ps and a compressibility of 4.5×10−5 bar. After this initial equilibration, productions runs were performed for a duration of 200 ns, with frames written to disk every 40 ps. Analysis Analysis was performed using a combination of the tools available from GROMACS, PyMOL [23] and auxiliary tools written for specific purposes. Principal component analysis (PCA) was performed using a new PyMOL module (manuscript in preparation). Post-hoc analysis of PCA results was performed in R [24]. The root mean square deviation To follow the relaxation over time, the root-mean-squared deviation (RMSD) was determined for each system. The RMSD was determined against the starting structure, as well as against the average structure from the last 100 ns and 20 ns of the simulations. The reason for this approach is that the RMSD with respect to the starting structure does not usually allow assessment of relaxation, as any plateau in RMSD corresponds to a volume of conformational space, with a size proportional to the RMSD. In other words, there may be drift in conformation, which is not reflected by a change in RMSD. Calculating the RMSD against the average structure does allow assessment of relaxation, because the average structure is an estimator of the mean structure, which is an ensemble property, and because the direct neighborhood of the average structure is associated with lower RMSD values and consequently smaller associated volume elements, resulting in a more direct relation between RMSD and conformational state. The RMSD of a given structure with respect to a reference structure rref, is calculated as RMSDðt Þ ¼
2 1=2 1 XN ref m r ð t Þ−r ; i i i¼1 i M
ð1Þ
where M=∑mi and ri(t) is the position of atom i at time t after least square fitting the structure to the reference structure.
J Mol Model (2014) 20:2212
Page 3 of 8, 2212
Distributions of Mg2+ For visualizing Mg2+ distributions around HpCheY1, first the molecular shaped unit cell representation was determined, which means that every particle is positioned in the periodic system where it is closest to a given reference group, in this case the protein [25]. Subsequently, all frames were fit to the WT starting structure, and the positions of the ions were projected around the reference structure.
Principal component analysis PCA in MD simulations starts with the construction of the positional covariance matrix. If x denotes a single structure consisting of n atoms represented as a single column vector of length 3n, then combining the structures sampled over time gives a data matrix X of dimensions 3n × m, with m being the number of frames. From this frame the covariance matrix can be calculated as S¼
1 −1 X I−ð j0 jÞ jj0 X0 : m−1
ð2Þ
Here I is the m × m identity matrix and j is a vector of ones with length m. The matrix (I−(j′j)−1jj′) is the idempotent centering matrix, and an apostrophe is used to indicate transposition. The covariance matrix S is decomposed into a matrix of eigenvectors P and a diagonal matrix of eigenvalues D to obtain the principal components: S ¼ PDP0 :
ð3Þ
Each principal component corresponds to a linear correlation between atomic displacements and can be regarded as a collective motion, although the real underlying motion may be non-linear and will then be dispersed over a number of components. We introduce here an alternative, hierarchical approach to use PCA more effectively to study conformation and dynamics in relation to mutations. The first stage consists of determining the average structures. These are used to determine the plane through conformational space in which the distinct states are most dispersed. The relaxation of the systems toward the respective end states then follows from the projections of the trajectories onto this plane. The second stage consists of PCA applied to each of the systems separately. This approach allows constructing a comprehensive framework to assess similarity and differences in structure and dynamics of different mutants.
PCA was applied to the last 100 ns of the simulations, as well as to the last 20 ns, allowing assessment of the robustness of the results in relation to the length of the trajectory.
Results and discussion The root mean square deviation To assess the stability of each structure, the RMSD was calculated for backbone atoms and for heavy atoms against the starting structure and against the mean structures obtained from the last 100 ns and the last 20 ns of the simulations as a function of time. These graphs are shown in Fig. 1. The RMSD against the starting structure, for both backbone (Fig. 1a) and all non-hydrogen atoms (Fig. 1b), increases over a time span of 15 ns to values of approximately 0.25 nm (backbone) or 0.35 nm (heavy atoms). The RMSD against the average obtained from the last 100 ns shows a steady decrease over 75 ns with little difference between backbone (Fig. 1c) and all heavy atoms (Fig. 1d). The T84A mutant reaches a stable state with little fluctuation, while the D53A mutant displays more flexibility. The wild type is the most dynamic, displaying a higher RMSD and showing larger fluctuations. A similar view is obtained from the RMSD against the average from the last 20 ns, which shows little difference with the other graphs for D53A and T84A, but displays significant fluctuations for wild type, temporarily reaching RMSD values against the average close to the value observed for the starting structure. Binding and distribution of magnesium The crystal structures of WT HpCheY1 and T84A each contain two magnesium ions. The D53A mutant contains only a single magnesium ion, as the binding site at D53 is removed. In Fig. 2, the positions of the magnesium ions over the whole length of the simulation are shown for all three systems to investigate the stability and distribution of Magnesium in and around HpCheY1. In the case of WT and T84A one ion remained bound for the whole length of the simulation, shown in red in Fig. 2. The other ion, as well as the single ion in the system containing the D53A mutant, diffused over the surface of the protein. This ion is shown in blue in Fig. 2. In the WT, there appear to be four sites on the protein surface with increased density, which are labeled B-E in Fig. 2. Removal of the D53 side chain causes a shift in density from site C to B, but also shifts the density from site E to an alternative site F near helix 5. In addition, a small cluster of density is observed above site B, which is labeled A, and which is located close to the binding site in the protein. Apparently, the loss of the internal magnesium ion causes a shift in the distributions in
2212, Page 4 of 8
J Mol Model (2014) 20:2212
Fig. 1 Root mean-square deviation (RMSD) over time. RMSD against starting structure for backbone atoms (a) and all heavy atoms (b). RMSD against average structures from the last 100 ns of each simulation for backbone atoms (c) and all heavy atoms (d). RMSD against average
structures from the last 20 ns of each simulation for backbone atoms (e) and all heavy atoms (f). Colors are black for wild type, magenta for D53A, and green for T84A
the direction of the binding site, suggesting a mechanism for catching and incorporating magnesium.
The three mean structures span a plane in conformational space. Two orthonormal vectors in this plane are easily obtained by performing PCA on those mean structures. Subsequently, the trajectories can be projected onto these vectors, which then allows assessment of the similarity or otherwise of the relaxation pathways and the ensembles. These projections are shown for the planes through the means from the last 100 ns and the last 20 ns in Fig. 4, where the starting structure and mean structures are indicated with large dots. From this Fig., it can be seen that the three structures relax in different directions and sample distinct areas in conformational space. Interestingly, the two mutants relax in opposite directions from the starting structure, while the WT relaxes in a direction almost orthogonal to the axis D53A-T84A. This is true for both views. Previously, it has been suggested that investigation of differences in structure and dynamics could be achieved by performing PCA on a combined trajectory [26]. However, that is only true if the dynamics orthogonal to the relaxation has much smaller variance than the dynamics in the plane of the
Conformations and dynamics To investigate the effects of the mutations on the conformation of HpCheY1, the mean structures from the last 100 ns and the last 20 ns of the simulations were determined, the latter of which are shown in Fig. 3. These structures show that there are small, but significant changes in the conformations, in particular in helix 4. In the T84A mutant, shown in green, this helix becomes more ordered compared to WT, while in the D53A mutant it becomes more disordered. In addition, the latter shows a shift in helices 1 and 5. Interestingly, the conformation of these helices in the D53A mutant is more similar to the crystal structures, which means that the helices are shifted in WT and T84A during the relaxation. The increased ordering of T84A is caused by E90, which is internalized to compensate for the loss of the threonine hydroxyl group. As a consequence, the loop K88-A89-E90-A91 becomes part of helix 4.
J Mol Model (2014) 20:2212
Fig. 2 Distributions of magnesium over time in and around HpCheY1 wildtype and mutants. WT HpCheY1 (a), D53A mutant (b), and T84A mutant (c). A total of six sites on the surface with increased Mg2+ density are observed, which are labeled A–F. The internal Mg2+ ion in WT (a) and T84A (c) is colored red and does not leave the binding site. In the WT, the main site on the surface is D. Removal of the internal binding site by mutation of D53 to alanine causes the density to shift toward B and A, as well as to a site F on the other side. Mutation of T84 to alanine retains the internal binding site, but also shifts the surface binding to site B
Fig. 3 Stereo image of HpCheY1 mean structures after equilibration. The image shows the mean structures from the last 20 ns of a 200 ns simulation of WT HpCheY1 (gray), D53A mutant (purple), and T84A mutant (green). The helices are numbered from N-terminus to C-
Page 5 of 8, 2212
relaxation, and the two are uncorrelated. In most cases, one of these conditions is not satisfied and the modes obtained mix relaxation and steady-state dynamics, making it difficult to comment on either. The procedure used here solves this issue and provides a comprehensive framework for investigating differences between states in terms of both structures and dynamics. From Fig. 4 it is clear that each system ends up in a distinct state. The overall picture does not change much when the last 100 ns is used for the analysis (Fig. 4a) or only the last 20 ns (Fig. 1b). In particular, the profiles for D53A and T84A hardly change, in agreement with the view obtained from the RMSD. For WT there is a difference in the projections of the trajectory around the mean. When using the mean structure from the last 100 ns, the projections are more dispersed and several clusters can be identified, which is not the case when using the average structure from the last 20 ns. This suggests that the last 20 ns alone are not sufficiently representative for the WT equilibrium ensemble, which is in accordance with the view obtained from the RMSD curves. For a first characterization of the equilibrium dynamics, ellipses are drawn around each mean structure, summarizing the in-plane dynamics of the last 100 (yellow) or 20 (orange) ns of each system. Again, D53A and T84A show significant overlap of the ellipses, suggesting that the collective motions in both regimes are similar, whereas the ellipses for WT show less overlap, indicating different mean projections and different main directions of motion. These ellipses also show that WT and D53A have differently oriented main radii, reflecting differences in the collective motions. In the T84A mutant, there seems to be little overall motion after relaxation, suggesting that the removal of the threonine side chains rigidifies the structure. Because the collective motions appear different in nature for each of the mutants, we decided to apply PCA to each system separately, of which the results are shown in Fig. 5. PCA was applied to the last 100 and last 20 ns from each
terminus. Residue E90 is highlighted with spheres. The image shows an extended helix 4 in T84A due to inward positioning of E90. In addition, helices 1 and 5 are shifted in D53A
2212, Page 6 of 8
J Mol Model (2014) 20:2212
Fig. 4 Projections of trajectories on the plane through mean structures. Projections on the plane determined from the mean structures from the last 100 ns (a) and on the plane determined from the last 20 ns (b). In both panels the starting structure is marked with s. The three 200 ns simulations are shown in gray/ black (wildtype), violet/purple (D53A mutant), and green/forest (T84A mutant). The darkest, thickest line marks the last 20 ns of each simulation, whereas the intermediate line marks the trajectory from 100 to 180 ns. The respective mean structures are indicated with a yellow dot, while orange ellipses are drawn to mark the principal directions of the motions in the plane in the last 100 ns, and yellow ellipses mark the principal directions in the plane in the last 20 ns
trajectory, corresponding to the thicker, darker parts of the trajectories drawn in Fig. 4. The traces of the covariance matrices from the last 100 ns and the last 20 ns are given in Table 1. The results from the last Fig. 5 Collective dynamics in WT HpCheY1 and mutants D53A and T84A. Visualization of the first principal component from simulations of HpCheY1. The loadings are drawn as rainbowcolored bicones, centered on the respective mean structures. The residues corresponding to D53 (magenta) and T84 (green) are highlighted and the side chain of E90 is shown as spheres. The images show the first principal component determined from the last 100 ns simulation for WT (a), D53A (b), and T84A (c), as well as the first component determined from the last 20 ns of simulation for WT (d), D53A (e), and T84A (f)
100 ns show that the wild type is the most dynamic, although the effect is more severe for T84A. If the PCA is performed on the last 20 ns only, the total variance for wild type is found to be lower than for D53A, but the values for D53A and T84A
J Mol Model (2014) 20:2212 Table 1 Total variance and first three eigenvalues
Page 7 of 8, 2212
Time
Structures
Total variance 2
Last 100 ns
Last 20 ns
WT D53A T84A WT D53A T84A
are unchanged. This actually shows that the part of the simulation, even when it is commonly considered to be in an equilibrium state, still has a large influence on the results from analysis of conformations and dynamics. The mean structure (Fig. 3) suggests that the decreased motility of T84A is caused by the extension and stabilization of helix 4, as discussed above. The first principal component obtained from T84A also only captures 8 % of the total variance, which suggests that the dynamics is mostly governed by thermal fluctuations, rather than by collective dynamics. For both WT and D53A, approximately 25–30 % of the variance is captured by the first component, suggesting a certain degree of collectivity. The collectivity can be assessed from Fig. 5. The top panel shows that in WT there is collectivity in the motions of helix 4 and the preceding loop. In particular, there seems to be a certain degree of rotational motion of the helix, related to the dynamics of the loop. This is observed in the analysis of the last 100 ns as well as of the last 20 ns, showing that, despite the differences, there is qualitative agreement between the regimes. In the D53A mutant, the rotational motions are largely absent. Interestingly, the region that shows collective dynamics in wild type, helix 4, forms part of the FliM interacting region. This suggests that the altered dynamics profiles of the D53A and T84A mutants could explain part of the differences observed in vivo, aside from effects on the phosphate binding. This applies in particular to T84A, which displays higher activity. In turn, this raises the suggestion that modulation of the dynamics by binding of phosphate plays a role in the binding of HpCheY1 and FliM, or even more broadly in binding of CheY to FliM in bacteria. This possibility will be further explored in a future study.
Conclusions We performed MD simulations to characterize the conformation and dynamics of apo WT HpCheY1, D53A mutant, and T84A mutant. The results suggest that the apo structures relax to distinct states, mainly characterized by changes in the FliM binding region helix 4. The T84A mutant is characterized by lengthening of the helix and stabilization due to internalization
Eigenvalue 1 2
Eigenvalue 2 2
Eigenvalue 3
nm
nm (%)
nm (%)
nm2 (%)
29.85 19.06 13.73 17.23 19.76 13.18
10.33 (35 %) 2.98 (16 %) 0.83 (6 %) 4.28 (25 %) 4.58 (23 %) 1.11 (8 %)
3.34 (11 %) 1.39 (7 %) 0.74 (5 %) 0.97 (6 %) 2.28 (12 %) 0.97 (7 %)
1.32 (4 1.23 (6 0.64 (5 0.90 (5 1.04 (5 0.70 (5
%) %) %) %) %) %)
of E90 to compensate the loss of the T84 hydroxyl group in the active site. The different states are also characterized by differences in the dynamics, with the wild type showing a certain degree of collectivity of backbone and side chain motions in helix 4 and the preceding loop. These motions are suppressed in the D53A mutant and diminished in the T84A mutant. Future studies should show whether such changes in dynamics are also invoked by phosphate activation and what role this dynamics plays in FliM binding and CheY mediated flagellar motor control.
References 1. Marshall BJ, Warren JR (1984) Unidentified curved bacilli in the stomach of patients with gastritis and peptic ulceration. Lancet 1: 1311–1315 2. Labigne A, de Reuse H (1996) Determinants of Helicobacter pylori pathogenicity. Infect Agents Dis 5:191–202 3. McColl KEL (1997) Helicobacter pylori: clinical aspects. J Infect 34: 7–13 4. Riegg SJ, Dunn BE, Blaser MJ (1995) Microbiology and pathogenesis of Helicobacter pylori. In: Blaser MJ et al (eds) Infections of the gastrointestinal tract. Raven, New York, pp 535–550 5. Covacci A, Telford JL, Del Giudice G, Parsonnet J, Rappuoli R (1999) Helicobacter pylori virulence and genetic geography. Science 284:1328–1333 6. Scharf BE, Fahrner KA, Turner L, Berg HC (1998) Control of direction of flagellar rotation in bacterial chemotaxis. Proc Natl Acad Sci 95:201–206 7. Foynes S, Dorrell N, Ward SJ, Stabler RA, McColm AA, Rycroft AN, Wren BW (2000) Helicobacter pylori possesses two CheY response regulators and a histidine kinase sensor, CheA, which are essential for chemotaxis and colonization of the gastric mucosa. Infect Immun 68:2016–2023 8. Berg HC, Brown DA (1972) Chemotaxis in E. coli analysed by threedimensional tracking. Nature 239:500–504 9. Schreiber S, Konradt M, Groll C, Scheid P, Hanauer G, Werling HO, Josenhans C, Suerbaum S (2004) The spatial orientation of Helicobacter pylori in the gastric mucus. Proc Natl Acad Sci 101: 5024–5029 10. Jimenez-Pearson MA, Delany I, Scarlato V, Beier D (2005) Phosphate flow in the chemotactic response system of Helicobacter pylori. Microbiology 151:3299–3311 11. Lam KH, Ling TK, Au SW (2010) Crystal structure of activated CheY1 from Helicobacter pylori. J Bacteriol 192:2324–2334
2212, Page 8 of 8 12. Lee SY, Cho HS, Pelton JG, Yan D, Henderson RK, King DS, Huang LS, Kustu S, Berry EA, Wemmer DE (2001) Crystal structure of an activated response regulator bound to its target. Nat Struct Mol Biol 81:52–56 13. Berendsen HJC, van der Spoel D, van Drunen R (1995) Gromacs: a message-passing parallel molecular dynamics implementation. Comput Phys Commun 91:43–46 14. Lindahl E, Hess B, van der Spoel D (2001) GROMACS 3.0: a package for molecular simulation and trajectory analysis. J Mol Model 7:306–317 15. Oostenbrink C, Villa A, Mark AE, van Gunsteren WF (2004) A biomolecular force field based on the free enthalpy of hydration and solvation: the GROMOS force-field parameter sets 53A5 and 53A6. J Comput Chem 25:1656–1676 16. Hess B, Bekker H, Berendsen HJC, Fraaije JGEM (1997) LINCS: a linear constraint solver for molecular simulations. J Comput Chem 18:1463–1472 17. Miyamoto S, Kollman PA (1992) Settle—an analytical version of the shake and rattle algorithm for rigid water models. J Comput Chem 13:952–962 18. Berendsen HJC, Postma JPM, Gunsteren WFV, Hermans J (1981) Interaction models for water in relation to protein hydration. In: Pullman B (ed) Intermolecular forces, vol 1. Reidel, Dordrecht, pp 331–342
J Mol Model (2014) 20:2212 19. Darden T, York D, Pedersen L (1993) Particle mesh Ewald: an N-log (N) method for Ewald sums in large systems. J Chem Phys 98: 10089–10092 20. Essmanm U, Perera L, Berkowitz ML, Darden T, Lee H, Pedersen LG (1995) A smooth particle mesh Ewald method. J Chem Phys 103: 8577–8593 21. Bussi G, Donadio D, Parrinello M (2007) Canonical sampling through velocity-rescaling. J Chem Phys 126:14101–14107 22. Berendsen HJC, Postma JPM, Gunsteren WFV, DiNola A, Haak JR (1984) Molecular dynamics with coupling to an external bath. J Chem Phys 81:3684–3690 23. DeLano WL, Bromberg S (2004) PyMOL user’s guide. DeLano Scientific LLC, San Francisco 24. R Development Core Team (2008) R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna 25. Wassenaar TA (2006) The molecular dynamics of sense and sensibility in processing and analysis of data. Faculty of Mathematics and Natural Sciences, University of Groningen, Groningen 26. Van Aalten DMF, Amadei A, Linsse ABM, Eijsink VGH, Vriend G, Berendsen HJC (1995) The essential dynamics of thermolysin: confirmation of the hinge‐bending motion and comparison of simulations in vacuum and water. Proteins 22:45–54
