PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
J Comput Chem. Author manuscript; available in PMC 2010 October 1.
Published in final edited form as:
PMCID: PMC2760447
NIHMSID: NIHMS147838

Exploring the Parameter Space of the Coarse-Grained UNRES Force Field by Random Search: Selecting a Transferable Medium-Resolution Force Field

Abstract

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 α-helical and a minimal β-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 α or α + β 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.

Keywords: protein folding, thermodynamic hypothesis, coarse-grained models, UNRES force field, force-field parameterization

Introduction

Protein folding is one of the most interesting problems that attracts much attention in many scientific fields. Although many years have passed, predicting protein structure and protein folding pathways based on Anfinsen’s thermodynamic hypothesis1 is still a great challenge. Knowledge-based methods, i.e., homology modeling, threading, and fragment methods, usually perform better in protein-structure prediction compared to physics-based methods.2-4 However, knowledge-based conformational-space search methods cannot be used in studying protein folding processes. On the other hand, knowledge-based statistical potentials have been very successful in simulating folding pathways5,6; however, in such applications,5,6 they are used in physics-based canonical or multicanonical simulations. Coarse-grained7 and all-atom8 statistical potentials can also locate native-like structures as the lowest in energy in unrestricted Monte Carlo simulations. The most commonly used fully physics-based methods are all-atom molecular dynamics or Monte Carlo methods.9 Although there are some successful examples of studying folding of small proteins by using this methodology,10,11 all-atom molecular dynamics with explicit solvent cannot be used for studying folding of large proteins, because of the high computational cost arising from the need to adopt very small time steps to simulate long-time folding processes. Larger proteins can be treated by all-atom molecular dynamics,12 global-optimization methods13,14 with implicit-solvent models or all-atom canonical or generalized-ensemble Monte Carlo methods15,16 with implicit solvent; however, the maximum length of polypeptide chains is about 60 amino-acid residues. Despite size-limitations of the all-atom approach, important conclusions can be drawn from such studies, the most recent example being the unveiling of a new folding mechanism in a study of the C-terminal fragment of Top7.16

In the past decade, we have developed a physics-based united-residue force field (UNRES) for off-lattice protein structure simulations.17-35 Because of the reduction of the number of interaction sites and variables, UNRES can be used in simulating the folding of large proteins in real time. As opposed to most united-residue force fields, which are largely knowledge-based, UNRES was carefully derived based on physics of interactions using a cluster-cumulant expansion of effective free energy of protein plus the surrounding solvent.20 A hierarchical procedure was designed to parameterize the force field.26-28 Recently,33 temperature dependence has been included in UNRES as a consequence of its definition as a restricted free-energy function; this modification of the force field enabled us to reproduce the thermodynamic characteristics of the folding transition. In addition, we have recently replaced the statistical potentials that account for the energetics of virtual-bond-angle bending34 and the side-chain rotamers35 with ones parameterized based on the results of quantum-mechanical calculations of model energy surfaces. We also carried out initial work on reparameterization of the temperature-dependent UNRES force field.33 However, the optimized parameters were good for α but not for β-proteins. Later,36 we included one α-(engralied homeodomain; PDB code: 1ENH) and one β-protein (FBP28 WW domain; PDB code: 1E0L) in optimization simultaneously and were able to reproduce the experimental structure and free-energy profiles of these proteins; however, the force field had only limited transferability. One of the reason for this could be the fact that we started force-field-parameter optimization from a single set of parameters, which results in searching only a limited range of energy-parameter space.

The aim of the work reported here was to explore the energy-parameter space to find good initial energy-parameter sets for optimization of the new UNRES force field. The primary motivation for this was our recent replacement of knowledge-based17 virtual-bond-angle bending and side-chain rotamer potentials with ones derived from quantum-mechanical calculations.34,35 Although these potentials control only very local properties of a polypeptide chain, their radical change led to loss of balance of the other UNRES energy terms, which are responsible for the overall fold of a polypeptide chain. Consequently, the energy-term weights found in our earlier studies27,28,33 are no longer good as initial approximations. We, therefore, carried out an extensive random search of energy-term weights and tested the resulting weights by performing replica-exchange (REMD) simulations15,37,38 on two peptides, namely the tryptophan cage39 (PDB code: 1L2Y) and the tryptophan zipper 240 (PDB code: 1LE1), which possess a minimum stable α-helix and β-hairpin fold, respectively. Having tested the parameter sets with these two peptides, we selected the parameter sets which produced native-like structures for both of them and used each of these sets to run multiplexed replica-exchange molecular dynamics (MREMD) simulations41 on two larger proteins 1ENH42 and 1E0L43 and selected the sets that found the native-like structures of these two proteins. The final parameter sets were tested with nine other proteins.

Methods

UNRES Force Field

In the UNRES model,17-35 a sequence of α-carbon atoms connected by virtual bonds with attached side chains is used to represent a polypeptide chain. Each amino-acid residue has two interaction sites, one located in the middle between two consecutive Cα atoms and the other one is in the mass center of the corresponding side chain. The Cα atoms serve only to define the geometry, as shown in Figure 1. The UNRES force field has been carefully derived as a restricted free energy (RFE) function of an all-atom polypeptide chain plus the surrounding solvent.19,20 The all-atom energy is averaged over the degrees of freedom which are lost when passing from the all-atom to the simplified representation. The RFE is further decomposed into factors representing the interactions within and between a given number of united interaction sites. Expansion of these factors into the Kubo generalized cumulant series44 enabled us to obtain approximate analytical expressions for the respective terms, including the multibody terms, which are derived in other force fields from structural databases or on a heuristic basis. The details of the theoretical basis have been described in our previous work.19,20 The energy function of the virtual-bond chain is expressed by eq. (1).

U(X,T)=wSCi<jUSCiSCj(X)+wSCpijUSCipj(X)+wppVDWi<j1UpipjVDW(X)+wppelf2(T)i<j1Upipjel(X)+wtorf2(T)iUtor(γi)+wtordf3(T)iUtord(γi,γi+1)+wbiUb(θi)+wrotiUrot(αSCi,βSCi)+wcorr(3)f3(T)Ucorr(3)(X)+wcorr(4)f4(T)Ucorr(4)(X)+wcorr(5)f5(T)Ucorr(5)(X)+wcorr(6)f6(T)Ucorr(6)(X)+wturn(3)f3(T)Uturn(3)(X)+wturn(4)f4(T)Uturn(4)(X)+wturn(6)f6(T)Uturn(6)(X)+wbondi=1nbondUbond(di)
(1)

where X denotes the set of coordinates of a conformation of the coarse-grained chain and T is the absolute temperature. The terms USCiSCj (X) are the potentials of mean force of the interaction of isolated side chains in water (with the contribution from the solvent already included). The terms USCipj (X) are the excluded-volume potentials of interactions between side chains and peptide groups. The terms UpipjVDW(X) are the Lennard-Jones potentials and Upipjel(X) are the averaged electrostatic-interaction potentials between peptide groups. The terms Utor(γi) and Utord(γi, γi+1) are the virtual-bond torsional and double-torsional potentials, respectively. The terms Ub(θi), Urot(αsci, βscj), and Ubond(di) denote the energy of virtual-bond angle bending, side-chain rotamers, and virtual-bond stretching, respectively. The terms Ucorr(m) represent the correlation or multibody contribution from the backbone-local and backbone-electrostatic interactions. The terms Uturn(m) denote the correlation contribution involving m consecutive peptide groups. The energy-term weights corresponding to second- and higher-order generalized cumulants are multiplied by the appropriate scaling factor, fn(T), which is defined by eq. (2), where n is the order of the cumulant.

fn(T)=ln[exp(1)+exp(1)]ln{exp[(TT0)n1]+exp[(TT0)n1]}
(2)

where T0 is the arbitrary reference temperature (we set T0 = 300 K), and T is the current temperature.

Figure 1
The UNRES model of the polypeptide chain. The united peptide groups (p) are shown as dark circles, and open circles represent the Cα atoms which serve only as geometric points. The side chains are represented by ellipsoids, SC’s; each ...

All parameters except the energy-term weights which were being determined in this work were taken from our earlier work.17,19,24,25,34,35 The parameters of the USCiSCj potentials [which have the Gay-Berne functional form45] were taken from ref. 17, where they were determined from protein-crystal data. The torsional (Utor) and double-torsional (Utord) potentials have been determined from the RFE surfaces of model system calculated by ab initio molecular quantum mechanics in ref. 24. Likewise, the peptide-group interaction (Upipjel and UpipjVDW) and correlation (Ucorr(m) and Uturn(m)) terms were calculated in ref. 25; the parameters of the correlation terms were further refined in ref. 27 by optimization of the energy landscape for selected training proteins. The virtual-bond-angle bending (Ub) terms were calculated using ab initio molecular quantum mechanics in ref. 34. The virtual-bond rotamer (Urot) and virtual-bond stretching (Ud) terms were calculated in ref. 35 from the energy surfaces of terminally blocked amino-acid residues calculated with the semiempirical AM1 method of molecular quantum mechanics. The USCipj terms were taken directly from ref. 19.

Search of Energy Parameter Space

We chose the energy-term weights [the ws of eq. (1)], because the force field is most sensitive to these parameters and they can be determined only by optimizing the force field on selected training proteins.21-23,28 In our previous work,21-23,28 we optimized the force field using old force field parameters as a starting point. This could result in too narrow a range of the search of the energy-parameter space and, consequently, reduced the performance of the force field. Because we recently modified the force field significantly by introducing temperature dependence33 and the new physics-based Ub34 and Urot35 potentials, we considered it worthwhile to perform a more extensive search of parameter space than optimization, started from a given set of parameters, can provide. The procedure was as follows:

  1. Generate a number of different sets of energy-term weights at random; in this work we generated 100 sets, each weight being selected from a uniform distribution from wxlowwxwxup, where wx denotes the weight of energy term x and the superscripts “low” and “up” mark the lower and upper boundary, respectively. There are 11 variable weights, which prohibits systematic grid search of the weight space.
  2. With each set, carry out REMD or MREMD runs on a peptide with a minimal α-helix fold (we chose 1L2Y) and a protein with a minimal β-sheet fold (we chose 1LE1 whose structure is a single β-hairpin).
  3. Select those sets of weights which result in the appearance of some native-like structures of the small test proteins in step 2.
  4. Test the sets of weights selected in step 3 by running REMD or MREMD on larger proteins containing α- and β-folds (we chose 1ENH and 1E0L, respectively, in this work).
  5. Select those sets of weights which result in the appearance of some native conformations of the test proteins in step 4.
  6. Test the sets of weights selected in step 5 on a larger set of proteins.

Following our earlier work,28 we considered only the correlation terms up to the fourth order and, consequently, set the weights of all correlation terms of order higher than 4 at zero. All other weights except that of the side-chain-interaction potentials, [USCiSCj; eq. (1)] which was fixed at wSC = 1 were variable. The value of wSC was fixed, because the USCiSCj potentials had been scaled in our earlier work17 to reproduce the experimental free energies of transfer of terminally blocked amino acids from water to n-octanol determined by Fauchere and Pliska;46 consequently, they have strict connection to the energetics of side-chain interaction. The boundaries of the weights are listed in Table 1. These boundaries were chosen based on the experience from our previous work21-23,28,33 in which the correct geometry of regular secondary-structure elements as well as reasonable ranges of UNRES energies were produced.

Table 1
Boundaries of Energy-Term Weights

Structure Comparison

We used the α-carbon root-mean-square deviation (Cα RMSD) to assess the overall similarity of the calculated and experimental structures. To assess the native-likeness of the two minimal-fold peptides: 1L2Y and 1LE1 and of fragments of larger proteins, we used the fraction of native contacts between the peptide groups as defined in ref. 27 and the quantity q defined by eq. (3) and introduced in our previous work,33 which is a measure of compatibility between the residue-residue contacts in a given conformation to those in the experimental conformation.

qi=11Ni{k{i}l{i}lkmexp[(dCkαClαdCkαClαnat)24(dCkαClαnat)2]+exp[(dSCkSCldSCkSClnat)24(dSCkSClnat)2]}
(3)

where {i} is the set of residues constituting fragment i (we note that the fragment does not have to be continuous; e.g., we can consider a β-sheet composed of two nonsequential strands as a fragment), Ni is the total number of distances (both between the Cα atoms and between the SC pseudoatoms) in fragment i, m is the minimum separation of residues in the amino-acid sequence at which the distances are calculated (we set m = 3), dCkαClα or dSCkSCl denotes the distance between the α-carbon atoms or side-chain centers of residues k and l, respectively, of the conformation under consideration, and the same quantities with the nat superscript denote the distances in the experimental structure. q = 0 if the conformation considered is the experimental structure and q > 0.5 indicates that the structure considered and the experimental structure are grossly dissimilar.

Simulation Details

All REMD and MREMD simulations were carried out at 20 temperatures T = 250, 260, 270, 280, 290, 300, 310, 320, 330, 340, 350, 360, 370, 380, 390, 400, 420, 440, 460, and 500 K. The frequency of replica exchange and the frequency of snapshot dumping was 20,000 MD steps. The Berendsen thermostat,47 implemented in UNRES/MD in our previous work,31 was used to maintain constant temperature; the coupling constant was τ = 0.489 fs, which was assessed in our previous work31 to provide a compromise between the extent of fluctuations of the potential and kinetic energy. The time increment in integrating the equations of motion was δt = 4.89 fs. For the two minimal-fold peptides: 1L2Y and 1LE1, one trajectory per temperature was run (a total of 20 trajectories; a REMD simulation) with the total of 5,000,000 MD steps per trajectory. For the other 11 systems 2 trajectories per temperature were run (a total of 40 trajectories; a MREMD simulation) with the total of 10,000,000 steps per trajectory. Thus, the total time per trajectory was 24.45 and 48.9 ns, respectively; it should be noted, however, that the UNRES time scale is stretched compared to the all-atom time scale because of averaging of the fast-varying degrees of freedom, especially those of the solvent.31

The weighted-histogram analysis method (WHAM),48 modified in our recent work33 to treat simulations with the UNRES force field, was implemented in the computation of heat-capacity as well as ensemble-averaged RMSD and q. We analyzed snapshots from the last 1,000,000 MD steps for 1L2Y and 1LE1 and from the last 2,000,0000 MD steps for the other systems. To dissect the conformational ensembles, we applied the procedure described in our recent work,33 in which a cluster analysis is carried out for the selected batch of conformations, and the probabilities of the clusters are computed by summing up the contributions from the constituent conformations at a selected temperature calculated with WHAM. The clusters are ranked following decreasing probability. The single-link algorithm49,50 was applied for clustering.

Results and Discussion

Initial Selection of Weight Sets by Tests with Minimal-Fold Peptides

As stated in the Methods section, we first randomly generated 100 sets of energy-term weights, each within the boundaries specified in Table 1. The sets of conformations of 1L2Y and 1LE1 obtained in the last 1,000,000 MD steps of simulations with each of the 100 sets of energy-term weights were analyzed for the presence of native-like conformations. For 1L2Y, we required that the Cα RMSD be < 3.0 Å and qAsn2-Leu10 < 0.35, which implies that the α-helical segment from residues Asn2 to Leu10 had a contact pattern characteristic of an α-helix. For 1LE1 we required that the Cα RMSD be < 2.5 Å and a contact pattern between peptide groups characteristic of a β-hairpin with a tolerance of ±3 residue shift with respect to the native β-hairpin. The native-like structures of both peptides were found in simulations carried out with eight sets of energy-term weights; these weights are summarized in Table 2. A summary of the results: the lowest RMSD values, the lowest ensemble-averaged RMSD values [computed by using the WHAM method as stated in ref. 33], and the ensemble-averaged q values at temperatures corresponding to those of the lowest-ensemble-averaged RMSD obtained with these eight sets of weights is presented in Table 3.

Table 2
Sets of Energy-Term Weights that Produced Native-Like Structures of Both 1L2Y and 1LE1 In MREMD Simulations
Table 3
Results of MREMD Simulation of 1L2Y and 1LE1 Obtained with Energy-Term Weights of Table 2

It can be seen that sets 1, 3, and 6 produced not only some native-like structures but the ensemble-averaged RMSD and q values also are low for both 1L2Y and 1LE1, which means that the folding transition occurs for structures which are close to the experimental structures of these two peptides. Figures Figures22 and and33 present graphs of heat capacity (Cv) and ensemble-averaged RMSD vs. temperature as well as the experimental structures and representative structures at T = 280 K (i.e., lower than the folding-transition temperature). It can be seen that the heat capacity curve of 1L2Y has a peak at T = 317 K (Fig. 2A), which corresponds to the inflection point of the ensemble-averaged RMSD curve (Fig. 2B). For 1LE1, the heat-capacity peak and the inflection point of the ensemble-averaged RMSD curve occur at T = 297 K. The experimental melting temperatures of these two peptides are Tm = 315 K39 and Tm = 323 K,40 respectively, which means that the calculated melting temperature of 1LE1 is significantly lower than the experimental melting temperature but the aim of our study was to find reasonable starting sets of energy-term weights and not to obtain exact correspondence between the calculated and experimental folding characteristics. Nevertheless, it is remarkable that the calculated melting temperatures, obtained not only with energy-term-weight set 1 but also with the remaining sets (Table 3), are within physiological range, a feature that occurred in our previous studies only when we explicitly requested that folding transition occur at a specified temperature.33,36 The force fields that we derived earlier to locate the native-like structures as the lowest minima of the UNRES energy function produced melting temperatures of the order of Tm = 1000 K or higher.38 That the melting temperatures are reasonable in our current study seems to result from our freezing the weight wSC at 1 and with the parameters of the USCiSCj potentials at the values determined in our previous work;17 although these parameters have been determined from protein-crystal data, they were subsequently scaled to reproduce the free energies of transfer of terminally blocked amino acids from water to n-octanol determined by Fauchere and Pliska.46 Therefore, the side-chain-interaction potentials are in strict correspondence with the thermodynamics of side chain-side chain interactions. The weights of the other UNRES energy terms had to accommodate to them in order to produce native-like structures of both peptides, and the thermodynamic parameters of folding are, therefore, in the range of the experimental ones even though no melting-temperature information was supplied as input to the procedure. In our previous optimizations, which did not consider thermodynamic data,27,28 the weight wSC and the parameters of the USCiSCj potentials were varied and the connection to thermodynamics was, therefore, lost.

Figure 2
Plots of (A) heat capacity, (B) ensemble-averaged RMSD from the native structure obtained in the MREMD simulation of 1L2Y with set 1 of energy-term weights, (C) Cα traces of the experimental structure, and (D) the representative structure of the ...
Figure 3
Plots of (A) heat capacity, (B) ensemble-averaged RMSD from the native structure obtained in the MREMD simulation of 1L2Y with set 2 of energy-term weights, (C) Cα traces of the experimental structure, and (D) the representative structure of the ...

Selection of Sets of Energy-Term Weights by Tests with 1ENH and 1E0L

The eight sets of energy-term weights were subsequently tested on two larger proteins: 1ENH (a three-helix bundle) and 1E0L (a three-stranded antiparallel β-sheet). Sets 1 and 2 produced native-like structures of both proteins; these results are collected in Table 4 and in Figure 4D and K. In Table 4 the probabilities and average RMSD values of the clusters containing native-like structures (detected based on topology rather than low ensemble-averaged RMSD values) are listed in addition to lowest RMSD values found for the conformations corresponding to the last 2,000,000 MD steps. However, as opposed to 1L2Y and 1LE1, the ensemble-averaged RMSD values were high at all temperatures below the folding-transition temperature. For the 1E0L run with set 1 of energy-term weights, the high ensemble-averaged RMSD value was caused by misalignment of the N- and the C-terminal fragments that cap the three-stranded β-sheet in the experimental structure, while conformations with the correct three-stranded antiparallel β-sheet form the dominant cluster. For the 1ENH run with both energy-term-weight sets, and for the 1E0L run with set 2, the clusters containing native-like structures have low probabilities (Table 4).

Figure 4
Most native-like structures obtained in MREMD simulations (left) and the experimental structures (right) of the proteins of the test set obtained with energy-term weights set 1. (A) 2HEP, (B) 1BDD, (C) 1GAB, (D) 1ENH, (E) 1LQ7, (F) 1POU, (G) 1CLB, (H) ...
Table 4
MREMD Results Obtained with Energy-Term-Weight Set 1 and 2 for all Test Proteins

The two sets of energy-term weights were subsequently tested on nine other proteins: 2HEP (an α-helical hairpin), 1BDD (a three-helix bundle), 1GAB (a three-helix bundle), 1LQ7 (a three-helix bundle), 1POU (α-helical structure), 1CLB (α-helical structure), 1FSD (α + β structure), 1E0G (α + β structure), and 1PGA (α + β structure), The results are also summarized in Table 4 and the best native-like structures obtained in the last 2,000,000 MD steps are shown in Figure 4. As shown in Table 4 and Figure 4, parameter set 1 produced medium-resolution native-like structures for most of these proteins except 1PGA, for which a mirror image of the native fold was obtained. Nevertheless, sets 1 and 2 of the energy-term weights produced quite transferable force fields with which to obtain medium-resolution structures. Set 1 appears better than set 2 as far as the resolution is concerned.

Conclusions

In this work, we have tested the performance of randomly generated sets of energy-term weights of the UNRES force field. Eight out of 100 sets were able to produce both α and β structures as assessed by tests with minimum-fold peptides: the tryptophan cage and the tryptophan zipper. Consequently, the space of reasonable UNRES parameters is not as restricted as it would appear from our previous studies33,36 in which we started force-field optimization from a given point in parameter space. Two of the energy-term-weight sets produced medium-resolution native-like structures of nine proteins with different secondary structure and size from 28 to 75 residues. Our results suggest that random search of UNRES force field parameter space is a good method to generate starting points for detailed optimization according to the procedure described in our previous work.33 More minimal-fold peptides or small proteins could be used as filters to select good parameters. This research is now being carried out in our laboratory.

Acknowledgments

This research was conducted by using the resources of (a) our 800-processor Beowulf cluster at the Baker Laboratory of Chemistry and Chemical Biology, Cornell University, (b) the National Science Foundation Terascale Computing System at the Pittsburgh Super-computer Center, (c) the John von Neumann Institute for Computing at the Central Institute for Applied Mathematics, Forschungszentrum Jülich, Germany, (d) the Beowulf cluster at the Department of Computer Science, Cornell University, (e) the resources of the Center for Computation and Technology at Louisiana State University, which is supported by funding from the Louisiana legislature, (f) our 45-processor Beowulf cluster at the Faculty of Chemistry, University of Gdańsk, (g) the Informatics Center of the Metropolitan Academic Network (IC MAN) in Gdańsk, and (h) the Interdisciplinary Center of Mathematical and Computer Modeling (ICM) at the University of Warsaw.

Contract/grant sponsor: National Science Foundation of China (NSCF); contract/grant number: 30525037

Contract/grant sponsor: Polish Ministry of Science and Higher Education; contract/grant numbers: BW-R300-5-0001-8, DS/8372-4-0-138-8

Contract/grant sponsor: U.S. National Institutes of Health; contract/grant number: GM-14312

Contract/grant sponsor: U.S. National Science Foundation; contract/grant number: MCB05-41633

References

1. Anfinsen CB. Science. 1973;181:223. [PubMed]
2. Bradley P, Chivian D, Meiler J, Misura KMS, Rohl CA, Schief WR, Wedemeyer WJ, Schueler-Furman O, Murphy P, Schonbrun J, Strauss CEM, Baker D. Proteins Struct Funct Genet. 2003;53:457. [PubMed]
3. Skolnick J, Zhang Y, Arakaki AK, Koliński A, Boniecki M, Szilagyi A, Kihara D. Proteins Struct Funct Genet. 2003;53:469. [PubMed]
4. Scheraga HA, Liwo A, Ołdziej S, Czaplewski C, Pillardy J, Ripoll DR, Vila JA, Kaźmierkiewicz R, Arnautova YA, Jagielska A, Chinchio M, Nanias M. Front Biosci. 2004;9:3296. [PubMed]
5. Kmiecik S, Kolinski A. Proc Natl Acad Sci USA. 2007;104:12330. [PubMed]
6. Kmiecik S, Kolinski A. Biophys J. 2008;94:726. [PubMed]
7. Kolinski A, Godzik A, Skolnick J. J Chem Phys. 1993;98:7420.
8. Yang JS, Chen WW, Skolnick J, Shakhnovich EI. Structure. 2007;115:53. [PubMed]
9. Zimmermann O, Hansmann UHE. Biochim Biophys Acta. 2008;1784:252. [PMC free article] [PubMed]
10. Duan Y, Kollman PA. Science. 1998;282:740. [PubMed]
11. Simmerling C, Strockbine B, Roitberg AE. J Am Chem Soc. 2002;124:11258. [PubMed]
12. Jang S, Kim E, Shin S, Pak Y. J Am Chem Soc. 2003;125:14841. [PubMed]
13. Vila JA, Ripoll DR, Scheraga HA. Proc Natl Acad Sci USA. 2003;100:14812. [PubMed]
14. Schug A, Wenzel W. Biophys J. 2006;90:4273. [PubMed]
15. Hansmann UHE. Chem Phys Lett. 1997;281:140.
16. Mohanty S, Meinke JH, Zimmermann O, Hansmann UHE. Proc Natl Acad Sci USA. 2008;105:8004. [PubMed]
17. Liwo A, Ołdziej S, Pincus MR, Wawak RJ, Rackovsky S, Scheraga HA. J Comput Chem. 1997;18:849.
18. Liwo A, Pincus MR, Wawak RJ, Rackovsky S, Ołdziej S, Scheraga HA. J Comput Chem. 1997;18:874.
19. Liwo A, Kaźmierkiewicz R, Czaplewski C, Groth M, Ołdziej S, Rackovsky RJ, Wawakand S, Pincus MR, Scheraga HA. J Comput Chem. 1998;19:259.
20. Liwo A, Czaplewski C, Pillardy J, Scheraga HA. J Chem Phys. 2001;115:2323.
21. Lee J, Ripoll DR, Czaplewski C, Pillardy J, Wedemeyer WJ, Scheraga HA. J Phys Chem B. 2001;105:7291.
22. Pillardy J, Czaplewski C, Liwo A, Wedemeyer WJ, Lee J, Ripoll DR, Arłukowicz P, Ołdziej S, Arnautova YA, Scheraga HA. J Phys Chem B. 2001;105:7299.
23. Liwo A, Arłukowicz P, Czaplewski C, Ołdziej S, Pillardy J, Scheraga HA. Proc Natl Acad Sci USA. 2002;99:1937. [PubMed]
24. Ołdziej S, Kozłowska U, Liwo A, Scheraga HA. J Phys Chem A. 2003;107:8035.
25. Liwo A, Ołdziej S, Czaplewski C, Kozłowska U, Scheraga HA. J Phys Chem B. 2004;108:9421.
26. Liwo A, Arłukowicz P, Ołdziej S, Czaplewski C, Makowski M, Scheraga HA. J Phys Chem B. 2004;108:16918.
27. Ołdziej S, Liwo A, Czaplewski C, Pillardy J, Scheraga HA. J Phys Chem B. 2004;108:16934.
28. Ołdziej S, Lagiewka J, Liwo A, Czaplewski C, Chinchio M, Nanias M, Scheraga HA. J Phys Chem B. 2004;108:16950.
29. Liwo A, Khalili M, Scheraga HA. Proc Natl Acad Sci USA. 2005;102:2362. [PubMed]
30. Khalili M, Liwo A, Rakowski F, Grochowski P, Scheraga HA. J Phys Chem B. 2005;109:13785. [PMC free article] [PubMed]
31. Khalili M, Liwo A, Jagielska A, Scheraga HA. J Phys Chem B. 2005;109:13798. [PMC free article] [PubMed]
32. Khalili M, Liwo A, Scheraga HA. J Mol Biol. 2006;355:536. [PubMed]
33. Liwo A, Khalili M, Czaplewski C, Kalinowski S, Ołdziej S, Wachucik K, Scheraga HA. J Phys Chem B. 2007;111:260. [PMC free article] [PubMed]
34. Kozłowska U, Liwo A, Scheraga HA. J. Phys Conden Matter. 2007;19:285203.
35. Kozłowska U, Maisuradze G, Liwo A, Scheraga HA. Protein Sci. 2007;16(Suppl 1):92. [PubMed]
36. Liwo A, Czaplewski C, Ołdziej S, Kozłowska U, Makowski M, Kalinowski S, Kaźmierkiewicz R, Shen H, Maisuradze G, Scheraga HA, John von Neumann Institute for Computing (NIC) In: Münster G, Wolf D, Kremer M, editors. NIC Series, NIC Symposium 2008; Jülich, Germany. 20-21 February 2008; 2008. 39. 63.
37. Sugita Y, Okamoto Y. Chem Phys Lett. 1999;314:141.
38. Nanias M, Czaplewski C, Scheraga HA. J Chem Theor Comput. 2006;2:513. [PMC free article] [PubMed]
39. Neidigh JW, Fesinmeyer RM, Andersen NH. Nat Struct Biol. 2002;9:425. [PubMed]
40. Cochran AG, Skelton NJ, Starovasnik MA. Proc Natl Acad Sci USA. 2001;98:5578. [PubMed]
41. Rhee YM, Pande VS. Biophys J. 2003;84:775. [PubMed]
42. Macias MJ, Gervais V, Civera C, Oschkinat H. Nat Struct Biol. 2000;7:375. [PubMed]
43. Clarke ND, Kissinger CR, Desjarlais J, Gilliland GL, Pabo CO. Protein Sci. 1994;3:1779. [PubMed]
44. Kubo RJ. Phys Soc Jpn. 1962;17:1100.
45. Gay JG, Berne BJ. J Chem Phys. 1981;74:3316.
46. Fauchere J-L, Pliška V. Eur J Med Chem. 1983;18:369.
47. Berendsen HJC, Postma JPM, van Gunsteren WF, DiNola A, Haak JR. J Chem Phys. 1984;81:3684.
48. Kumar S, Bouzida D, Swendsen RH, Kollman PA, Osenberg JMR. J Comput Chem. 1992;13:1011.
49. Murtagh F. Multidimensional Clustering Algorithms. Physica-Verlag; Vienna: 1985.
50. Murtagh F, Heck A. Multivariate Data Analysis. Kluwer Academic; Dordrecht: 1987.