Protein conformational landscapes: energy minimization and clustering of a long molecular dynamics trajectory

k-Means NANI: An Improved Clustering Algorithm for Molecular Dynamics Simulations

One of the key challenges of k-means clustering is the seed selection or the initial centroid estimation since the clustering result depends heavily on this choice. Alternatives such as k-means++ have mitigated this limitation by estimating the centroids using an empirical probability distribution. However, with high-dimensional and complex data sets such as those obtained from molecular simulation, k-means++ fails to partition the data in an optimal manner. Furthermore, stochastic elements in all flavors of k-means++ will lead to a lack of reproducibility. K-means N-Ary Natural Initiation (NANI) is presented as an alternative to tackle this challenge by using efficient n-ary comparisons to both identify high-density regions in the data and select a diverse set of initial conformations. Centroids generated from NANI are not only representative of the data and different from one another, helping k-means to partition the data accurately, but also deterministic, providing consistent cluster populations across replicates. From peptide and protein folding molecular simulations, NANI was able to create compact and well-separated clusters as well as accurately find the metastable states that agree with the literature. NANI can cluster diverse data sets and be used as a standalone tool or as part of our MDANCE clustering package.

k-Means NANI: an improved clustering algorithm for Molecular Dynamics simulations

One of the key challenges of k-means clustering is the seed selection or the initial centroid estimation since the clustering result depends heavily on this choice. Alternatives such as k-means++ have mitigated this limitation by estimating the centroids using an empirical probability distribution. However, with high-dimensional and complex datasets such as those obtained from molecular simulation, k-means++ fails to partition the data in an optimal manner. Furthermore, stochastic elements in all flavors of k-means++ will lead to a lack of reproducibility. K-means N-Ary Natural Initiation (NANI) is presented as an alternative to tackle this challenge by using efficient n-ary comparisons to both identify high-density regions in the data and select a diverse set of initial conformations. Centroids generated from NANI are not only representative of the data and different from one another, helping k-means to partition the data accurately, but also deterministic, providing consistent cluster populations across replicates. From peptide and protein folding molecular simulations, NANI was able to create compact and well-separated clusters as well as accurately find the metastable states that agree with the literature. NANI can cluster diverse datasets and be used as a standalone tool or as part of our MDANCE clustering package.

Embedding Nonrigid Solutes in an Averaged Environment: A Case Study on Rhodopsins

Many simulation methods concerning solvated molecules are based on the assumption that the solvated species and the solvent can be characterized by some representative structures of the solute and some embedding potential corresponding to this structure. While the averaging of the solvent configurations to obtain an embedding potential has been studied in great detail, this hinges on a single solute structure representation. This assumption is re-examined and generalized for conformationally flexible solutes and tested on 4 nonrigid systems. In this generalized approach, the solute is characterized by a set of representative structures and the corresponding embedding potentials. The representative structures are identified by means of subdividing the statistical ensemble, which in this work is generated by a constant-temperature molecular dynamics simulation. The embedding potential defined in the Frozen-Density Embedding Theory is used to characterize the average effect of the solvent in each subensemble. The numerical examples concern the vertical excitation energies of protonated retinal Schiff bases in protein environments. It is comprehensively shown that subensemble averaging leads to huge computational savings compared with explicit averaging of the excitation energies in the whole ensemble while introducing only minor errors in the case of the systems examined.

An Expedited Route to Optical and Electronic Properties at Finite Temperature via Unsupervised Learning

Electronic properties and absorption spectra are the grounds to investigate molecular electronic states and their interactions with the environment. Modeling and computations are required for the molecular understanding and design strategies of photo-active materials and sensors. However, the interpretation of such properties demands expensive computations and dealing with the interplay of electronic excited states with the conformational freedom of the chromophores in complex matrices (i.e., solvents, biomolecules, crystals) at finite temperature. Computational protocols combining time dependent density functional theory and ab initio molecular dynamics (MD) have become very powerful in this field, although they require still a large number of computations for a detailed reproduction of electronic properties, such as band shapes. Besides the ongoing research in more traditional computational chemistry fields, data analysis and machine learning methods have been increasingly employed as complementary approaches for efficient data exploration, prediction and model development, starting from the data resulting from MD simulations and electronic structure calculations. In this work, dataset reduction capabilities by unsupervised clustering techniques applied to MD trajectories are proposed and tested for the ab initio modeling of electronic absorption spectra of two challenging case studies: a non-covalent charge-transfer dimer and a ruthenium complex in solution at room temperature. The K-medoids clustering technique is applied and is proven to be able to reduce by ∼100 times the total cost of excited state calculations on an MD sampling with no loss in the accuracy and it also provides an easier understanding of the representative structures (medoids) to be analyzed on the molecular scale.

Pharmacophore Modeling of Targets Infested with Activity Cliffs via Molecular Dynamics Simulation Coupled with QSAR and Comparison with other Pharmacophore Generation Methods: KDR as Case Study

Activity cliffs (ACs) are defined as pairs of structurally similar compounds with large difference in their potencies against certain biotarget. We recently proposed that potent AC members induce significant entropically-driven conformational modifications of the target that unveil additional binding interactions, while their weakly-potent counterparts are enthalpically-driven binders with little influence on the protein target. We herein propose to extract pharmacophores for ACs-infested target(s) from molecular dynamics (MD) frames of purely "enthalpic" potent binder(s) complexed within the particular target. Genetic function algorithm/machine learning (GFA/ML) can then be employed to search for the best possible combination of MD pharmacophore(s) capable of explaining bioactivity variations within a list of inhibitors. We compared the performance of this approach with established ligand-based and structure-based methods. Kinase inserts domain receptor (KDR) was used as a case study. KDR plays a crucial role in angiogenic signaling and its inhibitors have been approved in cancer treatment. Interestingly, GFA/ML selected, MD-based, pharmacophores were of comparable performances to ligand-based and structure-based pharmacophores. The resulting pharmacophores and QSAR models were used to capture hits from the national cancer institute list of compounds. The most active hit showed anti-KDR IC50 of 2.76 µM.

From Data to Knowledge: Systematic Review of Tools for Automatic Analysis of Molecular Dynamics Output

An increasing number of crystal structures available on one side, and the boost of computational power available for computer-aided drug design tasks on the other, have caused that the structure-based drug design tools are intensively used in the drug development pipelines. Docking and molecular dynamics simulations, key representatives of the structure-based approaches, provide detailed information about the potential interaction of a ligand with a target receptor. However, at the same time, they require a three-dimensional structure of a protein and a relatively high amount of computational resources. Nowadays, as both docking and molecular dynamics are much more extensively used, the amount of data output from these procedures is also growing. Therefore, there are also more and more approaches that facilitate the analysis and interpretation of the results of structure-based tools. In this review, we will comprehensively summarize approaches for handling molecular dynamics simulations output. It will cover both statistical and machine-learning-based tools, as well as various forms of depiction of molecular dynamics output.

How to assess the structural dynamics of transcription factors by integrating sparse NMR and EPR constraints with molecular dynamics simulations

We review recent advances in modeling structural ensembles of transcription factors from nuclear magnetic resonance (NMR) and electron paramagnetic resonance (EPR) spectroscopic data, integrated with molecular dynamics (MD) simulations. We focus on approaches that confirm computed conformational ensembles by sparse constraints obtained from magnetic resonance. This combination enables the deduction of functional and structural protein models even if nuclear Overhauser effects (NOEs) are too scarce for conventional structure determination. We highlight recent insights into the folding-upon-DNA binding transitions of intrinsically disordered transcription factors that could be assessed using such integrative approaches.

Structural modelling and thermostability of a serine protease inhibitor belonging to the Kunitz-BPTI family from the Rhipicephalus microplus tick

rBmTI-A is a recombinant serine protease inhibitor that belongs to the Kunitz-BPTI family and that was cloned from Rhipicephalus microplus tick. rBmTI-A has inhibitory activities on bovine trypsin, human plasma kallikrein, human neutrophil elastase and plasmin with dissociation constants in nM range. It is characterized by two inhibitory domains and each domain presents six cysteines that form three disulfide bonds, which contribute to the high stability of its structure. Previous studies suggest that serine protease inhibitor rBmTI-A has a protective potential against pulmonary emphysema in mice and anti-inflammatory potential. Besides that, rBmTI-A presented a potent inhibitory activity against in vitro vessel formation. In this study, the tertiary structure of rBmTI-A was modeled. The structure stabilization was evaluated by molecular dynamics analysis. Circular dichroism spectroscopy data corroborated the secondary structure found by the homology modelling. Also, in circular dichroism data it was shown a thermostability of rBmTI-A until approximately 70 °C, corroborated by inhibitory assays toward trypsin.

Data-Driven Molecular Dynamics: A Multifaceted Challenge

The big data concept is currently revolutionizing several fields of science including drug discovery and development. While opening up new perspectives for better drug design and related strategies, big data analysis strongly challenges our current ability to manage and exploit an extraordinarily large and possibly diverse amount of information. The recent renewal of machine learning (ML)-based algorithms is key in providing the proper framework for addressing this issue. In this respect, the impact on the exploitation of molecular dynamics (MD) simulations, which have recently reached mainstream status in computational drug discovery, can be remarkable. Here, we review the recent progress in the use of ML methods coupled to biomolecular simulations with potentially relevant implications for drug design. Specifically, we show how different ML-based strategies can be applied to the outcome of MD simulations for gaining knowledge and enhancing sampling. Finally, we discuss how intrinsic limitations of MD in accurately modeling biomolecular systems can be alleviated by including information coming from experimental data.

Molecular enhanced sampling with autoencoders: On-the-fly collective variable discovery and accelerated free energy landscape exploration

Macromolecular and biomolecular folding landscapes typically contain high free energy barriers that impede efficient sampling of configurational space by standard molecular dynamics simulation. Biased sampling can artificially drive the simulation along prespecified collective variables (CVs), but success depends critically on the availability of good CVs associated with the important collective dynamical motions. Nonlinear machine learning techniques can identify such CVs but typically do not furnish an explicit relationship with the atomic coordinates necessary to perform biased sampling. In this work, we employ auto-associative artificial neural networks ("autoencoders") to learn nonlinear CVs that are explicit and differentiable functions of the atomic coordinates. Our approach offers substantial speedups in exploration of configurational space, and is distinguished from existing approaches by its capacity to simultaneously discover and directly accelerate along data-driven CVs. We demonstrate the approach in simulations of alanine dipeptide and Trp-cage, and have developed an open-source and freely available implementation within OpenMM. © 2018 Wiley Periodicals, Inc.

Collective variable discovery and enhanced sampling using autoencoders: Innovations in network architecture and error function design

Auto-associative neural networks ("autoencoders") present a powerful nonlinear dimensionality reduction technique to mine data-driven collective variables from molecular simulation trajectories. This technique furnishes explicit and differentiable expressions for the nonlinear collective variables, making it ideally suited for integration with enhanced sampling techniques for accelerated exploration of configurational space. In this work, we describe a number of sophistications of the neural network architectures to improve and generalize the process of interleaved collective variable discovery and enhanced sampling. We employ circular network nodes to accommodate periodicities in the collective variables, hierarchical network architectures to rank-order the collective variables, and generalized encoder-decoder architectures to support bespoke error functions for network training to incorporate prior knowledge. We demonstrate our approach in blind collective variable discovery and enhanced sampling of the configurational free energy landscapes of alanine dipeptide and Trp-cage using an open-source plugin developed for the OpenMM molecular simulation package.

Cluster analysis of molecular simulation trajectories for systems where both conformation and orientation of the sampled states are important

Clustering methods have been widely used to group together similar conformational states from molecular simulations of biomolecules in solution. For applications such as the interaction of a protein with a surface, the orientation of the protein relative to the surface is also an important clustering parameter because of its potential effect on adsorbed-state bioactivity. This study presents cluster analysis methods that are specifically designed for systems where both molecular orientation and conformation are important, and the methods are demonstrated using test cases of adsorbed proteins for validation. Additionally, because cluster analysis can be a very subjective process, an objective procedure for identifying both the optimal number of clusters and the best clustering algorithm to be applied to analyze a given dataset is presented. The method is demonstrated for several agglomerative hierarchical clustering algorithms used in conjunction with three cluster validation techniques. © 2016 Wiley Periodicals, Inc.

An Effective Approach for Clustering InhA Molecular Dynamics Trajectory Using Substrate-Binding Cavity Features

Protein receptor conformations, obtained from molecular dynamics (MD) simulations, have become a promising treatment of its explicit flexibility in molecular docking experiments applied to drug discovery and development. However, incorporating the entire ensemble of MD conformations in docking experiments to screen large candidate compound libraries is currently an unfeasible task. Clustering algorithms have been widely used as a means to reduce such ensembles to a manageable size. Most studies investigate different algorithms using pairwise Root-Mean Square Deviation (RMSD) values for all, or part of the MD conformations. Nevertheless, the RMSD only may not be the most appropriate gauge to cluster conformations when the target receptor has a plastic active site, since they are influenced by changes that occur on other parts of the structure. Hence, we have applied two partitioning methods (k-means and k-medoids) and four agglomerative hierarchical methods (Complete linkage, Ward's, Unweighted Pair Group Method and Weighted Pair Group Method) to analyze and compare the quality of partitions between a data set composed of properties from an enzyme receptor substrate-binding cavity and two data sets created using different RMSD approaches. Ensembles of representative MD conformations were generated by selecting a medoid of each group from all partitions analyzed. We investigated the performance of our new method for evaluating binding conformation of drug candidates to the InhA enzyme, which were performed by cross-docking experiments between a 20 ns MD trajectory and 20 different ligands. Statistical analyses showed that the novel ensemble, which is represented by only 0.48% of the MD conformations, was able to reproduce 75% of all dynamic behaviors within the binding cavity for the docking experiments performed. Moreover, this new approach not only outperforms the other two RMSD-clustering solutions, but it also shows to be a promising strategy to distill biologically relevant information from MD trajectories, especially for docking purposes.

Role of indirect readout mechanism in TATA box binding protein-DNA interaction

Gene expression generally initiates from recognition of TATA-box binding protein (TBP) to the minor groove of DNA of TATA box sequence where the DNA structure is significantly different from B-DNA. We have carried out molecular dynamics simulation studies of TBP-DNA system to understand how the DNA structure alters for efficient binding. We observed rigid nature of the protein while the DNA of TATA box sequence has an inherent flexibility in terms of bending and minor groove widening. The bending analysis of the free DNA and the TBP bound DNA systems indicate presence of some similar structures. Principal coordinate ordination analysis also indicates some structural features of the protein bound and free DNA are similar. Thus we suggest that the DNA of TATA box sequence regularly oscillates between several alternate structures and the one suitable for TBP binding is induced further by the protein for proper complex formation.

Visualizing energy landscapes with metric disconnectivity graphs

The visualization of multidimensional energy landscapes is important, providing insight into the kinetics and thermodynamics of a system, as well the range of structures a system can adopt. It is, however, highly nontrivial, with the number of dimensions required for a faithful reproduction of the landscape far higher than can be represented in two or three dimensions. Metric disconnectivity graphs provide a possible solution, incorporating the landscape connectivity information present in disconnectivity graphs with structural information in the form of a metric. In this study, we present a new software package, PyConnect, which is capable of producing both disconnectivity graphs and metric disconnectivity graphs in two or three dimensions. We present as a test case the analysis of the 69-bead BLN coarse-grained model protein and show that, by choosing appropriate order parameters, metric disconnectivity graphs can resolve correlations between structural features on the energy landscape with the landscapes energetic and kinetic properties.

Expected distributions of root-mean-square positional deviations in proteins

The atom positional root-mean-square deviation (RMSD) is a standard tool for comparing the similarity of two molecular structures. It is used to characterize the quality of biomolecular simulations, to cluster conformations, and as a reaction coordinate for conformational changes. This work presents an approximate analytic form for the expected distribution of RMSD values for a protein or polymer fluctuating about a stable native structure. The mean and maximum of the expected distribution are independent of chain length for long chains and linearly proportional to the average atom positional root-mean-square fluctuations (RMSF). To approximate the RMSD distribution for random-coil or unfolded ensembles, numerical distributions of RMSD were generated for ensembles of self-avoiding and non-self-avoiding random walks. In both cases, for all reference structures tested for chains more than three monomers long, the distributions have a maximum distant from the origin with a power-law dependence on chain length. The purely entropic nature of this result implies that care must be taken when interpreting stable high-RMSD regions of the free-energy landscape as "intermediates" or well-defined stable states.

Accurate prediction of the binding free energy and analysis of the mechanism of the interaction of replication protein A (RPA) with ssDNA

The eukaryotic replication protein A (RPA) has several pivotal functions in the cell metabolism, such as chromosomal replication, prevention of hairpin formation, DNA repair and recombination, and signaling after DNA damage. Moreover, RPA seems to have a crucial role in organizing the sequential assembly of DNA processing proteins along single stranded DNA (ssDNA). The strong RPA affinity for ssDNA, K(A) between 10(-9)-10(-10) M, is characterized by a low cooperativity with minor variation for changes on the nucleotide sequence. Recently, new data on RPA interactions was reported, including the binding free energy of the complex RPA70AB with dC(8) and dC(5), which has been estimated to be -10 ± 0.4 kcal mol(-1) and -7 ± 1 kcal mol(-1), respectively. In view of these results we performed a study based on molecular dynamics aimed to reproduce the absolute binding free energy of RPA70AB with the dC(5) and dC(8) oligonucleotides. We used several tools to analyze the binding free energy, rigidity, and time evolution of the complex. The results obtained by MM-PBSA method, with the use of ligand free geometry as a reference for the receptor in the separate trajectory approach, are in excellent agreement with the experimental data, with ±4 kcal mol(-1) error. This result shows that the MM-PB(GB)SA methods can provide accurate quantitative estimates of the binding free energy for interacting complexes when appropriate geometries are used for the receptor, ligand and complex. The decomposition of the MM-GBSA energy for each residue in the receptor allowed us to correlate the change of the affinity of the mutated protein with the ΔG(gas+sol) contribution of the residue considered in the mutation. The agreement with experiment is optimal and a strong change in the binding free energy can be considered as the dominant factor in the loss for the binding affinity resulting from mutation.

Binding selectivity of RecA to a single stranded DNA, a computational approach

Homologous recombination (HR) is the major DNA double strand break repair pathway which maintains the genomic integrity. It is fundamental for the survivability and functionality of all organisms. One of the initial steps in HR is the formation of the nucleoprotein filament composed by a single stranded DNA chain surrounded by the recombinases protein. The filament orchestrates the search for an undamaged homologue, as a template for the repair process. Our theoretical study was aimed at elucidating the selectivity of the interaction between a monomer of the recombinases enzyme in the Escherichia coli, EcRecA, the bacterial homologue of human Rad51, with a series of oligonucleotides of nine bases length. The complex, equilibrated for 20 ns with Langevian dynamics, was inserted in a periodic box with a 8 Å buffer of water molecules explicitly described by the TIP3P model. The absolute binding free energies are calculated in an implicit solvent using the Poisson-Boltzmann (PB) and the generalized Born (GB) solvent accessible surface area, using the MM-PB(GB)SA model. The solute entropic contribution is also calculated by normal mode analysis. The results underline how a significant contribution of the binding free energy is due to the interaction with the Arg196, a critical amino acid for the activity of the enzyme. The study revealed how the binding affinity of EcRecA is significantly higher toward dT₉ rather than dA₉, as expected from the experimental results.

Determination, prediction, and understanding of structures, using the energy landscapes of chemical systems – Part II

In the past decade, new theoretical approaches have been developed to determine, predict and understand the struc-ture of chemical compounds. The central element of these methods has been the investigation of the energy landscape of chemical systems. Applications range from extended crystalline and amorphous compounds over clusters and molecular crystals to proteins. In this review, we are going to give an introduction to energy landscapes and methods for their investigation, together with a number of examples. These include structure prediction of extended and mo-lecular crystals, structure prediction and folding of proteins, structure analysis of zeolites, and structure determination of crystals from powder diffraction data.

Characterizing rare-event property distributions via replicate molecular dynamics simulations of proteins

As computational resources increase, molecular dynamics simulations of biomolecules are becoming an increasingly informative complement to experimental studies. In particular, it has now become feasible to use multiple initial molecular configurations to generate an ensemble of replicate production-run simulations that allows for more complete characterization of rare events such as ligand-receptor unbinding. However, there are currently no explicit guidelines for selecting an ensemble of initial configurations for replicate simulations. Here, we use clustering analysis and steered molecular dynamics simulations to demonstrate that the configurational changes accessible in molecular dynamics simulations of biomolecules do not necessarily correlate with observed rare-event properties. This informs selection of a representative set of initial configurations. We also employ statistical analysis to identify the minimum number of replicate simulations required to sufficiently sample a given biomolecular property distribution. Together, these results suggest a general procedure for generating an ensemble of replicate simulations that will maximize accurate characterization of rare-event property distributions in biomolecules.

Modeling and simulation of chemomechanics at the cell-matrix interface

Chemomechanical characteristics of the extracellular materials with which cells interact can have a profound impact on cell adhesion and migration. To understand and modulate such complex multiscale processes, a detailed understanding of the feedback between a cell and the adjacent microenvironment is crucial. Here, we use computational modeling and simulation to examine the cell-matrix interaction at both the molecular and continuum lengthscales. Using steered molecular dynamics, we consider how extracellular matrix (ECM) stiffness and extracellular pH influence the interaction between cell surface adhesion receptors and extracellular matrix ligands, and we predict potential consequences for focal adhesion formation and dissolution. Using continuum level finite element simulations and analytical methods to model cell-induced ECM deformation as a function of ECM stiffness and thickness, we consider the implications toward design of synthetic substrata for cell biology experiments that intend to decouple chemical and mechanical cues.

Fast and reliable analysis of molecular motion using proximity relations and dimensionality reduction

The analysis of molecular motion starting from extensive sampling of molecular configurations remains an important and challenging task in computational biology. Existing methods require a significant amount of time to extract the most relevant motion information from such data sets. In this work, we provide a practical tool for molecular motion analysis. The proposed method builds upon the recent ScIMAP (Scalable Isomap) method, which, by using proximity relations and dimensionality reduction, has been shown to reliably extract from simulation data a few parameters that capture the main, linear and/or nonlinear, modes of motion of a molecular system. The results we present in the context of protein folding reveal that the proposed method characterizes the folding process essentially as well as ScIMAP. At the same time, by projecting the simulation data and computing proximity relations in a low-dimensional Euclidean space, it renders such analysis computationally practical. In many instances, the proposed method reduces the computational cost from several CPU months to just a few CPU hours, making it possible to analyze extensive simulation data in a matter of a few hours using only a single processor. These results establish the proposed method as a reliable and practical tool for analyzing motions of considerably large molecular systems and proteins with complex folding mechanisms.

Bayesian model based clustering analysis: application to a molecular dynamics trajectory of the HIV-1 integrase catalytic core

This work describes the application of a Bayesian method for clustering protein conformations sampled during a molecular dynamics simulation of the HIV-1 integrase catalytic core. A clustering analysis is carried out under the assumption of normal distribution without fixing the number of clusters in advance. Some performance measures, such as posterior probability and class cross entropy, are used to determine the most probable set of clusters. The Bayesian clustering method results in meaningful groups identifying transitions between conformational ensembles. The dihedral angles involved in such transitions are also examined in detail. The conformations in high dimensional space are projected into 3D space employing a multidimensional scaling technique to provide a visual inspection.

Quantifying Polypeptide Conformational Space: Sensitivity to Conformation and Ensemble Definition

Quantifying the density of conformations over phase space (the conformational distribution) is needed to model important macromolecular processes such as protein folding. In this work, we quantify the conformational distribution for a simple polypeptide (N-mer polyalanine) using the cumulative distribution function (CDF), which gives the probability that two randomly selected conformations are separated by less than a "conformational" distance and whose inverse gives conformation counts as a function of conformational radius. An important finding is that the conformation counts obtained by the CDF inverse depend critically on the assignment of a conformation's distance span and the ensemble (e.g., unfolded state model): varying ensemble and conformation definition (1 --> 2 A) varies the CDF-based conformation counts for Ala(50) from 10(11) to 10(69). In particular, relatively short molecular dynamics (MD) relaxation of Ala(50)'s random-walk ensemble reduces the number of conformers from 10(55) to 10(14) (using a 1 A root-mean-square-deviation radius conformation definition) pointing to potential disconnections in comparing the results from simplified models of unfolded proteins with those from all-atom MD simulations. Explicit waters are found to roughen the landscape considerably. Under some common conformation definitions, the results herein provide (i) an upper limit to the number of accessible conformations that compose unfolded states of proteins, (ii) the optimal clustering radius/conformation radius for counting conformations for a given energy and solvent model, (iii) a means of comparing various studies, and (iv) an assessment of the applicability of random search in protein folding.

A line integral reaction path approximation for large systems via nonlinear constrained optimization: Application to alanine dipeptide and the β hairpin of protein G

A variation of the line integral method of Elber with self-avoiding walk has been implemented using a state of the art nonlinear constrained optimization procedure. The new implementation appears to be robust in finding approximate reaction paths for small and large systems. Exact transition states and intermediates for the resulting paths can easily be pinpointed with subsequent application of the conjugate peak refinement method [S. Fischer and M. Karplus, Chem. Phys. Lett. 194, 252 (1992)] and unconstrained minimization, respectively. Unlike previous implementations utilizing a penalty function approach, the present implementation generates an exact solution of the underlying problem. Most importantly, this formulation does not require an initial guess for the path, which makes it particularly useful for studying complex molecular rearrangements. The method has been applied to conformational rearrangements of the alanine dipeptide in the gas phase and in water, and folding of the beta hairpin of protein G in water. In the latter case a procedure was developed to systematically sample the potential energy surface underlying folding and reconstruct folding pathways within the nearest-neighbor hopping approximation.

Effect of a Bound Non-Nucleoside RT Inhibitor on the Dynamics of Wild-Type and Mutant HIV-1 Reverse Transcriptase

HIV-1 reverse transcriptase (RT) is an important target for drugs used in the treatment of AIDS. Drugs known as non-nucleoside RT inhibitors (NNRTI) appear to alter the structural and dynamical properties of RT which in turn inhibit RT's ability to transcribe. Molecular dynamics (MD), principal component analysis (PCA), and binding free energy simulations are employed to explore the dynamics of RT and its interaction with the bound NNRTI nevirapine, for both wild-type and mutant (V106A, Y181C, Y188C) RT. These three mutations commonly arise in the presence of nevirapine and result in resistance to the drug. We show that a bound NNRTI hinders the motion of almost all RT amino acids. The mutations, located in the non-nucleoside RT inhibitor binding pocket, partially restore RT flexibility. The binding affinities calculated by molecular mechanics/Poisson-Boltzmann surface accessibility (MM-PBSA) show that nevirapine interacts stronger with wild-type RT than with mutant RT. The mutations cause a loss of van der Waals interactions between the drug and the binding pocket. The results from this study suggest that a good inhibitor should efficiently enter and maximally occupy the binding pocket, thereby interacting effectively with the amino acids around the binding pocket.

Dissociation of an antiviral compound from the internal pocket of human rhinovirus 14 capsid

WIN antiviral compounds bind human rhinovirus, as well as enterovirus and parechovirus, in an internal cavity located within the viral protein capsid. Access to the buried pocket necessitates deviation from the average viral protein structure identified by crystallography. We investigated the dissociation of WIN 52084 from the pocket in human rhinovirus 14 by using an adiabatic, biased molecular dynamics simulation method. Multiple dissociation trajectories are used to characterize the pathway. WIN 52084 exits between the polypeptide chain near the ends of betaC and betaH in a series of steps. Small, transient packing defects in the protein are sufficient for dissociation. A number of torsion-angle transitions of the antiviral compound are involved, which suggests that flexibility in antiviral compounds is important for binding. It is interesting to note that dissociation is associated with an increase in the conformational fluctuations of residues never in direct contact with WIN 52084 over the course of dissociation. These residues are N-terminal residues in the viral proteins VP3 and VP4 and are located in the interior of the capsid near the icosahedral 5-fold axis. The observed changes in dynamics may be relevant to structural changes associated with virion uncoating and its inhibition by antiviral compounds.

Functionally relevant protein motions: Extracting basin-specific collective coordinates from molecular dynamics trajectories

Functionally relevant motion of proteins has been associated with a number of atoms moving in a concerted fashion along so-called "collective coordinates." We present an approach to extract collective coordinates from conformations obtained from molecular dynamics simulations. The power of this technique for differentiating local structural fluctuations between classes of conformers obtained by clustering is illustrated by analyzing nanosecond-long trajectories for the response regulator protein Spo0F of Bacillus subtilis, generated both in vacuo and using an implicit-solvent representation. Conformational clustering is performed using automated histogram filtering of the inter-C(alpha) distances. Orthogonal (varimax) rotation of the vectors obtained by principal component analysis of these interresidue distances for the members of individual clusters is key to the interpretation of collective coordinates dominating each conformational class. The rotated loadings plots isolate significant variation in interresidue distances, and these are associated with entire mobile secondary structure elements. From this we infer concerted motions of these structural elements. For the Spo0F simulations employing an implicit-solvent representation, collective coordinates obtained in this fashion are consistent with the location of the protein's known active sites and experimentally determined mobile regions.

Distributions in protein conformation space: implications for structure prediction and entropy

By considering how polymer structures are distributed in conformation space, we show that it is possible to quantify the difficulty of structural prediction and to provide a measure of progress for prediction calculations. The critical issue is the probability that a conformation is found within a specified distance of another conformer. We address this question by constructing a cumulative distribution function (CDF) for the average probability from observations about its limiting behavior at small displacements and numerical simulations of polyalanine chains. We can use the CDF to estimate the likelihood that a structure prediction is better than random chance. For example, the chance of randomly predicting the native backbone structure of a 150-amino-acid protein to low resolution, say within 6 A, is 10(-14). A high-resolution structural prediction, say to 2 A, is immensely more difficult (10(-57)). With additional assumptions, the CDF yields the conformational entropy of protein folding from native-state coordinate variance. Or, using values of the conformational entropy change on folding, we can estimate the native state's conformational span. For example, for a 150-mer protein, equilibrium alpha-carbon displacements in the native ensemble would be 0.3-0.5 A based on T Delta S of 1.42 kcal/(mol residue).

Enhanced sampling by multiple molecular dynamics trajectories: carbonmonoxy myoglobin 10 μs A0 → A1–3 transition from ten 400 picosecond simulations

The utility of multiple trajectories to extend the time scale of molecular dynamics simulations is reported for the spectroscopic A-states of carbonmonoxy myoglobin (MbCO). Experimentally, the A0-->A(1-3) transition has been observed to be 10 micros at 300 K, which is beyond the time scale of standard molecular dynamics simulations. To simulate this transition, 10 short (400 ps) and two longer time (1.2 ns) molecular dynamics trajectories, starting from five different crystallographic and solution phase structures with random initial velocities centered in a 37 A radius sphere of water, have been used to sample the native-fold of MbCO. Analysis of the ensemble of structures gathered over the cumulative 5.6 ns reveals two biomolecular motions involving the side chains of His64 and Arg45 to explain the spectroscopic states of MbCO. The 10 micros A0-->A(1-3) transition involves the motion of His64, where distance between His64 and CO is found to vary up to 8.8 +/- 1.0 A during the transition of His64 from the ligand (A(1-3)) to bulk solvent (A0). The His64 motion occurs within a single trajectory only once, however the multiple trajectories populate the spectroscopic A-states fully. Consequently, multiple independent molecular dynamics simulations have been found to extend biomolecular motion from 5 ns of total simulation to experimental phenomena on the microsecond time scale.

Information content of molecular structures

For a completely enumerated set of conformers of a macromolecule or for exhaustive lattice walks of model polymers it is straightforward to use Shannon information theory to deduce the information content of the ensemble. It is also practicable to develop numerical measures of the information content of sets of exact distance constraints applied to specific conformational ensembles. We examine the effects of experimental uncertainties by considering "noisy" constraints. The introduction of noise requires additional assumptions about noise distribution and conformational clustering protocols that make the problem of measuring information content more complex. We make use of a standard concept in communication theory, the "noise sphere," to link uncertainty in measurements to information loss. Most of our numerical results are derived from two-dimensional lattice ensembles. Expressing results in terms of information per degree of freedom removes almost all of the chain length dependence. We also explore off-lattice polyalanine chains that yield surprisingly similar results.

A method for correlations analysis of coordinates: applications for molecular conformations

We describe a new method to analyze multiple correlations between subsets of coordinates that represent a sample. The correlation is established only between specific regions of interest at the coordinates. First, the region(s) of interest are selected at each molecular coordinate. Next, a correlation matrix is constructed for the selected regions. The matrix is subject to further analysis, illuminating the multidimensional structural characteristics that exist in the conformational space. The method's abilities are demonstrated in several examples: it is used to analyze the conformational space of complex molecules, it is successfully applied to compare related conformational spaces, and it is used to analyze a diverse set of protein folding trajectories.

Hamming distance geometry of a protein conformational space: application to the clustering of a 4-ns molecular dynamics trajectory of the HIV-1 integrase catalytic core

Protein structures can be encoded into binary sequences (Gabarro-Arpa et al., Comput Chem 2000;24:693-698) these are used to define a Hamming distance in conformational space: the distance between two different molecular conformations is the number of different bits in their sequences. Each bit in the sequence arises from a partition of conformational space in two halves. Thus, the information encoded in the binary sequences is also used to characterize the regions of conformational space visited by the system. We apply this distance and their associated geometric structures to the clustering and analysis of conformations sampled during a 4-ns molecular dynamics simulation of the HIV-1 integrase catalytic core. The cluster analysis of the simulation shows a division of the trajectory into two segments of 2.6 and 1.4 ns length, which are qualitatively different: the data points to the fact that equilibration is only reached at the end of the first segment. The Hamming distance is compared also to the r.m.s. deviation measure. The analysis of the cases studied so far shows that under the same conditions the two measures behave quite differently, and that the Hamming distance appears to be more robust than the r.m.s. deviation.

Top-down free-energy minimization on protein potential energy landscapes

The hierarchical properties of potential energy landscapes have been used to gain insight into thermodynamic and kinetic properties of protein ensembles. It also may be possible to use them to direct computational searches for thermodynamically stable macroscopic states, i.e., computational protein folding. To this end, we have developed a top-down search procedure in which conformation space is recursively dissected according to the intrinsic hierarchical structure of a landscape's effective-energy barriers. This procedure generates an inverted tree similar to the disconnectivity graphs generated by local minima-clustering methods, but it fundamentally differs in the manner in which the portion of the tree that is to be computationally explored is selected. A key ingredient is a branch-selection algorithm that takes advantage of statistically predictive properties of the landscape to guide searches down the tree branches that are most likely to lead to the physically relevant macroscopic states. Using the computational folding of a beta-hairpin-forming peptide as an example, we show that such predictive properties indeed exist and can be used for structure prediction by free-energy global minimization.

Application of self-organizing maps in conformational analysis of lipids

The characteristics of lipid assemblies are important for the functions of biological membranes. This has led to an increasing utilization of molecular dynamics simulations for the elucidation of the structural features of biomembranes. We have applied the self-organizing map (SOM) to the analysis of the complex conformational data from a 1-ns molecular dynamics simulation of PLPC phospholipids in a membrane assembly. Mapping of 1.44 million molecular conformations to a two-dimensional array of neurons revealed, without human intervention, the main conformational features in hours. Both the whole molecule and the characteristics of the unsaturated fatty acid chains were analyzed. All major structural features were easily distinguished, such as the orientational variability of the headgroup, the mainly trans state dihedral angles of the sn-1 chain, and both straight and bent conformations of the unsaturated sn-2 chain. Furthermore, presentation of the trajectory of an individual lipid molecule on the map provides information on conformational dynamics. The present results suggest that the SOM method provides a powerful tool for routinely gaining rapid insight to the main molecular conformations as well as to the conformational dynamics of any simulated molecular assembly without the requirement of a priori knowledge.