Free energy landscape of protein folding in water: explicit vs. implicit solvent

The Generalized Born (GB) continuum solvent model is arguably the most widely used implicit solvent model in protein folding and protein structure prediction simulations; however, it still remains an open question on how well the model behaves in these large-scale simulations. The current study uses the beta-hairpin from C-terminus of protein G as an example to explore the folding free energy landscape with various GB models, and the results are compared to the explicit solvent simulations and experiments. All free energy landscapes are obtained from extensive conformation space sampling with a highly parallel replica exchange method. Because solvation model parameters are strongly coupled with force fields, five different force field/solvation model combinations are examined and compared in this study, namely the explicit solvent model: OPLSAA/SPC model, and the implicit solvent models: OPLSAA/SGB (Surface GB), AMBER94/GBSA (GB with Solvent Accessible Surface Area), AMBER96/GBSA, and AMBER99/GBSA. Surprisingly, we find that the free energy landscapes from implicit solvent models are quite different from that of the explicit solvent model. Except for AMBER96/GBSA, all other implicit solvent models find the lowest free energy state not the native state. All implicit solvent models show erroneous salt-bridge effects between charged residues, particularly in OPLSAA/SGB model, where the overly strong salt-bridge effect results in an overweighting of a non-native structure with one hydrophobic residue F52 expelled from the hydrophobic core in order to make better salt bridges. On the other hand, both AMBER94/GBSA and AMBER99/GBSA models turn the beta-hairpin in to an alpha-helix, and the alpha-helical content is much higher than the previously reported alpha-helices in an explicit solvent simulation with AMBER94 (AMBER94/TIP3P). Only AMBER96/GBSA shows a reasonable free energy landscape with the lowest free energy structure the native one despite an erroneous salt-bridge between D47 and K50. Detailed results on free energy contour maps, lowest free energy structures, distribution of native contacts, alpha-helical content during the folding process, NOE comparison with NMR, and temperature dependences are reported and discussed for all five models.     

Ab initio construction of polypeptide fragments: efficient generation of accurate,representative ensembles

We describe a novel method to generate ensembles of conformations of the main-chain atoms [N, C(alpha), C, O, Cbeta] for a sequence of amino acids within the context of a fixed protein framework. Each conformation satisfies fundamental stereo-chemical restraints such as idealized geometry, favorable phi/psi angles, and excluded volume. The ensembles include conformations both near and far from the native structure. Algorithms for effective conformational sampling and constant time overlap detection permit the generation of thousands of distinct conformations in minutes. Unlike previous approaches, our method samples dihedral angles from fine-grained phi/psi state sets, which we demonstrate is superior to exhaustive enumeration from coarse phi/psi sets. Applied to a large set of loop structures, our method samples consistently near-native conformations, averaging 0.4, 1.1, and 2.2 A main-chain root-mean-square deviations for four, eight, and twelve residue long loops, respectively. The ensembles make ideal decoy sets to assess the discriminatory power of a selection method. Using these decoy sets, we conclude that quality of anchor geometry cannot reliably identify near-native conformations, though the selection results are comparable to previous loop prediction methods. In a subsequent study (de Bakker et al.: Proteins 2003;51:21-40), we demonstrate that the AMBER forcefield with the Generalized Born solvation model identifies near-native conformations significantly better than previous methods.     

Structural prediction of peptides bound to MHC class I

An ab initio structure prediction approach adapted to the peptide-major histocompatibility complex (MHC) class I system is presented. Based on structure comparisons of a large set of peptide-MHC class I complexes, a molecular dynamics protocol is proposed using simulated annealing (SA) cycles to sample the conformational space of the peptide in its fixed MHC environment. A set of 14 peptide-human leukocyte antigen (HLA) A0201 and 27 peptide-non-HLA A0201 complexes for which X-ray structures are available is used to test the accuracy of the prediction method. For each complex, 1000 peptide conformers are obtained from the SA sampling. A graph theory clustering algorithm based on heavy atom root-mean-square deviation (RMSD) values is applied to the sampled conformers. The clusters are ranked using cluster size, mean effective or conformational free energies, with solvation free energies computed using Generalized Born MV 2 (GB-MV2) and Poisson-Boltzmann (PB) continuum models. The final conformation is chosen as the center of the best-ranked cluster. With conformational free energies, the overall prediction success is 83% using a 1.00 Angstroms crystal RMSD criterion for main-chain atoms, and 76% using a 1.50 Angstroms RMSD criterion for heavy atoms. The prediction success is even higher for the set of 14 peptide-HLA A0201 complexes: 100% of the peptides have main-chain RMSD values < or =1.00 Angstroms and 93% of the peptides have heavy atom RMSD values < or =1.50 Angstroms. This structure prediction method can be applied to complexes of natural or modified antigenic peptides in their MHC environment with the aim to perform rational structure-based optimizations of tumor vaccines.     

Interplay of charge distribution and conformation in peptides: comparison of theory and experiment

We assessed the correlation between charge distribution and conformation of flexible peptides by comparing the theoretically calculated potentiometric-titration curves of two model peptides, Ac-Lys5-NHMe (a model of poly-L-lysine) and Ac-Lys-Ala11-Lys-Gly2-Tyr-NH2 (P1) in water and methanol, with the experimental curves. The calculation procedure consisted of three steps: (i) global conformational search of the peptide under study using the electrostatically driven Monte Carlo (EDMC) method with the empirical conformational energy program for peptides (ECEPP)/3 force field plus the surface-hydration (SRFOPT) or the generalized Born surface area (GBSA) solvation model as well as a molecular dynamics method with the assisted model building and energy refinement (AMBER)99/GBSA force field; (ii) reevaluation of the energy in the pH range considered by using the modified Poisson-Boltzmann approach and taking into account all possible protonation microstates of each conformation, and (iii) calculation of the average degree of protonation of the peptide at a given pH value by Boltzmann averaging over conformations. For Ac-Lys5-NHMe, the computed titration curve agrees qualitatively with the experimental curve of poly-L-lysine in 95% methanol. The experimental titration curves of peptide P1 in water and methanol indicate a remarkable downshift of the first pK(a) value compared to the values for reference compounds (n-butylamine and phenol, respectively), suggesting the presence of a hydrogen bond between the tyrosine hydroxyl oxygen and the H(epsilon) proton of a protonated lysine side chain. The theoretical titration curves agree well with the experimental curves, if conformations with such hydrogen bonds constitute a significant part of the ensemble; otherwise, the theory predicts too small a downward pH shift.     

Peptoid conformational free energy landscapes from implicit-solvent molecular simulations in AMBER

To test the accuracy of existing AMBER force field models in predicting peptoid conformation and dynamics, we simulated a set of model peptoid molecules recently examined by Butterfoss et al. (JACS 2009, 131, 16798-16807) using QM methods as well as three peptoid sequences with experimentally determined structures. We found that AMBER force fields, when used with a Generalized Born/Surface Area (GBSA) implicit solvation model, could accurately reproduce the peptoid torsional landscape as well as the major conformers of known peptoid structures. Enhanced sampling by replica exchange molecular dynamics (REMD) using temperatures from 300 to 800 K was used to sample over cis-trans isomerization barriers. Compared to (Nrch)5 and cyclo-octasarcosyl, the free energy of N-(2-nitro-3-hydroxyl phenyl)glycine-N-(phenyl)glycine has the most "foldable" free energy landscape, due to deep trans-amide minima dictated by N-aryl sidechains. For peptoids with (S)-N (1-phenylethyl) (Nspe) side chains, we observe a discrepancy in backbone dihedral propensities between molecular simulations and QM calculations, which may be due to force field effects or the inability to capture n --> n* interactions. For these residues, an empirical phi-angle biasing potential can "rescue" the backbone propensities seen in QM. This approach can serve as a general strategy for addressing force fields without resorting to a complete reparameterization. Overall, this study demonstrates the utility of implicit-solvent REMD simulations for efficient sampling to predict peptoid conformational landscapes, providing a potential tool for first-principles design of sequences with specific folding properties.     

Accurate and efficient loop selections by the DFIRE-based all-atom statistical potential

The conformations of loops are determined by the water-mediated interactions between amino acid residues. Energy functions that describe the interactions can be derived either from physical principles (physical-based energy function) or statistical analysis of known protein structures (knowledge-based statistical potentials). It is commonly believed that statistical potentials are appropriate for coarse-grained representation of proteins but are not as accurate as physical-based potentials when atomic resolution is required. Several recent applications of physical-based energy functions to loop selections appear to support this view. In this article, we apply a recently developed DFIRE-based statistical potential to three different loop decoy sets (RAPPER, Jacobson, and Forrest-Woolf sets). Together with a rotamer library for side-chain optimization, the performance of DFIRE-based potential in the RAPPER decoy set (385 loop targets) is comparable to that of AMBER/GBSA for short loops (two to eight residues). The DFIRE is more accurate for longer loops (9 to 12 residues). Similar trend is observed when comparing DFIRE with another physical-based OPLS/SGB-NP energy function in the large Jacobson decoy set (788 loop targets). In the Forrest-Woolf decoy set for the loops of membrane proteins, the DFIRE potential performs substantially better than the combination of the CHARMM force field with several solvation models. The results suggest that a single-term DFIRE-statistical energy function can provide an accurate loop prediction at a fraction of computing cost required for more complicate physical-based energy functions. A Web server for academic users is established for loop selection at the softwares/services section of the Web site http://theory.med.buffalo.edu/.     

Coupling between conformation and proton binding in proteins

Interest centers here on whether the use of a fixed charge distribution of a protein solute, or a treatment that considers proton-binding equilibria by solving the Poisson equation, is a better approach to discriminate native from non-native conformations of proteins. In this analysis of the charge distribution of 7 proteins, we estimate the solvation free energy contribution to the total free energy by exploring the 2(zeta) possible ionization states of the whole molecule, with zeta being the number of ionizable groups in the amino acid sequence, for every conformation in the ensembles of 7 proteins. As an additional consideration of the role of electrostatic interactions in determining the charge distribution of native folds, we carried out a comparison of alternative charge assignment models for the ionizable residues in a set of 21 native-like proteins. The results of this work indicate that (1) for 6 out of 7 proteins, estimation of solvent polarization based on the Generalized Born model with a fixed charge distribution provides the optimal trade-off between accuracy, with respect to the Poisson equation, and speed when compared to the accessible surface area model; for the seventh protein, consideration of all possible ionization states of the whole molecule appears to be crucial to discriminate the native from non-native conformations; (2) significant differences in the degree of ionization and hence the charge distribution for native folds are found between the different charge models examined; (3) the stability of the native state is determined by a delicate balance of all the energy components, and (4) conformational entropy, and hence the dynamics of folding, may play a crucial role for a successful ab initio protein folding prediction.     

Continuum solvent studies of the stability of RNA hairpin loops and helices

We apply continuum solvent models to investigate the relative stability of various conformational forms for two RNA sequences, GGAC(UUCG)GUCC and GGUG(UGAA)CACC. In the first part, we compare alternate hairpin conformations to explore the reliability of these models to discriminate between different local conformations. A second part looks at the hairpin-duplex conversion for the UUCG sequence, identifying major contributors to the thermodynamics of a much large scale transition. Structures were taken as snapshots from multi-nanosecond molecular dynamics simulations computed in a consistent fashion using explicit solvent and with long-range electrostatics accounted for using the Particle-Mesh Ewald procedure. The electrostatic contribution to solvation energies were computed using both a finite-difference Poisson-Boltzmann (PB) model and a pairwise Generalized Born model; non-electrostatic contributions were estimated with a surface-area dependent term. To these solvation free energies were added the mean solute internal energies (determined from a molecular mechanics potential) and estimates of the solute entropy (from a harmonic analysis). Consistent with experiment and with earlier solvated molecular dynamics simulations, the UUCG hairpin was found to prefer conformers close to a recent NMR structure determination in preference to those from an earlier NMR study. Similarly, results for the UGAA hairpin favored an NMR-derived structure over that to be expected for a generic GNRA hairpin loop. Experimental free energies are not known for the hairpin/duplex conversion, but must be close to zero since hairpins are seen in solution and duplexes in crystals; out calculations find a value near zero and illustrate the expected interplay of solvation, salt effects and entropy in affecting this equilibrium.     

Protein side chain conformation predictions with an MMGBSA energy function

The prediction of protein side chain conformations from backbone coordinates is an important task in structural biology, with applications in structure prediction and protein design. It is a difficult problem due to its combinatorial nature. We study the performance of an “MMGBSA” energy function, implemented in our protein design program Proteus, which combines molecular mechanics terms, a Generalized Born and Surface Area (GBSA) solvent model, with approximations that make the model pairwise additive. Proteus is not a competitor to specialized side chain prediction programs due to its cost, but it allows protein design applications, where side chain prediction is an important step and MMGBSA an effective energy model. We predict the side chain conformations for 18 proteins. The side chains are first predicted individually, with the rest of the protein in its crystallographic conformation. Next, all side chains are predicted together. The contributions of individual energy terms are evaluated and various parameterizations are compared. We find that the GB and SA terms, with an appropriate choice of the dielectric constant and surface energy coefficients, are beneficial for single side chain predictions. For the prediction of all side chains, however, errors due to the pairwise additive approximation overcome the improvement brought by these terms. We also show the crucial contribution of side chain minimization to alleviate the rigid rotamer approximation. Even without GB and SA terms, we obtain accuracies comparable to SCWRL4, a specialized side chain prediction program. In particular, we obtain a better RMSD than SCWRL4 for core residues (at a higher cost), despite our simpler rotamer library. Proteins 2016; 84:803–819. © 2016 Wiley Periodicals, Inc.     

Determination of a high-quality nuclear magnetic resonance solution structure of the bovine pancreatic trypsin inhibitor and comparison with three crystal structures.

A high-quality three-dimensional structure of the bovine pancreatic trypsin inhibitor (BPTI) in aqueous solution was determined by 1H nuclear magnetic resonance (n.m.r.) spectroscopy and compared to the three available high-resolution X-ray crystal structures. A newly collected input of 642 distance constraints derived from nuclear Overhauser effects and 115 dihedral angle constraints was used for the structure calculations with the program DIANA, followed by restrained energy minimization with the program AMBER. The BPTI solution structure is represented by a group of 20 conformers with an average root-mean-square deviation (RMSD) relative to the mean solution structure of 0.43 A for backbone atoms and 0.92 A for all heavy atoms of residues 2 to 56. The pairwise RMSD values of the three crystal structures relative to the mean solution structure are 0.76 to 0.85 A for the backbone atoms and 1.24 to 1.33 A for all heavy atoms of residues 2 to 56. Small local differences in backbone atom positions between the solution structure and the X-ray structures near residues 9, 25 to 27, 46 to 48 and 52 to 58, and conformational differences for individual amino acid side-chains were analyzed for possible correlations with intermolecular protein-protein contacts in the crystal lattices, using the pairwise RMSD values among the three crystal structures as a reference.     

Force field based conformational analysis of RNA structural motifs: GNRA tetraloops and their pyrimidine relatives

The protocol of conformational analysis applied here to ribonucleotide oligomers combines conformational search in the space of torsion angles and energy minimization using the AMBER4.1 force field with a continuum treatment of electrostatic solute-solvent interactions. RNA fragments with 5′-GGGCGNNAGCCU-3′ sequences commonly fold into hairpins with four-membered loops. The combinatorial search for acceptable conformations using the MC-SYM program was restricted to loop nucleotides and yielded roughly 1500 structures being compatible with a double-stranded stem. After energy minimization by the JUMNA program (without applying any experimental constraints), these structures converged into an ensemble of 74 different conformers including 26 structures which contained the sheared G-A base pair observed in experimental studies of GNRA tetraloops. Energetic analysis shows that inclusion of solvent electrostatic effects is critically important for the selection of conformers that agree with experimentally determined structures. The continuum model accounts for solvent polarization by means of the electrostatic reaction field. In the case of GNRA loop sequences, the contributions of the reaction field shift relative stabilities towards conformations showing most of the structural features derived from NMR studies. The agreement of computed conformations with the experimental structures of GAAA, GCAA, and GAGA tetraloops suggests that the continuum treatment of the solvent represents a definitive improvement over methods using simple damping models in electrostatic energy calculations. Application of the procedure described here to the evaluation of the relative stabilities of conformers resulting from searching the conformational space of RNA structural motifs provides some progress in (non-homology based) RNA 3D-structure prediction. Received: 20 January 1999 / Revised version: 4 June 1999 / Accepted: 10 June 1999     

Exploratory studies of ab initio protein structure prediction: multiple copy simulated annealing, AMBER energy functions, and a generalized born/solvent accessibility solvation model.

A theoretical and computational approach to ab initio structure prediction for polypeptides in water is described and applied to selected amino acid sequences for testing and preliminary validation. The method builds systematically on the extensive efforts applied to parameterization of molecular dynamics (MD) force fields, employs an empirically well-validated continuum dielectric model for solvation, and an eminently parallelizable approach to conformational search. The effective free energy of polypeptide chains is estimated from AMBER united atom potential functions, with internal degrees of freedom for both backbone and amino acid side chains explicitly treated. The hydration free energy of each structure is determined using the Generalized Born/Solvent Accessibility (GBSA) method, modified and reparameterized to include atom types consistent with the AMBER force field. The conformational search procedure employs a multiple copy, Monte Carlo simulated annealing (MCSA) protocol in full torsion angle space, applied iteratively on sets of structures of progressively lower free energy until a prediction of a structure with lowest effective free energy is obtained. Calibration tests for the effective energy function and search algorithm are performed on the alanine dipeptide, selected protein crystal structures, and united atom decoys on barnase, crambin, and six examples from the Rosetta set. Specific demonstration cases of the method are provided for the 8-mer sequence of Ala residues, a 12-residue peptide with longer side chains QLLKKLLQQLKQ, a de novo designed 16 residue peptide of sequence (AAQAA)3Y, a 15-residue sequence with a beta sheet motif, GEWTWDATKTFTVTE, and a 36 residue small protein, Villin headpiece. The Ala 8-mer readily formed an alpha-helix. An alpha-helix structure was predicted for the 16-mer, consistent with observed results from IR and CD spectroscopy and with the pattern in psi/straight phi angles of known protein structures. The predicted structure for the 12-mer, composed of a mix of helix and less regular elements of secondary structure, lies 2.65 A RMS from the observed crystal structure. Structure prediction for the 8-mer beta-motif resulted in form 4.50 A RMS from the crystal geometry. For Villin, the predicted native form is very close to the crystal structure, RMS values of 3.5 A (including sidechains), and 1.01 A (main chain only). The methodology permits a detailed analysis of the molecular forces which dominate various segments of the predicted folding trajectory. Analysis of the results in terms of internal torsional, electrostatic and van der Waals and the electrostatic and non-electrostatic contributions to hydration, including the hydrophobic effect, is presented.     

Effects of force fields on the conformational and dynamic properties of amyloid β(1‐40) dimer explored by replica exchange molecular dynamics simulations

The conformational space and structural ensembles of amyloid beta (Aβ) peptides and their oligomers in solution are inherently disordered and proven to be challenging to study. Optimum force field selection for molecular dynamics (MD) simulations and the biophysical relevance of results are still unknown. We compared the conformational space of the Aβ(1‐40) dimers by 300 ns replica exchange MD simulations at physiological temperature (310 K) using: the AMBER‐ff99sb‐ILDN, AMBER‐ff99sb*‐ILDN, AMBER‐ff99sb‐NMR, and CHARMM22* force fields. Statistical comparisons of simulation results to experimental data and previously published simulations utilizing the CHARMM22* and CHARMM36 force fields were performed. All force fields yield sampled ensembles of conformations with collision cross sectional areas for the dimer that are statistically significantly larger than experimental results. All force fields, with the exception of AMBER‐ff99sb‐ILDN (8.8 ± 6.4%) and CHARMM36 (2.7 ± 4.2%), tend to overestimate the α‐helical content compared to experimental CD (5.3 ± 5.2%). Using the AMBER‐ff99sb‐NMR force field resulted in the greatest degree of variance (41.3 ± 12.9%). Except for the AMBER‐ff99sb‐NMR force field, the others tended to under estimate the expected amount of β‐sheet and over estimate the amount of turn/bend/random coil conformations. All force fields, with the exception AMBER‐ff99sb‐NMR, reproduce a theoretically expected β‐sheet‐turn‐β‐sheet conformational motif, however, only the CHARMM22* and CHARMM36 force fields yield results compatible with collapse of the central and C‐terminal hydrophobic cores from residues 17‐21 and 30‐36. Although analyses of essential subspace sampling showed only minor variations between force fields, secondary structures of lowest energy conformers are different.     

A divide and conquer approach to fast loop modeling

We describe a fast ab initio method for modeling local segments in protein structures. The algorithm is based on a divide and conquer approach and uses a database of precalculated look-up tables, which represent a large set of possible conformations for loop segments of variable length. The target loop is recursively decomposed until the resulting conformations are small enough to be compiled analytically. The algorithm, which is not restricted to any specific loop length, generates a ranked set of loop conformations in 20-180 s on a desktop PC. The prediction quality is evaluated in terms of global RMSD. Depending on loop length the top prediction varies between 1.06 A RMSD for three-residue loops and 3.72 A RMSD for eight-residue loops. Due to its speed the method may also be useful to generate alternative starting conformations for complex simulations.     

Ab initio construction of all-atom loop conformations

In this study, a new ab initio method named CLOOP has been developed to build all-atom loop conformations. In this method, a loop main-chain conformation is generated by sampling main-chain dihedral angles from a restrained varphi/psi set, and the side-chain conformations are built randomly. The CHARMM all-atom force field was used to evaluate the loop conformations. Soft core potentials were used to treat the non-bond interactions, and a designed energy-minimization technique was used to close and optimize the loop conformations. It is shown that the two strategies improve the computational efficiency and the loop-closure rate substantially compared to normal minimization methods. CLOOP was used to construct the conformations of 4-, 8-, and 12-residue loops in Fiser's test set. The average main-chain root-mean-square deviations obtained in 1,000 trials for the 10 different loops of each size are 0.33, 1.27, and 2.77 A, respectively. CLOOP can build all-atom loop conformations with a sampling accuracy comparable with previous loop main-chain construction algorithms. [Figure: see text].     

Ab initio modeling of small, medium, and large loops in proteins.

This study presents different procedures for ab initio modeling of peptide loops of different sizes in proteins. Small loops (up to 8--12 residues) were generated by a straightforward procedure with subsequent "averaging" over all the low-energy conformers obtained. The averaged conformer fairly represents the entire set of low-energy conformers, root mean square deviation (RMSD) values being from 1.01 A for a 4-residue loop to 1.94 A for an 8-residue loop. Three-dimensional (3D) structures for several medium loops (20--30 residues) and for two large loops (54 and 61 residues) were predicted using residue-residue contact matrices divided into variable parts corresponding to the loops, and into a constant part corresponding to the known core of the protein. For each medium loop, a very limited number of sterically reasonable C(alpha) traces (from 1 to 3) was found; RMSD values ranged from 2.4 to 5.9 A. Single C(alpha) traces predicted for each of the large loops possessed RMSD values of 4.5 A. Generally, ab initio loop modeling presented in this work combines elements of computational procedures developed both for protein folding and for peptide conformational analysis.     

Effects of electric field on alamethicin bound at the lipid-water interface: a molecular mechanics study.

A systematic molecular mechanics study of the alamethicin molecule was made to determine a set of low-energy conformers in vacuo and in aqueous environment. The behavior of these conformers was investigated at the phase boundary which was modeled as a plane dividing two compartments with solvation properties of water and octanol with a constant electric field applied normal to the boundary. The calculations were performed with a molecular mechanics program for calculation of stable conformations at the phase boundary utilizing the Empiric Conformational Energy Program for Peptides force field and the Hopfinger-Scheraga solvation model. 371 minimum energy conformers of alamethicin, determined in vacuo with the build-up procedure, were used as starting conformations for energy minimization in aqueous environment and at the phase boundary. Only 49 interphase-bound structures were within 12 kcal/mol of the minima which was found. No helical structures having values close to the canonical parameters for an alpha- or 3(10)-helix were found despite the presence of eight alpha-methylalanine residues which favor the formation of these helices; four helix-like structures were found, having all negative phi, psi values. All the helical conformers have very high energies in water (approximately 14 kcal/mol), but are quite stable at the phase boundary (3.7-6.8 kcal/mol above the lowest minima found). The implications of these results for proposed mechanisms for membrane-binding and voltage-dependent gating are considered.     

High-resolution prediction of protein helix positions and orientations

We have developed a new method for predicting helix positions in globular proteins that is intended primarily for comparative modeling and other applications where high precision is required. Unlike helix packing algorithms designed for ab initio folding, we assume that knowledge is available about the qualitative placement of all helices. However, even among homologous proteins, the corresponding helices can demonstrate substantial differences in positions and orientations, and for this reason, improperly positioned helices can contribute significantly to the overall backbone root-mean-square deviation (RMSD) of comparative models. A helix packing algorithm for use in comparative modeling must obtain high precision to be useful, and for this reason we utilize an all-atom protein force field (OPLS) and a Generalized Born continuum solvent model. To reduce the computational expense associated with using a detailed, physics-based energy function, we have developed new hierarchical and multiscale algorithms for sampling the helices and flanking loops. We validate the method using a test suite of 33 cases, which are drawn from a diverse set of high-resolution crystal structures. The helix positions are reproduced with an average backbone RMSD of 0.6 A, while the average backbone RMSD of the complete loop-helix-loop region (i.e., the helix with the surrounding loops, which are also repredicted) is 1.3 A.