Srinivasa M.
Gopal
,
Alexander B.
Kuhn
and
Lars V.
Schäfer
*
Lehrstuhl für Theoretische Chemie, Ruhr-Universität Bochum, Germany. E-mail: lars.schaefer@rub.de; Fax: +49 234 3214045; Tel: +49 234 3221582
First published on 8th January 2015
In bundled SPC water models, the relative motion of groups of four water molecules is restrained by distance-dependent potentials. Bundled SPC models have been used in hybrid all-atom/coarse-grained (AA/CG) multiscale simulations, since they enable to couple atomistic SPC water with supra-molecular CG water models that effectively represent more than a single water molecule. In the present work, we systematically validated and critically tested bundled SPC water models as solvent for biomolecular simulations. To that aim, we investigated both thermodynamic and structural properties of various biomolecular systems through molecular dynamics (MD) simulations. Potentials of mean force of dimerization of pairs of amino acid side chains as well as hydration free energies of single side chains obtained with bundled SPC and standard (unrestrained) SPC water agree closely with each other and with experimental data. Decomposition of the hydration free energies into enthalpic and entropic contributions reveals that in bundled SPC, this favorable agreement of the free energies is due to a larger degree of error compensation between hydration enthalpy and entropy. The Ramachandran maps of Ala3, Ala5, and Ala7 peptides are similar in bundled and unrestrained SPC, whereas for the (GS)2 peptide, bundled water leads to a slight overpopulation of extended conformations. Analysis of the end-to-end distance autocorrelation times of the Ala5 and (GS)2 peptides shows that sampling in more viscous bundled SPC water is about two times slower. Pronounced differences between the water models were found for the structure of a coiled-coil dimer, which is instable in bundled SPC but not in standard SPC. In addition, the hydration of the active site of the serine protease α-chymotrypsin depends on the water model. Bundled SPC leads to an increased hydration of the active site region, more hydrogen bonds between water and catalytic triad residues, and a significantly slower exchange of water molecules between the active site and the bulk. Our results form a basis for assessing the accuracy that can be expected from bundled SPC water models. At the same time, this study also highlights the importance of evaluating beforehand the effects of water bundling on the biomolecular system of interest for a particular multiscale simulation application.
Multiscale simulations aim to achieve a balance between accuracy and efficiency by combining different levels of resolution. Multiscale methods can be classified into serial (or sequential) and parallel (or hybrid) approaches.6 In the serial approaches, the different models are used one after another: first, information from an atomistic simulation is used to parameterize a CG model. Then, back-mapping methods7 can be used to convert structures obtained from a CG simulation back into the underlying atomistic ensembles. In the parallel or hybrid schemes, by contrast, different resolution models are used for certain components or parts of the system. Hence, atomistic and coarse-grained representations are simultaneously present in the same simulation system, thus requiring direct interactions between them. One way to achieve this is to partition the simulation system into different spatial regimes. Often, the phenomena of interest are rather local, such as, e.g., the binding of a small molecule to a binding site or conformational changes within a (part of a) biological macromolecule. In such cases, AA/CG methods in which an AA model (for the small subsystem of interest) is coupled to a CG description for the remainder (mostly solvent) can be computationally very efficient.
In the adaptive resolution multiscale methods,8–12 the solute of interest is described at the all-atom level and embedded in a (usually spherical) shell of atomistic solvent molecules, which in turn is surrounded by CG solvent. A healing region of a particular width is introduced, in which the molecules gradually switch their resolution on-the-fly upon diffusing into or out of the atomistic zone. In this region, forces9 or potentials12 are scaled to enact smooth switching between the different levels of resolution. Adaptive resolution simulations in which a supra-molecular CG solvent is used require mapping of a group of atomistic solvent molecules to a single CG site located at the center of mass of these atomistic molecules. To map the atomistic coordinates in such a way that a low-energy CG configuration is obtained is impossible if the individual atomistic solvent molecules diffuse independently.13 To enable this mapping in multiscale simulations, atomistic bundled SPC water models14,15 have been developed, in which distance restraints are used to confine the relative motion of water molecules in groups of four. This choice was motivated by the corresponding four-to-one mapping of atomistic and CG water in the widely used CG-Martini force field,16–18 to which the bundled SPC models can thus be coupled. Recent applications have employed bundled SPC water in adaptive resolution simulations of solvents15,19 and a protein in water.20 These hybrid AA/CG simulations provided a computational speed-up of up to a factor 12 as compared to a fully atomistic system of the same size.15 Integrating the equations of motion in the AA and CG subsystems with different time steps would further boost efficiency. Of course, the computational efficiency gain that can be achieved for a particular simulation system depends on the system size and the ratio of atomistic and CG particles, though. The accuracy of such multiscale approaches depends critically on the properties of the bundled water, which is used as the inner shell solvent surrounding the solute of interest. However, thus far, the effects of bundling the water molecules on the structural and thermodynamic properties of the embedded biomolecular solutes have not been explored in much detail.
Here, we systematically evaluate and critically test bundled SPC water as solvent for various biomolecular systems. We did not perform AA/CG multiscale simulations, but focused on fully all-atom simulations in bundled SPC as compared to reference simulations in standard (unrestrained) water. This approach does not provide any computational speed-up, but it allowed us to study the effects of bundling the water molecules without possible additional influences of an AA/CG boundary. We calculated potentials of mean force (PMFs) of dimerization of selected pairs of amino acid side chains as well as free energies of hydration of amino acid side chains. These hydration free energies were decomposed into the enthalpic and entropic contributions. Furthermore, we investigated the conformational sampling of different polypeptides and a dimeric coiled-coil structure. Finally, we explored the effect of bundling the solvent molecules on the hydration of the active site of the protein α-chymotrypsin (α-CT). Our results show that for many of the investigated thermodynamic properties, the bundled SPC models yield results in agreement with unrestrained SPC and experiments. However, differences are observed for the structural sampling of the coiled-coil dimer and the arrangement of water molecules in the α-CT active site.
Model | k dr (kJ mol−1 nm−2) | C 12 (kJ mol−1 nm12) |
---|---|---|
Bundled SPC (MOD1) | 1000 | 3.25000 × 10−6 |
Bundled SPC (MOD2) | 4000 | 3.45000 × 10−6 |
SPC | n/a | 2.63413 × 10−6 |
Amino acid | Abbrev. | Analogue |
---|---|---|
Alanine | Ala (A) | Methane |
Asparagine | Asn (N) | Acetamide |
Cysteine | Cys (C) | Methanethiol |
Glutamine | Gln (Q) | Propionamide |
Leucine | Leu (L) | Isobutane |
Methionine | Met (M) | Methyl ethyl sulfide |
Backbone | NMA | n-Methylacetamide |
Phenylalanine | Phe (F) | Toluene |
Serine | Ser (S) | Methanol |
Threonine | Thr (T) | Ethanol |
Tryptophan | Trp (W) | 3-Methylindole |
Tyrosine | Tyr (Y) | p-Cresol |
Valine | Val (V) | Propane |
The constraint method, as outlined previously,36,37 was used for the PMF calculations. The constraint procedure involves a set of distance constraint simulations from which the PMF is calculated as , where 〈fc〉r is the average force on a constraint between the centers of mass of two solute molecules separated by distance r. The PMF at Rm, the maximum distance used, is set to zero. All systems were solvated with either 2400 SPC (or SPC/E) molecules or 600 bundled SPC water clusters (corresponding to 2400 water molecules) in a periodic rhombic dodecahedron simulation box with a size of 5.0 nm. The distance range (0.25 to 2.2 nm) and spacing (0.05 nm) is identical to the previous work.36 Additional distances were added in certain cases to explore in detail the local features of the PMF profiles. In total, ca. 40 constraint simulations were performed for each system. For each distance, the system was energy minimized (1000 steepest descent steps) and simulated for 11 ns, with constraint forces saved every 100 fs. For the PMF calculations, the first 1 ns of the trajectories was discarded. Statistical errors were estimated from the limiting values of the block averages,38 as implemented in the g_analyze tool of the Gromacs distribution. Our results in unrestrained SPC agree with the dimerization free energies obtained by de Jong and coworkers39 from counting the relative monomer/dimer populations in extended MD simulations.
Vsc(r) = (1 − λ)Vo([ασ6λ + r6](1/6)) |
To determine the hydration enthalpies, ΔHhyd, extended 110 ns MD simulations of the following systems were carried out: (i) the entire solute–solvent system, (ii) only solvent (comprising the same number of solvent molecules as in the entire system), and (iii) only the solute molecule in vacuo. The simulation parameters were identical to the ones described above, with the exception that instead of Langevin temperature coupling, the leap-frog integrator and the velocity rescaling29 thermostat (τT = 0.1 ps) were used. For the simulations of the solute in vacuo, no cut-offs were used for the non-bonded interactions, and overall translation and rotation were removed. Finally, the hydration enthalpy was calculated as ΔHhyd = 〈Usystem〉 − 〈Usolvent〉 − 〈Usolute〉. The pΔV contribution, which for the simulated systems is smaller than the statistical errors,40 was neglected. The entropy was obtained from the difference, TΔShyd = ΔHhyd − ΔGhyd. Statistical errors in ΔHhyd and TΔShyd are below 1.5 and 2.1 kJ mol−1, respectively, in all cases (see ESI†).
The PMF for the Phe–Phe pair is shown in Fig. 1A. The profile in unrestrained SPC (black) shows a strong contact minimum at 0.55 nm with a relative free energy of −3.0 kJ mol−1. A small barrier of 0.3 kJ mol−1 at 0.85 nm is followed by a solvent-separated shallow minimum at 1.0 nm. The results for SPC/E are similar to SPC, with a slightly less deep contact minimum of −2.7 kJ mol−1 (ESI†). Our Phe–Phe PMF is very similar to the benzene–benzene PMF obtained by Villa and coworkers46 with the Gromos 53a6 force field and SPC/E water. The PMF in bundled SPC MOD1 (red) is similar to unrestrained SPC, with a few notable differences. The contact minimum is deeper (−3.9 kJ mol−1), and the solvent-separated minimum is slightly shifted, lower in energy, and significantly broader. In MOD2 (blue), the first minimum is raised (to −3.1 kJ mol−1) compared to MOD1. The desolvation peak (0.5 kJ mol−1) is higher than in SPC, and a second shallow minimum occurs at 1.05 nm. The deviations of the PMF profiles in the bundled water models with respect to SPC are shown as inset in Fig. 1A. The overall root mean square deviation (RMSD) between the PMFs in bundled and unrestrained SPC is 0.5 and 0.1 kJ mol−1 for MOD1 and MOD2, respectively (Table 3). With respect to SPC/E, these RMSDs are significantly higher (3.2 and 2.9 kJ mol−1 for MOD1 and MOD2, respectively).
Fig. 1 PMFs in bundled SPC (red and blue curves for MOD1 and MOD2, respectively) as compared to the reference PMFs in unrestrained SPC (black). Insets show the deviation from the reference. |
System | MOD1 (kJ mol−1) | MOD2 (kJ mol−1) | ||
---|---|---|---|---|
SPC | SPC/E | SPC | SPC/E | |
Ser–Ser | 0.4 | 0.2 | 0.5 | 0.4 |
NMA–NMA | 0.3 | 0.3 | 0.4 | 0.5 |
Phe–Phe | 0.5 | 3.2 | 0.1 | 2.9 |
Val–Val | 0.3 | 1.7 | 0.4 | 1.7 |
Na+–Cl− | 1.2 | 0.9 | 1.9 | 1.2 |
Lys+–Glu− | 1.3 | 0.7 | 1.5 | 0.8 |
The PMF profile for Val–Val is similar to Phe–Phe, as shown in Fig. 1B. The contact minima in all water models are located at 0.5 nm. The dimers are more stable in MOD1 (−3.0 kJ mol−1) and MOD2 (−3.3 kJ mol−1) compared to SPC (−2.5 kJ mol−1) and SPC/E (−2.0 kJ mol−1). The desolvation barrier at 0.7 nm is similar, irrespective of the water model. As observed for Phe–Phe, the solvent-separated minimum is shallower and shifted by 0.15 nm in bundled water compared to SPC and SPC/E. The RMSD of the PMF profiles (inset) shows an oscillating behavior. The RMSD of the PMF profiles is very low, 0.3 and 0.4 kJ mol−1 for MOD1 and MOD2, respectively. The PMFs in SPC/E are provided in ESI.†
The Ser–Ser PMF in Fig. 1C shows that, in addition to a contact minimum at 0.45 nm, an additional sharp minimum exists at a shorter distance (0.33 nm). This peculiar feature differentiates this PMF profile from the others. In unrestrained SPC (black), the free energies of these two minima are −1.1 and −0.3 kJ mol−1, respectively. The relative energies in bundled SPC models MOD1 and MOD2 (+0.7 and +1.3 kJ mol−1, respectively) are higher for the sharp minimum. The minimum at 0.45 nm is similar to SPC for both MOD1 and MOD2. The solvent-separated minimum is shifted from 0.75 nm in SPC to 0.95 nm in bundled SPC, although the desolvation peak remains roughly unchanged.
Fig. 1D shows the PMF for a pair of NMA molecules. The global minimum is at 0.45 nm for all solvent models, with similar free energies of the contact minima (−2.3 kJ mol−1 in SPC compared to −2.7 kJ mol−1 in MOD1 and −2.2 kJ mol−1 in MOD2). The overall RMSDs between the PMFs in unrestrained SPC and bundled SPC are similarly low for both Ser–Ser and NMA–NMA (Table 3).
The PMFs of Na+–Cl− are shown in Fig. 1E. In unrestrained SPC (black), a minimum at 0.25 nm (−11.4 kJ mol−1) is separated by a barrier at 0.35 nm (+4.7 kJ mol−1) from the second minimum at 0.50 nm (−5.3 kJ mol−1). In bundled SPC, the locations of these points are similar, but the relative free energies are different. In particular, the first minimum is less deep (−7.6 and −5.8 kJ mol−1 for MOD1 and MOD2, respectively). In contrast to the side chain analogues, the position of the second minimum is not shifted in bundled SPC, as a result of the small size of the ions. Due to the free energy difference between the contact minima in unrestrained and bundled SPC, the overall RMSD between the PMFs is significant, 1.2 and 1.9 kJ mol−1 for MOD1 and MOD2, respectively (Table 3). In this respect, the bundled SPC models compare more favorably to unrestrained SPC/E, where the free energy of the contact minimum (−8.1 kJ mol−1) is similar. This agreement is most likely due to the more similar static dielectric permittivities of SPC/E (εr = 71) and bundled SPC19 (εr = 77), whereas εr = 66 for SPC. Indeed, the same effect is observed for the charged Lys+–Glu− pair (Fig. 1F), where the relative free energy of the contact minimum is −3.2 kJ mol−1 for unrestrained SPC/E and both bundled SPC models, as compared to −4.3 kJ mol−1 for unrestrained SPC.
In conclusion, the dimerization free energy profiles of the investigated molecules are overall very similar in bundled SPC water as compared to unrestrained SPC. One significant and recurring difference is that the solvent-separated minima are broader in bundled SPC. This may, at least to some extent, be a solvent packing effect, as reflected by the oxygen–oxygen RDF of bundled SPC itself that displays shifted and broadened minima.14 For the dimerization of larger molecules, however, we anticipate some stronger effects of the bundling potentials. For example, upon binding of two large (apolar) molecular surfaces, such as upon protein–protein encounter, bundling may well have an effect on the expulsion of the last hydration layers, since the water molecules have to get excluded from the contact interface in groups of four. Additional scenarios that, by design, cannot be expected to be captured with bundled water are files of single water molecules or interactions mediated by single water molecules, be it at protein–protein interfaces or within proteins. The performance of bundled water in simulations of a helix–helix dimer involving an apolar contact interface as well as the hydration of an enzyme active site will be discussed below.
Fig. 2 Calculated values of ΔGhyd for side chain analogues in unrestrained SPC (black triangles, data taken from ref. 40), MOD1 (red circles) and MOD2 (blue squares) compared to experimental values45 (black circles). |
Solvent | ΔGhyd | ΔHhyd | TΔShyd | |||
---|---|---|---|---|---|---|
Expt | SPC | Expt | SPC | Expt | SPC | |
SPC | 1.6 | — | 2.7 | — | 3.4 | — |
MOD1 | 2.4 | 1.0 | 5.4 | 3.7 | 4.4 | 3.1 |
MOD2 | 3.7 | 2.3 | 5.8 | 4.0 | 4.9 | 3.9 |
As opposed to ΔGhyd, hydration enthalpy and entropy are more sensitive to structural rearrangements of solvent molecules in the hydration shells around the solute. The hydration enthalpy for different solvent models is shown in Fig. 3A and ESI.† The MAE for ΔHhyd is significantly larger than for ΔGhyd (Table 4), suggesting that error compensation led to the more favorable agreement of ΔGhyd with experiment. This observation can be rationalized considering that hydration enthalpy is dominated by water–water interactions. The internal energy per water molecule is −33.8 and −31.6 kJ mol−1 for bundled water MOD1 and MOD2, respectively, which significantly differs from unrestrained SPC (−42.0 kJ mol−1). By design, the internal energy of bundled water has to be higher (less negative) due to the entropy loss as a result of the restraining potentials. Thus, when decomposing the free energy into enthalpic and entropic components, a larger degree of compensation is expected.
Fig. 3 The hydration enthalpy (A) and entropy (B) for side chain analogues in SPC (black triangles, data from ref. 40), MOD1 (red circles) and MOD2 (blue squares). Experimental values47 are shown as black filled circles. |
Hydration enthalpies are either overestimated (too positive), as in case of small hydrophobic residues, or underestimated (too negative), as observed for the polar side chains (Fig. 3A). The too positive hydration enthalpy of hydrophobic side chains in unrestrained SPC was attributed to the hydration volume and thermal expansion coefficient.40 To confirm this, we calculated the isobaric thermal expansion coefficients, αP, of MOD1 and MOD2. Indeed, the obtained values of 8.0 × 10−4 K−1 for MOD1 and 9.0 × 10−4 K−1 for MOD2 are higher than for unrestrained SPC (αP = 7.4 × 10−4 K−1), thus accounting for the observed too positive ΔHhyd. For polar side chains, additional effects play a role, a detailed investigation of which is beyond the scope of the present work.
The hydration entropy is shown in Fig. 3B. Both bundled SPC models either under- or overestimate TΔShyd. The overall MAE with respect to unrestrained SPC is 3.1 and 3.9 kJ mol−1 for MOD1 and MOD2, respectively. Like for the enthalpies, the hydration entropies of polar side chains are underestimated and those of the small hydrophobic side chains (Ala, Val and Leu) overestimated. As a consequence, error compensation leads to ΔGhyd that is close to experiment. One notable exception is Phe, where both ΔHhyd and TΔShyd closely match experimental values for both bundled SPC models.
As observed previously,48 the Gromos force field exhibits a broader minimum in the ppII region and a subdued C5 minimum. Additional features include two similar minima in the α+ basin. In addition, C7ax is not well-populated with the Gromos 53a6 force field. The transition region connecting αR and γ is slightly more prominent in the longer Ala7, along with a more pronounced γ region (around ψ = −75°, ϕ = −90°). The sampling in Ala3 is very similar in the bundled SPC models as depicted in Fig. 4A (middle and right panel). The major minima are well characterized. The only differences occur in the periphery of the dominant ϕ/ψ regions, which however we at least partly attribute to the larger statistical errors in these lower-populated regions at the rims of the wells. In addition, the γ region is stabilized in bundled SPC. A similar picture is obtained for Ala5 (Fig. 4B). For Ala7, by contrast, there are a few noticeable differences between the Ramachandran maps (Fig. 4C). In MOD2, the lower half of the α+ region is considerably more favorable compared to SPC, whereas the γ region is too high in free energy. In MOD1, the αL region is oversampled.
The ϕ/ψ sampling was explored further by comparing the populations of the various regions (Table 5). In spite of the overall similarity, there are differences in the ppII and β regions. The ppII region is the most favored region in all polypeptides, irrespective of the solvent model. In unrestrained SPC, this region accounts for about 33% of the overall population, which is slightly higher than for bundled models MOD1 and MOD2. Furthermore, the higher propensity of ppII compared to β is more evident in SPC, where the difference between the populations is 7.1%, 5.8%, and 7.9% for Ala3, Ala5, and Ala7 respectively. In bundled water, these differences are smaller (4.2%, 2.2%, and 3.5% in MOD1 and 3.1%, 2.8%, and 2.3% in MOD2). The α+ region is consistently more favored in MOD2 compared to the other solvents, with the largest difference observed for Ala7.
Solute | Solvent | % α+ | % αR | % αL | % ppII | % β | # hbonds |
---|---|---|---|---|---|---|---|
Ala3 | SPC | 13.6 | 7.0 | 0.9 | 32.7 | 25.6 | 8.1 |
MOD1 | 14.2 | 7.1 | 0.7 | 30.3 | 26.1 | 8.3 | |
MOD2 | 14.5 | 7.2 | 0.7 | 29.2 | 26.1 | 8.4 | |
Ala5 | SPC | 12.4 | 6.7 | 0.8 | 32.6 | 26.8 | 12.0 |
MOD1 | 11.8 | 6.0 | 1.0 | 30.1 | 27.9 | 12.3 | |
MOD2 | 13.4 | 7.0 | 1.0 | 29.6 | 26.8 | 12.4 | |
Ala7 | SPC | 11.8 | 6.7 | 1.0 | 33.5 | 25.6 | 15.8 |
MOD1 | 12.0 | 6.6 | 1.1 | 29.8 | 26.3 | 16.3 | |
MOD2 | 16.5 | 8.5 | 1.0 | 28.3 | 26.0 | 16.5 |
In addition to the alanine peptides, we investigated the conformational sampling of the (GS)2 peptide, which is frequently used in experiments as a model for flexible polypeptide linkers. Fig. 5 compares the end-to-end distance distributions in the different water models. In unrestrained SPC, the peptide explores a diverse range of compact and extended structures, with preference for the latter. In bundled SPC, these extended structures are even more favored. The populations of the compact structures (both minima at 0.35 and 0.6 nm combined) are 19%, 15%, and 14% for unrestrained SPC, MOD1 and MOD2, respectively.
Fig. 5 Distribution of end-to-end distance of (GS)2 in SPC (black), MOD1 (red), and MOD2 (blue). The end-to-end distance is calculated between the terminal N and C atoms. |
In the compact structures with 0.6 nm end-to-end distance, a single water molecule bridges the charged NH3+ and COO− termini. The bundling potentials make it more difficult for individual water molecules to adopt such a bridging position, hence reducing the population of this compact structure. The similar effect was also observed in an additional simulation of zwitterionic (instead of blocked) Ala7 (results not shown). In previous simulations of (GS)2 employing the CHARMM22 force field,49 the compact conformation was found to be more favorable than the extended conformation, at odds with our present findings obtained with the Gromos 53a6 force field in combination with any of the used water models, including unrestrained SPC. Taken together, our results suggest that the conformational sampling in the investigated peptides is largely determined by the protein force field and less so by the water model.
Fig. 6 The RDFs between peptides and water oxygen atoms. (A) CH3 group of central residue in Ala5. (B) Oγ of (GS)2 Ser residue. |
Fig. 7 Autocorrelation of the end to end distances of Ala5 (solid) and (GS)2 (dashed) peptides in SPC (black), MOD1 (red) and MOD2 (blue). |
Fig. 8 Native contacts of the coiled-coil dimer structure in the course of the MD simulations, using the first model in the NMR ensemble as a reference. Native contacts were calculated according to Best et al.51 The time traces in (A)–(C) show the three independent simulations (colored in cyan, green, and magenta) in SPC, MOD1, and MOD2, respectively. The inset in (A) shows the variation of the 20 structures in the NMR ensemble. (D) Probability distributions, obtained from the final 200 ns of each simulation. |
By contrast, in bundled SPC, large structural deviations were observed (Fig. 8B and C). These deviations involve partial unfolding of the helices and an increase in the solvent-accessible surface area. For both MOD1 and MOD2, a very large degree of unfolding was observed in one out of three simulations, although significant deviations were observed in all cases (Fig. 8D). These results show that the bundling restraints can have a dramatic effect in certain, sensitive cases in which several, structurally diverse states that are comparable in free energy are kinetically accessible during the simulations.
In all solvent models, the overall protein structure is stable on the 100 ns simulation time scale, as evident from the backbone RMSD (Fig. 10A). Stable protein structures in bundled SPC simulations have also been reported previously for lysozyme14 and ubiquitin.20Fig. 10B and C show the distance distributions of the hydrogen bonds between catalytic triad residues. In standard (non-bundled) SPC, the Ser-195–His-57 hydrogen bond is established for 60% of the simulation time (using cut-offs of 0.35 nm for the donor–acceptor distance and 30 degrees for the hydrogen–donor–acceptor angle, respectively). Likewise, the Asp-102–His-57 hydrogen bond is very stable in SPC, with one hydrogen bond present for 81% of the total simulation time, and even two hydrogen bonds (both carboxylic acid oxygens hydrogen-bonding to the His-57 NH) established in 36% of the frames. Both, the Ser-195–His-57 and the Asp-102–His-57 hydrogen bonds are significantly destabilized in bundled SPC. For bundled SPC models MOD1 and MOD2, the Ser-195–His-57 hydrogen bond occupancy drops to 20% and 10%, respectively. For the single and double Asp-102–His-57 hydrogen bonds, the percentages drop to 58%/35% and 21%/16% for MOD1 and MOD2, respectively.
Fig. 10 Comparison of structural parameters of α-CT in unrestrained SPC, MOD1, and MOD2. (A) Backbone RMSD. (B) Ser-195 Oγ–His-57 Nε distance distribution. (C) Asp-102 Cγ–His-57 Nδ distance distribution. The distances in the X-ray crystal structure are shown as dashed lines. The results shown here were obtained with the Gromos 54a7 force field; see ESI† for the corresponding Amber 99SB-ILDN results. |
To investigate the reasons for this differential stability of the hydrogen bonds between catalytic triad residues and to more thoroughly characterize the hydration of the active site, we analyzed the water molecules in the catalytic triad region. Upon setting up our simulations from the crystal structure, all crystallographic water molecules were removed, since it was not feasible to include them in the bundled SPC setup.‡ In the course of the MD simulations, the active site readily became rehydrated. First, we counted the number of water oxygen atoms within 0.6 nm of the center of mass of His-57, which is close to the geometric center of the active site. A water molecule was counted if its oxygen atom was within this cut-off. The analysis was repeated with different cut-offs (between 0.5 and 0.75 nm), with qualitatively similar results (data not shown). For unrestrained SPC, an average of 3.1 water molecules were found close to the catalytic triad, similar to the number of water molecules resolved in the crystal structure (Table 6). By contrast, for bundled SPC, an average of 4.5 and 10.9 water molecules (in MOD1 and MOD2, respectively) populate the active site region. This pronounced difference is also reflected in the number of hydrogen bonds that these water molecules form with the active site residues (Table 6). Due to their larger number, bundled SPC water molecules more successfully compete with intra-residue hydrogen bonds and form more hydrogen bonds with active site residues than standard SPC.
# water | # hbonds | ||||||
---|---|---|---|---|---|---|---|
Ser-195 | Asp-102 | His-57 | |||||
Backbone | Side chain | Back bone | Side chain | Back bone | Side chain | ||
X-ray42 | 4.0 | 1.0 | 1.0 | 1.0 | — | 1.0 | — |
SPC | 3.1 | 0.3 | 1.6 | 0.3 | 1.0 | 1.2 | 0.7 |
MOD1 | 4.5 | 0.2 | 1.8 | 0.0 | 1.8 | 1.2 | 1.9 |
MOD2 | 10.9 | 1.0 | 1.6 | 0.1 | 3.4 | 1.6 | 2.4 |
In addition, to characterize the dynamics of the active site hydration, we analyzed the residence times of the water molecules in the vicinity of the catalytic triad region. This was done by counting the time periods (τi) that every unique water molecule spend within the cut-off sphere, as defined above. These time periods were accumulated, i.e., also included rebinding of water molecules to the active site region. A mean residence time defined was as .53,54 In standard (non-bundled) SPC, more than 4300 unique water molecules (out of the total 8546 water molecules in the simulation box) visited the active site region during the 100 ns simulation, indicating that the active site hydration was well sampled on this time scale. Most of the active site water molecules stayed for less than 0.1 ns before exchanging with bulk waters, with a mean residence time of only 66 ps. The maximum residence time observed for any of the water molecules was 2.5 ns. By contrast, in the bundled SPC simulations, fewer unique water molecules visited the active site region (2609/8536 and 2815/8536 for MOD1 and MOD2, respectively), and the mean residence times were slightly longer (164 and 370 ps for MOD1 and MOD2, respectively). However, in this case, these mean residence times are less meaningful, since—unlike for standard SPC—the distributions are skewed and have a significant contribution of the long-time regime. Most strikingly, the maximum residence time found for a single water molecule belonging to a bundled 4-water cluster was 12 ns for MOD1, and even 87 ns for MOD2. Such long water residence times could overcompensate any speed-up of a hybrid AA/CG set-up, at least if one is interested in properties related to active site hydration. The slow down of water exchange kinetics may result from the fact that, to access or leave the partially buried active site region, the clusters of 4 water molecules have to deform, which is energetically penalized by the bundling distance restraints.
Significant differences between the water models were found for the larger biomolecular systems studied, such as the structure of the dimeric coiled-coil protein and the hydration of the active site of the enzyme α-chymotrypsin. For the coiled-coil dimer, severe structural distortions and loss of native contacts occurred in the bundled SPC simulations, but not in standard (unrestrained) SPC. These results suggest that, for sensitive systems in which the natively folded state is thermodynamically and kinetically only marginally more stable than alternative conformations, subtle solvation effects due to water bundling can have a pronounced influence on the structural ensemble sampled in the simulations. Likewise, we observed differences between the solvent models for the hydration of the active site of α-chymotrypsin, which is located close to the enzyme surface and accessible to solvent. In bundled water, a larger number of water molecules entered the active site region, leading to the formation of more hydrogen bonds between water and catalytic triad residues at the cost of intra-residue hydrogen bonds. Nevertheless, the overall structure of the enzyme was highly stable, irrespective of the solvent model used in the simulations. The α-chymotrypsin structure is less vulnerable to solvent effects, because in contrast to the coiled-coil dimer, it has a very stable fold including several disulphide bonds. In addition, the exchange of water molecules between the active site and the bulk was hindered for bundled water, with a maximum residence time of tens of ns, as opposed to 2 ns for standard SPC. This kinetic trap may impose severe sampling limitations, especially for the hydration dynamics of partially buried regions. In general, for the investigated thermodynamic, kinetic, as well as structural properties, out of the two bundled SPC models studied (MOD1 and MOD2), the model with the larger force constant of the bundling potential (MOD2) yielded slightly inferior results, although the major differences were observed between bundled and non-bundled (standard) water.
In conclusion, our results show that hybrid multiscale simulation approaches involving bundled water as inner shell solvent can, for some biomolecular systems, be comparably accurate as conventional fully atomistic models. At the same time, the shortcomings of the bundled water models become evident, for example in processes where the molecular nature of individual single waters plays a crucial role, such as the solvation of the α-chymotrypsin active site described here, single-file water molecules (e.g., in aquaporins or nanotubes), or interactions mediated by single waters (e.g., in protein–ligand binding). However, using bundled water may also not be recommendable in situations where its limitations are less obvious, as highlighted by our simulations of the coiled-coil dimer. Thus, it remains to be carefully tested for every simulation system beforehand whether it is recommendable to use bundled water as inner shell solvent in a multiscale setup.
Footnotes |
† Electronic supplementary information (ESI) available. See DOI: 10.1039/c4cp04784b |
‡ It would of course be possible – and in fact is usually recommendable – to retain the crystal waters in the simulations with standard (non-bundled) water, but for consistency we decided to remove them in these cases as well. |
This journal is © the Owner Societies 2015 |