QM/MM MD simulations of La(III)–phosphopeptide complexes

Abstract

Several bioanalytical enrichment techniques are based on the interactions of phosphopeptides with Ln(III) ions. In order to gain an improved understanding of these complexes and the respective ion–peptide interactions, hybrid quantum mechanics/molecular mechanics (QM/MM) molecular dynamics (MD) simulations of La(III) coordinating to the phosphopeptide VPQLEIVPNSpAEER were conducted. Simulations of di- as well as monoanionic phosphate groups were carried out. The La(III) ion and its first hydration layer, including the sidechain of the phosphoserine residue were treated quantum mechanically at RI-MP2/triple zeta level, whereas the remaining part of the system was treated with classical potentials. The simulation of the dianionic phosphopeptide revealed a 9-fold coordinated La(III) ion, with the phosphopeptide binding bi- as well as monodentate. The mean residence times (τ) of the first shell water molecules were 82 ps and 37 ps for the bi- and monodentate complexes, respectively, which is much higher compared to free La(III) in aqueous solution (τ = 17 ps). The simulation of the monoanionic La(III)–phosphopeptide complex revealed a bidentate coordination throughout the 80 ps sampling period. An intramolecular hydrogen bond between the hydrogen of the phosphate group and the backbone was observed and a τ value of 14 ps was obtained, which is much lower as for the dianionic complex.

---

1 Introduction

The reversible phosphorylation of phosphoproteins is an important regulatory mechanism for a number of cellular processes and therefore, the identification and characterization of phosphoproteins/phosphopeptides is an important issue.¹ Despite the fast development of mass spectrometric methods, the analysis of phosphopeptides remains challenging due to suppression effects and, therefore, selective enrichment steps prior to measurements are generally applied.² In recent years, several enrichment strategies based on the interactions of La(III) ions with phosphopeptides and phosphoproteins have been developed, showing exceptionally high selectivities.^3–9 Especially the co-precipitation of phosphoproteins^3,4 and phosphopeptides⁷ by lanthanoid(III) ions is a promising alternative to conventional approaches that are mainly based on stationary phases. However, structure and dynamics of La(III)–phosphopeptide complexes in aqueous solution are by now poorly understood.

Phosphopeptides and phosphoproteins have been studied by Fourier transform infrared spectroscopy (FTIR),^10,111H nuclear magnetic resonance^12,13 as well as molecular mechanics (MM) molecular dynamics (MD) simulations.¹⁴ Sanchez-Ruiz et al.¹⁰ reported a pK_a of about 6 for the phosphate monoester in O–phosphoserines by carrying out FTIR measurements. Fernández et al.¹¹ made similar observations for α-casein as a change in pH from 6.8 to 2.0 resulted in a dramatic reduction of the dianionic symmetric stretching peak at 976 cm⁻¹.

Hydrated lanthanum(III) was studied by EXAFS measurements¹⁵ as well as by a quantum mechanical charge field (QMCF) MD simulation¹⁶ and coordination numbers of around 9 were reported by both groups.

MD simulations are a powerful tool to study systems like La(III)–phosphopeptide complexes. However, the theoretical treatment of such highly charged systems requires quantum mechanical (QM) methods to achieve a sufficiently accurate description of the molecular interactions, since polarization and charge transfer effects play a crucial role.^17,18 Unfortunately, QM calculations are restricted to relatively small systems and hence combined quantum mechanical/molecular mechanical (QM/MM) techniques became a popular and highly successful approach.^19,20 Recently, the importance of QM/MM methods was underlined by the awarding of the Nobel prize in chemistry to Martin Karplus, Michael Levitt and Arieh Warshel.²¹ These methods are based on a partitioning of the molecular system into two region, where different levels of theory are applied: the chemically most relevant part is treated by computationally demanding quantum mechanical methods, whereas interactions in the surrounding region are treated with empirical potentials. Therefore, these methods enable a modeling of large molecular systems at reasonable cost while still providing high accuracy at the region of interest. The main challenge of this partitioning is the QM/MM boundary, especially if covalent bonds have to be cut. Two frequently applied approaches to this problem are the link-atom approach and the local self consistent field (LSCF) formalism, with each exhibiting its strengths and drawbacks.^22,23

In the present study, QM/MM MD simulations of the complex of La(III) with a phosphopeptide were conducted in order to gain an improved understanding of La(III) based bioanalytical enrichment techniques. The sequence of the studied phosphopeptide is VPQLEIVPNSpAEER. This peptide was chosen for this study as it is present in tryptic digests of α-casein, which are frequently used for bioanalytical studies.

2 Methods

2.1 Simulation protocol

2.1.1 Dianionic La(III)–phosphopeptide complex. As a starting structure, a cubic, periodic simulation box containing a La(III) ion, a phosphopeptide and 5312 water molecules has been set up. The density was fixed at 0.997 g cm⁻³ and periodic boundary conditions were employed. The aminoacid sequence of the phosphopeptide was VPQLEIVPNSpAEER (NVAL PRO GLN LEU GLH ILE VAL PRO ASN SER(P) ALA GLH GLH CARG), with the glutamate residues and the carboxyl terminus being protonated. Sp is the dianionic phosphoserine residue which coordinates bidentate to the La(III) ion. For the generation of the system's topology, the Amber leap software²⁴ was used.

The Berendsen algorithm²⁵ was applied with a relaxation time of 0.1 ps in order to keep the temperature at around 298 K. The velocity Verlet/RATTLE algorithm was used to integrate the equations of motions,²⁶ with constraints being applied to all bonds involving hydrogen atoms. The Coulombic cutoff was set to 15 Å and the long-range electrostatic interactions were considered via the reaction field method (ε = 78.36).²⁷ The FF99SB force field²⁸ was applied to all standard aminoacid residues and the potentials for phosphoserine²⁹ as well as for La(III)^30,31 were taken from literature. For water, the extended simple point charge (SPC/E) water model³² was used.

The La(III) ion has been placed in the center of the QM region and all water molecules within a radius of 4.2 Å of La(III) were treated quantum mechanically. To ensure a smooth transition between QM and MM region, a smoothing function³³ was applied to all particles located in the so called smoothing region (0.2 Å thick). The phosphoserine residue was further partitioned between C_α and C_β and only the resulting methylphosphate was treated quantum mechanically. The link-atom approach has been employed to satisfy unsaturated valence and the respective position was determined via:

r_CH = ρ(r_Cβ − r_Cα) (1)

To evaluate adequate settings for the link bond, QM/MM and full QM calculations of the model system ACE-SER(P)-NME coordinated to a La(III) ion were carried out (Fig. 1a). The corresponding QM setup is shown in Fig. 1b: due to the proximity of the C_α atom, it is excluded when computing the interaction between the QM system and the MM partial charges. Fig. 1c displays the energy obtained upon stretching of the QM/MM link bond for a full QM (black) and the QM/MM setup (red) at MP2 level, using an ideal factor ρ identified as 0.7135. If a properly adjusted MM potential is applied along the link bond (blue), the shape coincided well with the full QM treatment in the region close to the minimum.

Fig. 1 (a) Model system ACE-SER(P)-NME coordinated to a La(III) ion to determine the settings for the QM/MM link bond. (b) Corresponding QM setup: due to the proximity of the C_α atom, it is excluded when computing the interaction between the QM system and the MM partial charges. (c) Energy obtained upon stretching of the QM/MM link bond r_Cα–Cβ for a full QM (black) and the QM/MM setup (red) at MP2 level. If a properly adjusted MM potential is applied along the link bond (blue), the shape coincides well with the full QM treatment in the region close to the minimum.

The quantum mechanical calculations were conducted with the Turbomole program,³⁴ using the RI-MP2³⁵ method. The def-TZVP basis set³⁶ was applied to carbon, hydrogen, oxygen and phosphorus and the def2-TZVP³⁷ basis set with the def2-ECP^38,39 was used for La(III). The partial charges of the MM atoms were taken into account via a pertubational term in the core Hamiltonian, however C_α was not considered in order to avoid an unphysical polarization along the linker bond. To generate an appropriate starting structure, an MM simulation was carried out for 25 ns with the previously mentioned force field parameters. Prior to the 80 ps sampling period, the system was equilibrated for 20 ps. The simulation was carried out with the QMCF MD software package.^33,40,41

2.1.2 Monoanionic La(III)–phosphopeptide complex. A cubic simulation box with 5530 water molecules, a La(III) ion and a phosphopeptide with the sequence VPQLEIVPNSphAEER was taken as a starting structure. The glutamate residues and the carboxyl terminus is protonated and Sph is the monoanionic phosphoserine residue which coordinates bidentate to the La(III) ion. All other simulation settings are identical with the aforementioned settings.

2.1.3 Free La(III) in aqueous solution. A simulation box with one La(III) ion and 2000 water molecules was taken as starting structure and the identical simulation settings as chosen for the La(III)–phosphopeptide simulations were applied.

2.2 Analysis

The pair radial distribution function (RDF) gives the probability of finding an atom at a distance r from an other atom compared to the ideal gas distribution.⁴²ρ is the density of the ideal gas and N_bin are distance “bins” in which the number of certain particles are averaged over the entire simulation:where N_t is the number of timesteps and r_b and r_a are the inner and outer borders of these bins. In the present study, the RDFs were calculated between the La(III) ion and the oxygens of the water molecules.

The average coordination numbers (CN_av) were evaluated by a numerical integration via the trapezoidal rule⁴³ of all “bins” from 0 to the first minimum of the corresponding RDF (3.7 Å).

The angular distribution functions (ADF) are histograms based on ligand–ion–ligand angles within a certain radius. La(III) was considered as the center and all oxygens within the first minima of the corresponding RDF (3.7 Å) were considered as ligands.

Spectroscopic data was extracted by analyzing the autocorrelation of the particles' velocity vectors viaeqn (5), where N_t is the number of time origins in the time t_i utilized for the autocorrelation of all the velocity components v of a particle j while N is the number of particles under investigation.

The mean residence time (τ), indicating the average coordination time, was evaluated by the direct method:⁴⁴

where t_sim is the simulation time, N_av is the average number of particles in the respective shell and N^0.5_ex is the number of exchanges with a ligand displacement time of more than 0.5 ps, which is the mean lifetime of a hydrogen bond in pure water, determined by femtosecond laser pulse spectroscopy.⁴⁵ In the present study we considered solely the water molecules in the τ evaluations and defined the first shell in-between radii of 0.0 and 3.7 Å from the La(III) ion.

The Root-mean-square deviation was calculated via

where N is the number of atoms and δ are the distances between the positions of the aligned atoms at time t and time t₀ (starting structure). For the evaluation of the RMSD in the present study, solely the backbone atoms were considered.

3 Results and discussion

The pK_a value of the phosphate monoester in phosphoserines is about 6.¹⁰ Therefore, depending on the pH as well as on the La(III) concentration, phosphopeptides might be monoanionic or dianionic in aqueous solutions containing La(III) ions. As enrichment techniques are often applied at low pH, we considered the non-coordinating sites as protonated (carboxylate groups) and conducted simulations with La(III) coordinating towards a protonated as well as a deprotonated phosphoserine residue.

3.1 Dianionic La(III)–phosphopeptide

In the first 50 ps of the simulation the phosphate group coordinated bidentately to the La(III) ion (Fig. 2a). Subsequently, a monodentate complex (Fig. 2b) was formed and observed until the end of the sampling period. The mono- and bidentate structures have been evaluated separately and are denoted as La–pp_mono and La–pp_bi.

Fig. 2 Snapshot taken from the QM/MM MD simulation of (a) the dianionic bidentate coordinating La(III)–phosphopeptide complex and (b) the dianionic monodentate coordinating La(III)–phosphopeptide complex. Cyan: C; white: H; blue: N; red: O; yellow: P; gray: La(III).

Average coordination numbers (CN_av) of 8.9 and 9.0 were obtained for the bi- and monodentate complexes, respectively (Table 1), which are higher as for free La(III) in aqueous solution ([La(H₂O)_n]³⁺) where a CN_av of 8.5 was found. Fig. 3 depicts the La–O_wat RDFs for both La(III)–phosphopeptide complexes as well as for free La(III) in aqueous solution. Sharp peaks ranging from 2.3 to 3.2 Å indicate the first hydration layers. Compared to free La(III) in aqueous solution, the phosphopeptide complexes exhibit more narrow peaks, which are slightly shifted to larger distances. Furthermore, the peak of La–pp_bi is more narrow and slightly shifted to a shorter distance compared to the peak of La–pp_mono. The integrations of the first peaks of the RDFs are around 7 and 8 for La–pp_bi and La–pp_mono, respectively, which correlates with CN_av, considering the coordination sites of the phosphopeptide. Well defined peaks are observed at a distance of 3.9 to 5.5 Å, indicating the second hydration layers. Similar to the first hydration shell, slight differences were found between La–pp_mono, La–pp_bi and [La(H₂O)_n]³⁺. The water molecules of the second hydration layers of the La(III)–phosphopeptide complexes have a slightly larger distance from the center compared to the second shell water molecules of [La(H₂O)_n]³⁺.

Table 1 Mean residence times of first shell water molecules (τ), average number of coordinating oxygens (CN_av) and first peak maxima of the La–O_wat radial distribution functions of free La(III) in aqueous solution (La), dianionic bidentate coordinating La(III)–phosphopeptide complex (La–pp_bi), dianionic monodentate coordinating La(III)–phosphopeptide complex (La–pp_mono) and the monoanionic La(III)–phosphopeptide complex (La–Hpp) obtained from the QM/MM MD simulation

	τ (ps)	CNav	d La–Owat (Å)
La	17	8.5	2.58
La–ppbi	82	8.9	2.61
La–ppmono	37	9.0	2.62
La–Hpp	14	8.8	2.61

---

Fig. 3 La–O_wat radial distribution functions (solid lines) and corresponding integrations (dashed lines) of the monoanionic La(III)–phosphopeptide complex (red), free La(III) in aqueous solution (green), dianionic monodentate coordinating La(III)–phosphopeptide complex (blue) and dianionic bidentate coordinating La(III)–phosphopeptide complex (black) obtained from the QM/MM MD simulations.

Fig. 4 shows the O–La–O ADFs, with O being all oxygen atoms within a distance of 3.7 Å. The relatively narrow peaks with maxima around 71° and 137° indicate a quite rigid coordination structure, with geometries in-between a tri-capped trigonal prism and a capped square antiprism. The bidentate coordinating complex shows a shoulder at 59°, which results from the O_phosphate–La–O_phosphate angle. Consequently, the 7 water molecules in the first hydration layer of La–pp_bi exhibit larger O_wat–La–O_wat angles than the 8 water molecules in the first shell of La–pp_mono, which is indicated in the ADFs by a slight shift of the peak at 71° to larger angles (Fig. 4).

Fig. 4 O–La–O angular distribution functions, with O being all atoms within a distance of 3.7 Å to the La(III) ion, of the monoanionic La(III)–phosphopeptide complex (red), free La(III) in aqueous solution (green), dianionic monodentate coordinating La(III)–phosphopeptide complex (blue) and dianionic bidentate coordinating La(III)–phosphopeptide complex (black) obtained from the QM/MM MD simulation.

The average partial charge (Mulliken) of the La(III) ion is 2.03, 2.10, and 2.15 for La–pp_bi, La–pp_mono and [La(H₂O)_n]³⁺, respectively. The reduction of the La(III) partial charge is a result of the phosphate group-ion charge transfer. As can be expected this effect is more pronounced in the bidentate coordinating complex.

The τ values for all oxygens within a radius of 3.7 Å are 82 ps, 37 ps and 17 ps for La–pp_bi, La–pp_mono and [La(H₂O)_n]³⁺, respectively. It should be emphasized that these τ values, especially the higher ones, might have some statistical uncertainties. The exceptionally high τ of La–pp_bi indicates a more stable 9-fold coordinated structure compared to free La(III) in aqueous solution and this may be attributed to the relatively small O_phosphate–La–O_phosphate angle (Fig. 4), which results in a lower nuclei repulsion of the 7 first shell water molecules. The τ value for La–pp_mono is, although not as high as for La–pp_bi, higher as for [La(H₂O)_n]³⁺. In this complex, a further stabilizing effect might arise from the hydrogen bonds between the first shell water molecules and the oxygens from the phosphate group.

IR spectra were calculated via VACFs and Table 2 lists the most prominent P–O and La–O_wat bands. The La–O_wat bands of La–pp_mono and La–pp_bi differ only for about 10 cm⁻¹ and are slightly lower as for [La(H₂O)_n]³⁺. The P–O bands, which result from the asymmetric stretching mode, were found at 1090 cm⁻¹ and 1060 cm⁻¹ for La–pp_mono and La–pp_bi, respectively. The experimentally obtained asymmetric P–O stretching frequency of the dianionic phosphoserine is 1091 cm⁻¹.¹⁰

Table 2 The most prominent La–O_wat and P–O bands from the power spectra, calculated via velocity autocorrelation functions

	v La–Owat (cm^−1)	v P–O (cm^−1)
La	180
La–ppbi	150	1060
La–ppmono	160	1090
La–Hpp	160	1210

---

During the whole MM MD as well as QM/MM MD simulation, the La(III) ion interacted solely with the phosphate group of the studied phosphopeptide, whereas other aminoacid sidechains (e.g. glutamate) did not interact with the ion. This is in agreement with the experimental observation that Ln(III)-based enrichment techniques demonstrate high selectivities towards phosphorylated residues while non-phosphorylated peptides have been rarely detected using these approaches.^7,8 Therefore, interactions of La(III) ions with other monophosphorylated peptides might show similar binding properties. Multiphosphorylated peptides, however, could coordinate with more than one phosphate group toward the ion, possibly triggering a refolding of the peptide.

The root-mean-square deviation (RMSD) of the dianionic bidentate complex fluctuates about a constant value of 0.83 Å throughout the sampled 40 ps. The dianionic monodentate complex, however, shows a slight rearrangement after 15 ps, indicated by a small step in the RMSD plot (Fig. 5a). Such changes are in general a consequence of rotations of the backbone. The time series of the C_C(GLH₁₂)–N(GLH₁₃)–C_α(GLH₁₃)–C_C(GLH₁₃) dihedral angle (ϕ) and the running average is shown in Fig. 5(b) and indicate that the change in the RMSD plot results from a rotation around the N(GLH₁₃)–C_α(GLH₁₃) bond.

Fig. 5 (a) Root-mean-square deviation (RMSD) in Å and (b) time series of the C_C(GLH₁₂)–N(GLH₁₃)–C_α(GLH₁₃)–C_C(GLH₁₃) dihedral angle (ϕ) (black) as well as the corresponding running average of 1000 configurations (red) obtained from the QM/MM simulation of the bidentate dianionic La(III)–phosphopeptide complex.

3.2 Monoanionic La(III)–phosphopeptide complex

A bidentate coordination was observed throughout the sampling period and an intramolecular hydrogen bond between the hydrogen of the phosphate group and the backbone (Fig. 6) was formed. Fig. 7 displays the time series of the H_P–O_C distance (left) and the corresponding histogram representation (right), where H_P is the proton at the phosphate group and O_C the oxygen of the backbone's carbonyl group. The histogram displays the H_P–O_C distance distribution, with 1.85 Å being the maximum. The first 20 ps of sampling the H_P–O_C distance fluctuates about this value, indicating the formation of a hydrogen bond. In between 20 ps and 26 ps simulation time, H_P–O_C distances up to 5 Å were observed. During this period hydrogen bonds to water molecules from the bulk were formed, as shown in the left snapshot in Fig. 7. Similar events, although shorter, were observed in between 30 ps and 35 ps simulation time. Subsequently, the intramolecular hydrogen bond was more or less stable until a simulation time of 60 ps. In between 60 ps and 72 ps sampling, hydrogen bonds to bulk water molecules were observed again (right snapshot in Fig. 7). At a simulation time of 72 ps, the intramolecular hydrogen bond was formed and remained stable until the end of the simulation.

Fig. 6 Snapshot taken from the QM/MM MD simulation of the monoanionic La(III)–phosphopeptide complex. Cyan: C; white: H; blue: N; red: O; yellow: P; gray: La(III).

Fig. 7 Time series of H_P–O_C distance (left) and corresponding histogram representation (right) obtained from the QM/MM simulation of the monoanionic La(III)–phosphopeptide.

The La–O_wat RDF of the protonated La(III)–phosphopeptide complex (La–Hpp) in Fig. 3 shows a slightly smaller first peak and the CN_av is marginally lower (Table 1) compared to La–pp_bi. The second peak of the La–O_wat RDF is shifted to larger distances, compared to La–pp_bi. The reduced number of water molecules at around 6 Å can be attributed to the presence of a valine residue. The partial charge of the La(III) ion is 2.05, which is similar as for the La–pp_bi complex.

The τ value for the coordinating oxygen atoms is 14 ps, which is similar to the τ of free La(III) in aqueous solution (Table 1). Although the 9-fold coordinated structure is stabilized by the small O_phosphate–La–O_phosphate angle as discussed above, much more water exchange processes were observed as for the La–pp_bi complex. One reason therefore could be the reduced distance to the valine residue. Phosphopeptides with other aminoacid residues at this position (e.g. aspartate) could possibly have the adverse effect by forming hydrogen bonds and hence stabilizing the hydration structure.

An asymmetric stretching frequency of 1210 cm⁻¹ was found which is considerably higher as for the dianionic complex. The experimentally obtained asymmetric stretching frequency of O–phosphoserine is 1194 cm⁻¹.¹⁰

The RMSD plot obtained from the simulation of the monoanionic La(III)–phosphopeptide complex is shown in Fig. 8. Significant changes of the positions of the backbone atoms were observed in between 20 ps and 30 ps as well as in between 40 ps and 50 ps simulation time. These changes could result from rotations around the N(ILE₆)–C_α(ILE₆) bond, as the C_C(GLH₅)–N(ILE₆)–C_α(ILE₆)–C_C(ILE₆) dihedral angle (ϕ) changes its value almost 100° at 25 ps and 45 ps sampling time (Fig. 8b). Furthermore, the rearrangement in between 20 ps and 30 ps simulation time could be related to the first break of the intramolecular hydrogen bond, which takes place at the same sampling time (Fig. 7).

Fig. 8 (a) Root-mean-square deviation (RMSD) in Å and (b) time series of the C_C(GLH₅)–N(ILE₆)–C_α(ILE₆)–C_C(ILE₆) dihedral angle (ϕ) obtained from the QM/MM simulation of the monoanionic La(III)–phosphopeptide complex.

4 Conclusion

QM/MM simulations of La(III)–phosphopeptide complexes revealed detailed structural as well as dynamical data. The definition of a QM region in combination with the link-atom approach has proven successful in studying the solvation properties of such systems. Due to the constant increase in computational power, the simulations could be conducted at RI-MP2/triple zeta level. The dianionic phosphopeptide was found to coordinate mono- as well as bidentate to the La(III) ion. The exchange rate of the first shell water molecules is much lower as for free La(III) in aqueous solution, indicating an energetically favorable structure which could possibly explain the high selectivity of La(III) based bioanalytical approaches. The observed intramolecular hydrogen in the simulation of the monoanionic La(III)–phosphopeptide complex could result in a stabilization of the acid state and a lowered pK_a of phosphoserine.

Acknowledgements

This work was supported by the Austrian Ministry of Science BMWF UniInfrastrukturprogramm as part of the Research Focal Point Scientific Computing at the University of Innsbruck.

References

1. G. Manning, G. D. Plowman, T. Hunter and S. Sudarsanam, Trends Biochem. Sci., 2002, 27, 514–520.
2. M. R. Larsen, T. E. Thingholm, O. N. Jensen, P. Roepstorff and T. J. Jørgensen, Mol. Cell. Proteomics, 2005, 4, 873–886.
3. M. Pink, N. Verma, F. Polato, G. K. Bonn, H. A. Baba, A. W. Rettenmeier and S. Schmitz-Spanke, J. Proteomics, 2011, 75, 375–383.
4. Y. Güzel, M. Rainer, M. R. Mirza and G. K. Bonn, Anal. Bioanal. Chem., 2012, 403, 1323–1331.
5. Y. Güzel, M. Rainer, M. R. Mirza, C. B. Messner and G. K. Bonn, Analyst, 2013, 138, 2897–2905.
6. M. Pink, C. Stein, N. Verma, A. W. Rettenmeier and S. Schmitz-Spanke, Electrophoresis, 2013, 34, 501–504.
7. M. R. Mirza, M. Rainer, Y. Güzel, I. M. Choudhary and G. K. Bonn, Anal. Bioanal. Chem., 2012, 404, 853–862.
8. M. R. Mirza, M. Rainer, C. B. Messner, Y. Guezel, D. Schemeth, T. Stasyk, M. I. Choudhary, L. A. Huber, B. M. Rode and G. K. Bonn, Analyst, 2013, 138, 2995–3004.
9. A. Saeed, M. Najam-ul Haq, F. Jabeen and F. Svec, Anal. Chem., 2013, 85, 8979–8986.
10. J. M. Sanchez-Ruiz and M. Martinez-Carrion, Biochemistry, 1988, 27, 3338–3342.
11. C. Fernández, S. F. Ausar, R. G. Badini, L. F. Castagna, I. D. Bianco and D. M. Beltramo, Int. Dairy J., 2003, 13, 897–901.
12. A. Tholey, A. Lindemann, V. Kinzel and J. Reed, Biophys. J., 1999, 76, 76–87.
13. J.-T. Du, Y.-M. Li, W. Wei, G.-S. Wu, Y.-F. Zhao, K. Kanazawa, T. Nemoto and H. Nakanishi, J. Am. Chem. Soc., 2005, 127, 16350–16351.
14. S. E. Wong, K. Bernacki and M. Jacobson, J. Phys. Chem. B, 2005, 109, 5249–5258.
15. I. Persson, P. D'Angelo, S. De Panfilis, M. Sandström and L. Eriksson, Chem. – Eur. J., 2008, 14, 3056–3066.
16. O. Lutz, T. S. Hofer, B. R. Randolf and B. M. Rode, Chem. Phys. Lett., 2012, 536, 50–54.
17. T. S. Hofer, B. R. Randolf and B. M. Rode, Phys. Chem. Chem. Phys., 2005, 7, 1382–1387.
18. T. S. Hofer, A. K. H. Weiss, B. R. Randolf and B. M. Rode, Chem. Phys. Lett., 2011, 512, 139–145.
19. M. J. Field, P. A. Bash and M. Karplus, J. Comput. Chem., 1990, 11, 700–733.
20. D. Bakowies and W. Thiel, J. Phys. Chem., 1996, 100, 10580–10594.
21. The Nobel Prize in Chemistry 2013 - Press Release http://www.nobelprize.org/nobel_prizes/chemistry/laureates/2013/press.html, accessed July 2014.
22. N. Reuter, A. Dejaegere, B. Maigret and M. Karplus, J. Phys. Chem. A, 2000, 104, 1720–1735.
23. H. Lin and D. G. Truhlar, Theor. Chem. Acc., 2007, 117, 185–199.
24. D. A. Case, T. A. Darden, T. E. Cheatham, C. L. Simmerling, J. Wang, R. E. Duke, R. Luo, R. Walker, W. Zhang, K. M. Merz, B. Roberts, S. Hayik, A. Roitberg, G. Seabra, J. Swails, A. W. Goetz, I. Kolossvry, K. F. Wong, F. Paesani, J. Vanicek, R. M. Wolf, J. Liu, X. Wu, W. R. Brozell, T. Steinbrecher, H. Gohlke, Q. Cai, X. Ye, J. Wang, M.-J. Hsieh, G. Cui, D. R. Roe, D. H. Mathews, M. G. Seetin, R. Salomon-Ferrer, C. Sagui, V. Babin, T. Luchko, S. Gusarov, A. Kovalenko and P. A. Kollman, University of California, San Francisco, 2012.
25. H. J. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola and J. Haak, J. Chem. Phys., 1984, 81, 3684.
26. H. C. Andersen, J. Comput. Phys., 1983, 52, 24–34.
27. D. J. Adams, E. M. Adams and G. J. Hills, Mol. Phys., 1979, 38, 387–400.
28. V. Hornak, R. Abel, A. Okur, B. Strockbine, A. Roitberg and C. Simmerling, Proteins: Struct., Funct., Bioinf., 2006, 65, 712–725.
29. N. Homeyer, A. H. Horn, H. Lanig and H. Sticht, J. Mol. Model., 2006, 12, 281–289.
30. M. Baaden, M. Burgard, C. Boehme and G. Wipff, Phys. Chem. Chem. Phys., 2001, 3, 1317–1325.
31. F. C. van Veggel and D. N. Reinhoudt, Chem. – Eur. J., 1999, 5, 90–95.
32. H. Berendsen, J. Grigera and T. Straatsma, J. Phys. Chem., 1987, 91, 6269–6271.
33. B. M. Rode, T. S. Hofer, B. R. Randolf, C. F. Schwenk, D. Xenides and V. Vchirawongkwin, Theor. Chem. Acc., 2006, 115, 77–85.
34. R. Ahlrichs, M. Bär, M. Häser, H. Horn and C. Kölmel, Chem. Phys. Lett., 1989, 162, 165–169.
35. F. Weigend, M. Häser, H. Patzelt and R. Ahlrichs, Chem. Phys. Lett., 1998, 294, 143–152.
36. A. Schäfer, C. Huber and R. Ahlrichs, J. Chem. Phys., 1994, 100, 5829.
37. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
38. M. Dolg, H. Stoll and H. Preuss, Theor. Chim. Acta, 1993, 85, 441–450.
39. M. Dolg, H. Stoll, A. Savin and H. Preuss, Theor. Chim. Acta, 1989, 75, 173–194.
40. B. M. Rode and T. S. Hofer, et al. , Pure Appl. Chem., 2006, 78, 525–539.
41. T. S. Hofer, A. B. Pribil and B. R. Randolf, et al. , Adv. Inorg. Chem., 2010, 62, 143–175.
42. A. R. Leach, Molecular modelling: principles and applications, Pearson Education, 2001.
43. K. E. Atkinson, An introduction to numerical analysis, John Wiley & Sons, 2008.
44. T. S. Hofer, H. T. Tran, C. F. Schwenk and B. M. Rode, J. Comput. Chem., 2004, 25, 211–217.
45. A. J. Lock, S. Woutersen and H. J. Bakker, Vibrational relaxation and energy equilibration in liquid water, in Femtochemistry and Femtobiology, ed. A. Douhal and J. Santamaria, World Scientific, Singapore, 2002, pp. 234–239.

---