A solvation-free-energy functional: A reference-modified density functional formulation

RISMiCal: A software package to perform fast RISM/3D-RISM calculations

Solvent plays an essential role in a variety of chemical, physical, and biological processes that occur in the solution phase. The reference interaction site model (RISM) and its three-dimensional extension (3D-RISM) serve as powerful computational tools for modeling solvation effects in chemical reactions, biological functions, and structure formations. We present the RISM integrated calculator (RISMiCal) program package, which is based on RISM and 3D-RISM theories with fast GPU code. RISMiCal has been developed as an integrated RISM/3D-RISM program that has interfaces with external programs such as Gaussian16, GAMESS, and Tinker. Fast 3D-RISM programs for single- and multi-GPU codes written in CUDA would enhance the availability of these hybrid methods because they require the performance of many computationally expensive 3D-RISM calculations. We expect that our package can be widely applied for chemical and biological processes in solvent. The RISMiCal package is available at https://rismical-dev.github.io.

Effect of Main and Side Chains on the Folding Mechanism of the Trp-Cage Miniprotein

Proteins that do not fold into their functional native state have been linked to diseases. In this study, the influence of the main and side chains of individual amino acids on the folding of the tryptophan cage (Trp-cage), a designed 20-residue miniprotein, was analyzed. For this purpose, we calculated the solvation free energy (SFE) contributions of individual atoms by using the 3D-reference interaction site model with the atomic decomposition method. The mechanism by which the Trp-cage is stabilized during the folding process was examined by calculating the total energy, which is the sum of the conformational energy and SFE. The folding process of the Trp-cage resulted in a stable native state, with a total energy that was 62.4 kcal/mol lower than that of the unfolded state. The solvation entropy, which is considered to be responsible for the hydrophobic effect, contributed 31.3 kcal/mol to structural stabilization. In other words, the contribution of the solvation entropy accounted for approximately half of the total contribution to Trp-cage folding. The hydrophobic core centered on Trp6 contributed 15.6 kcal/mol to the total energy, whereas the solvation entropy contribution was 6.3 kcal/mol. The salt bridge formed by the hydrophilic side chains of Asp9 and Arg16 contributed 10.9 and 5.0 kcal/mol, respectively. This indicates that not only the hydrophobic core but also the salt bridge of the hydrophilic side chains gain solvation entropy and contribute to stabilizing the native structure of the Trp-cage.

Analysis of Protein Folding Simulation with Moving Root Mean Square Deviation

We apply moving root-mean-square deviation (mRMSD), which does not require a reference structure, as a method for analyzing protein dynamics. This method can be used to calculate the root-mean-square deviation (RMSD) of structure between two specified time points and to analyze protein dynamics behavior through time series analysis. We applied this method to the Trp-cage trajectory calculated by the Anton supercomputer and found that it shows regions of stable states as well as the conventional RMSD. In addition, we extracted a characteristic structure in which the side chains of Asp1 and Arg16 form hydrogen bonds near the most stable structure of the Trp-cage. We also determined that ≥20 ns is an appropriate time interval to investigate protein dynamics using mRMSD. Applying this method to NuG2 protein, we found that mRMSD can be used to detect regions of metastable states in addition to the stable state. This method can be applied to molecular dynamics simulations of proteins whose stable structures are unknown.

Solvation entropy, enthalpy and free energy prediction using a multi-task deep learning functional in 1D-RISM

Simultaneous calculation of entropies, enthalpies and free energies has been a long-standing challenge in computational chemistry, partly because of the difficulty in obtaining estimates of all three properties from a single consistent simulation methodology. This has been particularly true for methods from the Integral Equation Theory of Molecular Liquids such as the Reference Interaction Site Model which have traditionally given large errors in solvation thermodynamics. Recently, we presented pyRISM-CNN, a combination of the 1 Dimensional Reference Interaction Site Model (1D-RISM) solver, pyRISM, with a deep learning based free energy functional, as a method of predicting solvation free energy (SFE). With this approach, a 40-fold improvement in prediction accuracy was delivered for a multi-solvent, multi-temperature dataset when compared to the standard 1D-RISM theory [Fowles et al., Digital Discovery, 2023, 2, 177-188]. Here, we report three further developments to the pyRISM-CNN methodology. Firstly, solvation free energies have been introduced for organic molecular ions in methanol or water solvent systems at 298 K, with errors below 4 kcal mol^-1 obtained without the need for corrections or additional descriptors. Secondly, the number of solvents in the training data has been expanded from carbon tetrachloride, water and chloroform to now also include methanol. For neutral solutes, prediction errors nearing or below 1 kcal mol^-1 are obtained for each organic solvent system at 298 K and water solvent systems at 273-373 K. Lastly, pyRISM-CNN was successfully applied to the simultaneous prediction of solvation enthalpy, entropy and free energy through a multi-task learning approach, with errors of 1.04, 0.98 and 0.47 kcal mol^-1, respectively, for water solvent systems at 298 K.

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.

Structural Stability Analysis of Proteins Using End-to-End Distance: A 3D-RISM Approach

The stability of a protein is determined from its properties and surrounding solvent. In our previous study, the total energy as a sum of the conformational and solvation free energies was demonstrated to be an appropriate energy function for evaluating the stability of a protein in a protein folding system. We plotted the various energies against the root mean square deviation, required as a reference structure. Herein, we replotted the various energies against the end-to-end distance between the N- and C-termini, which is not a required reference and is experimentally measurable. The solvation free energies for all proteins tend to be low as the end-to-end distance increases, whereas the conformational energies tend to be low as the end-to-end distance decreases. The end-to-end distance is one of interesting measures to study the behavior of proteins.

Common and unique strategies of myoglobin evolution for deep-sea adaptation of diving mammals

Myoglobin (Mb) is highly concentrated in the myocytes of diving mammals such as whales and seals, in comparison with land animals, and its molecular evolution has played a crucial role in their deep-sea adaptation. We previously resurrected ancestral whale Mbs and demonstrated the evolutional strategies for higher solubility under macromolecular crowding conditions. Pinnipeds, such as seals and sea lions, are also expert diving mammals with Mb-rich muscles. In the present study, we resurrected ancestral pinniped Mbs and investigated their biochemical and structural properties. Comparisons between pinniped and whale Mbs revealed the common and distinctive strategies for the deep-sea adaptation. The overall evolution processes, gaining precipitant tolerance and improving thermodynamic stability, were commonly observed. However, the strategies for improving the folding stability differed, and the pinniped Mbs exploited the shielding of hydrophobic surfaces more effectively than the whale Mbs.

Water-mediated interactions destabilize proteins

Proteins are folded to avoid exposure of the nonpolar groups to water because water‐mediated interactions between nonpolar groups are a promising factor in the thermodynamic stabilities of proteins—which is a well‐accepted view as one of the unique effects of hydrophobic interactions. This article poses a critical question for this classical view by conducting an accurate solvation free‐energy calculation for a thermodynamic cycle of a protein folding using a liquid‐state density functional theory. Here, the solvation‐free energy for a leucine zipper formation was examined in the coiled‐coil protein GCN4‐p1, a typical model for hydrophobic interactions, which demonstrated that water‐mediated interactions were unfavorable for the association of nonpolar groups in the native state, while the dispersion forces between them were, instead, responsible for the association. Furthermore, the present analysis well predicted the isolated helical state stabilized by pressure, which was previously observed in an experiment. We reviewed the problems in the classical concept and semiempirical presumption that the energetic cost of the hydration of nonpolar groups is a driving force of folding.

Mutation-induced change in chignolin stability from π-turn to α-turn

A mutation from threonine to proline at the eighth residue in chignolin changes π-turn to α-turn.

Analysis of molecular dynamics simulations of 10-residue peptide, chignolin, using statistical mechanics: Relaxation mode analysis and three-dimensional reference interaction site model theory

Molecular dynamics simulation is a fruitful tool for investigating the structural stability, dynamics, and functions of biopolymers at an atomic level. In recent years, simulations can be performed on time scales of the order of milliseconds using special purpose systems. Since the most stable structure, as well as meta-stable structures and intermediate structures, is included in trajectories in long simulations, it is necessary to develop analysis methods for extracting them from trajectories of simulations. For these structures, methods for evaluating the stabilities, including the solvent effect, are also needed. We have developed relaxation mode analysis to investigate dynamics and kinetics of simulations based on statistical mechanics. We have also applied the three-dimensional reference interaction site model theory to investigate stabilities with solvent effects. In this paper, we review the results for designing amino-acid substitution of the 10-residue peptide, chignolin, to stabilize the misfolded structure using these developed analysis methods.

Theoretical analysis on thermodynamic stability of chignolin

Understanding the dominant factor in thermodynamic stability of proteins remains an open challenge. Kauzmann's hydrophobic interaction hypothesis, which considers hydrophobic interactions between nonpolar groups as the dominant factor, has been widely accepted for about sixty years and attracted many scientists. The hypothesis, however, has not been verified or disproved because it is difficult, both theoretically and experimentally, to quantify the solvent effects on the free energy change in protein folding. Here, we developed a computational method for extracting the dominant factor behind thermodynamic stability of proteins and applied it to a small, designed protein, chignolin. The resulting free energy profile quantitatively agreed with the molecular dynamics simulations. Decomposition of the free energy profile indicated that intramolecular interactions predominantly stabilized collapsed conformations, whereas solvent-induced interactions, including hydrophobic ones, destabilized them. These results obtained for chignolin were consistent with the site-directed mutagenesis and calorimetry experiments for globular proteins with hydrophobic interior cores.

Tracing whale myoglobin evolution by resurrecting ancient proteins

Extant cetaceans, such as sperm whale, acquired the great ability to dive into the ocean depths during the evolution from their terrestrial ancestor that lived about 50 million years ago. Myoglobin (Mb) is highly concentrated in the myocytes of diving animals, in comparison with those of land animals, and is thought to play a crucial role in their adaptation as the molecular aqualung. Here, we resurrected ancestral whale Mbs, which are from the common ancestor between toothed and baleen whales (Basilosaurus), and from a further common quadrupedal ancestor between whale and hippopotamus (Pakicetus). The experimental and theoretical analyses demonstrated that whale Mb adopted two distinguished strategies to increase the protein concentration in vivo along the evolutionary history of deep sea adaptation; gaining precipitant tolerance in the early phase of the evolution, and increase of folding stability in the late phase.

Are the Sublimation Thermodynamics of Organic Molecules Predictable?

We compare a range of computational methods for the prediction of sublimation thermodynamics (enthalpy, entropy, and free energy of sublimation). These include a model from theoretical chemistry that utilizes crystal lattice energy minimization (with the DMACRYS program) and quantitative structure property relationship (QSPR) models generated by both machine learning (random forest and support vector machines) and regression (partial least squares) methods. Using these methods we investigate the predictability of the enthalpy, entropy and free energy of sublimation, with consideration of whether such a method may be able to improve solubility prediction schemes. Previous work has suggested that the major source of error in solubility prediction schemes involving a thermodynamic cycle via the solid state is in the modeling of the free energy change away from the solid state. Yet contrary to this conclusion other work has found that the inclusion of terms such as the enthalpy of sublimation in QSPR methods does not improve the predictions of solubility. We suggest the use of theoretical chemistry terms, detailed explicitly in the Methods section, as descriptors for the prediction of the enthalpy and free energy of sublimation. A data set of 158 molecules with experimental sublimation thermodynamics values and some CSD refcodes has been collected from the literature and is provided with their original source references.

Accuracy enhancement in the estimation of molecular hydration free energies by implementing the intramolecular hydrogen bond effects

Background

The formation of intramolecular hydrogen bonds (IHBs) may induce the remarkable changes in molecular physicochemical properties. Within the framework of the extended solvent-contact model, we investigate the effect of implementing the IHB interactions on the accuracy in estimating the molecular hydration free energies.

Results

The performances of hydration free energy functions including and excluding the IHB parameters are compared using the molecules distributed for SAMPL4 blind prediction challenge and those in Free Solvation Database (FSD). The calculated hydration free energies with IHB effects are found to be in considerably better agreement with the experimental data than those without them. For example, the root mean square error of the estimation decreases from 2.56 to 1.66 and from 1.73 to 1.54 kcal/mol for SAMPL4 and FSD molecules, respectively, due to the extension of atomic parameter space to cope with IHBs.

Conclusions

These improvements are made possible by reducing the overestimation of attractive interactions between water and the solute molecules involving IHBs. The modified hydration free energy function is thus anticipated to be useful for estimating the desolvation cost for various organic molecules.

Electronic supplementary material

The online version of this article (doi:10.1186/s13321-015-0106-2) contains supplementary material, which is available to authorized users.