Energetically unfavorable protein angles: Exploration of a conserved dihedral angle in triosephosphate isomerase.

Over the past 3.5 billion years of evolution, enzymes have adopted a myriad of conformations to suit life on earth. However, torsional angles of proteins have settled into limited zones of energetically favorable dihedrals observed in Ramachandran plots. Areas outside said zones are believed to be disallowed to all amino acids, except glycine, due to steric hindrance. Triosephosphate isomerase (TIM), a homodimer with a catalytic rate approaching the diffusion limit, contains an active site lysine residue (K13) with dihedrals within the fourth quadrant (Φ = +51/Ψ = -143). Both the amino acid and the dihedral angles are conserved across all species of TIM and known crystal structures regardless of ligand. Only crystal structures of the engineered monomeric version (1MSS) show accepted ß-sheet dihedral values of Φ = -135/Ψ = +170 but experiments show a 1000-fold loss in activity. Based on these results, we hypothesized that adopting the unfavorable torsion angle for K13 contributes to catalysis. Using both, computational and experimental approaches, four residues that interact with K13 (N11, M14, E97, and Q64) were mutated to alanine. In silico molecular dynamics (MD) simulations were performed using 2JK2 unliganded human TIM as a starting structure. Ramachandran plots, containing K13 dihedral values reveal full or partial loss of disallowed zone angles. N11A showed no detectable catalytic activity and lost the unfavorable K13 dihedral angles across four separate force fields during simulation while all other mutants plus wild type retained activity and retained the conserved K13 dihedral angles.

A Comprehensive Review of Cholinesterase Modeling and Simulation.

The present article reviews published efforts to study acetylcholinesterase and butyrylcholinesterase structure and function using computer-based modeling and simulation techniques. Structures and models of both enzymes from various organisms, including rays, mice, and humans, are discussed to highlight key structural similarities in the active site gorges of the two enzymes, such as flexibility, binding site location, and function, as well as differences, such as gorge volume and binding site residue composition. Catalytic studies are also described, with an emphasis on the mechanism of acetylcholine hydrolysis by each enzyme and novel mutants that increase catalytic efficiency. The inhibitory activities of myriad compounds have been computationally assessed, primarily through Monte Carlo-based docking calculations and molecular dynamics simulations. Pharmaceutical compounds examined herein include FDA-approved therapeutics and their derivatives, as well as several other prescription drug derivatives. Cholinesterase interactions with both narcotics and organophosphate compounds are discussed, with the latter focusing primarily on molecular recognition studies of potential therapeutic value and on improving our understanding of the reactivation of cholinesterases that are bound to toxins. This review also explores the inhibitory properties of several other organic and biological moieties, as well as advancements in virtual screening methodologies with respect to these enzymes.

Overcoming the Heuristic Nature of k-Means Clustering: Identification and Characterization of Binding Modes from Simulations of Molecular Recognition Complexes.

The accurate and reproducible detection and description of thermodynamic states in computational data is a nontrivial problem, particularly when the number of states is unknown a priori and for large, flexible chemical systems and complexes. To this end, we report a novel clustering protocol that combines high-resolution structural representation, brute-force repeat clustering, and optimization of clustering statistics to reproducibly identify the number of clusters present in a data set (k) for simulated ensembles of butyrylcholinesterase in complex with two previously studied organophosphate inhibitors. Each structure within our simulated ensembles was depicted as a high-dimensionality vector with components defined by specific protein-inhibitor contacts at the chemical group level and the magnitudes of these components defined by their respective extents of pair-wise atomic contact, thus allowing for algorithmic differentiation between varying degrees of interaction. These surface-weighted interaction fingerprints were tabulated for each of over 1 million structures from more than 100 µs of all-atom molecular dynamics simulation per complex and used as the input for repetitive k-means clustering. Minimization of cluster population variance and range afforded accurate and reproducible identification of k, thereby allowing for the characterization of discrete binding modes from molecular simulation data in the form of contact tables that concisely encapsulate the observed intermolecular contact motifs. While the protocol presented herein to determine k and achieve non-heuristic clustering is demonstrated on data from massive atomistic simulation, our approach is generalizable to other data types and clustering algorithms, and is tractable with limited computational resources.

Understanding the enzyme-ligand complex: insights from all-atom simulations of butyrylcholinesterase inhibition.

All-atom molecular dynamics simulations of butyrylcholinesterase (BChE) sans inhibitor and in complex with each of 15 dialkyl phenyl phosphate derivatives were conducted to characterize inhibitor binding modes and strengths. Each system was sampled on the 250 ns timescale in explicit ionic solvent, for a total of over 4 µs of simulation time. A K-means algorithm was used to cluster the resulting structures into distinct binding modes, which were further characterized based on atomic-level contacts between inhibitor chemical groups and active site residues. Comparison of experimentally observed inhibition constants (KI) with the resulting contact tables provides structural explanations for relative binding coefficients and highlights several notable interaction motifs. These include ubiquitous contact between glycines in the oxyanion hole and the inhibitor phosphate group; a sterically driven binding preference for positional isomers that extend aromaticity; a stereochemical binding preference for choline-containing inhibitors, which mimic natural BChE substrates; and the mechanically induced opening of the omega loop region to fully expose the active site gorge in the presence of choline-containing inhibitors. Taken together, these observations can greatly inform future design of BChE inhibitors, and the approach reported herein is generalizable to other enzyme-inhibitor systems and similar complexes that depend on non-covalent molecular recognition.Communicated by Ramaswamy H. Sarma.

Ensemble Molecular Dynamics of a Protein-Ligand Complex: Residual Inhibitor Entropy Enhances Drug Potency in Butyrylcholinesterase.

Butyrylcholinesterase is a key enzyme that catalyzes the hydrolysis of the neurotransmitter acetylcholine and shows an increased activity in patients suffering from Alzheimer's disease (AD), making this enzyme a primary target in treating AD. Central to this problem, and to similar scenarios involving biomolecular recognition, is our understanding of the nature of the protein-ligand complex. The butyrylcholinesterase enzyme was studied via all-atom, explicit solvent, ensemble molecular dynamics simulations sans inhibitor and in the presence of three dialkyl phenyl phosphate inhibitors of known potency to a cumulative sampling of over 40 µs. Following the relaxation of these ensembles to conformational equilibria, binding modes for each inhibitor were identified. While classical models, which assume significant reduction in protein and ligand conformational entropies, continue to be favored in contemporary studies, our observations contradict those assumptions: bound ligands occupy many conformational states, thereby stabilizing the complex, while also promoting protein flexibility.

Synthesis, biochemical evaluation, and molecular modeling studies of aryl and arylalkyl di-n-butyl phosphates, effective butyrylcholinesterase inhibitors.

A series of dialkyl aryl phosphates and dialkyl arylalkyl phosphates were synthesized. Their inhibitory activities were evaluated against acetylcholinesterase (AChE) and butyrylcholinesterase (BChE). The di-n-butyl phosphate series consistently displayed selective inhibition of BChE over AChE. The most potent inhibitors of butyrylcholinesterase were di-n-butyl-3,5-dimethylphenyl phosphate (4b) [KI=1.0±0.4µM] and di-n-butyl 2-naphthyl phosphate (5b) [KI=1.9±0.4µM]. Molecular modeling was used to uncover three subsites within the active site gorge that accommodate the three substituents attached to the phosphate group. Phosphates 4b and 5b were found to bind to these three subsites in analogous fashion with the aromatic groups in both analogs being accommodated by the "lower region," while the lone pairs on the PO oxygen atoms were oriented towards the oxyanion hole. In contrast, di-n-butyl-3,4-dimethylphenyl phosphate (4a) [KI=9±1µM], an isomer of 4b, was found to orient its aromatic group in the "upper left region" subsite as placement of this group in the "lower region" resulted in significant steric hindrance by a ridge-like region in this subsite. Future studies will be designed to exploit these features in an effort to develop inhibitors of higher inhibitory strength against butyrylcholinesterase.

Ensemble simulations: folding, unfolding and misfolding of a high-efficiency frameshifting RNA pseudoknot.

Massive all-atom molecular dynamics simulations were conducted across a distributed computing network to study the folding, unfolding, misfolding and conformational plasticity of the high-efficiency frameshifting double mutant of the 26 nt potato leaf roll virus RNA pseudoknot. Our robust sampling, which included over 40 starting structures spanning the spectrum from the extended unfolded state to the native fold, yielded nearly 120 µs of cumulative sampling time. Conformational microstate transitions on the 1.0 ns to 10.0 µs timescales were observed, with post-equilibration sampling providing detailed representations of the conformational free energy landscape and the complex folding mechanism inherent to the pseudoknot motif. Herein, we identify and characterize two alternative native structures, three intermediate states, and numerous misfolded states, the latter of which have not previously been characterized via atomistic simulation techniques. While in line with previous thermodynamics-based models of a general RNA folding mechanism, our observations indicate that stem-strand-sequence-separation may serve as an alternative predictor of the order of stem formation during pseudoknot folding. Our results contradict a model of frameshifting based on structural rigidity and resistance to mechanical unfolding, and instead strongly support more recent studies in which conformational plasticity is identified as a determining factor in frameshifting efficiency.

The extent of pyrene excimer fluorescence emission is a reflector of distance and flexibility: analysis of the segment linking the LDL receptor-binding and tetramerization domains of apolipoprotein E3.

Pyrene is a spatially sensitive probe that displays an ensemble of monomeric fluorescence emission peaks (375-405 nm) and an additional band (called excimer) at ~460 nm when two fluorophores are spatially proximal. We examined if there is a correlation between distance between two pyrenes on an α-helical structure and excimer/monomer (e/m) ratio. Using structure-guided design, pyrene maleimide was attached to pairs of Cys residues separated by ~5 Å increments on helix 2 of the N-terminal domain of apolipoprotein E3 (apoE3). Fluorescence spectral analysis revealed an intense excimer band when the probes were ~5 Å from each other with an e/m ratio of ~3.0, which decreased to ~1.0 at 20 Å. An inverse correlation between e/m ratio and the distance between pyrenes was observed, with the probe and helix flexibility also contributing to the extent of excimer formation. We verified this approach by estimating the distance between T57C and C112 (located on helices 2 and 3, respectively) to be 5.2 Å (4.9 Å from NMR and 5.7 Å from the X-ray structure). Excimer formation was also noted to a significant extent with probes located in the linker segment, suggesting spatial proximity (10-15 Å) to corresponding sites on neighboring molecules in the tetrameric configuration of apoE. We infer that oligomerization via the C-terminal domain juxtaposes the linker segments from neighboring apoE molecules. This study offers new insights into the conformation of tetrameric apoE and presents the use of pyrene as a powerful probe for studying protein spatial organization.

Evaluating molecular mechanical potentials for helical peptides and proteins.

Multiple variants of the AMBER all-atom force field were quantitatively evaluated with respect to their ability to accurately characterize helix-coil equilibria in explicit solvent simulations. Using a global distributed computing network, absolute conformational convergence was achieved for large ensembles of the capped A(21) and F(s) helical peptides. Further assessment of these AMBER variants was conducted via simulations of a flexible 164-residue five-helix-bundle protein, apolipophorin-III, on the 100 ns timescale. Of the contemporary potentials that had not been assessed previously, the AMBER-99SB force field showed significant helix-destabilizing tendencies, with beta bridge formation occurring in helical peptides, and unfolding of apolipophorin-III occurring on the tens of nanoseconds timescale. The AMBER-03 force field, while showing adequate helical propensities for both peptides and stabilizing apolipophorin-III, (i) predicts an unexpected decrease in helicity with ALA-->ARG(+) substitution, (ii) lacks experimentally observed 3(10) helical content, and (iii) deviates strongly from average apolipophorin-III NMR structural properties. As is observed for AMBER-99SB, AMBER-03 significantly overweighs the contribution of extended and polyproline backbone configurations to the conformational equilibrium. In contrast, the AMBER-99phi force field, which was previously shown to best reproduce experimental measurements of the helix-coil transition in model helical peptides, adequately stabilizes apolipophorin-III and yields both an average gyration radius and polar solvent exposed surface area that are in excellent agreement with the NMR ensemble.

Equilibrium conformational dynamics in an RNA tetraloop from massively parallel molecular dynamics.

Conformational equilibrium within the ubiquitous GNRA tetraloop motif was simulated at the ensemble level, including 10 000 independent all-atom molecular dynamics trajectories totaling over 110 micros of simulation time. This robust sampling reveals a highly dynamic structure comprised of 15 conformational microstates. We assemble a Markov model that includes transitions ranging from the nanosecond to microsecond timescales and is dominated by six key loop conformations that contribute to fluctuations around the native state. Mining of the Protein Data Bank provides an abundance of structures in which GNRA tetraloops participate in tertiary contact formation. Most predominantly observed in the experimental data are interactions of the native loop structure within the minor groove of adjacent helical regions. Additionally, a second trend is observed in which the tetraloop assumes non-native conformations while participating in multiple tertiary contacts, in some cases involving multiple possible loop conformations. This tetraloop flexibility can act to counterbalance the energetic penalty associated with assuming non-native loop structures in forming tertiary contacts. The GNRA motif has thus evolved not only to readily participate in simple tertiary interactions involving native loop structure, but also to easily adapt tetraloop secondary conformation in order to participate in larger, more complex tertiary interactions.

Molecular simulation of multistate peptide dynamics: a comparison between microsecond timescale sampling and multiple shorter trajectories.

Molecular dynamics simulations of the RN24 peptide, which includes a diverse set of structurally heterogeneous states, are carried out in explicit solvent. Two approaches are employed and compared directly under identical simulation conditions. Specifically, we examine sampling by two individual long trajectories (microsecond timescale) and many shorter (MS) uncoupled trajectories. Statistical analysis of the structural properties indicates a qualitative agreement between these approaches. Microsecond timescale sampling gives large uncertainties on most structural metrics, while the shorter timescale of MS simulations results in slight structural memory for beta-structure starting states. Additionally, MS sampling detects numerous transitions on a relatively short timescale that are not observed in microsecond sampling, while long simulations allow for detection of a few transitions on significantly longer timescales. A correlation between the complex free energy landscape and the kinetics of the equilibrium is highlighted by principal component analysis on both simulation sets. This report highlights the increased precision of the MS approach when studying the kinetics of complex conformational change, while revealing the complementary insight and qualitative agreement offered by far fewer individual simulations on significantly longer timescales.

Nanotube confinement denatures protein helices.

In striking contrast to simple polymer physics theory, which does not account for solvent effects, we find that physical confinement of solvated biopolymers decreases solvent entropy, which in turn leads to a reduction in the organized structural content of the polymer. Since our theory is based on a fundamental property of water-protein statistical mechanics, we expect it to have broad implications in many biological and material science contexts.

The solvation interface is a determining factor in peptide conformational preferences.

The 21 residue polyalanine-based F(s) peptide was studied using thousands of long, explicit solvent, atomistic molecular dynamics simulations that reached equilibrium at the ensemble level. Peptide conformational preference as a function of hydrophobicity was examined using a spectrum of explicit solvent models, and the peptide length-dependence of the hydrophilic and hydrophobic components of solvent-accessible surface area for several ideal conformational types was considered. Our results demonstrate how the character of the solvation interface induces several conformational preferences, including a decrease in mean helical content with increased hydrophilicity, which occurs predominantly through reduced nucleation tendency and, to a lesser extent, destabilization of helical propagation. Interestingly, an opposing effect occurs through increased propensity for 3(10)-helix conformations, as well as increased polyproline structure. Our observations provide a framework for understanding previous reports of conformational preferences in polyalanine-based peptides including (i) terminal 3(10)-helix prominence, (ii) low pi-helix propensity, (iii) increased polyproline conformations in short and unfolded peptides, and (iv) membrane helix stability in the presence and absence of water. These observations provide physical insight into the role of water in peptide conformational equilibria at the atomic level, and expand our view of the complexity of even the most "simple" of biopolymers. Whereas previous studies have focused predominantly on hydrophobic effects with respect to tertiary structure, this work highlights the need for consideration of such effects at the secondary structural level.

Direct calculation of the binding free energies of FKBP ligands.

Direct calculations of the absolute free energies of binding for eight ligands to FKBP protein were performed using the Fujitsu BioServer massively parallel computer. Using the latest version of the general assisted model building with energy refinement (AMBER) force field for ligand model parameters and the Bennett acceptance ratio for computing free-energy differences, we obtained an excellent linear fit between the calculated and experimental binding free energies. The rms error from a linear fit is 0.4 kcal/mol for eight ligand complexes. In comparison with a previous study of the binding energies of these same eight ligand complexes, these results suggest that the use of improved model parameters can lead to more predictive binding estimates, and that these estimates can be obtained with significantly less computer time than previously thought. These findings make such direct methods more attractive for use in rational drug design.

Unusual compactness of a polyproline type II structure.

Polyproline type II (PPII) helix has emerged recently as the dominant paradigm for describing the conformation of unfolded polypeptides. However, most experimental observables used to characterize unfolded proteins typically provide only short-range, sequence-local structural information that is both time- and ensemble-averaged, giving limited detail about the long-range structure of the chain. Here, we report a study of a long-range property: the radius of gyration of an alanine-based peptide, Ace-(diaminobutyric acid)2-(Ala)7-(ornithine)2-NH2. This molecule has previously been studied as a model for the unfolded state of proteins under folding conditions and is believed to adopt a PPII fold based on short-range techniques such as NMR and CD. By using synchrotron radiation and small-angle x-ray scattering, we have determined the radius of gyration of this peptide to be 7.4 +/- 0.5 angstroms, which is significantly less than the value expected from an ideal PPII helix in solution (13.1 angstroms). To further study this contradiction, we have used molecular dynamics simulations using six variants of the AMBER force field and the GROMOS 53A6 force field. However, in all cases, the simulated ensembles underestimate the PPII content while overestimating the experimental radius of gyration. The conformational model that we propose, based on our small angle x-ray scattering results and what is known about this molecule from before, is that of a very flexible, fluctuating structure that on the level of individual residues explores a wide basin around the ideal PPII geometry but is never, or only rarely, in the ideal extended PPII helical conformation.

How well can simulation predict protein folding kinetics and thermodynamics?

Simulation of protein folding has come a long way in five years. Notably, new quantitative comparisons with experiments for small, rapidly folding proteins have become possible. As the only way to validate simulation methodology, this achievement marks a significant advance. Here, we detail these recent achievements and ask whether simulations have indeed rendered quantitative predictions in several areas, including protein folding kinetics, thermodynamics, and physics-based methods for structure prediction. We conclude by looking to the future of such comparisons between simulations and experiments.

Empirical force-field assessment: The interplay between backbone torsions and noncovalent term scaling.

The kinetic and thermodynamic aspects of the helix-coil transition in polyalanine-based peptides have been studied at the ensemble level using a distributed computing network. This study builds on a previous report, which critically assessed the performance of several contemporary force fields in reproducing experimental measurements and elucidated the complex nature of helix-coil systems. Here we consider the effects of modifying backbone torsions and the scaling of noncovalent interactions. Although these elements determine the potential of mean force between atoms separated by three covalent bonds (and thus largely determine the local conformational distributions observed in simulation), we demonstrate that the interplay between these factors is both complex and force field dependent. We quantitatively assess the heliophilicity of several helix-stabilizing potentials as well as the changes in heliophilicity resulting from such modifications, which can "make or break" the accuracy of a given force field, and our findings suggests that future force field development may need to better consider effect that vary with peptide length. This report also serves as an example of the utility of distributed computing in analyzing and improving upon contemporary force fields at the level of absolute ensemble equilibrium, the next step in force field development.

Does water play a structural role in the folding of small nucleic acids?

Nucleic acid structure and dynamics are known to be closely coupled to local environmental conditions and, in particular, to the ionic character of the solvent. Here we consider what role the discrete properties of water and ions play in the collapse and folding of small nucleic acids. We study the folding of an experimentally well-characterized RNA hairpin-loop motif (sequence 5'-GGGC[GCAA]GCCU-3') via ensemble molecular dynamics simulation and, with nearly 500 micros of aggregate simulation time using an explicit representation of the ionic solvent, report successful ensemble folding simulations with a predicted folding time of 8.8(+/-2.0) micros, in agreement with experimental measurements of approximately 10 micros. Comparing our results to previous folding simulations using the GB/SA continuum solvent model shows that accounting for water-mediated interactions is necessary to accurately characterize the free energy surface and stochastic nature of folding. The formation of the secondary structure appears to be more rapid than the fastest ionic degrees of freedom, and counterions do not participate discretely in observed folding events. We find that hydrophobic collapse follows a predominantly expulsive mechanism in which a diffusion-search of early structural compaction is followed by the final formation of native structure that occurs in tandem with solvent evacuation.

Exploring the helix-coil transition via all-atom equilibrium ensemble simulations.

The ensemble folding of two 21-residue alpha-helical peptides has been studied using all-atom simulations under several variants of the AMBER potential in explicit solvent using a global distributed computing network. Our extensive sampling, orders of magnitude greater than the experimental folding time, results in complete convergence to ensemble equilibrium. This allows for a quantitative assessment of these potentials, including a new variant of the AMBER-99 force field, denoted AMBER-99 phi, which shows improved agreement with experimental kinetic and thermodynamic measurements. From bulk analysis of the simulated AMBER-99 phi equilibrium, we find that the folding landscape is pseudo-two-state, with complexity arising from the broad, shallow character of the "native" and "unfolded" regions of the phase space. Each of these macrostates allows for configurational diffusion among a diverse ensemble of conformational microstates with greatly varying helical content and molecular size. Indeed, the observed structural dynamics are better represented as a conformational diffusion than as a simple exponential process, and equilibrium transition rates spanning several orders of magnitude are reported. After multiple nucleation steps, on average, helix formation proceeds via a kinetic "alignment" phase in which two or more short, low-entropy helical segments form a more ideal, single-helix structure.

Simulations of the role of water in the protein-folding mechanism.

There are many unresolved questions regarding the role of water in protein folding. Does water merely induce hydrophobic forces, or does the discrete nature of water play a structural role in folding? Are the nonadditive aspects of water important in determining the folding mechanism? To help to address these questions, we have performed simulations of the folding of a model protein (BBA5) in explicit solvent. Starting 10,000 independent trajectories from a fully unfolded conformation, we have observed numerous folding events, making this work a comprehensive study of the kinetics of protein folding starting from the unfolded state and reaching the folded state and with an explicit solvation model and experimentally validated rates. Indeed, both the raw TIP3P folding rate (4.5 +/- 2.5 micros) and the diffusion-constant corrected rate (7.5 +/- 4.2 micros) are in strong agreement with the experimentally observed rate of 7.5 +/- 3.5 micros. To address the role of water in folding, the mechanism is compared with that predicted from implicit solvation simulations. An examination of solvent density near hydrophobic groups during folding suggests that in the case of BBA5, there are water-induced effects not captured by implicit solvation models, including signs of a "concurrent mechanism" of core collapse and desolvation.