Determination of side-chain-rotamer and side-chain and backbone virtual-bond-stretching potentials of mean force from AM1 energy surfaces of terminally-blocked amino-acid residues, for coarse-grained simulations of protein structure and folding. II. Results, comparison with statistical potentials, and implementation in the UNRES force field

Using the harmonic-approximation approach of the accompanying article and AM1 energy surfaces of terminally blocked amino-acid residues, we determined physics-based side-chain rotamer potentials and the side-chain virtual-bond-deformation potentials of 19 natural amino-acid residues with side chains. The potentials were approximated by analytical formulas and implemented in the UNRES mesoscopic dynamics program. For comparison, the corresponding statistical potentials were determined from 19,682 high-resolution protein structures. The low free-energy region of both the AM1-derived and the statistical potentials is determined by the valence geometry and the L-chirality, and its size increases with side-chain flexibility and decreases with increasing virtual-bond-angle theta. The differences between the free energies of rotamers are greater for the AM1-derived potentials compared with the statistical potentials and, for alanine and other residues with small side chains, a region corresponding to the C(ax)(7) conformation has remarkably low free-energy for the AM1-derived potentials, as opposed to the statistical potentials. These differences probably result from the interactions between neighboring residues and indicate the need for introduction of cooperative terms accounting for the coupling between side-chain rotamer and backbone interactions. Both AM1-derived and statistical virtual-bond-deformation potentials are multimodal for flexible side chains and are topologically similar; however, the regions of minima of the statistical potentials are much narrower, which probably results from imposing restraints in structure determination. The force field with the new potentials was preliminarily optimized using the FBP WW domain (1E0L) and the engrailed homeodomain (1ENH) as training proteins and assessed to be reasonably transferable.

Investigation of protein folding by coarse-grained molecular dynamics with the UNRES force field

Coarse-grained molecular dynamics simulations offer a dramatic extension of the time-scale of simulations compared to all-atom approaches. In this article, we describe the use of the physics-based united-residue (UNRES) force field, developed in our laboratory, in protein-structure simulations. We demonstrate that this force field offers about a 4000-times extension of the simulation time scale; this feature arises both from averaging out the fast-moving degrees of freedom and reduction of the cost of energy and force calculations compared to all-atom approaches with explicit solvent. With massively parallel computers, microsecond folding simulation times of proteins containing about 1000 residues can be obtained in days. A straightforward application of canonical UNRES/MD simulations, demonstrated with the example of the N-terminal part of the B-domain of staphylococcal protein A (PDB code: 1BDD, a three-alpha-helix bundle), discerns the folding mechanism and determines kinetic parameters by parallel simulations of several hundred or more trajectories. Use of generalized-ensemble techniques, of which the multiplexed replica exchange method proved to be the most effective, enables us to compute thermodynamics of folding and carry out fully physics-based prediction of protein structure, in which the predicted structure is determined as a mean over the most populated ensemble below the folding-transition temperature. By using principal component analysis of the UNRES folding trajectories of the formin-binding protein WW domain (PDB code: 1E0L; a three-stranded antiparallel beta-sheet) and 1BDD, we identified representative structures along the folding pathways and demonstrated that only a few (low-indexed) principal components can capture the main structural features of a protein-folding trajectory; the potentials of mean force calculated along these essential modes exhibit multiple minima, as opposed to those along the remaining modes that are unimodal. In addition, a comparison between the structures that are representative of the minima in the free-energy profile along the essential collective coordinates of protein folding (computed by principal component analysis) and the free-energy profile projected along the virtual-bond dihedral angles gamma of the backbone revealed the key residues involved in the transitions between the different basins of the folding free-energy profile, in agreement with existing experimental data for 1E0L .

Mechanism of formation of the C-terminal beta-hairpin of the B3 domain of the immunoglobulin binding protein G from Streptococcus. III. Dynamics of long-range hydrophobic interactions

A 20-residue peptide, IG(42-61), derived from the C-terminal beta-hairpin of the B3 domain of the immunoglobulin binding protein G from Streptoccocus was studied using circular dichroism, nuclear magnetic resonance (NMR) spectroscopy at various temperatures and by differential scanning calorimetry (DSC). Unlike other related peptides studied so far, this peptide displays two heat capacity peaks in DSC measurements (at a scanning rate of 1.5 deg/min at a peptide concentration of 0.07 mM), which suggests a three-state folding/unfolding process. The results from DSC and NMR measurements suggest the formation of a dynamic network of hydrophobic interactions stabilizing the structure, which resembles a beta-hairpin shape over a wide range of temperatures (283-313 K). Our results show that IG (42-61) possesses a well-organized three-dimensional structure stabilized by long-range hydrophobic interactions (Tyr50 ... Phe57 and Trp48 ... Val59) at T = 283 K and (Trp48 ... Val59) at 305 and 313 K. The mechanism of beta-hairpin folding and unfolding, as well as the influence of peptide length on its conformational properties, are also discussed.

Potential of mean force of association of large hydrophobic particles: toward the nanoscale limit

The potentials of mean force (PMFs) were determined, in both water with the TIP3P water model and in vacuo, for systems involving formation of nonpolar dimers composed of bicyclooctane, adamantane (both an all-atom model and a sphere with the radius of 3.4 A representing adamantane), and fullerene, respectively. A series of umbrella-sampling molecular dynamics simulations with the AMBER force field were carried out for each pair under both environmental conditions. The PMFs were calculated by using the weighted histogram analysis method. The results were compared with our previously determined PMF for neopentane. The shape of the PMFs for dimers of all four nonpolar molecules is characteristic of hydrophobic interactions with contact and solvent-separated minima and desolvation maxima. The positions of all these minima and maxima change with the size of the nonpolar molecule; for larger molecules they shift toward larger distances. Comparison of the PMFs of the bicyclooctane, adamantane, and fullerene dimers in water and in vacuo shows that hydrophobic interactions in each dimer are different from that for the dimer of neopentane. Interactions in the bicyclooctane, adamantane, and fullerene dimers are stronger in vacuo than in water. These dimers cannot be treated as classical, spherical, hydrophobic objects. The solvent contribution to the PMF was also computed by subtracting the PMF determined in vacuo from that in explicit solvent. The solvent contribution to the PMFs of bicyclooctane, adamantane, and fullerene is positive, as opposed to that of neopentane. The water molecules in the first solvation sphere of both adamantane and neopentane dimers are more ordered as compared to bulk water, with their dipole moments pointing away from the surface of the dimers. The average number of hydrogen bonds per water molecule in the first hydration shell of adamantane is smaller compared to that in bulk water, but this shell is thicker for all-atom adamantane than for neopentane or a spherical model of adamantane. In the second hydration shell, the average number of hydrogen bonds is greater compared to that in bulk water only for neopentane and a spherical model of adamantane but not for the all-atom model. The strength of the hydrophobic interactions shows a linear dependence on the number of carbon atoms both in water and in vacuo. Smaller nonpolar particles interact more strongly in water than in vacuo. For larger molecules, such as bicyclooctane, adamantane and fullerene, the reversed tendency is observed.

Implementation of molecular dynamics and its extensions with the coarse-grained UNRES force field on massively parallel systems; towards millisecond-scale simulations of protein structure, dynamics, and thermodynamics

We report the implementation of our united-residue UNRES force field for simulations of protein structure and dynamics with massively parallel architectures. In addition to coarse-grained parallelism already implemented in our previous work, in which each conformation was treated by a different task, we introduce a fine-grained level in which energy and gradient evaluation are split between several tasks. The Message Passing Interface (MPI) libraries have been utilized to construct the parallel code. The parallel performance of the code has been tested on a professional Beowulf cluster (Xeon Quad Core), a Cray XT3 supercomputer, and two IBM BlueGene/P supercomputers with canonical and replica-exchange molecular dynamics. With IBM BlueGene/P, about 50 % efficiency and 120-fold speed-up of the fine-grained part was achieved for a single trajectory of a 767-residue protein with use of 256 processors/trajectory. Because of averaging over the fast degrees of freedom, UNRES provides an effective 1000-fold speed-up compared to the experimental time scale and, therefore, enables us to effectively carry out millisecond-scale simulations of proteins with 500 and more amino-acid residues in days of wall-clock time.

Relation between free energy landscapes of proteins and dynamics

By examining the molecular dynamics (MD) of protein folding trajectories, generated with the coarse-grained UNRES force field, for the B-domain of staphylococcal protein A and the triple β-strand WW domain from the Formin binding protein 28 (FBP), by principal component analysis (PCA), it is demonstrated how different free energy landscapes (FELs) and folding pathways of trajectories can be, even though they appear to be very similar by visual inspection of the time-dependence of the root-mean-square deviation (rmsd). Approaches to determine the minimal dimensionality of FELs for a correct description of protein folding dynamics are discussed. The correlation between the amplitude of the fluctuations of proteins and the dimensionality of the FELs is shown. The advantage of internal coordinate PCA over Cartesian PCA for small proteins is also illustrated.

Mechanism of formation of the C-terminal beta-hairpin of the B3 domain of the immunoglobulin-binding protein G from Streptococcus. IV. Implication for the mechanism of folding of the parent protein

A 34-residue alpha/beta peptide [IG(28-61)], derived from the C-terminal part of the B3 domain of the immunoglobulin binding protein G from Streptoccocus, was studied using CD and NMR spectroscopy at various temperatures and by differential scanning calorimetry. It was found that the C-terminal part (a 16-residue-long fragment) of this peptide, which corresponds to the sequence of the beta-hairpin in the native structure, forms structure similar to the beta-hairpin only at T = 313 K, and the structure is stabilized by non-native long-range hydrophobic interactions (Val47-Val59). On the other hand, the N-terminal part of IG(28-61), which corresponds to the middle alpha-helix in the native structure, is unstructured at low temperature (283 K) and forms an alpha-helix-like structure at 305 K, and only one helical turn is observed at 313 K. At all temperatures at which NMR experiments were performed (283, 305, and 313 K), we do not observe any long-range connectivities which would have supported packing between the C-terminal (beta-hairpin) and the N-terminal (alpha-helix) parts of the sequence. Such interactions are absent, in contrast to the folding pathway of the B domain of protein G, proposed recently by Kmiecik and Kolinski (Biophys J 2008, 94, 726-736), based on Monte-Carlo dynamics studies. Alternative folding mechanisms are proposed and discussed.

Exploring the parameter space of the coarse-grained UNRES force field by random search: selecting a transferable medium-resolution force field

We explored the energy-parameter space of our coarse-grained UNRES force field for large-scale ab initio simulations of protein folding, to obtain good initial approximations for hierarchical optimization of the force field with new virtual-bond-angle bending and side-chain-rotamer potentials which we recently introduced to replace the statistical potentials. 100 sets of energy-term weights were generated randomly, and good sets were selected by carrying out replica-exchange molecular dynamics simulations of two peptides with a minimal alpha-helical and a minimal beta-hairpin fold, respectively: the tryptophan cage (PDB code: 1L2Y) and tryptophan zipper (PDB code: 1LE1). Eight sets of parameters produced native-like structures of these two peptides. These eight sets were tested on two larger proteins: the engrailed homeodomain (PDB code: 1ENH) and FBP WW domain (PDB code: 1E0L); two sets were found to produce native-like conformations of these proteins. These two sets were tested further on a larger set of nine proteins with alpha or alpha + beta structure and found to locate native-like structures of most of them. These results demonstrate that, in addition to finding reasonable initial starting points for optimization, an extensive search of parameter space is a powerful method to produce a transferable force field.

Mechanism of formation of the C-terminal beta-hairpin of the B3 domain of the immunoglobulin binding protein G from Streptococcus. I. Importance of hydrophobic interactions in stabilization of beta-hairpin structure

We previously studied a 16-amino acid-residue fragment of the C-terminal beta-hairpin of the B3 domain (residues 46-61), [IG(46-61)] of the immunoglobulin binding protein G from Streptoccocus, and found that hydrophobic interactions and the turn region play an important role in stabilizing the structure. Based on these results, we carried out systematic structural studies of peptides derived from the sequence of IG (46-61) by systematically shortening the peptide by one residue at a time from both the C- and the N-terminus. To determine the structure and stability of two resulting 12- and 14-amino acid-residue peptides, IG(48-59) and IG(47-60), respectively, we carried out circular dichroism, NMR, and calorimetric studies of these peptides in pure water. Our results show that IG(48-59) possesses organized three-dimensional structure stabilized by hydrophobic interactions (Tyr50-Phe57 and Trp48-Val59) at T = 283 and 305 K. At T = 313 K, the structure breaks down because of increased chain entropy, but the turn region is preserved in the same position observed for the structure of the whole protein. The breakdown of structure occurs near the melting temperature of this peptide (T(m) = 310 K) measured by differential scanning calorimetry (DSC). The melting temperature of IG(47-60) determined by DSC is T(m) = 330 K and its structure is similar to that of the native beta-hairpin at all (lower) temperatures examined (283-313 K). Both of these truncated sequences are conserved in all known amino acid sequences of the B domains of the immunoglobulin binding protein G from bacteria. Thus, this study contributes to an understanding of the mechanism of folding of this whole family of proteins, and provides information about the mechanism of formation and stabilization of a beta-hairpin structural element.

Mechanism of formation of the C-terminal beta-hairpin of the B3 domain of the immunoglobulin binding protein G from Streptococcus. II. Interplay of local backbone conformational dynamics and long-range hydrophobic interactions in hairpin formation

Two peptides, corresponding to the turn region of the C-terminal beta-hairpin of the B3 domain of the immunoglobulin binding protein G from Streptococcus, consisting of residues 51-56 [IG(51-56)] and 50-57 [IG(50-57)], respectively, were studied by circular dichroism and NMR spectroscopy at various temperatures and by differential scanning calorimetry. Our results show that the part of the sequence corresponding to the beta-turn in the native structure (DDATKT) of the B3 domain forms bent conformations similar to those observed in the native protein. The formation of a turn is observed for both peptides in a broad range of temperatures (T = 283-323 K), which confirms the conclusion drawn from our previous studies of longer sequences from the C-terminal beta-hairpin of the B3 domain of the immunoglobulin binding protein G (16, 14, and 12 residues), that the DDATKT sequence forms a nucleation site for formation of the beta-hairpin structure of peptides corresponding to the C-terminal part of all the B domains of the immunoglobulin binding protein G. We also show and discuss the role of long-range hydrophobic interactions as well as local conformational properties of polypeptide chains in the mechanism of formation of the beta-hairpin structure.

Conformational studies of the C-terminal 16-amino-acid-residue fragment of the B3 domain of the immunoglobulin binding protein G from Streptococcus

The structure and stability of the 16-amino-acid-residue fragment [IG(46-61)] corresponding to the C-terminal beta-hairpin of the B3 domain of the immunoglobulin binding protein G from Streptococcus was investigated by means of CD and NMR spectroscopy and by differential scanning calorimetry. The CD and 2D NMR experiments were carried out (i) in water at different temperatures and (ii) at one temperature (305 K), with only CD, at different TFE concentrations. Our results show that the IG(46-61) peptide possesses organized three-dimensional structure at all investigated temperatures. The three-dimensional structure of the IG(46-61) peptide resembles the general shape of a beta-hairpin that is also observed for this peptide in the experimental structure of the B3 domain in the whole G protein; the structure is stabilized by hydrophobic interactions between nonpolar side chains. Our study shows that the melting temperature of the IG(46-61) peptide is about 320 K which supports the hypothesis that the investigated peptide can serve as a folding initiation site of the B3 domain of the immunoglobulin binding protein G.

An improved functional form for the temperature scaling factors of the components of the mesoscopic UNRES force field for simulations of protein structure and dynamics

Coarse-grained or mesoscopic models of proteins and the corresponding force fields are of great importance because they enable us to reduce the folding simulation time by several orders of magnitude compared to the all-atom approach and, consequently, reach the millisecond time scale of simulations. In the coarse-grained UNRES model for simulations of protein structure and dynamics, developed by our group, each amino acid residue is represented by a united side chain and a united peptide group located in the middle between the two neighboring alpha-carbon atoms, which assist only in the definition of the geometry. The prototype of the UNRES force field has been defined as a potential of mean force or restricted free-energy function corresponding to averaging out the degrees of freedom not present in the coarse-grained representation, which has further been approximated by a truncated Kubo cumulant series to enable us to derive analytical expressions for the corresponding terms. This force field should depend on temperature, and in its simplest form, a term corresponding to the cumulant of order n should be multiplied by f(n) = 1/T(n-1). The temperature dependence has been introduced in recent work ( J. Phys. Chem. B , 2007 , 111 , 260 - 285 ), and in order to prevent too steep a variation with temperature, the factors at the nth order cumulant terms were assumed to have a form f(n) = ln[exp(1) + exp(-1)]/ln{exp[(T/T(0))(n-1)] + exp[-(T/T(0))(n-1)]}, where T(0) = 300 K is the reference temperature. In this work, we have introduced a modified scaling factor f(n) = ln[exp(c) + exp(-c)]/ln{exp[c(T/T(0))(n-1)] + exp[-c(T/T(0))(n-1)]}, where c is an adjustable parameter, and determined c by fitting the analytical approximation of the temperature dependence of the virtual bond torsional term corresponding to rotation about the C(alpha)...C(alpha) virtual bond in terminally blocked dialanine to the respective potential of mean force calculated from the MP2/6-31G(d, p) ab initio energy surfaces of terminally blocked alanine (Ac-Ala-NHMe) and, independently, by optimizing it to obtain a sharp heat capacity curve and the lowest ensemble-averaged root-mean-square deviation over the C(alpha) atoms of 1GAB used as a training protein. Both approaches gave consistent results, and c = 1.4 has been selected as the optimal value of this parameter. The force field with the new temperature scaling factors has been optimized using 1GAB as the training protein. The new force field has been tested on a series of medium size alpha-helical proteins and found to perform better than that with the original temperature scaling factors.

Towards temperature-dependent coarse-grained potentials of side-chain interactions for protein folding simulations. I: molecular dynamics study of a pair of methane molecules in water at various temperatures

By means of molecular dynamics simulations of a pair of methane molecules in a TIP3P periodic water box with the NVT scheme at six temperatures and, additionally, the NPT scheme at three temperatures ranging from T = 283 to 373 K, we determined the potential of mean force (PMF) of pairs of interacting methane molecules in water as functions of distance between the methane molecules. The PMFs converge to a single baseline only for r> 11 A at all temperatures. The curves of the dimensionless PMF obtained at different temperatures with the NVT scheme overlap almost perfectly in the region of the contact minimum and still very well in the regions of the desolvation maximum and the solvent-separated minimum, which suggests that the temperature-dependent hydrophobic interaction potentials at constant volume in united-residue force fields can be obtained by scaling the respective dimensionless potentials by RT, R being the universal gas constant. For the dimensionless potentials of mean force obtained with the NPT scheme, the depth of the contact minimum increases, whereas the height of the desolvation maximum and the depth of the solvent-separated minimum decrease with temperature, in agreement with results reported in the literature.

How adequate are one- and two-dimensional free energy landscapes for protein folding dynamics?

The molecular dynamics trajectories of protein folding or unfolding, generated with the coarse-grained united-residue force field for the B domain of staphylococcal protein A, were analyzed by principal component analysis (PCA). The folding or unfolding process was examined by using free-energy landscapes (FELs) in PC space. By introducing a novel multidimensional FEL, it was shown that the low-dimensional FELs are not always sufficient for the description of folding or unfolding processes. Similarities between the topographies of FELs along low- and high-indexed principal components were observed.

Application of Multiplexed Replica Exchange Molecular Dynamics to the UNRES Force Field: Tests with alpha and alpha+beta Proteins

The replica exchange (RE) method is increasingly used to improve sampling in molecular dynamics (MD) simulations of biomolecular systems. Recently, we implemented the united-residue UNRES force field for mesoscopic MD. Initial results from UNRES MD simulations show that we are able to simulate folding events that take place in a microsecond or even a millisecond time scale. To speed up the search further, we applied the multiplexing replica exchange molecular dynamics (MREMD) method. The multiplexed variant (MREMD) of the RE method, developed by Rhee and Pande, differs from the original RE method in that several trajectories are run at a given temperature. Each set of trajectories run at a different temperature constitutes a layer. Exchanges are attempted not only within a single layer but also between layers. The code has been parallelized and scales up to 4000 processors. We present a comparison of canonical MD, REMD, and MREMD simulations of protein folding with the UNRES force-field. We demonstrate that the multiplexed procedure increases the power of replica exchange MD considerably and convergence of the thermodynamic quantities is achieved much faster.

Acidic-basic properties of three alanine-based peptides containing acidic and basic side chains: comparison between theory and experiment

The purpose of this work was to evaluate the effect of the nature of the ionizable end groups, and the solvent, on their acid-base properties in alanine-based peptides. Hence, the acid-base properties of three alanine-based peptides: Ac-KK-(A)(7)-KK-NH(2) (KAK), Ac-OO-(A)(7)-DD-NH(2) (OAD), Ac-KK-(A)(7)-EE-NH(2) (KAE), where A, D, E, K, and O denote alanine, aspartic acid, glutamic acid, lysine, and ornithine, respectively, were determined in water and in methanol by potentiometry. With the availability of these data, the ability of two theoretical methods to simulate pH-metric titration of those peptides was assessed: (i) the electrostatically driven Monte Carlo method with the ECEPP/3 force field and the Poisson-Boltzmann approach to compute solvation energy (EDMC/PB/pH), and (ii) the molecular dynamics method with the AMBER force field and the Generalized Born model (MD/GB/pH). For OAD and KAE, pK(a1) and pK(a2) correspond to the acidic side chains. For all three compounds in both solvents, the pK(a1) value is remarkably lower than the pK(a) of a compound modeling the respective isolated side chain, which can be explained by the influence of the electrostatic field from positively charged ornithine or lysine side chains. The experimental titration curves are reproduced well by the MD/GB/pH approach, the agreement being better if restraints derived from NMR measurements are incorporated in the conformational search. Poorer agreement is achieved by the EDMC/PB/pH method.

Conformational studies of the alpha-helical 28-43 fragment of the B3 domain of the immunoglobulin binding protein G from Streptococcus

To determine whether the alpha-helix in the B3 immunoglobulin binding domain of protein G from group G Streptococcus has conformational stability as an isolated fragment, we carried out a CD and NMR study of the 16-residue peptide in solution corresponding to this alpha-helix. Based on two-dimensional H-NMR spectra recorded at three different temperatures (283, 305, and 313 K), it was found that this peptide is mostly unstructured in water at these temperatures. Weak signals corresponding to i,i+3 or i,i+4 interactions, which are characteristic of formation of turn-like structures, were observed in the ROE spectra at all temperatures. The absence of a stable three-dimensional structure of the investigated peptide supports an earlier study (Blanco and Serrano, Eur J Biochem 1995, 230, 634-649) of a possible mechanism for folding of other (B1 and B2) immunoglobulin binding domains of Protein G.

Principal component analysis for protein folding dynamics

Protein folding is considered here by studying the dynamics of the folding of the triple beta-strand WW domain from the Formin-binding protein 28. Starting from the unfolded state and ending either in the native or nonnative conformational states, trajectories are generated with the coarse-grained united residue (UNRES) force field. The effectiveness of principal components analysis (PCA), an already established mathematical technique for finding global, correlated motions in atomic simulations of proteins, is evaluated here for coarse-grained trajectories. The problems related to PCA and their solutions are discussed. The folding and nonfolding of proteins are examined with free-energy landscapes. Detailed analyses of many folding and nonfolding trajectories at different temperatures show that PCA is very efficient for characterizing the general folding and nonfolding features of proteins. It is shown that the first principal component captures and describes in detail the dynamics of a system. Anomalous diffusion in the folding/nonfolding dynamics is examined by the mean-square displacement (MSD) and the fractional diffusion and fractional kinetic equations. The collisionless (or ballistic) behavior of a polypeptide undergoing Brownian motion along the first few principal components is accounted for.

Implementations of Nosé-Hoover and Nosé-Poincaré thermostats in mesoscopic dynamic simulations with the united-residue model of a polypeptide chain

Molecular dynamics (MD) simulations generate a canonical ensemble only when integration of the equations of motion is coupled to a thermostat. Three extended phase space thermostats, one version of Nose-Hoover and two versions of Nose-Poincare, are compared with each other and with the Berendsen thermostat and Langevin stochastic dynamics. Implementation of extended phase space thermostats was first tested on a model Lennard-Jones fluid system; subsequently, they were implemented with our physics-based protein united-residue (UNRES) force field MD. The thermostats were also implemented and tested for the multiple-time-step reversible reference system propagator (RESPA). The velocity and temperature distributions were analyzed to confirm that the proper canonical distribution is generated by each simulation. The value of the artificial mass constant, Q, of the thermostat has a large influence on the distribution of the temperatures sampled during UNRES simulations (the velocity distributions were affected only slightly). The numerical stabilities of all three algorithms were compared with each other and with that of microcanonical MD. Both Nose-Poincare thermostats, which are symplectic, were not very stable for both the Lennard-Jones fluid and UNRES MD simulations started from nonequilibrated structures which implies major changes of the potential energy throughout a trajectory. Even though the Nose-Hoover thermostat does not have a canonical symplectic structure, it is the most stable algorithm for UNRES MD simulations. For UNRES with RESPA, the "extended system inside-reference system propagator algorithm" of the RESPA implementation of the Nose-Hoover thermostat was the only stable algorithm, and enabled us to increase the integration time step.

Simple physics-based analytical formulas for the potentials of mean force for the interaction of amino acid side chains in water. IV. Pairs of different hydrophobic side chains

The potentials of mean force of 21 heterodimers of the molecules modeling hydrophobic amino acid side chains: ethane (for alanine), propane (for proline), isobutane (for valine), isopentane (for leucine and isoleucine), ethylbenzene (for phenylalanine), methyl propyl sulfide (for methionine), and indole (for tryptophane) were determined by umbrella-sampling molecular dynamics simulations in explicit water as functions of distance and orientation. Analytical expressions consisting of the Gay-Berne term to represent effective van der Waals interactions and the cavity term proposed in our earlier work were fitted to the potentials of mean force. The positions and depths of the contact minima and the positions and heights of the desolvation maxima, including their dependence on the orientation of the molecules, are well represented by the analytical expressions for all systems; large deviations between the MD-determined PMF and the analytical approximations are observed for pairs involving the least spheroidal solutes: ethylbenzene, indole, and methyl propyl sulfide at short distances at which the PMF is high and, consequently, these regions are rarely visited. When data from the PMF within only 10 kcal/mol above the global minimum are considered, the standard deviation between the MD-determined and the fitted PMF is from 0.25 to 0.55 kcal/mol (the relative standard deviation being from 4% to 8%); it is larger for pairs involving nonspherical solute molecules. The free energies of contact computed from the PMF surfaces are well correlated with those determined from protein-crystal data with a slope close to that relating the free energies of transfer of amino acids (from water to n-octanol) to the average contact free energies determined from protein-crystal data. These observations justify future use of the determined potentials in coarse-grained protein-folding simulations.