Explicit Solvent Calculations Research Articles

Tuning Potential Functions to Host-Guest Binding Data.

Software to more rapidly and accurately predict protein-ligand binding affinities is of high interest for early-stage drug discovery, and physics-based methods are among the most widely used technologies for this purpose. The accuracy of these methods depends critically on the accuracy of the potential functions that they use. Potential functions are typically trained against a combination of quantum chemical and experimental data. However, although binding affinities are among the most important quantities to predict, experimental binding affinities have not to date been integrated into the experimental data set used to train potential functions. In recent years, the use of host-guest complexes as simple and tractable models of binding thermodynamics has gained popularity due to their small size and simplicity, relative to protein-ligand systems. Host-guest complexes can also avoid ambiguities that arise in protein-ligand systems such as uncertain protonation states. Thus, experimental host-guest binding data are an appealing additional data type to integrate into the experimental data set used to optimize potential functions. Here, we report the extension of the Open Force Field Evaluator framework to enable the systematic calculation of host-guest binding free energies and their gradients with respect to force field parameters, coupled with the curation of 126 host-guest complexes with available experimental binding free energies. As an initial application of this novel infrastructure, we optimized generalized Born (GB) cavity radii for the OBC2 GB implicit solvent model against experimental data for 36 host-guest systems. This refitting led to a dramatic improvement in accuracy for both the training set and a separate test set with 90 additional host-guest systems. The optimized radii also showed encouraging transferability from host-guest systems to 59 protein-ligand systems. However, the new radii are significantly smaller than the baseline radii and lead to excessively favorable hydration free energies (HFEs). Thus, users of the OBC2 GB model currently may choose between GB cavity radii that yield more accurate binding affinities and GB cavity radii that yield more accurate HFEs. We suspect that achieving good accuracy on both will require more far-reaching adjustments to the GB model. We note that binding free-energy calculations using the OBC2 model in OpenMM gain about a 10× speedup relative to corresponding explicit solvent calculations, suggesting a future role for implicit solvent absolute binding free-energy (ABFE) calculations in virtual compound screening. This study proves the principle of using host-guest systems to train potential functions that are transferrable to protein-ligand systems and provides an infrastructure that enables a range of applications.

Identifying Systematic Force Field Errors Using a 3D-RISM Element Counting Correction.

Hydration free energies of small molecules are commonly used as benchmarks for solvation models. However, errors in predicting hydration free energies are partially due to the force fields used and not just the solvation model. To address this, we have used the 3D reference interaction site model (3D-RISM) of molecular solvation and existing benchmark explicit solvent calculations with a simple element count correction (ECC) to identify problems with the non-bond parameters in the general AMBER force field (GAFF). 3D-RISM was used to calculate hydration free energies of all 642 molecules in the FreeSolv database, and a partial molar volume correction (PMVC), ECC, and their combination (PMVECC) were applied to the results. The PMVECC produced a mean unsigned error of 1.01±0.04kcal/mol and root mean squared error of 1.44±0.07kcal/mol, better than the benchmark explicit solvent calculations from FreeSolv, and required less than 15 s of computing time per molecule on a single CPU core. Importantly, parameters for PMVECC showed systematic errors for molecules containing Cl, Br, I, and P. Applying ECC to the explicit solvent hydration free energies found the same systematic errors. The results strongly suggest that some small adjustments to the Lennard-Jones parameters for GAFF will lead to improved hydration free energy calculations for all solvent models.

Explicit and Hybrid Solvent Models for Estimates of Parameters Relevant to the Reduction Potential of Ethylene Carbonate.

Using ethylene carbonate as a sample solvent, we investigated two molecular parameters used to estimate the reduction potential of the solvent: electron affinity, and the energy of the lowest unoccupied molecular orbital (LUMO). The results showed that the values of these parameters are inconsistent for a single ethylene carbonate molecule in vacuum calculations and in the continuous effective solvent. We performed a series of calculations employing explicit or hybrid (explicit/continuous) solvent models for aggregates of solvent molecules or solvated salt ions. In the hybrid solvent model, values of the two estimates extrapolated to an infinite system size converged to one common value, whereas the difference of 1 eV was calculated in the purely explicit solvent. The values of the gap between the highest occupied molecular orbital (HOMO) and the LUMO obtained in the hybrid model were significantly larger than those resulting from the explicit solvent calculations. We related these differences to the differences in frontier orbitals and changes of electron density obtained in the two solvent models. In the hybrid solvent model, the location of the additional electron in the reduced system usually corresponds to the LUMO orbital of the oxidized system. The presence of salt ions in the solvent affects the extrapolated values of the electron affinity and LUMO energy.

Solution Structural Studies of Pre-amyloid Oligomer States of the Biofilm Protein Aap

Staphylococcus epidermidis is a commensal bacterium on human skin that is also the leading cause of medical device-related infections. The accumulation-associated protein (Aap) from S. epidermidis is a critical factor for infection via its ability to mediate biofilm formation. The B-repeat superdomain of Aap is composed of 5 to 17 Zn2+-binding B-repeats, which undergo rapid, reversible assembly to form dimer and tetramer species. The tetramer can then undergo a conformational change and nucleate highly stable functional amyloid fibrils. In this study, multiple techniques including analytical ultracentrifugation (AUC) and small-angle X-ray scattering (SAXS) are used to probe a panel of B-repeat mutant constructs that assemble to distinct oligomeric states to define the structural characteristics of B-repeat dimer and tetramer species. The B-repeat region from Aap forms an extremely elongated conformation that presents several challenges for standard SAXS analyses. Specialized approaches, such as cross-sectional analyses, allowed for in-depth interpretation of data, while explicit-solvent calculations via WAXSiS allowed for accurate evaluation of atomistic models. The resulting models suggest mechanisms by which Aap functional amyloid fibrils form, illuminating an important contributing factor to recurrent staphylococcal infections.

Ursolic acid as a potential inhibitor of Mycobacterium tuberculosis cytochrome bc1 oxidase-a molecular modelling perspective.

The escalating burden of tuberculosis disease and drastic effects of current medicine has stimulated a search for alternative drugs. A medicinal plant Warburgia salutaris has been reported to possess inhibitory properties against M. tuberculosis. In this study, we apply computational methods to investigate the probability of W. salutaris compounds as potential inhibitors of M. tuberculosis QcrB protein. We performed molecular docking, molecular dynamics simulations, radius of gyration, principal component analysis (PCA), and molecular mechanics-generalized born surface area (MM-GBSA) binding-free energy calculations in explicit solvent to achieve our objective. The results suggested that ursolic acid (UA) and ursolic acid acetate (UAA) could serve as preferred potential inhibitors of mycobacterial QcrB compared to lansoprazole sulphide (LSPZ) and telacebec (Q203)-UA and UAA have a higher binding affinity to QcrB compared to LSPZ and Q203 drugs. UA binding affinity is attributed to hydrogen bond formation with Val120, Arg364 and Arg366, and largely resonated from van der Waals forces resulting from UA interactions with hydrophobic amino acids in its vicinity. UAA binds to the porphyrin ring binding site with higher binding affinity compared to LSPZ. The binding affinity results primarily from van der Waals forces between UAA and hydrophobic residues of QcrB in the porphyrin ring binding site where UAA binds competitively. UA and UAA formed stable complexes with the protein with reduced overall residue mobility, consequently supporting the magnitude of binding affinity of the respective ligands. UAA could potentially compete with the porphyrin ring for the binding site and deprive the mycobacterial cell from oxygen, consequently disturbing mycobacterial oxygen-dependent metabolic processes. Therefore, discovery of a compound that competes with porphyrin ring for the binding site may be useful in QcrB pharmocological studies. UA proved to be a superior compound, although its estimated toxicity profile revealed UA to be hepatotoxic within acceptable parameters. Although preliminary findings of this report still warrant experimental validation, they could serve as a baseline for the development of new anti-tubercular drugs from natural resources that target QcrB.

Quantum Mechanical Modeling of Reaction Rate Acceleration in Microdroplets.

Organic reactions in microdroplets can be orders of magnitude faster than their bulk counterparts. We hypothesize that solvation energy differences between bulk and interface play a key role in the intrinsic rate constant increase and test the hypothesis with explicit solvent calculations. We demonstrate for both the protonated phenylhydrazine reagent and the hydrazone transition state (TSB) that molecular orientations which place the charge sites at the surface confer high energy. A pathway in which this high-energy form transforms into a fully solvated TSB has a lower activation energy than bulk by some 59 kJ/mol, a result that is consistent with experimental rate acceleration studies.

Combined Small-Angle X-ray and Neutron Scattering Restraints in Molecular Dynamics Simulations.

Small-angle X-ray and small-angle neutron scattering (SAXS/SANS) provide unique structural information on biomolecules and their complexes in solution. SANS may provide multiple independent data sets by means of contrast variation experiments, that is, by measuring at different D2O concentrations and different perdeuteration conditions of the biomolecular complex. However, even the combined data from multiple SAXS/SANS sets is by far insufficient to define all degrees of freedom of a complex, leading to a significant risk of overfitting when refining biomolecular structures against SAXS/SANS data. Hence, to control against overfitting, the low-information SAXS/SANS data must be complemented by accurate physical models, and, if possible, refined models should be cross-validated against independent data not used during the refinement. We present a method for refining atomic biomolecular structures against multiple sets of SAXS and SANS data using all-atom molecular dynamics simulations. Using the protein citrate synthase and the protein/RNA complex Sxl-Unr-msl2 mRNA as test cases, we demonstrate how multiple SAXS and SANS sets may be used for refinement and cross-validation, thereby excluding overfitting during refinement. For the Sxl-Unr-msl2 complex, we find that perdeuteration of the Unr domain leads to a unique, slightly compacted conformation, whereas other perdeuteration conditions lead to similar solution conformations compared to the nondeuterated state. In line with our previous method for predicting SAXS curves, SANS curves were predicted with explicit-solvent calculations, taking atomic models for both the hydration layer and the excluded solvent into account, thereby avoiding the use of solvent-related fitting parameters and solvent-reduced neutron scattering lengths. We expect the method to be useful for deriving and validating solution structures of biomolecules and soft-matter complexes, and for critically assessing whether multiple SAXS and SANS sets are mutually compatible.

Comparison of Implicit and Explicit Solvent Models for the Calculation of Solvation Free Energy in Organic Solvents.

2015, 55, 1192-1201) and experimental data. The comparison indicates that both PB and GB give poor agreement with explicit solvent calculations and even worse agreement with experiments (root-mean-square deviation ≈ 15 kJ/mol). The main problem seems to be the prediction of the apolar contribution, which should include the solvent entropy. The quantum mechanical-based SMD model gives significantly better agreement with experimental data than do PB or GB, but it is not as good as explicit solvent calculation results. The dielectric constant ε of the solvent is found to be a powerful predictor for the polar contribution to the free energy in implicit models; however, the Onsager relation may not hold for realistic solvent, as suggested by explicit solvent and SMD calculations. From the comparison, we also find that with an optimization of the apolar contribution, the PB model gives slightly better agreement with experiments than the SMD model, whereas the correlation between the optimized GB models and experiments remains poor. Further optimization of the apolar contribution is needed for GB models to be able to treat solvents other than water.

Accurate Atom-by-Atom Predictions of Solvation Electrostatics Using a Hydration-Shell Poisson-Boltzmann Model

Traditional Poisson-Boltzmann models of molecular electrostatics fail to model the behavior of the hydration-shell solvent, which deviates significantly from bulk. Here we present a hydration-shell Poisson-Boltzmann (HSPB) model that accounts for these deviations by perturbing the displacement-field boundary condition. The HSPB model correctly predicts charging free energies for Born ions, amino acids in multiple protonation states, and challenging test cases for charge-hydration asymmetry. Remarkably, the model achieves this accuracy using standard Lennard-Jones radii from molecular-dynamics force fields, and thus it has an order of magnitude fewer fitting parameters than traditional PB theories. Free-energy calculations in explicit solvent establish that HSPB also accurately predicts the charging free energies of individual atoms in biological molecules.In contrast, calculations using the traditional Poisson--Boltzmann model indicate substantial errors in intramolecular screening. The HSPB model's accuracy in predicting free-energy changes recommend it for drug design and protein engineering applications, where electrostatic contributions to binding are often assessed relative to hydrophobic isosteres.

NMR Restrained Protein Structure Calculations in Implicit Water/Membrane Environments

The aqueous and lipid membrane environments that make up the intra- and extra-cellular compartments of biological organisms are critical for supporting the structural and functional integrity of both soluble proteins and membrane proteins. NMR spectroscopy enables characterization of the structures of both types of proteins in their respective, near-native environments. However, most NMR structures are not calculated in physically realistic environments and de novo structure calculations in explicit solvent or explicit lipids are computationally expensive. Traditionally, NMR structures are calculated with a simplified repulsive term to prevent atom clashing. Such treatment accelerates the calculation but sacrifices non-bonded interactions. To facilitate NMR structure calculations in natural environments we are developing a computationally practical implicit solvent and membrane potential for Xplor-NIH. Here we show that this potential affords significant improvements in the conformational quality, accuracy and precision of the calculated structures. Further, it provides correct embedding of membrane proteins in lipid bilayers, as well as physically meaningful insights about residue-residue and protein-membrane interactions. We also describe extensions of the potential to bicelle and micelle sample environments as well as various membrane systems.This research was supported by grants from the National Institutes of Health (GM 110658; GM 100265; EB002031; CA 179087; CA 030199).

Quantifying the Energy Landscape of Ribosome Function

The vast amount of structural and kinetic data available for the ribosome makes it the ideal system for investigating the physical-chemical properties that enable biological dynamics. Through the use of simplified models and explicit-solvent simulations, we are probing the balance between structural flexibility and energetics during large-scale conformational transitions. Explicit-solvent simulations (200ns - 1.3 microseconds) are used to obtain diffusion coefficients, which help provide a quantitative relationship between tRNA movement, large-scale rotary motions and the energetic barriers encountered during function. To complement explicit-solvent calculations, all-atom simulations that employ simplified energetic descriptions allow us to determine the relative contributions of configurational entropy and sterics. Together, these calculations provide a theoretical foundation that quantitatively bridges experimental kinetics, single-molecule measurements, structural/mutational data and theoretical calculations.

Generalized Born and Explicit Solvent Models for Free Energy Calculations in Organic Solvents: Cyclodextrin Dimerization.

Evaluation of solvation (binding) free energies with implicit solvent models in different dielectric environments for biological simulations as well as high throughput ligand screening remain challenging endeavors. In order to address how well implicit solvent models approximate explicit ones we examined four generalized Born models (GB(Still), GB(HCT), GB(OBC)I, and GB(OBC)II) for determining the dimerization free energy (ΔG(0)) of β-cyclodextrin monomers in 17 implicit solvents with dielectric constants (D) ranging from 5 to 80 and compared the results to previous free energy calculations with explicit solvents ( Zhang et al. J. Phys. Chem. B 2012 , 116 , 12684 - 12693 ). The comparison indicates that neglecting the environmental dependence of Born radii appears acceptable for such calculations involving cyclodextrin and that the GB(Still) and GB(OBC)I models yield a reasonable estimation of ΔG(0), although the details of binding are quite different from explicit solvents. Large discrepancies between implicit and explicit solvent models occur in high-dielectric media with strong hydrogen bond (HB) interruption properties. ΔG(0) with the GB models is shown to correlate strongly to 2(D-1)/(2D+1) (R(2) ∼ 0.90) in line with the Onsager reaction field ( Onsager J. Am. Chem. Soc. 1936 , 58 , 1486 - 1493 ) but to be very sensitive to D (D < 10) as well. Both high-dielectric environments where hydrogen bonds are of interest and low-dielectric media such as protein binding pockets and membrane interiors therefore need to be considered with caution in GB-based calculations. Finally, a literature analysis of Gibbs energy of solvation of small molecules in organic liquids shows that the Onsager relation does not hold for real molecules since the correlation between ΔG(0) and 2(D-1)/(2D+1) is low for most solutes. Interestingly, explicit solvent calculations of the solvation free energy ( Zhang et al. J. Chem. Inf. Model . 2015 , 55 , 1192 - 1201 ) reproduce the weak experimental correlations with 2(D-1)/(2D+1) very well.

NMR-Restrained Structure Calculations of Membrane Proteins in Implicit Lipid Bilayer Membranes

The biological functions and molecular structures of proteins are highly dependent on the physical and chemical properties of the surrounding environment. Just as water is essential for supporting the native states of soluble proteins, the lipid bilayer is critical for preserving the functional and structural integrity of membrane proteins.The principal advantage of NMR spectroscopy as a method for structure determination is its ability to examine proteins in samples that are very close to their functional environments. NMR is well suited for studying soluble proteins in water and membrane proteins in lipid or detergent environments.However, NMR structure calculations in explicit solvent or explicit lipids are computationally expensive and many NMR structures are calculated with a simplified repulsive term to prevent atom clashing. This accelerates the calculation but sacrifices analysis of the non-bonded interactions that play important roles in structure and function. To facilitate NMR structure calculations in a physically realistic environment we are developing a computationally efficient implicit solvent and implicit membrane potential for the Xplor-NIH, NMR structure refinement package.Here we show that the potential provides significant improvements both in the quality and precision of the calculated structures, provides correct embedding of membrane proteins in lipids, and provides physically meaningful views of residue-residue interactions and residue-membrane interactions.

Simulating Protein and Nucleic Acid Dynamics on the Microsecond to Millisecond Timescale

Molecular dynamics simulation is an invaluable tool for studying biomolecular structure, dynamics and function. Currently this approach suffers from high computational cost, restricting direct simulation of dynamics to the microsecond timescale. Implicit treatment of solvent can offers orders of magnitude increases in sampling, due to both reduced viscosity and high computational efficiency on consumer graphics processors (GPUs). In the past, these fast implicit solvent models typically had significant weaknesses, and thus the advantages in speed were offset by poor accuracy. We present the application of our newly developed implicit solvent model and physics-based force field, with application to unraveling puzzles in protein and nucleic acid folding. Specifically, we show that within weeks of simulation in implicit solvent, we were able to fold a benchmark series of protein with a variety of sizes and topologies. Comparable calculations in explicit solvent would take years and are largely intractable. Our solvent model is also extended to highly charged (and more challenging) systems such as DNA and RNA.

Small molecule solvation changes due to the presence of salt are governed by the cost of solvent cavity formation and dispersion.

We are interested in the free energies of transferring nonpolar solutes into aqueous NaCl solutions with salt concentrations upwards of 2 M, the Hofmeister regime. We use the semi-explicit assembly (SEA) computational model to represent these electrolyte solutions. We find good agreement with experiments (Setschenow coefficients) on 43 nonpolar and polar solutes and with TIP3P explicit-solvent simulations. Besides being much faster than explicit solvent calculations, SEA is more accurate than the PB models we tested, successfully capturing even subtle salt effects in both the polar and nonpolar components of solvation. We find that the salt effects are mainly due to changes in the cost of forming nonpolar cavities in aqueous NaCl solutions, and not mainly due to solute-ion electrostatic interactions.

Solvatochromic probe in molecular solvents: implicit versus explicit solvent model

Solvent-induced shifts in the absorption spectrum of N,N-diethyl-4-nitroaniline were studied by quantum-chemical methods in water, dimethylsulfoxide, acetonitrile and acetone. TDDFT methodology and semiempirical ZINDO/S and PM6-CIS approaches were used to calculate excitation energies. Solvent effect was modeled in implicit solvent model by different variants of the PCM approach. Classical molecular dynamics was applied to obtain solute–solvent geometries used in explicit solvent modeling. Most implicit solvent models fail to reproduce the sequence of solvatochromic shifts for four studied solvents, usually yielding too small effect for water. The best result of the PCM method was obtained with SMD atomic radii. Semiempirical quantum-chemical methods in explicit solvent model did not provide satisfactory description of solvatochromic shifts with the largest disagreement to experiment observed for water. TDDFT explicit solvent calculations performed the best in modeling of spectral shifts. Problems with reproduction of experimental data were attributed to specific interactions.

Box size effects are negligible for solvation free energies of neutral solutes.

Hydration free energy calculations in explicit solvent have become an integral part of binding free energy calculations and a valuable test of force fields. Most of these simulations follow the conventional norm of keeping edge length of the periodic solvent box larger than twice the Lennard-Jones (LJ) cutoff distance, with the rationale that this should be sufficient to keep the interactions between copies of the solute to a minimum. However, for charged solutes, hydration free energies can exhibit substantial box size-dependence even at typical box sizes. Here, we examine whether similar size-dependence affects hydration of neutral molecules. Thus, we focused on two strongly polar molecules with large dipole moments, where any size-dependence should be most pronounced, and determined how their hydration free energies vary as a function of simulation box size. In addition to testing a variety of simulation box sizes, we also tested two LJ cut-off distances, 0.65 and 1.0nm. We show from these simulations that the calculated hydration free energy is independent of the box-size as well as the LJ cut-off distance, suggesting that typical hydration free energy calculations of neutral compounds indeed need not be particularly concerned with finite-size effects as long as standard good practices are followed.

Long wavelength absorbing carbostyrils as test cases for different TDDFT procedures and solvent models.

The vertical excitation energies of 3,4-dicyano-6-methoxy and 3,4-dicyano-6,7-dimethoxy carbostyril have been computed with different approximations for the time-dependent density functional theory (TD-DFT) procedure and with different implementations of the continuum solvation model COSMO. Different DFT functionals were tested in TD-DFT and Tamm-Dancoff approximations (TDA) for the excitation energies in the gas phase. TDA-B3LYP showed the best agreement with the experimental data. Then TDA-B3LYP computations were performed combined with the COSMO model of solvation comparing a linear response (LR) and a post-configuration interaction (CI) implementation of the fast solvent reorganization. The post-CI solvent model overestimates the π→π* transitions and strongly underestimates the n→π* transition. The TDA approximation in combination with the linear response implementation of the COSMO solvation model perfectly computes the experimental results. TDA-LR is the most reliable method for the computation of the vertical excitation energies in a solvent. Comparison with explicit solvent calculations shows there is only a minor effect on the energies of the electronic interaction of the solute with the solvent.

Reduced Model Captures Mg2+-RNA Interaction Free Energy of Riboswitches

The stability of RNA tertiary structures depends heavily on Mg2+. The Mg2+-RNA interaction free energy that stabilizes an RNA structure can be computed experimentally through fluorescence-based assays that measure Γ2+, the number of excess Mg2+ associated with an RNA molecule. Previous explicit-solvent simulations predict that the majority of excess Mg2+ ions interact closely and strongly with the RNA, unlike monovalent ions such as K+, suggesting that an explicit treatment of Mg2+ is important for capturing RNA dynamics. Here we present a reduced model that accurately reproduces the thermodynamics of Mg2+-RNA interactions. This model is able to characterize long-timescale RNA dynamics coupled to Mg2+ through the explicit representation of Mg2+ ions. KCl is described by Debye-Hückel screening and a Manning condensation parameter, which represents condensed K+ and models its competition with condensed Mg2+. The model contains one fitted parameter, the number of condensed K+ ions in the absence of Mg2+. Values of Γ2+ computed from molecular dynamics simulations using the model show excellent agreement with both experimental data on the adenine riboswitch and previous explicit-solvent simulations of the SAM-I riboswitch. This agreement confirms the thermodynamic accuracy of the model via the direct relation of Γ2+ to the Mg2+-RNA interaction free energy, and provides further support for the predictions from explicit-solvent calculations. This reduced model will be useful for future studies of the interplay between Mg2+ and RNA dynamics.

NMR-Restrained Protein Structure Calculations in Implicit Environment

The aqueous and lipid membrane environments that make up the intra- and extra-cellular compartments of biological organisms are critical for supporting the structural and functional integrity of both soluble proteins and membrane proteins. NMR spectroscopy is extremely versatile and adept at characterizing the structures and functions of both types of proteins in samples that closely resemble the native environment. However, structure calculations in explicit solvent or explicit lipids are computationally expensive. Therefore, even when NMR structural restraints are measured in a native-like protein environment, protein structures are typically calculated using simple repulsive potentials without physically meaningful terms for electrostatics, non-covalent interactions or atom solvation. To facilitate NMR structure calculations in the proper environment we are developing a computationally less demanding implicit solvent potential for the XPLOR-NIH NMR structure refinement package. Here we show that the potential provides significant improvements both in the quality and precision of the calculated structures; it improves the accuracy of structures determined with sparse restraints; it provides numerous physically meaningful hydrogen bonding connections; and it provides correct embedding of membrane proteins in lipid bilayer membranes.This research was supported by grants from the National Institutes of Health (AI074805; GM100265; GM094727; EB002031).