Structure prediction of cyclic peptides by molecular dynamics + machine learning

Recent computational methods have made strides in discovering well-structured cyclic peptides that preferentially populate a single conformation. However, many successful cyclic-peptide therapeutics adopt multiple conformations in solution. In fact, the chameleonic properties of some cyclic peptides are likely responsible for their high cell membrane permeability. Thus, we require the ability to predict complete structural ensembles for cyclic peptides, including the majority of cyclic peptides that have broad structural ensembles, to significantly improve our ability to rationally design cyclic-peptide therapeutics. Here, we introduce the idea of using molecular dynamics simulation results to train machine learning models to enable efficient structure prediction for cyclic peptides. Using molecular dynamics simulation results for several hundred cyclic pentapeptides as the training datasets, we developed machine-learning models that can provide molecular dynamics simulation-quality predictions of structural ensembles for all the hundreds of thousands of sequences in the entire sequence space. The prediction for each individual cyclic peptide can be made using less than 1 second of computation time. Even for the most challenging classes of poorly structured cyclic peptides with broad conformational ensembles, our predictions were similar to those one would normally obtain only after running multiple days of explicit-solvent molecular dynamics simulations. The resulting method, termed StrEAMM (Structural Ensembles Achieved by Molecular Dynamics and Machine Learning), is the first technique capable of efficiently predicting complete structural ensembles of cyclic peptides without relying on additional molecular dynamics simulations, constituting a seven-order-of-magnitude improvement in speed while retaining the same accuracy as explicit-solvent simulations.

The StrEAMM method enables predicting the structural ensembles of cyclic peptides that adopt multiple conformations in solution.     

Molecular activity prediction by means of supervised subspace projection based ensembles of classifiers

This paper proposes a method for molecular activity prediction in QSAR studies using ensembles of classifiers constructed by means of two supervised subspace projection methods, namely nonparametric discriminant analysis (NDA) and hybrid discriminant analysis (HDA). We studied the performance of the proposed ensembles compared to classical ensemble methods using four molecular datasets and eight different models for the representation of the molecular structure. Using several measures and statistical tests for classifier comparison, we observe that our proposal improves the classification results with respect to classical ensemble methods. Therefore, we show that ensembles constructed using supervised subspace projections offer an effective way of creating classifiers in cheminformatics.     

RNACluster: An integrated tool for RNA secondary structure comparison and clustering

RNA structure comparison is a fundamental problem in structural biology, structural chemistry, and bioinformatics. It can be used for analysis of RNA energy landscapes, conformational switches, and facilitating RNA structure prediction. The purpose of our integrated tool RNACluster is twofold: to provide a platform for computing and comparison of different distances between RNA secondary structures, and to perform cluster identification to derive useful information of RNA structure ensembles, using a minimum spanning tree (MST) based clustering algorithm. RNACluster employs a cluster identification approach based on a MST representation of the RNA ensemble data and currently supports six distance measures between RNA secondary structures. RNACluster provides a user-friendly graphical interface to allow a user to compare different structural distances, analyze the structure ensembles, and visualize predicted structural clusters.     

Computation of Free Energy

Many quantities that are standardly used to characterize a chemical system are related to free‐energy differences between particular states of the system. By statistical mechanics, free‐energy differences may be expressed in terms of averages over ensembles of atomic configurations for the molecular system of interest. Here, we review the most useful formulae to calculate free‐energy differences from ensembles generated by molecular simulation, illustrate a number of recent developments, and highlight practical aspects of such calculations with examples selected from the literature.     

Simulation studies of the fidelity of biomolecular structure ensemble recreation

We examine the ability of Bayesian methods to recreate structural ensembles for partially folded molecules from averaged data. Specifically we test the ability of various algorithms to recreate different transition state ensembles for folding proteins using a multiple replica simulation algorithm using input from "gold standard" reference ensembles that were first generated with a Go-like Hamiltonian having nonpairwise additive terms. A set of low resolution data, which function as the "experimental" phi values, were first constructed from this reference ensemble. The resulting phi values were then treated as one would treat laboratory experimental data and were used as input in the replica reconstruction algorithm. The resulting ensembles of structures obtained by the replica algorithm were compared to the gold standard reference ensemble, from which those "data" were, in fact, obtained. It is found that for a unimodal transition state ensemble with a low barrier, the multiple replica algorithm does recreate the reference ensemble fairly successfully when no experimental error is assumed. The Kolmogorov-Smirnov test as well as principal component analysis show that the overlap of the recovered and reference ensembles is significantly enhanced when multiple replicas are used. Reduction of the multiple replica ensembles by clustering successfully yields subensembles with close similarity to the reference ensembles. On the other hand, for a high barrier transition state with two distinct transition state ensembles, the single replica algorithm only samples a few structures of one of the reference ensemble basins. This is due to the fact that the phi values are intrinsically ensemble averaged quantities. The replica algorithm with multiple copies does sample both reference ensemble basins. In contrast to the single replica case, the multiple replicas are constrained to reproduce the average phi values, but allow fluctuations in phi for each individual copy. These fluctuations facilitate a more faithful sampling of the reference ensemble basins. Finally, we test how robustly the reconstruction algorithm can function by introducing errors in phi comparable in magnitude to those suggested by some authors. In this circumstance we observe that the chances of ensemble recovery with the replica algorithm are poor using a single replica, but are improved when multiple copies are used. A multimodal transition state ensemble, however, turns out to be more sensitive to large errors in phi (if appropriately gauged) and attempts at successful recreation of the reference ensemble with simple replica algorithms can fall short.     

The Effect of Pulling and Twisting Forces on Chameleon Sequence Peptides**

Chameleon sequences are amino acid sequences found in several distinct configurations in experiment. They challenge our understanding of the link between sequence and structure, and provide insight into structural competition in proteins. Here, we study the energy landscapes for three such sequences, and interrogate how pulling and twisting forces impact the available structural ensembles. Chameleon sequences do not necessarily exhibit multiple structural ensembles on a multifunnel energy landscape when we consider them in isolation. The application of even small forces leads to drastic changes in the energy landscapes. For pulling forces, we observe transitions from helical to extended structures in a very small span of forces. For twisting forces, the picture is much more complex, and highly dependent on the magnitude and handedness of the applied force as well as the reference angle for the twist. Depending on these parameters, more complex and more simplistic energy landscapes are observed alongside more and less diverse structural ensembles. The impact of even small forces is significant, confirming their likely role in folding events. In addition, small forces exerted by the remaining scaffold of a protein may be sufficient to lead to the adoption of a specific structural ensemble by a chameleon sequence.     

Protonation/deprotonation effects on the stability of the Trp-cage miniprotein

The Trp-cage miniprotein is a 20 amino acid peptide that exhibits many of the properties of globular proteins. In this protein, the hydrophobic core is formed by a buried Trp side chain. The folded state is stabilized by an ion pair between aspartic acid and an arginine side chain. The effect of protonating the aspartic acid on the Trp-cage miniprotein folding/unfolding equilibrium is studied by explicit solvent molecular dynamics simulations of the protein in the charged and protonated Asp9 states. Unbiased Replica Exchange Molecular Dynamics (REMD) simulations, spanning a wide temperature range, are carried out to the microsecond time scale, using the AMBER99SB forcefield in explicit TIP3P water. The protein structural ensembles are studied in terms of various order parameters that differentiate the folded and unfolded states. We observe that in the folded state the root mean square distance (rmsd) from the backbone of the NMR structure shows two highly populated basins close to the native state with peaks at 0.06 nm and 0.16 nm, which are consistent with previous simulations using the same forcefield. The fraction of folded replicas shows a drastic decrease because of the absence of the salt bridge. However, significant populations of conformations with the arginine side chain exposed to the solvent, but within the folded basin, are found. This shows the possibility to reach the folded state without formation of the ion pair. We also characterize changes in the unfolded state. The equilibrium populations of the folded and unfolded states are used to characterize the thermodynamics of the system. We find that the change in free energy difference due to the protonation of the Asp amino acid is 3 kJ mol(-1) at 297 K, favoring the charged state, and resulting in ΔpK(1) = 0.5 units for Asp9. We also study the differences in the unfolded state ensembles for the two charge states and find significant changes at low temperature, where the protonated Asp side chain makes multiple hydrogen bonds to the protein backbone.     

Comparison of chemical and thermal protein denaturation by combination of computational and experimental approaches. II

Chemical and thermal denaturation methods have been widely used to investigate folding processes of proteins in vitro. However, a molecular understanding of the relationship between these two perturbation methods is lacking. Here, we combined computational and experimental approaches to investigate denaturing effects on three structurally different proteins. We derived a linear relationship between thermal denaturation at temperature T(b) and chemical denaturation at another temperature T(u) using the stability change of a protein (ΔG). For this, we related the dependence of ΔG on temperature, in the Gibbs-Helmholtz equation, to that of ΔG on urea concentration in the linear extrapolation method, assuming that there is a temperature pair from the urea (T(u)) and the aqueous (T(b)) ensembles that produces the same protein structures. We tested this relationship on apoazurin, cytochrome c, and apoflavodoxin using coarse-grained molecular simulations. We found a linear correlation between the temperature for a particular structural ensemble in the absence of urea, T(b), and the temperature of the same structural ensemble at a specific urea concentration, T(u). The in silico results agreed with in vitro far-UV circular dichroism data on apoazurin and cytochrome c. We conclude that chemical and thermal unfolding processes correlate in terms of thermodynamics and structural ensembles at most conditions; however, deviations were found at high concentrations of denaturant.     

General method to simulate molecular processes involving complex interactions between combining subsystems

We consider the transformation process of one molecular subsystem into another (for example, structural isomer-isomer transformation) under the condition of a group of close levels in the first subsystem that have the energy on average coincident with the mean energy of the second subsystem (quasi-degeneration). It is shown that, similar to the previously discussed resonance between two levels of two subsystems, it is also possible in this case to compose an oscillating wave packet leading to a resonant transition from one subsystem to another. The calculation procedure is described that can be applied to atomic ensembles of any complexity with any number of quasi-resonant levels.     

Spontaneous transmission of chirality through multiple length scales

The hierarchical transfer of chirality in nature, from the nano‐, to meso‐, to macroscopic length scales, is very complex, and as of yet, not well understood. The advent of scanning probes has allowed chirality to be monitored at the single molecule or monolayer level and has opened up the possibility to track enantiospecific interactions and chiral self‐assembly with molecular‐scale detail. This paper describes the self‐assembly of a simple, model molecule (naphtho[2,3‐a]pyrene) that is achiral in the gas phase, but becomes chiral when adsorbed on a surface. This polyaromatic hydrocarbon forms a stable and reversibly ordered system on Cu(111) in which the transmission of chirality from single surface‐bound molecules to complex 2D chiral architectures can be monitored as a function of molecular packing density and surface temperature. In addition to the point chirality of the surface‐bound molecule, the unit cell of the molecular domains was also found to be chiral due to the incommensurate alignment of the molecular rows with respect to the underlying metal lattice. These molecular domains always aggregated in groups of three, all of the same chirality, but with different rotational orientations, forming homochiral “tri‐lobe” ensembles. At a larger length scale, these tri‐lobe ensembles associated with nearest‐neighbor tri‐lobe units of opposite chirality at lower packing densities before forming an extended array of homochiral tri‐lobe ensembles at higher converges. This system displayed chirality at a variety of size scales from the molecular (≈1 nm) and domain (≈5 nm) to the tri‐lobe ensemble (≈10 nm) and extended array (>25 nm) levels. The chirality of the tri‐lobe ensembles dictated how the overall surface packing occurred and both homo‐ and heterochiral arrays could be reproducibly and reversibly formed and interchanged as a function of surface coverage. Finally, these chirally templated surfaces displayed remarkable enantiospecificity for naphtho[2,3‐a]pyrene molecules adsorbed in the second layer. Given its simplicity, reversibility, and rich degree of order, this system represents an ideal test bed for the investigation of symmetry breaking and the hierarchical transmission of chirality through multiple length scales.     

Electronic coupling between base pair dimers of LNA:DNA oligomers

We calculated ab initio electronic coupling elements between neighboring base-pair dimers in a set of LNA:DNA oligomers with different numbers of locked nucleotides and compared them by averaging the values over ensembles of snapshots from molecular dynamics trajectories. Averaging was based on coupling elements for various ensembles comprising of 33,000 structures. The known pronounced variations of coupling elements on the nanosecond timescale due to thermal fluctuations of the DNA structure were confirmed. We found significant differences in electronic coupling at the dimer level between a non-modified DNA:DNA duplex and the corresponding duplex containing one fully LNA-substituted strand. We rationalized these differences by very dissimilar overlap in the pi-stack as a consequence of the LNA-modified system approximating an A-DNA-type helix. The calculated coupling elements for the non-modified reference duplex were similar to those of standard B-DNA and those for the fully modified oligomer resembled the matrix elements estimated for standard A-DNA.     

MC(JBW): Simple but smart Monte Carlo algorithm for free energy simulations of multiconformational molecules

Many of the most common molecular simulation methods, including Monte Carlo (MC) and molecular or stochastic dynamics (MD or SD), have significant difficulties in sampling the space of molecular potential energy surfaces characterized by multiple conformational minima and significant energy barriers. In such cases improved sampling can be obtained by special techniques that lower such barriers or somehow direct search steps toward different low energy regions of space. We recently described a hybrid MC/SD algorithm [MC(JBW)/SD] incorporating such a technique that directed MC moves of selected torsion and bond angles toward known low energy regions of conformational space. Exploration of other degrees of freedom was left to the SD part of the hybrid algorithm. In the work described here, we develop a related but simpler simulation algorithm that uses only MC to sample all degrees of freedom (e.g., stretch, bend, and torsion). We term this algorithm MC(JBW). Using simulations on various model potential energy surfaces and on simple molecular systems (n-pentane, n-butane, and cyclohexane), MC(JBW) is shown to generate ensembles of states that are indistinguishable from the canonical ensembles generated by classical Metropolis MC in the limit of very long simulations. We further demonstrate the utility of MC(JBW) by evaluating the room temperature free energy differences between conformers of various substituted cyclohexanes and the larger ring hydrocarbons cycloheptane, cyclooctane, cyclononane, and cyclodecane. The results compare favorably with available experimental data and results from previously reported MC(JBW)/SD conformational free energy calculations. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1736–1745, 1998     

Recognition of RNA by amide modified backbone nucleic acids: molecular dynamics simulations of DNA-RNA hybrids in aqueous solution

Thermodynamic and structural properties of a chemically modified DNA-RNA hybrid in which a phosphodiester linkage is replaced by a neutral amide-3 linkage (3'-CH(2)-CONH-5') were investigated using UV melting experiments, molecular dynamics simulations in explicit water, and continuum solvent models. van't Hoff analysis of the experimental UV melting curves suggests that the significant increase of the thermodynamic stability of a 15-mer DNA-RNA with seven alternated amide-3 modifications (+11 degrees C) is mainly due to an increased binding enthalpy. To further evaluate the origin in the observed affinities differences, the electrostatic contribution to the binding free energy was calculated by solving the Poisson-Boltzmann equation numerically. The nonelectrostatic contribution was estimated as the product of a hydrophobic surface tension coefficient and the surface area that is buried upon double strand formation. Structures were taken from 10 ns molecular dynamics simulations computed in a consistent fashion using explicit solvent, counterions, and the particle-mesh Ewald procedure. The present preliminary thermodynamic study suggests that the favorable binding free energy of the amide-3 DNA single strand to the complementary RNA is equally driven by electrostatic and nonpolar contributions to the binding compared to their natural analogues. In addition, molecular dynamics simulations in explicit water were performed on an amide-3 DNA single strand and the corresponding natural DNA. Results from the conformations cluster analysis of the simulated amide-3 DNA single strand ensembles suggest that the 25% of the population sampled within 10 ns has a pre-organized conformation where the sugar C3' endo pucker is favored at the 3'-flanking nucleotides. These structural and thermodynamic features contribute to the understanding of the observed increased affinities of the amide-3 DNA-RNA hybrids at the microscopic level.     

Protein Structural Ensembles Visualized by Solvent Paramagnetic Relaxation Enhancement

A protein can be in different conformations when fulfilling its function. Yet depiction of protein structural ensembles remains difficult. Here we show that the accurate measurement of solvent paramagnetic relaxation enhancement (sPRE) in the presence of an inert paramagnetic cosolute allows the assessment of protein dynamics. Demonstrated with two multi‐domain proteins, we present a method to characterize protein microsecond–millisecond dynamics based on the analysis of the sPRE. Provided with the known structures of a protein, our method uncovers an ensemble of structures that fully accounts for the observed sPRE. In conjunction with molecular dynamics simulations, our method can identify protein alternative conformation that has only been theorized before. Together, our method expands the application of sPRE beyond structural characterization of rigid proteins and complements the established PRE NMR technique.     

Multiple pathways in conformational transitions of the alanine dipeptide: an application of dynamic importance sampling

We present multiple dynamic transition pathways on the two-dimensional dihedral plane between conformational states of the alanine dipeptide. The method used in this study is dynamic importance sampling (DIMS). To perform DIMS, unbiased molecular dynamic simulations are used to generate equilibrium ensembles for the alanine dipeptide within different states. Free energy surfaces on the dihedral plane are calculated from the equilibrium simulations, and four energy minima defined from the surface are used as the starting and ending points for DIMS dynamics. The DIMS method represents an important step towards finding multiple transition pathways within complex biomolecular systems.     

Ensembles of the gallium,indium, and scandium complexes with the macrocycles

Specific features of coordination chemistry of three-charge cations were considered for the formation of ensembles of crown ether 1,4,7,10,13,16-hexaoxocyclooctadecane (18C6) or its complexes M[18C6]⁺ (M = Na, K) with the gallium, indium, and scandium chloride and thiocyanate complexes. These features were geometric isomerism, disproportionation of molecular forms to ionic ones, and the influence of solvate molecules on the structural features of the supramolecular ensembles. The conditions for the recognition and stabilization of the complexes and the formation of heterometallic polymers were revealed.     

Hybrid Hamiltonian replica exchange molecular dynamics simulation method employing the Poisson-Boltzmann model

A hybrid Hamiltonian replica exchange molecular dynamics simulation scheme based on explicit water model hybrided with Poisson-Boltzmann model is brought out. In this method the motions of atoms are governed by potential energy obtained from explicit water model. However, the exchanges between different replicas under different temperatures are controlled by the solvation energies of the solute calculated using the Poisson-Boltzmann model. In order to get the correct canonical ensembles, the van der Waals radii, which are used to define the dielectric boundary, have to be optimized. The conformational spaces of three distinct pentapeptides, Met-enkephalin, alanine 5, and glycine 5, are explored. We find that with the optimized radii the structural ensembles are nearly identical to those obtained by standard replica exchange simulations while the number of replica needed is reduced greatly.     

Some aspects of dynamic chaos in small Lennard-Jones clusters

For the purpose of studying stochastic phenomena in the clusters due to an instability of their trajectories in phase space, an evolution of ensembles of the clusters with the representative points initially located in small regions of phase space is investigated using the molecular dynamics method. A path of a cloud of the representative points through the subspace of phase space available for the clusters, and cloud growth in size and in volume are considered. Specific results are presented for Ar₄. A “cold” cluster (subjected to neither an isomerization nor an evaporation) is examined in details. A brief comparison with a “hot” cluster (both the phenomena are available) is given.