Explicit Treatment Of Solvent Research Articles

Implicit Solvents for the Polarizable Atomic Multipole AMOEBA Force Field.

Computational protein design, ab initio protein/RNA folding, and protein-ligand screening can be too computationally demanding for explicit treatment of solvent. For these applications, implicit solvent offers a compelling alternative, which we describe here for the polarizable atomic multipole AMOEBA force field based on three treatments of continuum electrostatics: numerical solutions to the nonlinear and linearized versions of the Poisson-Boltzmann equation (PBE), the domain-decomposition conductor-like screening model (ddCOSMO) approximation to the PBE, and the analytic generalized Kirkwood (GK) approximation. The continuum electrostatics models are combined with a nonpolar estimator based on novel cavitation and dispersion terms. Electrostatic model parameters are numerically optimized using a least-squares style target function based on a library of 103 small-molecule solvation free energy differences. Mean signed errors for the adaptive Poisson-Boltzmann solver (APBS), ddCOSMO, and GK models are 0.05, 0.00, and 0.00 kcal/mol, respectively, while the mean unsigned errors are 0.70, 0.63, and 0.58 kcal/mol, respectively. Validation of the electrostatic response of the resulting implicit solvents, which are available in the Tinker (or Tinker-HP), OpenMM, and Force Field X software packages, is based on comparisons to explicit solvent simulations for a series of proteins and nucleic acids. Overall, the emergence of performative implicit solvent models for polarizable force fields opens the door to their use for folding and design applications.

Host-guest interactions between p-sulfonatocalix[4]arene and p-sulfonatothiacalix[4]arene and group IA, IIA and f-block metal cations: a DFT/SMD study.

The molecular recognition in aqueous solution is extremely important because most biological processes occur in aqueous solution. Water-soluble members of the calix[n]arene family (e.g., p-sulfonato substituted) can serve as model systems for studying the nature and manner of interactions between biological receptors and small ions. The complex formation behavior of water-soluble p-sulfonatocalix[4]arene and thiacalix[4]arene and group IA, IIA and f-block metal cations has been investigated computationally by means of density functional theory computations in the gas phase and in aqueous environment. The calculated Gibbs free energy values of the complex formation reaction of these ligands with the bare metal cations suggest a spontaneous and energy-favorable process for all metal cations in the gas phase and only for Na+, Mg2+, Lu3+ cations in water environment. For one of the studied cations (La3+) a supramolecular approach with explicit solvent treatment has been applied in the study of the effect of metal hydration on the complexation process. The La3+ binding to the p-sulfonatocalix[4]arene host molecule (now in the metal’s second coordination shell) is still exergonic as evidenced by the negative Gibbs free energy values (ΔG1 and ΔG78). The combination of implicit/explicit solvent treatment seems useful in the modeling of the p-sulfonatocalix[4]arene (and thiacalix[4]arene) complexes with metal cations and in the prediction of the thermodynamic parameters of the complex formation reactions.

An Accurate and Efficient, Computational Method for the Hydration Free Energy of Large and Complex Molecules

The hydration free energy (HFE) is an essential factor to characterize various chemical processes in solvent. It is then needed to construct computational methods of the HFE. In the computational point of view, not only accuracy but efficiency is necessary when large and complex molecules such as proteins are analyzed. For example, the explicit treatment of solvent is effective for accuracy, however, the treatment for protein solutions sacrifices computational efficiency. To construct a computational method satisfying both accuracy and efficiency, we decompose the HFE into two parts called the direct and indirect parts: one depends on the solute-solvent interaction, and the other is the remaining determined by the geometry of solute. The former is computed by the all-atom treatment. The latter is computed by a hybrid method of the energy representation (ER) method and the morphometric approach (MA). The ER method provides the same accuracy as standard methods such as the thermodynamic integration and free energy perturbation methods and achieves efficient computation by the end-point calculation. The ER method has computational advantageous, however, its protocols can be reduced by the hybrid method. The MA provides thermodynamic quantities primarily governed by the geometry of solute and achieves rapid computation by the fitting function constructed by geometric measures and its corresponding coefficients. Since the indirect part of the HFE depends on the geometry of solute, the MA can be applied to the indirect part and can provide same accuracy as the ER method when the coefficients of the MA are determined by fitting to values calculated by the ER method. We show the accuracy and efficiency of the hybrid method, and demonstrate performance in the discrimination of the native structure from decoys.

Evaluation of (95)Mo Nuclear Shielding and Chemical Shift of [Mo6X14](2-) Clusters in the Liquid Phase.

[Mo6X14](2-) octahedral molybdenum clusters are the main building blocks of a large range of materials. Although (95)Mo nuclear magnetic resonance was proposed to be a powerful tool to characterize their structural and dynamical properties in solution, these measurements have never been complemented by theoretical studies which can limit their interpretation for complex systems. In this Article, we use quantum chemical calculations to evaluate the (95)Mo chemical shift of three clusters: [Mo6Cl14](2-), [Mo6Br14](2-), and [Mo6I14](2-). In particular, we test various computational parameters influencing the quality of the results: size of the basis set, treatment of relativistic and solvent effects. Furthermore, to provide quantum chemical calculations that are directly comparable with experimental data, we evaluate for the first time the (95)Mo nuclear magnetic shielding of the experimental reference, namely, MoO4(2-) in aqueous solution. This is achieved by combining ab initio molecular dynamics simulations with a periodic approach to evaluate the (95)Mo nuclear shieldings. The results demonstrate that, despite the difficulty to obtain accurate (95)Mo chemical shifts, relative values for a cluster series can be fairly well-reproduced by DFT calculations. We also show that performing an explicit solvent treatment for the reference compound improves by ∼50 ppm the agreement between theory and experiment. Finally, the standard deviation of ∼70 ppm that we calculate for the (95)Mo nuclear shielding of the reference provides an estimation of the accuracy we can achieve for the calculation of the (95)Mo chemical shifts using a static approach. These results demonstrate the growing ability of quantum chemical calculations to complement and interpret complex experimental measurements.

A Reparametrized Implicit Solvent Model for Accurate Computation of Hydration Free Energies

While solvation plays an important role in many chemical and biological processes, detailed solvent behavior is typically of much less interest than the effect of the solvent on the solute. Although explicit solvent treatment provides a very accurate description of the considered system, it is not necessary for many applications and can be replaced by implicit solvent models. A typical application is the computation of hydration free energies for small molecules. However a recent study [1] showed that there is still a gap between the accuracy of implicit models typically used in biomolecular simulations and explicit solvent treatment. To address this issue we have reparametrized an implicit solvent model consisting of a generalized Born term to model polar solvation and an extended nonpolar term that models cavity formation and Lennard-Jones interactions with water. When fitted to experimental hydration free energies of a large set of 499 small molecules, the model is able to reproduce hydration free energies with a similar error than those derived from explicit solvent computation, while still offering some degrees of freedom in the fit parameters. To test the transferability of the resulting model, we have applied it to the prediction of hydration free energies of Glycine-X-Glycine tri-peptides. This model will be part of the SIMONA [2] simulation package.[1] Knight, J. L. & Brooks III, C. L. Surveying implicit solvent models for estimating small molecule absolute hydration free energies. Journal of Computational Chemistry 32, 2909-2923 (2011).[2] Strunk, T. et al. SIMONA 1.0: An efficient and versatile framework for stochastic simulations of molecular and nanoscale systems. Journal of Computational Chemistry (2012).doi:10.1002/jcc.23089

Introducing Charge Hydration Asymmetry in the Realm of Continuum Solvation

Charge hydration asymmetry (CHA) manifests itself in the experimentally observed strong dependence of free energy of ion hydration on the sign of the ion charge. This asymmetry is not consistently accounted for by popular models of solvation; its magnitude varies greatly between the models. While it is clear that CHA is somehow related to charge distribution within a water molecule, the exact nature of this relationship is unknown. We propose a simple, yet general and rigorous criterion that relates rotational and charge inversion properties of a water molecule's charge distribution with its ability to cause CHA. We show which electric multipole components of a water molecule are key to explain its ability for asymmetric charge hydration. We then test several popular water models and explain why specific models show none, little, or strong CHA in simulations. We use the gained insight to derive an analogue of the Born equation that includes the missing physics necessary to account for CHA and does not rely on redefining the continuum dielectric boundary. The proposed formula is as simple as the original, does not contain any fitting parameters, and predicts hydration free energies and entropies of spherical cations and anions within experimental uncertainty. We further incorporate CHA into the framework of Generalized Born continuum solvation model; our findings suggest that the gap between the practical continuum electrostatics framework and the more fundamental explicit solvent treatment can be reduced considerably by explicitly introducing CHA into the existing framework.

Charge Hydration Asymmetry: The Basic Principle and How to Use It to Test and Improve Water Models

Charge hydration asymmetry (CHA) manifests itself in the experimentally observed strong dependence of free energy of ion hydration on the sign of the ion charge. This asymmetry is not consistently accounted for by popular models of solvation; its magnitude varies greatly between the models. While it is clear that CHA is somehow related to charge distribution within a water molecule, the exact nature of this relationship is unknown. We propose a simple, yet general and rigorous criterion that relates rotational and charge inversion properties of a water molecule's charge distribution with its ability to cause CHA. We show which electric multipole components of a water molecule are key to explain its ability for asymmetric charge hydration. We then test several popular water models and explain why specific models show none, little, or strong CHA in simulations. We use the gained insight to derive an analogue of the Born equation that includes the missing physics necessary to account for CHA and does not rely on redefining the continuum dielectric boundary. The proposed formula is as simple as the original, does not contain any fitting parameters, and predicts hydration free energies and entropies of spherical cations and anions within experimental uncertainty. Our findings suggest that the gap between the practical continuum electrostatics framework and the more fundamental explicit solvent treatment may be reduced considerably by explicitly introducing CHA into the existing continuum framework.

Efficient Computation of the Total Solvation Energy of Small Molecules via the R6 Generalized Born Model.

Efficient and accurate methodologies to compute solvation free energies of small molecules are relevant for many biological and industrial research areas including rational drug design. In this work we test the performance of a recently developed generalized Born method, GB_NSR6 (Aguilar et al. J. Chem. Theory Comput. 2010, 6, 3613-3639) on a common benchmark set of 504 small molecules. The computed solvation energies are compared with those obtained previously by explicit solvent models and experiment. The dominant polar component of the solvation energy is computed by GB_NSR6 with no adjustable parameters, producing a root mean square deviation (RMSD) of 0.89 kcal/mol with respect to explicit solvent (TIP3P). The relatively small nonpolar contribution is estimated using the Gallicchio et al. (J. Comput. Chem. 2005, 25, 479-499) approach. Our results show that GB_NSR6 offers a reasonable balance between efficiency and accuracy: the RMSD from the experiment of computed solvation energies is 1.2 kcal/mol, which is essentially the same as the accuracy of the much more computationally expensive explicit solvent treatment. The average computational time needed to compute the total solvation energy per molecule via GB_NSR6 is only tens of milliseconds on a commodity PC for a typical molecule of about 20 atoms. All of the software developed in this work is freely available from http://people.cs.vt.edu/onufriev/software.php .

Parallel Generalized Born Implicit Solvent Calculations with NAMD.

Accurate electrostatic descriptions of aqueous solvent are critical for simulation studies of bio-molecules, but the computational cost of explicit treatment of solvent is very high. A computationally more feasible alternative is a generalized Born implicit solvent description which models polar solvent as a dielectric continuum. Unfortunately, the attainable simulation speedup does not transfer to the massive parallel computers often employed for simulation of large structures. Longer cutoff distances, spatially heterogenous distribution of atoms and the necessary three-fold iteration over atom-pairs in each timestep combine to challenge efficient parallel performance of generalized Born implicit solvent algorithms. Here we report how NAMD, a parallel molecular dynamics program, meets the challenge through a unique parallelization strategy. NAMD now permits efficient simulation of large systems whose slow conformational motions benefit most from implicit solvent descriptions due to the inherent low viscosity. NAMD's implicit solvent performance is benchmarked and then illustrated in simulating the ratcheting Escherichia coli ribosome involving ~250,000 atoms.

Theoretical Study of O–O Single Bond Formation in the Oxidation of Water by the Ruthenium Blue Dimer

The first key step in the oxidation of water to O(2) by the oxidized species [(bpy)(2)(O)Ru(V)ORu(V)(O)(bpy)(2)](4+) of the Ru blue dimer is studied using density functional theory (DFT) and an explicit solvent treatment. In the model reaction system [L(2)(O)Ru(V)ORu(V)(O)L(2)](4+)·(H(2)O)(4)·W(76), the surrounding water solvent molecules W are described classically while the inner core reaction system is described quantum mechanically using smaller model ligands (L). The reaction path found for the O--O single bond formation involves a proton relay chain: direct participation of two water molecules in two proton transfers to yield the product [L(2)(HOO)Ru(IV)ORu(IV)(OH)L(2)](4+)·(H(2)O)(3)·W(76). The calculated ∼3 kcal/mol reaction free energy and ∼15 kcal/mol activation free energy barrier at 298 K are consistent with experiment. Structural changes and charge flow along the intrinsic reaction coordinate, the solvent's role in the reaction barrier, and their significance for water oxidation catalysis are examined in detail.

Structure of the Alanine Hydration Shell as Probed by NMR Chemical Shifts and Indirect Spin−Spin Coupling

The structure of the alanine hydration shell was modeled by Carr-Parinello molecular dynamics (CPMD) to explain subtle differences in NMR chemical shifts and indirect spin-spin coupling constants of the neutral (zwitterionic), cationic, and anionic forms of this amino acid. In comparison with classical molecular dynamics (MD), the quantum mechanical CPMD approach revealed a more structured solvent and significant differences in the radial and angular distributions of the water molecules around the solute. In particular, the solvent was predicted to be organized around the uncharged COOH and NH(2) residues to a similar degree as that for the charged ones. This was not the case with MD. For snapshot CPMD configurations, the NMR parameters were computed by density functional theory (DFT) and averaged. Obtained values were significantly closer to experimental parameters known for (15)N and (13)C isotopically labeled alanine than those calculated by the conventional implicit dielectric solvent model. The NMR results also quantitatively reflect a superiority of the CPMD over the MD explicit solvent treatment. A further improvement of the computed spin-spin coupling constants could be achieved by taking into account vibrational averaging beyond the harmonic approximation. Differently positioned water molecules in the clusters cause an unexpectedly large scattering of the NMR parameters. About 10-15 dynamics snapshots were required for a satisfactory convergence of the shifts and couplings. The NMR chemical shift was found to be much more sensitive to the molecular hydration than the coupling. The results thus indicate a large potential of the NMR spectroscopy and quantum simulations to probe not only the structure of molecules but also their interactions with the environment.

Cluster Hydration Model for Binding Energy Calculations of Protein−Ligand Complexes

We present a new model for predicting the binding affinity in protein-ligand complexes based on an explicit solvent treatment. The model is referred to as the cluster hydration model. Test calculations were performed on complexes of FK506-binding protein (FKBP) and its six ligands. The calculated binding energies had a good correlation with experimental binding affinities; the correlation coefficient r(2) was 0.93. Moreover, we examined the competition between the residue-ligand interactions and the desolvation energies of residues. The results suggested that, for the stabilization of the complexes in solution, the contributions of nonpolar residues at the binding site are larger than those of charged and polar ones because the electrostatic interaction energies between the residues and the ligand were mostly canceled by the desolvation penalties.

Mapping Electron Paramagnetic Resonance Spin Label Conformations by the Simulated Scaling Method

In order to efficiently simulate spin label behavior when attached to the protein backbone we developed a novel approach that enhances local conformational sampling. The simulated scaling (SS) approach (Li, H., et al. J. Chem. Phys. 2007, 126, 24106) couples the random walk of a potential scaling parameter and molecular dynamics in the framework of hybrid Monte Carlo. This approach allows efficient barrier crossings between conformations. The method retains the thermodynamic detailed balance allowing for determination of relative free energies between various conformations. The accuracy of our method was validated by comparison with the recently resolved X-ray crystal structure of a spin labeled T4 lysozyme in which the spin label was in the interior of the protein. Consistent potentials of mean force (PMF) are obtained for the spin label torsion angles to illustrate their behavior in various protein environments: surface, semiburied, and buried. These PMFs reflect the experimentally observed trends and provide the rationale for the spin label dynamics. We have used this method to compare an implicit and explicit solvent model in spin label modeling. The implicit model, which is computationally faster, was found to be in excellent agreement with the explicit solvent treatment. Based on this collection of results, we believe that the presented approach has great potential in the general strategy of describing the behavior of the spin label using molecular modeling and using this information in the interpretation of EPR measurements in terms of protein conformation and dynamics.

Theoretical Determination of the Standard Reduction Potentials of Pheophytin-ainN,N-Dimethyl Formamide and Membrane

Quantum mechanical/molecular mechanics (QM/MM) calculations were performed on the neutral, anionic, and dianionic forms of Pheophytin-a (Pheo-a) in N,N-dimethyl formamide (DMF) in order to calculate the absolute free energy of reduction of Pheo-a in solution. The geometry of the solvated species was optimized by restricted open-shell density functional treatment (ROB3LYP) using the 6-31G(d) basis set for the molecular species while the primary solvent shell consisting of 45 DMF molecules was treated by the MM method using the universal force field (UFF). Electronic energies of the neutral, anionic, and dianionic species were obtained by carrying out single point density functional theory (DFT) calculations using the 6-311+G(2d,2p) basis set on the respective ONIOM optimized geometries. The CHARMM27 force field was used to account for the dynamical nature of the primary solvation shell of 45 DMF molecules. In the calculations using solvent shells, the required atomic charges for each solvent molecule were obtained from ROB3LYP/6-31G(d) calculation on a single isolated DMF molecule. Randomly sampled configurations generated by the Monte Carlo (MC) technique were used to determine the contribution of the primary shell to the free energy of solvation of the three species. The dynamical nature of the primary shell significantly corrects the free energy of solvation. Frequency calculations at the ROB3LYP/6-31G(d) level were carried out on the optimized geometries of truncated 47-atom models of the neutral and ionic species in vacuum so as to determine the differences in thermal energy and molecular entropy. The Born energy of ion-dielectric interaction, the Onsager energy of dipole-dielectric interaction, and the Debye-Hückel energy of ion-ionic cloud interaction for the pheophytin-solvent aggregate were added as perturbative corrections. The Born interaction also makes a large contribution to the absolute free energy of reduction. An implicit solvation model (DPCM) was also employed for the calculation of standard reduction potentials in DMF. Both the models were successful in reproducing the standard reduction potentials. An explicit solvent treatment(QM/MM/MC + Born + Onsager + Debye corrections) yielded the one electron reduction potential of Pheo-a as -0.92 +/- 0.27 V and the two electron reduction potential as -1.34 +/- 0.25 V at 298.15 K, while the implicit solvent treatment yielded the corresponding values as -1.03 +/- 0.17 and -1.30 +/- 0.17 V, respectively. The calculated values more or less agree with the experimental midpoint potentials of -0.90 and -1.25 V, respectively. Moreover, a numerical finite difference Poisson-Boltzmann solver (FDPB) along with the DPCM methodology was employed to obtain the reduction potential of pheophytin in the thylakoid membrane. The calculated reduction potential value of -0.58 V is in excellent agreement with the reported value -0.61 V.

Oxidation of Hydroxylamine by Nitrous and Nitric Acids. Model Development from First Principle SCRF Calculations

Ab initio molecular orbital calculations have been performed to develop an elementary reaction mechanism for the autocatalytic and scavenging reactions of hydroxylamine in an aqueous nitric acid medium. An improved understanding of the titled reactions is needed to determine the "stability boundary of hydroxylamine" for safe operations of the plutonium-uranium reduction extraction (PUREX) process. Under the operating conditions of the PUREX process, namely, 6 M nitric acid, the reactive forms of hydroxylamine are NH2OH, NH3OH+, and the complex NH3OH.NO3, and those of nitrous acid are NO+, H2ONO+, N2O4, N2O3, NO2, and NO. High-level CBSQB3/IEFPCM and CBSQB3/COSMO calculations were performed using GAUSSIAN03 to investigate the energy landscape and to explore a large number of possible ion-ion, ion-radical, ion-molecule, radical-radical, radical-molecule, and molecule-molecule pathways available to the reactive forms of the reactants in solution. It was found that in solution the autocatalytic generation of nitrous acid proceeds through free radical pathways at low-hydroxylamine concentrations from unprotonated NH2OH via hydrogen abstraction. At high [NH3OH+], we suggest a possible involvement of the NH3ONO+ intermediate via the reaction NH2ONO + NO2 --> HNO + HONO + NO. The NH3ONO+ intermediate, in turn, is formed favorably via the ion-ion reactions of NH3OH+ with NO+ and/or the reaction between NO+ and hydroxylammonium nitrate (HAN). The intermediates involved in the scavenging reaction of nitrous acid by hydroxylamine are NH3ONO+, NH2ONO, NH2(NO)O, NH(NO)OH, and HONNOH and the rate-determining step is the 1,2-NO migration in NH2ONO leading to NH2(NO)O. Reactions NH2ONO --> NH2(NO)O and NH2(NO)O --> NH(NO)OH were studied with two explicit water molecules and the results are discussed in the context of the importance of the explicit treatment of solvent in the determination of the energetics and mechanism of these processes. The rate constants for the reactions were estimated using transition-state theory and other traditional techniques. The kinetic parameters obtained at the B3LYP/CBSB7/IEFPCM level are in reasonable agreement with the limited experimental value. IEFPCM results on free energy of undelocalized polar ions such as NO3-, NO2-, and NH3OH+ are not very accurate and have difficulties in predicting the right direction of acid dissociation equilibrium of HONO2, HONO, and NH3OH+. Explicit incorporation of a solvation shell to these ions improves the theoretical descriptions of acid ionization equilibria as it captures some of the nonlocal effects of these ions. Additional work is needed to correctly describe the solvation shell and to introduce consistency in the theoretical treatment involving explicit solvent molecules. Nevertheless, this systematic exploration of reactions in solution and mechanism development for a solution phase process based on self-consistent reaction field (SCRF) results is likely to be one of the first of its kind.

Conformational studies of the two anomeric methyl glycosides of α-D-Manp-(1→2)-D-Glcp by molecular simulations and NMR 1H,1H T-ROESY experiments

A conformational study has been performed of the two disaccharides α-D-Manp-(1→2)-α-D-Glcp-OMe (1) and α-D-Manp-(1→2)-β-D-Glcp-OMe (2) differing only in the anomeric configuration at the methyl glucoside. Molecular simulations were carried out with five different potential energy functions ranging from pure in vacuo to explicit solvent treatment. For comparison, 1H,1H T-ROESY and heteronuclear trans-glycosidic correlation NMR experiments were performed which resulted in proton–proton distances and 3JC,H values, respectively. In addition, interpretation of 13C NMR glycosylation shifts indicated slightly different average conformations of the two disaccharides. From the molecular simulations selected parameters were calculated for comparison to experiments. The combined analysis revealed the conformational region of the disaccharides and showed that the average ψ torsion angle in 2 was more eclipsed than in 1. Moreover, the translational diffusion coefficients of the disaccharides were determined at 25 °C in D2O by the Longitudinal Eddy-Delayed experiment which resulted in Dt = 0.37 × 10−9 m2 s−1. These values were compared to those obtained from molecular dynamics simulations with explicit water as solvent. The relative differences in translational diffusion between the disaccharides and water were similar for experiment and simulation. However, it is well-known that the TIP3P water model overestimates Dt, so also in these simulations.