The nature in which the protecting osmolyte trimethylamine-N-oxide (TMAO) and the denaturing osmolyte urea affect protein stability is investigated simulating a deca-alanine peptide model in multiple conformations of the denatured ensemble. Binary solutions of both osmolytes and mixed osmolyte solutions at physiologically-relevant concentrations of 2:1 (urea:TMAO) are studied using standard molecular dynamics simulations and solvation free energy calculations. Component analysis reveals the differences in the importance of the van der Waals (vdW) and electrostatic interactions for protecting and denaturing osmolytes. We find that urea denaturation governed by transfer free energy differences is dominated by vdW attractions, whereas TMAO exerts its effect by causing unfavorable electrostatic interactions both in the binary solution and mixed osmolyte solution. Analysis of the results showed no evidence in the ternary solution of disruption of the correlations among the peptide and osmolytes, nor of significant changes in the strength of the water hydrogen bond network.

The study of organic osmolytes has been pivotal in demonstrating the role of solvent effects on the protein backbone in the folding process. Whereas a thermodynamic description of the interactions between the protein backbone and osmolyte has been well defined, the structural analysis of the effect of osmolyte on the protein backbone has been incomplete. Therefore, we have carried out simulations of a peptide backbone model, glycine 15 in protecting osmolyte TMAO solution in order to determine the effect of the solution structure on the conformation of the peptide backbone. We demonstrate that the models chosen show that the ensemble of backbone structures shifts towards a more collapsed state in TMAO solution as compared to pure water solution. The collapse is consistent with preferential exclusion of the osmolyte caused by unfavorable interactions between osmolyte and peptide backbone. The exclusion is due to strong triplet correlations of osmolyte, water, and peptide backbone. This provides a clear mechanism demonstrating that even a modest concentration of TMAO forces the protein backbone to adopt a more collapsed structure in the absence of sidechain effects.

We performed molecular dynamics simulations of urea solutions at different concentrations with two urea models (OPLS and KBFF) to examine the structures responsible for the thermodynamic solution properties. Our simulation results showed that hydrogen-bonding properties such as the average number of hydrogen bonds and their lifetime distributions were nearly constant at all concentrations between infinite dilution and the solubility limit. This implies that the characterization of urea-water solutions in the molarity concentration scale as nearly ideal is a result of facile local hydrogen bonding rather than a global property. Thus, urea concentration does not influence the local propensity for hydrogen bonds, only how they are satisfied. By comparison, the KBFF model of urea donated fewer hydrogen bonds than OPLS. We found that the KBFF urea model in TIP3P water better reproduced the experimental density and diffusion constant data. Preferential solvation analysis showed that there were weak urea-urea and water-water associations in OPLS solution at short distances, but there were no strong associations. We divided urea molecules into large, medium, and small clusters to examine fluctuation properties and found that any particular urea molecule did not stay in the same cluster for a long time. We found neither persistent nor large clusters.

The transfer model implying additivity of the peptide backbone free energy of transfer is computationally tested. Molecular dynamics simulations are used to determine the extent of change in transfer free energy (DG tr ) with increase in chain length of oligoglycine with capped end groups. Solvation free energies of oligoglycine models of varying lengths in pure water and in the osmolyte solutions, 2M urea and 2M trimethylamine N-oxide (TMAO), were calculated from simulations of all atom models, and DG tr values for peptide backbone transfer from water to the osmolyte solutions were determined. The results show that the transfer free energies change linearly with increasing chain length, demonstrating the principle of additivity, and provide values in reasonable agreement with experiment. The peptide backbone transfer free energy contributions arise from van der Waals interactions in the case of transfer to urea, but from electrostatics on transfer to TMAO solution. The simulations used here allow for the calculation of the solvation and transfer free energy of longer oligoglycine models to be evaluated than is currently possible through experiment. The peptide backbone unit computed transfer free energy of 254 cal/mol/M compares quite favorably with 243 cal/mol/M determined experimentally.

The binding mode prediction is of great importance to structure-based drug design. The discrimination of various binding poses of ligand generated by docking is a great challenge not only to docking score functions but also to the relatively expensive free energy calculation methods. Here we systematically analyzed the stability of various ligand poses under molecular dynamics (MD) simulation. First, a data set of 120 complexes was built based on the typical physicochemical properties of drug-like ligands. Three potential binding poses (one correct pose and two decoys) were selected for each ligand from self-docking in addition to the experimental pose. Then, five independent MD simulations for each pose were performed with different initial velocities for the statistical analysis. Finally, the stabilities of ligand poses under MD were evaluated and compared with the native one from crystal structure. We found that about 94% of the native poses were maintained stable during the simulations, which suggests that MD simulations are accurate enough to judge most experimental binding poses as stable properly. Interestingly, incorrect decoy poses were maintained much less and 38-44% of decoys could be excluded just by performing equilibrium MD simulations, though 56-62% of decoys were stable. The computationally-heavy binding free energy calculation can be performed only for these survived poses.

We have developed a prediction method for the binding structures of ligands with proteins. Our method consists of three steps. First, replica-exchange umbrella sampling simulations are performed along the distance between a putative binding site of a protein and a ligand as the reaction coordinate. Second, we obtain the potential of mean force (PMF) of the unbiased system using the weighted histogram analysis method and determine the distance that corresponds to the global minimum of PMF. Third, structures that have this global-minimum distance and energy values around the average potential energy are collected and analyzed using the principal component analysis. We predict the binding structure as the global-minimum free energy state on the free energy landscapes along the two major principal component axes. As test cases, we applied our method to five protein-ligand complex systems. Starting from the configuration in which the protein and the ligand are far away from each other in each system, our method predicted the ligand binding structures in excellent agreement with the experimental data from Protein Data Bank.

Activity coefficients of urea solutions are calculated to explore the mechanism of its solution properties, which form the basis for its well-known use as a strong protein denaturant. We perform free energy simulations of urea solutions in different urea concentrations using two urea models (OPLS and KBFF models) to calculate and decompose the activity coefficients. For the case of urea, we clarify the concept of the ideal solution in different concentration scales and standard states and its effect on our subsequent analysis. The analytical form of activity coefficients depends on the concentration units and standard states. For both models studied, urea displays a weak concentration dependence for excess chemical potential. However, for the OPLS force-field model, this results from contributions that are independent of concentration to the van der Waals and electrostatic components whereas for the KBFF model those components are nontrivial but oppose each other. The strong ideality of urea solutions in some concentration scales (incidentally implying a lack of water perturbation) is discussed in terms of recent data and ideas on the mechanism of urea denaturation of proteins.

The electrostatic (ΔGel), van der Waals cavity-formation (ΔGvdw), and total (ΔG) solvation free energies for 10 alanine peptides ranging in length (n) from 1 to 10 monomers were calculated. The free energies were computed both with fixed, extended conformations of the peptides and again for some of the peptides without constraints. The solvation free energies, ΔGel, and components ΔGvdw, and ΔG, were found to be linear in n, with the slopes of the best-fit lines being γel, γvdw, and γ, respectively. Both γel and γ were negative for fixed and flexible peptides, and γvdw was negative for fixed peptides. That γvdw was negative was surprising, as experimental data on alkanes, theoretical models, and MD computations on small molecules and model systems generally suggest that γvdw should be positive. A negative γvdw seemingly contradicts the notion that ΔGvdw drives the initial collapse of the protein when it folds by favoring conformations with small surface areas. When we computed ΔGvdw for the flexible peptides, thereby allowing the peptides to assume natural ensembles of more compact conformations, γvdw was positive. Because most proteins do not assume extended conformations, a ΔGvdw that increases with increasing surface area may be typical for globular proteins. An alternative hypothesis is that the collapse is driven by intramolecular interactions. We find few intramolecular h-bonds but show that the intramolecular van der Waal’s interaction energy is more favorable for the flexible than for the extended peptides, seemingly favoring this hypothesis. The large fluctuations in the vdw energy may make attributing the collapse of the peptide to this intramolecular energy difficult.