Application of torsion angle molecular dynamics for efficient sampling of protein conformations

The SCR-17 and SCR-18 glycans in human complement Factor H enhance its regulatory function

Human complement factor H (CFH) plays a central role in regulating activated C3b to protect host cells. CFH contain 20 short complement regulator (SCR) domains and eight N-glycosylation sites. The N-terminal SCR domains mediate C3b degradation while the C-terminal CFH domains bind to host cell surfaces to protect these. Our earlier study of Pichia-generated CFH fragments indicated a self-association site at SCR-17/18 that comprises a dimerization site for human factor H. Two N-linked glycans are located on SCR-17 and SCR-18. Here, when we expressed SCR-17/18 without glycans in an E. coli system, analytical ultracentrifugation showed that no dimers were now formed. To investigate this novel finding, full-length CFH and its C-terminal fragments were purified from human plasma and Pichia pastoris respectively, and their glycans were enzymatically removed using PNGase F. Using size-exclusion chromatography, mass spectrometry, and analytical ultracentrifugation, SCR-17/18 from Pichia showed notably less dimer formation without its glycans, confirming that the glycans are necessary for the formation of SCR-17/18 dimers. By surface plasmon resonance, affinity analyses interaction showed decreased binding of deglycosylated full-length CFH to immobilised C3b, showing that CFH glycosylation enhances the key CFH regulation of C3b. We conclude that our study revealed a significant new aspect of CFH regulation based on its glycosylation and its resulting dimerisation.

The solution structure of the heavy chain-only C5-Fc nanobody reveals exposed variable regions that are optimal for COVID-19 antigen interactions

Heavy chain-only antibodies can offer advantages of higher binding affinities, reduced sizes, and higher stabilities than conventional antibodies. To address the challenge of SARS-CoV-2 coronavirus, a llama-derived single-domain nanobody C5 was developed previously that has high COVID-19 virus neutralization potency. The fusion protein C5-Fc comprises two C5 domains attached to a glycosylated Fc region of a human IgG1 antibody and shows therapeutic efficacy in vivo. Here, we have characterized the solution arrangement of the molecule. Two 1443 Da N-linked glycans seen in the mass spectra of C5-Fc were removed and the glycosylated and deglycosylated structures were evaluated. Reduction of C5-Fc with 2-mercaptoethylamine indicated three interchain Cys-Cys disulfide bridges within the hinge. The X-ray and neutron Guinier RG values, which provide information about structural elongation, were similar at 4.1 to 4.2 nm for glycosylated and deglycosylated C5-Fc. To explain these RG values, atomistic scattering modeling based on Monte Carlo simulations resulted in 72,737 and 56,749 physically realistic trial X-ray and neutron structures, respectively. From these, the top 100 best-fit X-ray and neutron models were identified as representative asymmetric solution structures, similar to that of human IgG1, with good R-factors below 2.00%. Both C5 domains were solvent exposed, consistent with the functional effectiveness of C5-Fc. Greater disorder occurred in the Fc region after deglycosylation. Our results clarify the importance of variable and exposed C5 conformations in the therapeutic function of C5-Fc, while the glycans in the Fc region are key for conformational stability in C5-Fc.

Conformational Search for the Building Block of Proteins Based on the Gradient Gravitational Search Algorithm (ConfGGS) Using Force Fields: CHARMM, AMBER, and OPLS-AA

Proteins are linear polymers built from a repertoire of 20 different amino acids, which are considered building blocks of proteins. The diversity and versatility of these 20 building blocks with regard to their conformations are key to adopting three-dimensional structures that facilitate proteins to undergo important mechanistic biological processes in living systems. The present investigation reports a conformational search of 20 different amino acids, building blocks of proteins, using three different force fields, CHARMM, AMBER, and OPLS-AA, implemented in the gradient gravitational search algorithm. The search technique (ConfGGS) includes the contribution from both bonded and nonbonded terms using Cartesian coordinates. The efficiency of such conformational searches has also been compared with other optimization algorithms: DE/Best, DE/Rand, and PSO algorithms with respect to computational time and accuracy based on the minimum number of iteration steps and computed lowest mean absolute error (MAE) and mean standard deviation (MSD) values for dihedral angles of respective near-optimal structures. Moreover, the ConfGGS technique has also been extended to an ordered protein fragment (PQITL) extracted from HIV-1 protease (PDB ID: 1YTH), an intrinsically disordered protein fragment, i.e., an amyloid-forming segment (AVVTGVTAV), from the NAC domain of Parkinson's disease protein α-synuclein, residues 69-77 (PDB ID: 4RIK), the experimental NMR atomic-resolution structure of α-synuclein fibrils (PDB ID: 2N0A), and a disulfide bond-containing protein fragment sequence (PCYGWPVCY), residues 59-67 (PDB ID: 6Y4F) toward structure prediction as a close homologue compared with experimental accuracy, using the CHARMM force field. The MolProbity validation results for the protein fragment (PQITL) obtained by ConfGGS/CHARMM are in better agreement with the native protein fragment structure of HIV-1 protease (PDB ID: 1YTH). Furthermore, the computed results have also been compared with the coordinates obtained from the AlphaFold network.

Characterization of an Extensive Interface on Vitronectin for Binding to Plasminogen Activator Inhibitor-1: Adoption of Structure in an Intrinsically Disordered Region

Small-angle neutron scattering (SANS) measurements were pursued to study human vitronectin, a protein found in tissues and the circulation that regulates cell adhesion/migration and proteolytic cascades that govern hemostasis and pericellular proteolysis. Many of these functions occur via interactions with its binding partner, plasminogen activator inhibitor-1 (PAI-1), the chief inhibitor of proteases that lyse and activate plasminogen. We focused on a region of vitronectin that remains uncharacterized from previous X-ray scattering, nuclear magnetic resonance, and computational modeling approaches and which we propose is involved in binding to PAI-1. This region, which bridges the N-terminal somatomedin B (SMB) domain with a large central β-propeller domain of vitronectin, appears unstructured and has characteristics of an intrinsically disordered domain (IDD). The effect of osmolytes was evaluated using circular dichroism and SANS to explore the potential of the IDD to undergo a disorder-to-order transition. The results suggest that the IDD favors a more ordered structure under osmotic pressure; SANS shows a smaller radius of gyration (R_g) and a more compact fold of the IDD upon addition of osmolytes. To test whether PAI-1 binding is also coupled to folding within the IDD structure, a set of SANS experiments with contrast variation were performed on the complex of PAI-1 with a vitronectin fragment corresponding to the N-terminal 130 amino acids (denoted the SMB-IDD because it contains the SMB domain and IDD in linear sequence). Analysis of the SANS data using the Ensemble Optimization Method confirms that the SMB-IDD adopts a more compact configuration when bound to PAI-1. Calculated structures for the PAI-1:SMB-IDD complex suggest that the IDD provides an interaction surface outside of the primary PAI-1-binding site located within the SMB domain; this binding is proposed to lead to the assembly of higher-order structures of vitronectin and PAI-1 commonly found in tissues.

Efficient minimization of multipole electrostatic potentials in torsion space

The development of models of macromolecular electrostatics capable of delivering improved fidelity to quantum mechanical calculations is an active field of research in computational chemistry. Most molecular force field development takes place in the context of models with full Cartesian coordinate degrees of freedom. Nevertheless, a number of macromolecular modeling programs use a reduced set of conformational variables limited to rotatable bonds. Efficient algorithms for minimizing the energies of macromolecular systems with torsional degrees of freedom have been developed with the assumption that all atom-atom interaction potentials are isotropic. We describe novel modifications to address the anisotropy of higher order multipole terms while retaining the efficiency of these approaches. In addition, we present a treatment for obtaining derivatives of atom-centered tensors with respect to torsional degrees of freedom. We apply these results to enable minimization of the Amoeba multipole electrostatics potential in a system with torsional degrees of freedom, and validate the correctness of the gradients by comparison to finite difference approximations. In the interest of enabling a complete model of electrostatics with implicit treatment of solvent-mediated effects, we also derive expressions for the derivative of solvent accessible surface area with respect to torsional degrees of freedom.

Hamiltonian Monte Carlo with Constrained Molecular Dynamics as Gibbs Sampling

Compared to fully flexible molecular dynamics, simulations of constrained systems can use larger time steps and focus kinetic energy on soft degrees of freedom. Achieving ergodic sampling from the Boltzmann distribution, however, has proven challenging. Using recent generalizations of the equipartition principle and Fixman potential, here we implement Hamiltonian Monte Carlo based on constrained molecular dynamics as a Gibbs sampling move. By mixing Hamiltonian Monte Carlo based on fully flexible and torsional dynamics, we are able to reproduce free energy landscapes of simple model systems and enhance sampling of macrocycles.

Computational protein structure refinement: Almost there, yet still so far to go

Protein structures are essential in modern biology yet experimental methods are far from being able to catch up with the rapid increase in available genomic data. Computational protein structure prediction methods aim to fill the gap while the role of protein structure refinement is to take approximate initial template-based models and bring them closer to the true native structure. Current methods for computational structure refinement rely on molecular dynamics simulations, related sampling methods, or iterative structure optimization protocols. The best methods are able to achieve moderate degrees of refinement but consistent refinement that can reach near-experimental accuracy remains elusive. Key issues revolve around the accuracy of the energy function, the inability to reliably rank multiple models, and the use of restraints that keep sampling close to the native state but also limit the degree of possible refinement. A different aspect is the question of what exactly the target of high-resolution refinement should be as experimental structures are affected by experimental conditions and different biological questions require varying levels of accuracy. While improvement of the global protein structure is a difficult problem, high-resolution refinement methods that improves local structural quality such as favorable stereochemistry and the avoidance of atomic clashes are much more successful.

Atomistic modelling of scattering data in the Collaborative Computational Project for Small Angle Scattering (CCP-SAS)

The capabilities of current computer simulations provide a unique opportunity to model small-angle scattering (SAS) data at the atomistic level, and to include other structural constraints ranging from molecular and atomistic energetics to crystallography, electron microscopy and NMR. This extends the capabilities of solution scattering and provides deeper insights into the physics and chemistry of the systems studied. Realizing this potential, however, requires integrating the experimental data with a new generation of modelling software. To achieve this, the CCP-SAS collaboration (http://www.ccpsas.org/) is developing open-source, high-throughput and user-friendly software for the atomistic and coarse-grained molecular modelling of scattering data. Robust state-of-the-art molecular simulation engines and molecular dynamics and Monte Carlo force fields provide constraints to the solution structure inferred from the small-angle scattering data, which incorporates the known physical chemistry of the system. The implementation of this software suite involves a tiered approach in which GenApp provides the deployment infrastructure for running applications on both standard and high-performance computing hardware, and SASSIE provides a workflow framework into which modules can be plugged to prepare structures, carry out simulations, calculate theoretical scattering data and compare results with experimental data. GenApp produces the accessible web-based front end termed SASSIE-web, and GenApp and SASSIE also make community SAS codes available. Applications are illustrated by case studies: (i) inter-domain flexibility in two- to six-domain proteins as exemplified by HIV-1 Gag, MASP and ubiquitin; (ii) the hinge conformation in human IgG2 and IgA1 antibodies; (iii) the complex formed between a hexameric protein Hfq and mRNA; and (iv) synthetic 'bottlebrush' polymers.

Assessment of the Sampling Performance of Multiple-Copy Dynamics versus a Unique Trajectory

The goal of the present study was to ascertain the differential performance of a long molecular dynamics trajectory versus several shorter ones starting from different points in the phase space and covering the same sampling time. For this purpose, we selected the 16-mer peptide Bak₁₆BH3 as a model for study and carried out several samplings in explicit solvent. These samplings included an 8 μs trajectory (sampling S1); two 4 μs trajectories (sampling S2); four 2 μs trajectories (sampling S3); eight 1 μs trajectories (sampling S4); 16 0.5 μs trajectories (sampling S5), and 80 0.1 μs trajectories (sampling S6). Moreover, the 8 μs trajectory was further extended to 16 μs to have reference values of the diverse properties measured. The diverse samplings were compared qualitatively and quantitatively. Among the former, we carried out a comparison of the conformational profiles of the peptide using cluster analysis. Moreover, we also gained insight into the interchange among these structures along the sampling process. Among the latter, we computed the number of new conformational patterns sampled with time using strings defined from the conformations attained by each of the residues in the peptide. We also compared the locations and depths of the obtained minima on the free energy surface using principal component analysis. Finally, we also compared the helical profiles per residue at the end of the sampling process. The results suggest that a few short molecular dynamics trajectories may provide better sampling than one unique trajectory. Moreover, this procedure can also be advantageous to avoid getting trapped in a local minimum. However, caution should be exercised since short trajectories need to be long enough to overcome local barriers surrounding the starting point and the required sampling time depends on the number of degrees of freedom of the system under study. An effective way to gain insight into the minimum MD trajectory length is to monitor the convergence of different structural features, as shown in the present work.

Implementation and application of helix-helix distance and crossing angle restraint potentials

Based on the definition of helix-helix distance and crossing angle introduced by Chothia et al. (J Mol Biol 1981, 145, 215), we have developed the restraint potentials by which the distance and crossing angle of two selected helices can be maintained around target values during molecular dynamics simulations. A series of assessments show that calculated restraint forces are numerically accurate. Since the restraint forces are only exerted on atoms which define the helical principal axes, each helix can rotate along its helical axis, depending on the helix-helix intermolecular interactions. Such a restraint potential enables us to characterize the helix-helix interactions at atomic details by sampling their conformational space around specific distance and crossing angle with (restraint) force-dependent fluctuations. Its efficacy is illustrated by calculating the potential of mean force as a function of helix-helix distance between two transmembrane helical peptides in an implicit membrane model.

Overcoming potential energy distortions in constrained internal coordinate molecular dynamics simulations

The Internal Coordinate Molecular Dynamics (ICMD) method is an attractive molecular dynamics (MD) method for studying the dynamics of bonded systems such as proteins and polymers. It offers a simple venue for coarsening the dynamics model of a system at multiple hierarchical levels. For example, large scale protein dynamics can be studied using torsional dynamics, where large domains or helical structures can be treated as rigid bodies and the loops connecting them as flexible torsions. ICMD with such a dynamic model of the protein, combined with enhanced conformational sampling method such as temperature replica exchange, allows the sampling of large scale domain motion involving high energy barrier transitions. Once these large scale conformational transitions are sampled, all-torsion, or even all-atom, MD simulations can be carried out for the low energy conformations sampled via coarse grained ICMD to calculate the energetics of distinct conformations. Such hierarchical MD simulations can be carried out with standard all-atom forcefields without the need for compromising on the accuracy of the forces. Using constraints to treat bond lengths and bond angles as rigid can, however, distort the potential energy landscape of the system and reduce the number of dihedral transitions as well as conformational sampling. We present here a two-part solution to overcome such distortions of the potential energy landscape with ICMD models. To alleviate the intrinsic distortion that stems from the reduced phase space in torsional MD, we use the Fixman compensating potential. To additionally alleviate the extrinsic distortion that arises from the coupling between the dihedral angles and bond angles within a force field, we propose a hybrid ICMD method that allows the selective relaxing of bond angles. This hybrid ICMD method bridges the gap between all-atom MD and torsional MD. We demonstrate with examples that these methods together offer a solution to eliminate the potential energy distortions encountered in constrained ICMD simulations of peptide molecules.

Deuterium Labeling Together with Contrast Variation Small-Angle Neutron Scattering Suggests How Skp Captures and Releases Unfolded Outer Membrane Proteins

In Gram-negative bacteria, the chaperone protein Skp forms specific and stable complexes with membrane proteins while they are transported across the periplasm to the outer membrane. The jellyfish-like architecture of Skp is similar to the eukaryotic and archaeal prefoldins and the mitochondrial Tim chaperones, that is the α-helical "tentacles" extend from a β-strand "body" to create an internal cavity. Contrast variation small-angle neutron scattering (SANS) experiments on Skp alone in solution and bound in two different complexes to unfolded outer membrane proteins (uOMPs), OmpA and OmpW, demonstrate that the helical tentacles of Skp bind their substrate in a clamp-like mechanism in a conformation similar to that previously observed in the apo crystal structure of Skp. Deuteration of the uOMP component combined with contrast variation analysis allowed the shapes of Skp and uOMP as well as the location of uOMP with respect to Skp to be determined in both complexes. This represents unique information that could not be obtained without deuterium labeling of the uOMPs. The data yield the first direct structural evidence that the α-helical Skp tentacles move closer together on binding its substrate and that the structure of Skp is different when binding different uOMPs. This work presents, by example, a tutorial on performing SANS experiments using both deuterium labeling and contrast variation, including SANS theory, sample preparation, data collection, sample quality validation, data analysis, and structure modeling.

Internal coordinate molecular dynamics: a foundation for multiscale dynamics

Internal coordinates such as bond lengths, bond angles, and torsion angles (BAT) are natural coordinates for describing a bonded molecular system. However, the molecular dynamics (MD) simulation methods that are widely used for proteins, DNA, and polymers are based on Cartesian coordinates owing to the mathematical simplicity of the equations of motion. However, constraints are often needed with Cartesian MD simulations to enhance the conformational sampling. This makes the equations of motion in the Cartesian coordinates differential-algebraic, which adversely impacts the complexity and the robustness of the simulations. On the other hand, constraints can be easily placed in BAT coordinates by removing the degrees of freedom that need to be constrained. Thus, the internal coordinate MD (ICMD) offers an attractive alternative to Cartesian coordinate MD for developing multiscale MD method. The torsional MD method is a special adaptation of the ICMD method, where all the bond lengths and bond angles are kept rigid. The advantages of ICMD simulation methods are the longer time step size afforded by freezing high frequency degrees of freedom and performing a conformational search in the more important low frequency torsional degrees of freedom. However, the advancements in the ICMD simulations have been slow and stifled by long-standing mathematical bottlenecks. In this review, we summarize the recent mathematical advancements we have made based on spatial operator algebra, in developing a robust long time scale ICMD simulation toolkit useful for various applications. We also present the applications of ICMD simulations to study conformational changes in proteins and protein structure refinement. We review the advantages of the ICMD simulations over the Cartesian simulations when used with enhanced sampling methods and project the future use of ICMD simulations in protein dynamics.

GneimoSim: a modular internal coordinates molecular dynamics simulation package

The generalized Newton-Euler inverse mass operator (GNEIMO) method is an advanced method for internal coordinates molecular dynamics (ICMD). GNEIMO includes several theoretical and algorithmic advancements that address longstanding challenges with ICMD simulations. In this article, we describe the GneimoSim ICMD software package that implements the GNEIMO method. We believe that GneimoSim is the first software package to include advanced features such as the equipartition principle derived for internal coordinates, and a method for including the Fixman potential to eliminate systematic statistical biases introduced by the use of hard constraints. Moreover, by design, GneimoSim is extensible and can be easily interfaced with third party force field packages for ICMD simulations. Currently, GneimoSim includes interfaces to LAMMPS, OpenMM, and Rosetta force field calculation packages. The availability of a comprehensive Python interface to the underlying C++ classes and their methods provides a powerful and versatile mechanism for users to develop simulation scripts to configure the simulation and control the simulation flow. GneimoSim has been used extensively for studying the dynamics of protein structures, refinement of protein homology models, and for simulating large scale protein conformational changes with enhanced sampling methods. GneimoSim is not limited to proteins and can also be used for the simulation of polymeric materials.

A simple molecular mechanics integrator in mixed rigid body and dihedral angle space

We propose a numerical scheme to integrate equations of motion in a mixed space of rigid-body and dihedral angle coordinates. The focus of the presentation is biomolecular systems and the framework is applicable to polymers with tree-like topology. By approximating the effective mass matrix as diagonal and lumping all bias torques into the time dependencies of the diagonal elements, we take advantage of the formal decoupling of individual equations of motion. We impose energy conservation independently for every degree of freedom and this is used to derive a numerical integration scheme. The cost of all auxiliary operations is linear in the number of atoms. By coupling the scheme to one of two popular thermostats, we extend the method to sample constant temperature ensembles. We demonstrate that the integrator of choice yields satisfactory stability and is free of mass-metric tensor artifacts, which is expected by construction of the algorithm. Two fundamentally different systems, viz., liquid water and an α-helical peptide in a continuum solvent are used to establish the applicability of our method to a wide range of problems. The resultant constant temperature ensembles are shown to be thermodynamically accurate. The latter relies on detailed, quantitative comparisons to data from reference sampling schemes operating on exactly the same sets of degrees of freedom.

Protein structure refinement of CASP target proteins using GNEIMO torsional dynamics method

A longstanding challenge in using computational methods for protein structure prediction is the refinement of low-resolution structural models derived from comparative modeling methods into highly accurate atomistic models useful for detailed structural studies. Previously, we have developed and demonstrated the utility of the internal coordinate molecular dynamics (MD) technique, generalized Newton-Euler inverse mass operator (GNEIMO), for refinement of small proteins. Using GNEIMO, the high-frequency degrees of freedom are frozen and the protein is modeled as a collection of rigid clusters connected by torsional hinges. This physical model allows larger integration time steps and focuses the conformational search in the low frequency torsional degrees of freedom. Here, we have applied GNEIMO with temperature replica exchange to refine low-resolution protein models of 30 proteins taken from the continuous assessment of structure prediction (CASP) competition. We have shown that GNEIMO torsional MD method leads to refinement of up to 1.3 Å in the root-mean-square deviation in coordinates for 30 CASP target proteins without using any experimental data as restraints in performing the GNEIMO simulations. This is in contrast with the unconstrained all-atom Cartesian MD method performed under the same conditions, where refinement requires the use of restraints during the simulations.

QM/MM investigations of organic chemistry oriented questions

About 35 years after its first suggestion, QM/MM became the standard theoretical approach to investigate enzymatic structures and processes. The success is due to the ability of QM/MM to provide an accurate atomistic picture of enzymes and related processes. This picture can even be turned into a movie if nuclei-dynamics is taken into account to describe enzymatic processes. In the field of organic chemistry, QM/MM methods are used to a much lesser extent although almost all relevant processes happen in condensed matter or are influenced by complicated interactions between substrate and catalyst. There is less importance for theoretical organic chemistry since the influence of nonpolar solvents is rather weak and the effect of polar solvents can often be accurately described by continuum approaches. Catalytic processes (homogeneous and heterogeneous) can often be reduced to truncated model systems, which are so small that pure quantum-mechanical approaches can be employed. However, since QM/MM becomes more and more efficient due to the success in software and hardware developments, it is more and more used in theoretical organic chemistry to study effects which result from the molecular nature of the environment. It is shown by many examples discussed in this review that the influence can be tremendous, even for nonpolar reactions. The importance of environmental effects in theoretical spectroscopy was already known. Due to its benefits, QM/MM can be expected to experience ongoing growth for the next decade.In the present chapter we give an overview of QM/MM developments and their importance in theoretical organic chemistry, and review applications which give impressions of the possibilities and the importance of the relevant effects. Since there is already a bunch of excellent reviews dealing with QM/MM, we will discuss fundamental ingredients and developments of QM/MM very briefly with a focus on very recent progress. For the applications we follow a similar strategy.

Fixman compensating potential for general branched molecules

The technique of constraining high frequency modes of molecular motion is an effective way to increase simulation time scale and improve conformational sampling in molecular dynamics simulations. However, it has been shown that constraints on higher frequency modes such as bond lengths and bond angles stiffen the molecular model, thereby introducing systematic biases in the statistical behavior of the simulations. Fixman proposed a compensating potential to remove such biases in the thermodynamic and kinetic properties calculated from dynamics simulations. Previous implementations of the Fixman potential have been limited to only short serial chain systems. In this paper, we present a spatial operator algebra based algorithm to calculate the Fixman potential and its gradient within constrained dynamics simulations for branched topology molecules of any size. Our numerical studies on molecules of increasing complexity validate our algorithm by demonstrating recovery of the dihedral angle probability distribution function for systems that range in complexity from serial chains to protein molecules. We observe that the Fixman compensating potential recovers the free energy surface of a serial chain polymer, thus annulling the biases caused by constraining the bond lengths and bond angles. The inclusion of Fixman potential entails only a modest increase in the computational cost in these simulations. We believe that this work represents the first instance where the Fixman potential has been used for general branched systems, and establishes the viability for its use in constrained dynamics simulations of proteins and other macromolecules.

Mapping conformational dynamics of proteins using torsional dynamics simulations

All-atom molecular dynamics simulations are widely used to study the flexibility of protein conformations. However, enhanced sampling techniques are required for simulating protein dynamics that occur on the millisecond timescale. In this work, we show that torsional molecular dynamics simulations enhance protein conformational sampling by performing conformational search in the low-frequency torsional degrees of freedom. In this article, we use our recently developed torsional-dynamics method called Generalized Newton-Euler Inverse Mass Operator (GNEIMO) to study the conformational dynamics of four proteins. We investigate the use of the GNEIMO method in simulations of the conformationally flexible proteins fasciculin and calmodulin, as well as the less flexible crambin and bovine pancreatic trypsin inhibitor. For the latter two proteins, the GNEIMO simulations with an implicit-solvent model reproduced the average protein structural fluctuations and sample conformations similar to those from Cartesian simulations with explicit solvent. The application of GNEIMO with replica exchange to the study of fasciculin conformational dynamics produced sampling of two of this protein's experimentally established conformational substates. Conformational transition of calmodulin from the Ca(2+)-bound to the Ca(2+)-free conformation occurred readily with GNEIMO simulations. Moreover, the GNEIMO method generated an ensemble of conformations that satisfy about half of both short- and long-range interresidue distances obtained from NMR structures of holo to apo transitions in calmodulin. Although unconstrained all-atom Cartesian simulations have failed to sample transitions between the substates of fasciculin and calmodulin, GNEIMO simulations show the transitions in both systems. The relatively short simulation times required to capture these long-timescale conformational dynamics indicate that GNEIMO is a promising molecular-dynamics technique for studying domain motion in proteins.

Statistical torsion angle potential energy functions for protein structure modeling: a bicubic interpolation approach

A set of grid type knowledge-based energy functions is introduced for ϕ-χ₁ , ψ-χ₁ , ϕ-ψ, and χ₁ -χ₂ torsion angle combinations. Boltzmann distribution is assumed for the torsion angle populations from protein X-ray structures, and the functions are named as statistical torsion angle potential energy functions. The grid points around periodic boundaries are duplicated to force periodicity, and the remedy relieves the derivative discontinuity problem. The devised functions rapidly improve the quality of model structures. The potential bias in the functions and the usefulness of additional secondary structure information are also investigated. The proposed guiding functions are expected to facilitate protein structure modeling, such as protein structure prediction, protein design, and structure refinement.

Advanced techniques for constrained internal coordinate molecular dynamics

Internal coordinate molecular dynamics (ICMD) methods provide a more natural description of a protein by using bond, angle, and torsional coordinates instead of a Cartesian coordinate representation. Freezing high-frequency bonds and angles in the ICMD model gives rise to constrained ICMD (CICMD) models. There are several theoretical aspects that need to be developed to make the CICMD method robust and widely usable. In this article, we have designed a new framework for (1) initializing velocities for nonindependent CICMD coordinates, (2) efficient computation of center of mass velocity during CICMD simulations, (3) using advanced integrators such as Runge-Kutta, Lobatto, and adaptive CVODE for CICMD simulations, and (4) cancelling out the "flying ice cube effect" that sometimes arises in Nosé-Hoover dynamics. The Generalized Newton-Euler Inverse Mass Operator (GNEIMO) method is an implementation of a CICMD method that we have developed to study protein dynamics. GNEIMO allows for a hierarchy of coarse-grained simulation models based on the ability to rigidly constrain any group of atoms. In this article, we perform tests on the Lobatto and Runge-Kutta integrators to determine optimal simulation parameters. We also implement an adaptive coarse-graining tool using the GNEIMO Python interface. This tool enables the secondary structure-guided "freezing and thawing" of degrees of freedom in the molecule on the fly during molecular dynamics simulations and is shown to fold four proteins to their native topologies. With these advancements, we envision the use of the GNEIMO method in protein structure prediction, structure refinement, and in studying domain motion.

Equipartition Principle for Internal Coordinate Molecular Dynamics

The principle of equipartition of (kinetic) energy for all-atom Cartesian molecular dynamics states that each momentum phase space coordinate on the average has ½kT of kinetic energy in a canonical ensemble. This principle is used in molecular dynamics simulations to initialize velocities, and to calculate statistical properties such as entropy. Internal coordinate molecular dynamics (ICMD) models differ from Cartesian models in that the overall kinetic energy depends on the generalized coordinates and includes cross-terms. Due to this coupled structure, no such equipartition principle holds for ICMD models. In this paper we introduce non-canonical modal coordinates to recover some of the structural simplicity of Cartesian models and develop a new equipartition principle for ICMD models. We derive low-order recursive computational algorithms for transforming between the modal and physical coordinates. The equipartition principle in modal coordinates provides a rigorous method for initializing velocities in ICMD simulations thus replacing the ad hoc methods used until now. It also sets the basis for calculating conformational entropy using internal coordinates.

Computational models of protein kinematics and dynamics: beyond simulation

Physics-based simulation represents a powerful method for investigating the time-varying behavior of dynamic protein systems at high spatial and temporal resolution. Such simulations, however, can be prohibitively difficult or lengthy for large proteins or when probing the lower-resolution, long-timescale behaviors of proteins generally. Importantly, not all questions about a protein system require full space and time resolution to produce an informative answer. For instance, by avoiding the simulation of uncorrelated, high-frequency atomic movements, a larger, domain-level picture of protein dynamics can be revealed. The purpose of this review is to highlight the growing body of complementary work that goes beyond simulation. In particular, this review focuses on methods that address kinematics and dynamics, as well as those that address larger organizational questions and can quickly yield useful information about the long-timescale behavior of a protein.

Structure refinement of protein low resolution models using the GNEIMO constrained dynamics method

The challenge in protein structure prediction using homology modeling is the lack of reliable methods to refine the low resolution homology models. Unconstrained all-atom molecular dynamics (MD) does not serve well for structure refinement due to its limited conformational search. We have developed and tested the constrained MD method, based on the generalized Newton-Euler inverse mass operator (GNEIMO) algorithm for protein structure refinement. In this method, the high-frequency degrees of freedom are replaced with hard holonomic constraints and a protein is modeled as a collection of rigid body clusters connected by flexible torsional hinges. This allows larger integration time steps and enhances the conformational search space. In this work, we have demonstrated the use of torsional GNEIMO method without using any experimental data as constraints, for protein structure refinement starting from low-resolution decoy sets derived from homology methods. In the eight proteins with three decoys for each, we observed an improvement of ~2 Å in the rmsd in coordinates to the known experimental structures of these proteins. The GNEIMO trajectories also showed enrichment in the population density of native-like conformations. In addition, we demonstrated structural refinement using a "freeze and thaw" clustering scheme with the GNEIMO framework as a viable tool for enhancing localized conformational search. We have derived a robust protocol based on the GNEIMO replica exchange method for protein structure refinement that can be readily extended to other proteins and possibly applicable for high throughput protein structure refinement.

50 Years of Lifson-Roig Models: Application to Molecular Simulation Data

Simple helix-coil transition theories have been indispensable tools in the analysis of data reporting on the reversible folding of α-helical polypeptides. They provide a transferable means to not only characterize different systems but to also compare different techniques, viz., experimental probes monitoring helix-coil transitions in vitro or biomolecular force fields in silico. This article addresses several issues with the application of Lifson-Roig theory to helix-coil transition data. We use computer simulation to generate two sets of ensembles for the temperature-controlled, reversible folding of the 21-residue, alanine-rich FS peptide. Ensembles differ in the rigidity of backbone bond angles and are analyzed using two distinct descriptors of helicity. The analysis unmasks an underlying phase diagram that is surprisingly complex. The complexities give rise to fitted nucleation and propagation parameters that are difficult to interpret and that are inconsistent with the distribution of isolated residues in the α-helical basin. We show that enthalpies of helix formation are more robustly determined using van't Hoff analysis of simple measures of helicity rather than fitted propagation parameters. To overcome some of these issues, we design a simple variant of the Lifson-Roig model that recovers physical interpretability of the obtained parameters by allowing bundle formation to be described in simple fashion. The relevance of our results is discussed in relation to the applicability of Lifson-Roig models to both in silico and in vitro data.

Solid-state NMR ensemble dynamics as a mediator between experiment and simulation

Solid-state NMR (SSNMR) is a powerful technique to describe the orientations of membrane proteins and peptides in their native membrane bilayer environments. The deuterium ((2)H) quadrupolar splitting (DQS), one of the SSNMR observables, has been used to characterize the orientations of various single-pass transmembrane (TM) helices using a semistatic rigid-body model such as the geometric analysis of labeled alanine (GALA) method. However, dynamic information of these TM helices, which could be related to important biological function, can be missing or misinterpreted with the semistatic model. We have investigated the orientation of WALP23 in an implicit membrane of dimyristoylglycerophosphocholine by determining an ensemble of structures using multiple conformer models with a DQS restraint potential. When a single conformer is used, the resulting helix orientation (tilt angle (τ) of 5.6 ± 3.2° and rotation angle (ρ) of 141.8 ± 40.6°) is similar to that determined by the GALA method. However, as the number of conformers is increased, the tilt angles of WALP23 ensemble structures become larger (26.9 ± 6.7°), which agrees well with previous molecular dynamics simulation results. In addition, the ensemble structure distribution shows excellent agreement with the two-dimensional free energy surface as a function of WALP23's τ and ρ. These results demonstrate that SSNMR ensemble dynamics provides a means to extract orientational and dynamic information of TM helices from their SSNMR observables and to explain the discrepancy between molecular dynamics simulation and GALA-based interpretation of DQS data.

Transmembrane helix orientation and dynamics: insights from ensemble dynamics with solid-state NMR observables

As the major component of membrane proteins, transmembrane helices embedded in anisotropic bilayer environments adopt preferential orientations that are characteristic or related to their functional states. Recent developments in solid-state nuclear magnetic resonance (SSNMR) spectroscopy have made it possible to measure NMR observables that can be used to determine such orientations in a native bilayer environment. A quasistatic single conformer model is frequently used to interpret the SSNMR observables, but important motional information can be missing or misinterpreted in the model. In this work, we have investigated the orientation of the single-pass transmembrane domain of viral protein "u" (VpuTM) from HIV-1 by determining an ensemble of structures using multiple conformer models based on the SSNMR ensemble dynamics technique. The resulting structure ensemble shows significantly larger orientational fluctuations while the ensemble-averaged orientation is compatible with the orientation based on the quasistatic model. This observation is further corroborated by comparison with the VpuTM orientation from comparative molecular dynamics simulations in explicit bilayer membranes. SSNMR ensemble dynamics not only reveals the importance of transmembrane helix dynamics in interpretation of SSNMR observables, but also provides a means to simultaneously extract both transmembrane helix orientation and dynamics information from the SSNMR measurements.

Folding of small proteins using constrained molecular dynamics

The focus of this paper is to examine whether conformational search using constrained molecular dynamics (MD) method is more enhanced and enriched toward "native-like" structures compared to all-atom MD for the protein folding as a model problem. Constrained MD methods provide an alternate MD tool for protein structure prediction and structure refinement. It is computationally expensive to perform all-atom simulations of protein folding because the processes occur on a time scale of microseconds. Compared to the all-atom MD simulation, constrained MD methods have the advantage that stable dynamics can be achieved for larger time steps and the number of degrees of freedom is an order of magnitude smaller, leading to a decrease in computational cost. We have developed a generalized constrained MD method that allows the user to "freeze and thaw" torsional degrees of freedom as fit for the problem studied. We have used this method to perform all-torsion constrained MD in implicit solvent coupled with the replica exchange method to study folding of small proteins with various secondary structural motifs such as, α-helix (polyalanine, WALP16), β-turn (1E0Q), and a mixed motif protein (Trp-cage). We demonstrate that constrained MD replica exchange method exhibits a wider conformational search than all-atom MD with increased enrichment of near-native structures. "Hierarchical" constrained MD simulations, where the partially formed helical regions in the initial stretch of the all-torsion folding simulation trajectory of Trp-cage were frozen, showed a better sampling of near-native structures than all-torsion constrained MD simulations. This is in agreement with the zipping-and-assembly folding model put forth by Dill and co-workers for folding proteins. The use of hierarchical "freeze and thaw" clustering schemes in constrained MD simulation can be used to sample conformations that contribute significantly to folding of proteins.

Development of a new physics-based internal coordinate mechanics force field and its application to protein loop modeling

We report the development of internal coordinate mechanics force field (ICMFF), new force field parameterized using a combination of experimental data for crystals of small molecules and quantum mechanics calculations. The main features of ICMFF include: (a) parameterization for the dielectric constant relevant to the condensed state (ε = 2) instead of vacuum, (b) an improved description of hydrogen-bond interactions using duplicate sets of van der Waals parameters for heavy atom-hydrogen interactions, and (c) improved backbone covalent geometry and energetics achieved using novel backbone torsional potentials and inclusion of the bond angles at the C(α) atoms into the internal variable set. The performance of ICMFF was evaluated through loop modeling simulations for 4-13 residue loops. ICMFF was combined with a solvent-accessible surface area solvation model optimized using a large set of loop decoys. Conformational sampling was carried out using the biased probability Monte Carlo method. Average/median backbone root-mean-square deviations of the lowest energy conformations from the native structures were 0.25/0.21 Å for four residues loops, 0.84/0.46 Å for eight residue loops, and 1.16/0.73 Å for 12 residue loops. To our knowledge, these results are significantly better than or comparable with those reported to date for any loop modeling method that does not take crystal packing into account. Moreover, the accuracy of our method is on par with the best previously reported results obtained considering the crystal environment. We attribute this success to the high accuracy of the new ICM force field achieved by meticulous parameterization, to the optimized solvent model, and the efficiency of the search method.

Using NMR chemical shifts as structural restraints in molecular dynamics simulations of proteins

We introduce a procedure to determine the structures of proteins by incorporating NMR chemical shifts as structural restraints in molecular dynamics simulations. In this approach, the chemical shifts are expressed as differentiable functions of the atomic coordinates and used to compute forces to generate trajectories that lead to the reduction of the differences between experimental and calculated chemical shifts. We show that this strategy enables the folding of a set of proteins with representative topologies starting from partially denatured initial conformations without the use of additional experimental information. This method also enables the straightforward combination of chemical shifts with other standard NMR restraints, including those derived from NOE, J-coupling, and residual dipolar coupling measurements. We illustrate this aspect by calculating the structure of a transiently populated excited state conformation from chemical shift and residual dipolar coupling data measured by relaxation dispersion NMR experiments.

Bioactive conformational generation of small molecules: a comparative analysis between force-field and multiple empirical criteria based methods

Background

Conformational sampling for small molecules plays an essential role in drug discovery research pipeline. Based on multi-objective evolution algorithm (MOEA), we have developed a conformational generation method called Cyndi in the previous study. In this work, in addition to Tripos force field in the previous version, Cyndi was updated by incorporation of MMFF94 force field to assess the conformational energy more rationally. With two force fields against a larger dataset of 742 bioactive conformations of small ligands extracted from PDB, a comparative analysis was performed between pure force field based method (FFBM) and multiple empirical criteria based method (MECBM) hybrided with different force fields.

Results

Our analysis reveals that incorporating multiple empirical rules can significantly improve the accuracy of conformational generation. MECBM, which takes both empirical and force field criteria as the objective functions, can reproduce about 54% (within 1Å RMSD) of the bioactive conformations in the 742-molecule testset, much higher than that of pure force field method (FFBM, about 37%). On the other hand, MECBM achieved a more complete and efficient sampling of the conformational space because the average size of unique conformations ensemble per molecule is about 6 times larger than that of FFBM, while the time scale for conformational generation is nearly the same as FFBM. Furthermore, as a complementary comparison study between the methods with and without empirical biases, we also tested the performance of the three conformational generation methods in MacroModel in combination with different force fields. Compared with the methods in MacroModel, MECBM is more competitive in retrieving the bioactive conformations in light of accuracy but has much lower computational cost.

Conclusions

By incorporating different energy terms with several empirical criteria, the MECBM method can produce more reasonable conformational ensemble with high accuracy but approximately the same computational cost in comparison with FFBM method. Our analysis also reveals that the performance of conformational generation is irrelevant to the types of force field adopted in characterization of conformational accessibility. Moreover, post energy minimization is not necessary and may even undermine the diversity of conformational ensemble. All the results guide us to explore more empirical criteria like geometric restraints during the conformational process, which may improve the performance of conformational generation in combination with energetic accessibility, regardless of force field types adopted.

Prediction of protein-ligand binding affinity by free energy simulations: assumptions, pitfalls and expectations

Many limitations of current computer-aided drug design arise from the difficulty of reliably predicting the binding affinity of a small molecule to a biological target. There is thus a strong interest in novel computational methodologies that claim predictions of greater accuracy than current scoring functions, and at a throughput compatible with the rapid pace of drug discovery in the pharmaceutical industry. Notably, computational methodologies firmly rooted in statistical thermodynamics have received particular attention in recent years. Yet free energy calculations can be daunting to learn for a novice user because of numerous technical issues and various approaches advocated by experts in the field. The purpose of this article is to provide an overview of the current capabilities of free energy calculations and to discuss the applicability of this technology to drug discovery.

Implementation of the CHARMM Force Field in GROMACS: Analysis of Protein Stability Effects from Correction Maps, Virtual Interaction Sites, and Water Models

CHARMM27 is a widespread and popular force field for biomolecular simulation, and several recent algorithms such as implicit solvent models have been developed specifically for it. We have here implemented the CHARMM force field and all necessary extended functional forms in the GROMACS molecular simulation package, to make CHARMM-specific features available and to test them in combination with techniques for extended time steps, to make all major force fields available for comparison studies in GROMACS, and to test various solvent model optimizations, in particular the effect of Lennard-Jones interactions on hydrogens. The implementation has full support both for CHARMM-specific features such as multiple potentials over the same dihedral angle and the grid-based energy correction map on the ϕ, ψ protein backbone dihedrals, as well as all GROMACS features such as virtual hydrogen interaction sites that enable 5 fs time steps. The medium-to-long time effects of both the correction maps and virtual sites have been tested by performing a series of 100 ns simulations using different models for water representation, including comparisons between CHARMM and traditional TIP3P. Including the correction maps improves sampling of near native-state conformations in our systems, and to some extent it is even able to refine distorted protein conformations. Finally, we show that this accuracy is largely maintained with a new implicit solvent implementation that works with virtual interaction sites, which enables performance in excess of 250 ns/day for a 900-atom protein on a quad-core desktop computer.

An Evaluation of Explicit Receptor Flexibility in Molecular Docking Using Molecular Dynamics and Torsion Angle Molecular Dynamics

Incorporating receptor flexibility into molecular docking should improve results for flexible proteins. However, the incorporation of explicit all-atom flexibility with molecular dynamics for the entire protein chain may also introduce significant error and "noise" that could decrease docking accuracy and deteriorate the ability of a scoring function to rank native-like poses. We address this apparent paradox by comparing the success of several flexible receptor models in cross-docking and multiple receptor ensemble docking for p38α mitogen-activated protein (MAP) kinase. Explicit all-atom receptor flexibility has been incorporated into a CHARMM-based molecular docking method (CDOCKER) using both molecular dynamics (MD) and torsion angle molecular dynamics (TAMD) for the refinement of predicted protein-ligand binding geometries. These flexible receptor models have been evaluated, and the accuracy and efficiency of TAMD sampling is directly compared to MD sampling. Several flexible receptor models are compared, encompassing flexible side chains, flexible loops, multiple flexible backbone segments, and treatment of the entire chain as flexible. We find that although including side chain and some backbone flexibility is required for improved docking accuracy as expected, docking accuracy also diminishes as additional and unnecessary receptor flexibility is included into the conformational search space. Ensemble docking results demonstrate that including protein flexibility leads to to improved agreement with binding data for 227 active compounds. This comparison also demonstrates that a flexible receptor model enriches high affinity compound identification without significantly increasing the number of false positives from low affinity compounds.

Computational multiscale modeling in protein--ligand docking

In biological systems, the binding of small molecule ligands to proteins is a crucial process for almost every aspect of biochemistry and molecular biology. Enzymes are proteins that function by catalyzing specific biochemical reactions that convert reactants into products. Complex organisms are typically composed of cells in which thousands of enzymes participate in complex and interconnected biochemical pathways. Some enzymes serve as sequential steps in specific pathways (such as energy metabolism), while others function to regulate entire pathways and cellular functions [1]. Small molecule ligands can be designed to bind to a specific enzyme and inhibit the biochemical reaction. Inhibiting the activity of key enzymes may result in the entire biochemical pathways being turned on or off [2], [3]. Many small molecule drugs marketed today function in this generic way as enzyme inhibitors. If research identifies a specific enzyme as being crucial to the progress of disease, then this enzyme may be targeted with an inhibitor, which may slow down or reverse the progress of disease. In this way, enzymes are targeted from specific pathogens (e.g., virus, bacteria, fungi) for infectious diseases [4], [5], and human enzymes are targeted for noninfectious diseases such as cardiovascular disease, cancer, diabetes, and neurodegenerative diseases [6].

Beta-hairpin restraint potentials for calculations of potentials of mean force as a function of beta-hairpin tilt, rotation, and distance

We have developed a set of restraint potentials for beta-hairpin tilt relative to the membrane normal, beta-hairpin rotation around the beta-hairpin axis, and hairpin-hairpin distance. Such restraint potentials enable us to characterize the molecular basis of specific beta-hairpin tilt and rotation in membranes and hairpin-hairpin interactions at the atomic level by sampling their conformational space along these degrees of freedom, i.e., reaction coordinates, during molecular dynamics simulations. We illustrate the efficacy of the beta-hairpin restraint potentials by calculating the potentials of mean force (PMFs) as a function of tilt and rotation angles of protegrin-1 (PG-1), a beta-hairpin antimicrobial peptide, in an implicit membrane model. The peptide association in the membrane is also examined by calculating the PMFs as a function of distance between two PG-1 peptides in various dimer interfaces. These novel restraint potentials are found to perform well in each of these cases and are expected to be a useful means to study the microscopic driving forces of insertion, tilting, and rotation of beta-hairpin peptides in membranes as well as their association in aqueous solvent or membrane environments particularly when combined with explicit solvent models.

Molecular dynamics-based approaches for enhanced sampling of long-time, large-scale conformational changes in biomolecules

The rugged energy landscape of biomolecules together with shortcomings of traditional molecular dynamics (MD) simulations require specialized methods for capturing large-scale, long-time configurational changes along with chemical dynamics behavior. In this report, MD-based methods for biomolecules are surveyed, involving modification of the potential, simulation protocol, or algorithm as well as global reformulations. While many of these methods are successful at probing the thermally accessible configuration space at the expense of altered kinetics, more sophisticated approaches like transition path sampling or Markov chain models are required to obtain mechanistic information, reaction pathways, and/or reaction rates. Divide-and-conquer methods for sampling and for piecing together reaction rate information are especially suitable for readily available computer cluster networks. Successful applications to biomolecules remain a challenge.

Monte Carlo, harmonic approximation, and coarse-graining approaches for enhanced sampling of biomolecular structure

The rugged energy landscape of biomolecules and associated large-scale conformational changes have triggered the development of many innovative enhanced sampling methods, either based or not based on molecular dynamics (MD) simulations. Surveyed here are methods in the latter class - including Monte Carlo methods, harmonic approximations, and coarse graining - many of which yield valuable conformational insights into biomolecular structure and flexibility, despite altered kinetics. MD-based methods are surveyed in an upcoming issue of F1000 Biology Reports.

Principles of flexible protein-protein docking

Treating flexibility in molecular docking is a major challenge in cell biology research. Here we describe the background and the principles of existing flexible protein-protein docking methods, focusing on the algorithms and their rational. We describe how protein flexibility is treated in different stages of the docking process: in the preprocessing stage, rigid and flexible parts are identified and their possible conformations are modeled. This preprocessing provides information for the subsequent docking and refinement stages. In the docking stage, an ensemble of pre-generated conformations or the identified rigid domains may be docked separately. In the refinement stage, small-scale movements of the backbone and side-chains are modeled and the binding orientation is improved by rigid-body adjustments. For clarity of presentation, we divide the different methods into categories. This should allow the reader to focus on the most suitable method for a particular docking problem.

Application of solid-state NMR restraint potentials in membrane protein modeling

We have developed a set of orientational restraint potentials for solid-state NMR observables including (15)N chemical shift and (15)N-(1)H dipolar coupling. Torsion angle molecular dynamics simulations with available experimental (15)N chemical shift and (15)N-(1)H dipolar coupling as target values have been performed to determine orientational information of four membrane proteins and to model the structures of some of these systems in oligomer states. The results suggest that incorporation of the orientational restraint potentials into molecular dynamics provides an efficient means to the determination of structures that optimally satisfy the experimental observables without an extensive geometrical search.

Computational techniques for efficient conformational sampling of proteins

In this review, we summarize the computational methods for sampling the conformational space of biomacromolecules. We discuss the methods applicable to find only lowest energy conformations (global minimization of the potential-energy function) and to generate canonical ensembles (canonical Monte Carlo method and canonical molecular dynamics method and their extensions). Special attention is devoted to the use of coarse-grained models that enable simulations to be enhanced by several orders of magnitude.

CAESAR: a new conformer generation algorithm based on recursive buildup and local rotational symmetry consideration

A highly efficient conformer search algorithm based on a divide-and-conquer and recursive conformer build-up approach is presented in this paper. This approach is combined with consideration of local rotational symmetry so that conformer duplicates due to topological symmetry in the systematic search can be efficiently eliminated. This new algorithm, termed CAESAR (Conformer Algorithm based on Energy Screening and Recursive Buildup), has been implemented in Discovery Studio 1.7 as part of the Catalyst Component Collection. CAESAR has been validated by comparing the conformer models generated by the new method and Catalyst/FAST. CAESAR is consistently 5-20 times faster than Catalyst/FAST for all data sets investigated. The speedup is even more dramatic for molecules with high topological symmetry or for molecules that require a large number of conformers to be sampled. The quality of the conformer models generated by CAESAR has been validated by assessing the ability to reproduce the receptor-bound X-ray conformation of ligands extracted for the Protein Data Bank (PDB) and assessing the ability to adequately cover the pharmacophore space. It is shown that CAESAR is able to reproduce the receptor-bound conformation slightly better than the Catalyst/FAST method for a data set of 918 ligands retrieved from the PDB. In addition, it is shown that CEASAR covers the pharmacophore space as well or better than Catalyst/FAST.