Implicit Solvent Methods Research Articles

Predicting pKa in Implicit Solvents: Current Status and Future Directions

Computational prediction of condensed phase acidity is a topic of much interest in the field today. We introduce the methods available for predicting gas phase acidity and pKas in aqueous and non-aqueous solvents including high-level electronic structure methods, empirical linear free energy relationships (LFERs), implicit solvent methods, explicit solvent statistical free energy methods, and hybrid implicit–explicit approaches. The focus of this paper is on implicit solvent methods, and we review recent developments including new electronic structure methods, cluster-continuum schemes for calculating ionic solvation free energies, as well as address issues relating to the choice of proton solvation free energy to use with implicit solvation models, and whether thermodynamic cycles are necessary for the computation of pKas. A comparison of the scope and accuracy of implicit solvent methods with ab initio molecular dynamics free energy methods is also presented. The present status of the theory and future directions are outlined.

Molecular Dynamics Folding Simulation of Amyloid A4 Peptide in Implicit Solvent

The objectives of this study is to investigate the pathway of Amyloid A4 peptide (Pdb id: 1AML) folding and to study the interactions involved during the process. The folding process was performed using molecular dynamics (MD) simulation in implicit solvent method. This protein was simulated up to 350 ns. The trajectories of the folded protein were analyzed based on RMSD, hydrogen bond formation and secondary structure evolution. The results showed that there were a few crucial interactions involved in the folding process and important residues that stabilized the folded peptide were also identified.

Predicting pKafor proteins using COSMO-RS

We have used the COSMO-RS implicit solvation method to calculate the equilibrium constants, pKa, for deprotonation of the acidic residues of the ovomucoid inhibitor protein, OMTKY3. The root mean square error for comparison with experimental data is only 0.5 pH units and the maximum error 0.8 pH units. The results show that the accuracy of pKa prediction using COSMO-RS is as good for large biomolecules as it is for smaller inorganic and organic acids and that the method compares very well to previous pKa predictions of the OMTKY3 protein using Quantum Mechanics/Molecular Mechanics. Our approach works well for systems of about 1000 atoms or less, which makes it useful for small proteins as well as for investigating portions of larger proteins such as active sites in enzymes.

Could MM-GBSA be accurate enough for calculation of absolute protein/ligand binding free energies?

Implicit solvation methods such as MM-GBSA, when applied to evaluating protein/ligand binding free energies, are widely believed to be accurate only for the estimation of relative binding free energies for a congeneric series of ligands. In this work, we show that the MM-GBSA flavor of Prime 3.0, VSGB-2.0, with a variable dielectric model and a novel energy function, could be approaching the accuracy required for evaluating absolute binding free energies, albeit, through a linear regression fit. The data-set used for validation includes 106 protein–ligand complexes that were carefully selected to control for variability in the affinity data as well as error in the modeled complexes. Through systematic analysis, we also quantify the degradation in the R2 of fit between experimental and calculated values with either greater variability in the affinity data or an increase in error in the modeled protein/ligand complexes. Limitations for its application in drug discovery are discussed along with the identification of areas for future development.

Is Uracil Aromatic? The Enthalpies of Hydrogenation in the Gaseous and Crystalline Phases, and in Aqueous Solution, as Tools to Obtain an Answer

The enthalpy of hydrogenation of uracil was derived from the experimental enthalpies of formation, in the gaseous phase, of uracil and 5,6-dihydrouracil, in order to analyze its aromaticity. The enthalpy of formation of 5,6-dihydrouracil was obtained from combustion calorimetry, Knudsen effusion technique and Calvet microcalorimetry results. High-level computational methods were tested for the enthalpy of hydrogenation of uracil, but only with G3 was possible to obtain results in agreement with the experimental ones. It was found that uracil possesses 30.0% of aromatic character in the gaseous phase. Using both implicit, explicit, and hybrid solvation methods, it was possible to obtain a reference value for the enthalpy of hydrogenation of uracil in the aqueous solution and the effect of polarity and hydrogen bonds on the aromaticity of uracil was analyzed. The value of the hydrogenation enthalpy of uracil in aqueous solution was compared with the experimental value in the crystal phase, also dominated by polarity and hydrogen bonds, derived from combustion calorimetry results. The supramolecular effects on the crystal lattice were explored by the computational simulation of π-π staking dimers and hydrogen bonded dimers.

Assessing the Accuracy and Performance of Implicit Solvent Models for Drug Molecules: Conformational Ensemble Approaches

The accuracy and performance of implicit solvent methods for solvation free energy calculations were assessed on a set of 20 neutral drug molecules. Molecular dynamics (MD) provided ensembles of conformations in water and water-saturated octanol. The solvation free energies were calculated by popular implicit solvent models based on quantum mechanical (QM) electronic densities (COSMO-RS, MST, SMD) as well as on molecular mechanical (MM) point-charge models (GB, PB). The performance of the implicit models was tested by a comparison with experimental water-octanol transfer free energies (ΔG(ow)) by using single- and multiconformation approaches. MD simulations revealed difficulties in a priori estimation of the flexibility features of the solutes from simple structural descriptors, such as the number of rotatable bonds. An increasing accuracy of the calculated ΔG(ow) was observed in the following order: GB1 ~ PB < GB7 ≪ MST < SMD ~ COSMO-RS with a clear distinction identified between MM- and QM-based models, although for the set excluding three largest molecules, the differences among COSMO-RS, MST, and SMD were negligible. It was shown that the single-conformation approach applied to crystal geometries provides a rather accurate estimate of ΔG(ow) for rigid molecules yet fails completely for the flexible ones. The multiconformation approaches improved the performance, but only when the deformation contribution was ignored. It was revealed that for large-scale calculations on small molecules a recent GB model, GB7, provided a reasonable accuracy/speed ratio. In conclusion, the study contributes to the understanding of solvation free energy calculations for physical and medicinal chemistry applications.

Solvation of complex surfaces via molecular density functional theory

We show that classical molecular density functional theory, here in the homogeneous reference fluid approximation in which the functional is inferred from the properties of the bulk solvent, is a powerful new tool to study, at a fully molecular level, the solvation of complex surfaces and interfaces by polar solvents. This implicit solvent method allows for the determination of structural, orientational, and energetic solvation properties that are on a par with all-atom molecular simulations performed for the same system, while reducing the computer time by two orders of magnitude. This is illustrated by the study of an atomistically-resolved clay surface composed of over a thousand atoms wetted by a molecular dipolar solvent. The high numerical efficiency of the method is exploited to carry a systematic analysis of the electrostatic and non-electrostatic components of the surface-solvent interaction within the popular Clay Force Field (CLAYFF). Solvent energetics and structure are found to depend weakly upon the atomic charges distribution of the clay surface, even for a rather polar solvent. We conclude on the consequences of such findings for force-field development.

DFTr Studies of Five‐ and Six‐Residue Cyclic‐β(1→4) Cellulosic Molecules

Density functional (DFT) conformational in vacuo studies of cellobiose have shown that ϕ(H) -anti conformations are low in energy relative to the syn forms, while the ψ(H) -anti forms are higher in energy. Further, as the cellulosic fragments became larger than a disaccharide and new hydrogen bonding interactions between multiple residues become available, stable low energy ϕ(H) -anti, and ψ(H) -anti cellulosic structures became possible. To test the stability of cyclic anti-conformations, a number of β-linked five- and six-residue molecules were created and then energy optimized in solvent (water, n-heptane) using the implicit solvation method COSMO at the B3LYP level of theory. The created symmetric cyclic structures were without distortion. Upon optimization some cyclic conformations were found to be of low energy when compared with linear five- and six-residue chains, after correcting the energy for the exclusion of a water molecule upon cyclization. It was also obvious from the hydrogen bonding network formed above and below the plane of the cyclic structure that these structures could exhibit strong synergistic tendencies. The conformational energy preferences for clockwise "c" and counter-clockwise "r" hydroxyl groups and preference for the hydroxymethyl rotamers is described. Because these structures contain energetically unfavorable flipped conformations in water, that is, dihedral angles of ∼180°/0° or ∼0°/180° in ϕ(H) /ψ(H) , it is clear that the synthesis of these compounds will be challenging.

DFT Energy Optimization of a Large Carbohydrate: Cyclomaltohexaicosaose (CA-26)

CA-26 is the largest cyclodextrin (546 atoms) for which refined X-ray structural data is available. Because of its size, 26 D-glucose residues, it is beyond the scope of study of most ab initio or density functional methods and to date has only been computationally examined using empirical force fields. The crystal structure of CA-26 is folded like a figure "8" into two 10 D-glucoses long antiparallel left-handed V (Verkleisterung)-type helices with a "band-flip" and "kink" at the top and bottom of the helices. DFTr methods were applied to CA-26 to determine if a carbohydrate molecule of this size could be geometry optimized and if it would show structural variances from application of dispersion and/or solvation. The DFTr reduced basis set method developed by the authors uses 4-31G on the carbon atoms of the glucose rings and 6-31+G* on all other atoms. B3LYP is the density functional used to successfully optimize CA-26, and other density functionals were then applied, including a self-consistent charge density functional tight binding (SCC-DFTB) method and the B97D (dispersion-corrected) and B97D-PCM (dispersion + implicit solvent) methods. Heavy atom coordinates were taken from one X-ray structure, fitted with hydrogen atoms, and geometry optimized using PM3 followed by B3LYP/6-31+G*/4-31G optimization. After optimization, the heavy atom rms deviation of the optimized DFTr (B3LYP) structure to the crystal structure was 0.89 Å, the rmsd of the B97D optimization was 1.38 Å, that for B97D-PCM was 0.95 Å, and that for SCC-DFTB was 0.94 Å. These results are very good considering that no explicit water molecules were included in the computational analysis and there were ~32-38 water molecules around each CA-26 molecule in the crystal structure. Tables of internal coordinates and puckering parameters were compared to the X-ray structures, and close correspondence was found.

Pseudo-random number generation for Brownian Dynamics and Dissipative Particle Dynamics simulations on GPU devices

Brownian Dynamics (BD), also known as Langevin Dynamics, and Dissipative Particle Dynamics (DPD) are implicit solvent methods commonly used in models of soft matter and biomolecular systems. The interaction of the numerous solvent particles with larger particles is coarse-grained as a Langevin thermostat is applied to individual particles or to particle pairs. The Langevin thermostat requires a pseudo-random number generator (PRNG) to generate the stochastic force applied to each particle or pair of neighboring particles during each time step in the integration of Newton’s equations of motion. In a Single-Instruction-Multiple-Thread (SIMT) GPU parallel computing environment, small batches of random numbers must be generated over thousands of threads and millions of kernel calls. In this communication we introduce a one-PRNG-per-kernel-call-per-thread scheme, in which a micro-stream of pseudorandom numbers is generated in each thread and kernel call. These high quality, statistically robust micro-streams require no global memory for state storage, are more computationally efficient than other PRNG schemes in memory-bound kernels, and uniquely enable the DPD simulation method without requiring communication between threads.

Predicting Fragment Binding Poses Using a Combined MCSS MM-GBSA Approach

Improved methods are required to predict the position and orientation (pose) of binding to the target protein of low molecular weight compounds identified in fragment screening campaigns. This is particularly important to guide initial chemistry to generate structure-activity relationships for the cases where a high resolution structure cannot be obtained. We have assessed the benefit of an implicit solvent method for assessment of fragment binding poses generated by the Multiple Copy Simultaneous Search (MCSS) method in CHARMm. Additionally, the effect of using multiple receptor structures for a flexible receptor is investigated. The original MCSS performance -50% of fragment positions accurately predicted and scored - was increased up to 67% by scoring MCSS energy minima with a Molecular Mechanics Generalized Born approach with molecular volume integration and Surface Area model (MM-GBSA). The same increase in performance (but occasionally for different targets) was observed when using the docking program GOLD followed by MM-GBSA rescoring. The combined results from both methods resulted in a higher success rate emphasizing that a comparison of different docking methods can increase the correct identification of binding poses. For a receptor where multiple structures are available, Hsp90, the average performance on randomly adding receptor structures was also investigated. The results suggest that predictions using these docking methods can be used with some confidence to guide chemical optimization, if the structure of the target either remains relatively fixed on ligand binding, or if a number of crystal structures are available with diverse ligands bound and there is information on the positions of key water molecules in the binding site.

Rapidly calculated DFT relaxed iso-potential ϕ/ψ maps: β-cellobiose

New cellobiose ϕH/ψH maps are generated using a mixed basis set DFTr method, found to achieve a high level of confidence while reducing computer resources dramatically. Relaxed iso-potential maps are created for different conformational states of cellobiose, showing how glycosidic bond dihedral angles vary as different sets of hydroxymethyl rotamers and hydroxyl directions are examined. These maps are generated, fixing the dihedral ϕH and ψH values at ten degree intervals and energy optimizing the remaining geometry using the B3LYP/6-31+G* functional for all atoms except carbon atoms, where the functional B3LYP was used with the mixed basis set, 4-31G. Mapping results are compared between in vacuo structures using the mixed basis set, in vacuo using the full basis set, and those in which the implicit solvent method, COSMO, is included with the mixed basis set. Results show significant changes in position of energy minima with variation in hydroxyl rotamers and with application of solvent. Unique to this study is the mapping of the hydration energy at each ϕH/ψH point on the map using the energy derived at each point by applying COSMO. Using hydration gradients as a guide one observes directional solvent driven changes in the minimum energy positions. Interesting internal coordinate variances are described.

DFTMD studies of β-cellobiose: conformational preference using implicit solvent

Previous DFT in vacuo studies on the conformational preferences for cellobiose showed that upon optimization the φ H -anti conformations were of lower energy than the syn forms. Upon optimization using an implicit solvation method, COSMO, the syn or observed form was still not predicted to be of lower energy than the φ H -anti form, even though optimization after addition of several explicit water molecules did show a relative energy difference favoring the syn form. In order to examine the predictive ability of COSMO on this carbohydrate, constant energy dynamics, DFTMD, simulations were carried out on low energy syn and φ H -anti conformations with and without COSMO included during the dynamics. The resulting analysis confirmed that when COSMO is included in the dynamics, the syn conformations become energetically favored over the φ H -anti forms suggesting that both solvent and entropy play roles in dictating the solution conformation of cellobiose. Analysis of the dynamic runs includes distributions of selected dihedral angles versus time, conformational transitions, and populations of some quasi-planar, boat, skew forms during the simulations.

Conformational preferences and pKa value of selenocysteine residue

The conformational preferences of the L-selenocysteine (Sec) dipeptides with selenol and selenolate groups (Ac-Sec-NHMe and Ac-Sec(-) -NHMe, respectively) and the apparent (i.e., macroscopic) pK(a) value of the Sec residue have been studied using the dispersion-corrected density functionals M06-2X and B2PLYP-D with the implicit solvation method in the gas phase and in water. In the gas phase, the backbone-to-backbone and/or side chain-to-backbone hydrogen bonds are found to contribute in stabilizing the most preferred conformations for the Sec and Sec(-) residues, as seen for the Cys and Cys(-) residues. However, the polyproline II-like conformations prevail over the conformations with the backbone-to-backbone hydrogen bonds in water because of the weakened hydrogen bonds by the favorable direct interactions between the backbone CO and HN groups and water molecules. The Sec and Sec(-) residues are found to adopt more various conformations than the Cys and Cys(-) residues in water, although the most preferred conformations of the neutral and/or anionic forms of the two residues are similar each other in the gas phase and in water. Using the statistically weighted free energies of the Sec and Sec(-) dipeptides in the gas phase and their solvation free energies, the pK(a) value of the Sec residue is estimated to be 5.47 at 25°C, which is in good agreement with the experimental value of 5.43 ± 0.02. It is found that the lower pK(a) value of the selenol side chain for the Sec residue by ∼3 units than the thiol side chain for the Cys residue is ascribed to the higher gas-phase acidity of the Sec residue.

Desolvation Costs of Salt Bridges across Protein Binding Interfaces: Similarities and Differences between Implicit and Explicit Solvent Models.

The prevalence of salt bridges across protein binding interfaces is surprising given the significant costs of desolvating the two charged groups upon binding. These desolvation costs, which are difficult to examine using laboratory experiments, have been computed in previous studies using the Poisson−Boltzmann (PB) implicit solvent model. Here, for the first time, we directly compare the PB implicit solvent model with several explicit water models in computing the desolvation penalties of salt bridges across protein−protein interfaces. We report both overall agreement as well as significant differences between the implicit and explicit solvent results. These differences highlight challenges to be faced in the application of implicit solvent methods.

Application of a novel in silico high-throughput screen to identify selective inhibitors for protein–protein interactions

Increasing numbers of target protein structures available for computational studies makes the structure-based screening paradigm more attractive for initial hit indentification. We have developed a novel in silico screening methodology incorporating Molecular Mechanics (MM)/implicit solvent methods to evaluate binding free energies and applied this technology to the identification of inhibitors of the TLR4/MD-2 interaction.

27 ps DFT molecular dynamics simulation of α‐maltose: A reduced basis set study

DFT molecular dynamics simulations are time intensive when carried out on carbohydrates such as alpha-maltose. In a recent publication (Momany et al., J. Mol. Struct. THEOCHEM, submitted) forces for dynamics were generated from B3LYP/6-31+G* electronic structure calculations. The implicit solvent method COSMO was applied to simulate the solution environment. Here we present a modification of the DFT method that keeps the critical aspects of the larger basis set (B3LYP/6-31+G*) while allowing the less-essential atom interactions to be calculated using a smaller basis set, thus allowing for faster completion without sacrificing the interactions dictating the hydrogen bonding networks in alpha-maltose. In previous studies, the gg'-gg-c solvated form quickly converged to the "r" form during a 5 ps dynamics run. This important conformational transition is tested by carrying out a long 27 ps simulation. The trend for the "r" conformer to be most stable during dynamics when fully solvated, is confirmed, resulting in approximately 20/80% c/r population. Further, the study shows that considerable molecular end effects are important, the reducing end being fairly stable, the O6-H pointing at the O5, while the nonreducing end moves freely to take on different conformations. Some "kink" and transition state forms are populated during the simulation. The average H1'...H4 distance of 2.28 A confirms that the syn form is the primary glycosidic conformation, while the average C1'-O1'-C4 bond angle was 118.8 degrees , in excellent agreement with experimental values. The length of this simulation allowed the evaluation of vibrational frequencies by Fourier transform of the velocity correlation function, taken from different time segments along the simulation path.

Molecular Dynamics Studies of Human β-Glucuronidase

Problem statement: The enzyme β-glucuronidase is being used as a reporter molecule in the area of genetic engineering, as a component of prodrug therapy in cancer treatment and in the scouring process of cotton fabrics. However, a detailed understanding of the factors responsible for the stability and the activity of this enzyme is still not available. Molecular Dynamics (MD) simulations provide an estimate of equilibrium and dynamic properties of enzyme systems that cannot be calculated analytically. With this perspective, molecular dynamics simulations of human β-glucuronidase (GUS) have been carried out to determine the behavior of this enzyme in vacuum and solvent environments at a defined temperature. Approach: CHARMM force field along with distance dependent dielectric model was used to represent the solvent environment in the MD simulations. The parameters employed in various stages of MD simulations had been selected based on repeated trials under various conditions as a method of choosing the optimum parameters for each stage. Results: It was found that simulations in vacuum caused the backbone of GUS to have smaller fluctuations from their mean values compared with the fluctuation in implicit solvent simulations, due to the fact that vacuum environment did not provide for the electrostatic interactions affecting the backbone of GUS that may otherwise exist in a solvent environment. Conclusion: Inclusion of solvent effects in MD simulations is crucial in understanding structural flexibility and stability of β-glucuronidase. Implicit solvent method can provide a realistic inclusion of backbone flexibility and structural compactness of GUS, which will have profound influence on the stability and activity of the enzymes, with a marginal increase in computational time.

DFT molecular dynamics (DFTMD) simulations of carbohydrates: COSMO solvated α-maltose

Density functional molecular dynamics (DFTMD) is carried out on low-energy conformations of α-maltose. Finite temperature molecular dynamics trajectories are generated with forces obtained from B3LYP/6-31+G∗ electronic structure calculations. The implicit solvent method COSMO is applied to simulate the solution environment. Each simulation is carried out for ∼5 ps, starting from low-energy optimized geometries, including different hydroxymethyl rotamers and hydroxyl clockwise, ‘ c’, and counterclockwise, ‘ r’, orientations. The gg′-gg- r solvated form is of lowest relative energy by ∼0.6 kcal/mol relative to the solvated gg′-gg- c form, the latter conformation tending to converge to the ‘ r’ form during dynamics. Conformational transitions and conformers residing as ‘kink’ forms were observed during vacuum runs. In one case, the syn gt’-gt- r + COSMO conformer moved during dynamics into a ‘kink’ conformation, remaining there for most of the 5 ps simulation. However, when this same conformer was started from a minimum energy ‘kink’ form, it rapidly reverted into the normal syn conformation in which the H1′⋯H4 hydrogen atoms across the glycosidic bond were ∼2.15 Å on average. Other solvated structures showed a perfunctory preference for the ‘kink’ conformation during dynamics, even though in previous optimization studies the ‘kink’ conformations were of higher energy than the syn conformers. Similarly, ‘band-flip’ conformational studies showed that the ‘ c’ and ‘ r’ forms differed, ‘ c’ undergoing transitions to the syn form, the ‘ r’ form staying in the band-flip conformation over the 5 ps simulation. The trend for the ‘ r’ conformers to be more stable than the ‘ c’ conformers when fully solvated, appears to confirm optimization studies, although transitions to stable partial ‘ c’ forms produces some confusing conformational effects. For example, in one case a hydroxymethyl O6 H⋯O6′ hydrogen bond locked into a glycosidic ‘kink’ structure, allowing the O3 H and O2′ H’ hydroxyl groups between rings to rotate more freely because of the enlarged distances between the groups across the glycosidic bridge.

Molecular dynamics study of interaction and substrate channeling between neuron‐specific enolase and B‐type phosphoglycerate mutase

Phosphoglycerate mutase (PGM) and enolase are consecutive enzymes in the glycolytic pathway. We used molecular dynamics simulation to examine the interaction of human B-type PGM (dPGM-B) and neuron-specific enolase (NSE). Specifically, we studied the interactions of 31 orientations of these enzymes by means of the effective energy function implicit solvation method. Interactions between active regions of the enzymes occurred preferentially, although the strongest interactions appeared to be between the back side of NSE and the active regions of dPGM-B. Cleavage of 2PG from dPGM-B was investigated, and the Ser(14)-Leu(30) loop of dPGM-B is suggested as a cleavage site and, likely, another entrance site of a ligand. Substrate channeling between the enzymes was observed when NSE with its active regions Leu(11)-Asn(16), Arg(49)-Lys(59), and Gly(155)-Ala(158) covered the Ser(14)-Leu(30) loop of dPGM-B. Analyses of the results make us believe that the channeling between PGM and enolase "benefits" from weak interaction. The probability of formation of channeling favorable complex is estimated to be up to 5%, while functional interaction between NSE and dPGM-B might be as high as 20%. NSE and dPGM-B functional interaction seems not to be isotype specific.