A Local Rigid Body Framework for Global Optimization of Biomolecules

Energy landscapes for clusters of hexapeptides

We present the results for energy landscapes of hexapeptides obtained using interfaces to the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) program. We have used basin-hopping global optimization and discrete path sampling to explore the landscapes of hexapeptide monomers, dimers, and oligomers containing 10, 100, and 200 monomers modeled using a residue-level coarse-grained potential, Mpipi, implemented in LAMMPS. We find that the dimers of peptides containing amino acid residues that are better at promoting phase separation, such as tyrosine and arginine, have melting peaks at higher temperature in their heat capacity compared to phenylalanine and lysine, respectively. This observation correlates with previous work on the same uncapped hexapeptide monomers modeled using atomistic potential. For oligomers, we compare the variation in monomer conformations with radial distance and observe trends for selected angles calculated for each monomer. The LAMMPS interfaces to the GMIN and OPTIM programs for landscape exploration offer new opportunities to investigate larger systems and provide access to the coarse-grained potentials implemented within LAMMPS.

Energy Landscapes and Heat Capacity Signatures for Monomers and Dimers of Amyloid-Forming Hexapeptides

Amyloid formation is a hallmark of various neurodegenerative disorders. In this contribution, energy landscapes are explored for various hexapeptides that are known to form amyloids. Heat capacity (CV) analysis at low temperature for these hexapeptides reveals that the low energy structures contributing to the first heat capacity feature above a threshold temperature exhibit a variety of backbone conformations for amyloid-forming monomers. The corresponding control sequences do not exhibit such structural polymorphism, as diagnosed via end-to-end distance and a dihedral angle defined for the monomer. A similar heat capacity analysis for dimer conformations obtained using basin-hopping global optimisation shows clear features in end-to-end distance versus dihedral correlation plots, where amyloid-forming sequences exhibit a preference for larger end-to-end distances and larger positive dihedrals. These results hold true for sequences taken from tau, amylin, insulin A chain, a de novo designed peptide, and various control sequences. While there is a little overall correlation between the aggregation propensity and the temperature at which the low-temperature CV feature occurs, further analysis suggests that the amyloid-forming sequences exhibit the key CV feature at a lower temperature compared to control sequences derived from the same protein.

RNA Modeling with the Computational Energy Landscape Framework

The recent advances in computational abilities, such as the enormous speed-ups provided by GPU computing, allow for large scale computational studies of RNA molecules at an atomic level of detail. As RNA molecules are known to adopt multiple conformations with comparable energies, but different two-dimensional structures, all-atom models are necessary to better describe the structural ensembles for RNA molecules. This point is important because different conformations can exhibit different functions, and their regulation or mis-regulation is linked to a number of diseases. Problematically, the energy barriers between different conformational ensembles are high, resulting in long time scales for interensemble transitions. The computational potential energy landscape framework was designed to overcome this problem of broken ergodicity by use of geometry optimization. Here, we describe the algorithms used in the energy landscape explorations with the OPTIM and PATHSAMPLE programs, and how they are used in biomolecular simulations. We present a recent case study of the 5'-hairpin of RNA 7SK to illustrate how the method can be applied to interpret experimental results, and to obtain a detailed description of molecular properties.

Crystal Structure Prediction for Benzene Using Basin-Hopping Global Optimization

Organic molecules can be stable in distinct crystalline forms, known as polymorphs, which have significant consequences for industrial applications. Here, we predict the polymorphs of crystalline benzene computationally for an accurate anisotropic model parametrized to reproduce electronic structure calculations. We adapt the basin-hopping global optimization procedure to the case of crystalline unit cells, simultaneously optimizing the molecular coordinates and unit cell parameters to locate multiple low-energy structures from a variety of crystal space groups. We rapidly locate all the well-established experimental polymorphs of benzene, each of which corresponds to a single local energy minimum of the model. Our results show that basin-hopping can be both an efficient and effective tool for polymorphic crystal structure prediction, requiring no a priori experimental knowledge of cell parameters or symmetry.

Atlasing of Assembly Landscapes using Distance Geometry and Graph Rigidity

This Article describes a novel geometric methodology for analyzing free energy and kinetics of assembly driven by short-range pair-potentials in an implicit solvent and provides a proof-of-concept illustration of its unique capabilities. An atlas is a labeled partition of the assembly landscape into a roadmap of maximal, contiguous, nearly-equipotential-energy conformational regions or macrostates, together with their neighborhood relationships. The new methodology decouples the roadmap generation from sampling and produces: (1) a queryable atlas of local potential energy minima, their basin structure, energy barriers, and neighboring basins; (2) paths between a specified pair of basins, each path being a sequence of conformational regions or macrostates below a desired energy threshold; and (3) approximations of relative path lengths, basin volumes (configurational entropy), and path probabilities. Results demonstrating the core algorithm's capabilities and high computational efficiency have been generated by a resource-light, curated open source software implementation EASAL (Efficient Atlasing and Search of Assembly Landscapes, ACM Trans. Math. Softw. 2018 44, 1-48. 10.1145/3204472; see software, Efficient Atlasing and Search of Assembly Landscapes, 2016. https://bitbucket.org/geoplexity/easal; video, Video Illustrating the opensource software EASAL, 2016. https://cise.ufl.edu/~sitharam/EASALvideo.mpeg; and user guide, EASAL software user guide, 2016. https://bitbucket.org/geoplexity/easal/src/master/CompleteUserGuide.pdf). Running on a laptop with Intel(R) Core(TM) i7-7700@3.60 GHz CPU with 16GB of RAM, EASAL atlases several hundred thousand conformational regions or macrostates in minutes using a single compute core. Subsequent path and basin computations each take seconds. A parallelized EASAL version running on the same laptop with 4 cores gives a 3× speedup for atlas generation. The core algorithm's correctness, time complexity, and efficiency-accuracy trade-offs are formally guaranteed using modern distance geometry, geometric constraint systems and combinatorial rigidity. The methodology further links the shape of the input assembling units to a type of intuitive and queryable bar-code of the output atlas, which in turn determine stable assembled structures and kinetics. This succinct input-output relationship facilitates reverse analysis and control toward design. A novel feature that is crucial to both the high sampling efficiency and decoupling of roadmap generation from sampling is a recently developed theory of convex Cayley (distance-based) custom parametrizations specific to assembly, as opposed to folding. Representing microstates with macrostate-specific Cayley parameters, to generate microstate samples, avoids gradient-descent search used by all prevailing methods. Further, these parametrizations convexify conformational regions or macrostates. This ratchets up sampling efficiency, significantly reducing number of repeated and discarded samples. These features of the new stand-alone methodology can also be used to complement the strengths of prevailing methodologies including Molecular Dynamics, Monte Carlo, and Fast Fourier Transform based methods.

Entropy of Proteins Using Multiscale Cell Correlation

A new multiscale method is presented to calculate the entropy of proteins from molecular dynamics simulations. Termed Multiscale Cell Correlation (MCC), the method decomposes the protein into sets of rigid-body units based on their covalent-bond connectivity at three levels of hierarchy: molecule, residue, and united atom. It evaluates the vibrational and topographical entropy from forces, torques, and dihedrals at each level, taking into account correlations between sets of constituent units that together make up a larger unit at the coarser length scale. MCC gives entropies in close agreement with normal-mode analysis and smaller than those using quasiharmonic analysis as well as providing much faster convergence. Moreover, MCC provides an insightful decomposition of entropy at each length scale and for each type of amino acid according to their solvent exposure and whether they are terminal residues. While the residue entropy depends weakly on solvent exposure, there is greater variation in entropy components for larger, more polar amino acids, which have increased conformational entropy but reduced vibrational entropy with greater solvent exposure.

Sphere Encapsulated Monte Carlo: Obtaining Minimum Energy Configurations of Large Aromatic Systems

We introduce a simple global optimization approach that is able to find minimum energy configurations of clusters containing aromatic molecules. The translational and rotational perturbations required in Monte Carlo-based methods often lead to unrealistic configurations within which two or more molecular rings intersect, causing many of the computational steps to be rejected and the optimization process to be inefficient. Here we develop a modification of the basin-hopping global optimization procedure tailored to tackle problems with intersecting molecular rings. Termed the Sphere Encapsulated Monte Carlo (SEMC) method, this method introduces sphere-based rearrangement and minimization steps at each iteration, and its performance is shown through the exploration of potential energy landscapes of polycyclic aromatic hydrocarbon (PAH) clusters, systems of interest in combustion and astrophysics research. The SEMC method provides clusters that are accurate to 5% mean difference of the minimum energy at a 10-fold speed up compared to previous work using advanced molecular dynamics simulations. Importantly, the SEMC method captures key structural characteristics and molecular size partitioning trends as measured by the molecular radial distances and coordination numbers. The advantages of the SEMC method are further highlighted in its application to previously unstudied heterogeneous PAH clusters.

Mutational Basin-Hopping: Combined Structure and Sequence Optimization for Biomolecules

The study of energy landscapes has led to a good understanding of how and why proteins and nucleic acids adopt their native structure. Through evolution, sequences have adapted until they exhibit a strongly funneled energy landscape, stabilizing the native fold. Design of artificial biomolecules faces the challenge of creating similar stable, minimally frustrated, and functional sequences. Here we present a biminimization approach, mutational basin-hopping, in which we simultaneously use global optimization to optimize the energy and a target function describing a desired property of the system. This optimization of structure and sequence is a generalized basin-hopping method and produces an efficient design process, which can target properties such as binding affinity or solubility.

Exploring Energy Landscapes

Recent advances in the potential energy landscapes approach are highlighted, including both theoretical and computational contributions. Treating the high dimensionality of molecular and condensed matter systems of contemporary interest is important for understanding how emergent properties are encoded in the landscape and for calculating these properties while faithfully representing barriers between different morphologies. The pathways characterized in full dimensionality, which are used to construct kinetic transition networks, may prove useful in guiding such calculations. The energy landscape perspective has also produced new procedures for structure prediction and analysis of thermodynamic properties. Basin-hopping global optimization, with alternative acceptance criteria and generalizations to multiple metric spaces, has been used to treat systems ranging from biomolecules to nanoalloy clusters and condensed matter. This review also illustrates how all this methodology, developed in the context of chemical physics, can be transferred to landscapes defined by cost functions associated with machine learning.

Atomic clusters with addressable complexity

A general formulation for constructing addressable atomic clusters is introduced, based on one or more reference structures. By modifying the well depths in a given interatomic potential in favour of nearest-neighbour interactions that are defined in the reference(s), the potential energy landscape can be biased to make a particular permutational isomer the global minimum. The magnitude of the bias changes the resulting potential energy landscape systematically, providing a framework to produce clusters that should self-organise efficiently into the target structure. These features are illustrated for small systems, where all the relevant local minima and transition states can be identified, and for the low-energy regions of the landscape for larger clusters. For a 55-particle cluster, it is possible to design a target structure from a transition state of the original potential and to retain this structure in a doubly addressable landscape. Disconnectivity graphs based on local minima that have no direct connections to a lower minimum provide a helpful way to visualise the larger databases. These minima correspond to the termini of monotonic sequences, which always proceed downhill in terms of potential energy, and we identify them as a class of biminimum. Multiple copies of the target cluster are treated by adding a repulsive term between particles with the same address to maintain distinguishable targets upon aggregation. By tuning the magnitude of this term, it is possible to create assemblies of the target cluster corresponding to a variety of structures, including rings and chains.

Exploring biomolecular energy landscapes

The potential energy landscape perspective provides both a conceptual and a computational framework for predicting, understanding and designing molecular properties. In this Feature Article, we highlight some recent advances that greatly facilitate structure prediction and analysis of global thermodynamics and kinetics in proteins and nucleic acids. The geometry optimisation procedures, on which these calculations are based, can be accelerated significantly using local rigidification of selected degrees of freedom, and through implementations on graphics processing units. Results of progressive local rigidification are first summarised for trpzip1, including a systematic analysis of the heat capacity and rearrangement rates. Benchmarks for all the essential optimisation procedures are then provided for a variety of proteins. Applications are then illustrated from a study of how mutation affects the energy landscape for a coiled-coil protein, and for transitions in helix morphology for a DNA duplex. Both systems exhibit an intrinsically multifunnel landscape, with the potential to act as biomolecular switches.

DynaDom: structure-based prediction of T cell receptor inter-domain and T cell receptor-peptide-MHC (class I) association angles

Background

T cell receptor (TCR) molecules are involved in the adaptive immune response as they distinguish between self- and foreign-peptides, presented in major histocompatibility complex molecules (pMHC). Former studies showed that the association angles of the TCR variable domains (Vα/Vβ) can differ significantly and change upon binding to the pMHC complex. These changes can be described as a rotation of the domains around a general Center of Rotation, characterized by the interaction of two highly conserved glutamine residues.

Methods

We developed a computational method, DynaDom, for the prediction of TCR Vα/Vβ inter-domain and TCR/pMHC orientations in TCRpMHC complexes, which allows predicting the orientation of multiple protein-domains. In addition, we implemented a new approach to predict the correct orientation of the carboxamide endgroups in glutamine and asparagine residues, which can also be used as an external, independent tool.

Results

The approach was evaluated for the remodeling of 75 and 53 experimental structures of TCR and TCRpMHC (class I) complexes, respectively. We show that the DynaDom method predicts the correct orientation of the TCR Vα/Vβ angles in 96 and 89% of the cases, for the poses with the best RMSD and best interaction energy, respectively. For the concurrent prediction of the TCR Vα/Vβ and pMHC orientations, the respective rates reached 74 and 72%. Through an exhaustive analysis, we could show that the pMHC placement can be further improved by a straightforward, yet very time intensive extension of the current approach.

Conclusions

The results obtained in the present remodeling study prove the suitability of our approach for interdomain-angle optimization. In addition, the high prediction rate obtained specifically for the energetically highest ranked poses further demonstrates that our method is a powerful candidate for blind prediction. Therefore it should be well suited as part of any accurate atomistic modeling pipeline for TCRpMHC complexes and potentially other large molecular assemblies.

Electronic supplementary material

The online version of this article (doi:10.1186/s12900-016-0071-7) contains supplementary material, which is available to authorized users.

Geometry optimization made simple with translation and rotation coordinates

The effective description of molecular geometry is important for theoretical studies of intermolecular interactions. Here we introduce a new translation-rotation-internal coordinate (TRIC) system which explicitly includes the collective translations and rotations of molecules, or parts of molecules such as monomers or ligands, as degrees of freedom. The translations are described as the centroid position and the orientations are represented with the exponential map parameterization of quaternions. When TRIC is incorporated into geometry optimization calculations, the performance is consistently superior to existing coordinate systems for a diverse set of systems including water clusters, organic semiconductor donor-acceptor complexes, and small proteins, all of which are characterized by nontrivial intermolecular interactions. The method also introduces a new way to scan the molecular orientations while allowing orthogonal degrees of freedom to relax. Our findings indicate that an explicit description of molecular translation and rotation is a natural way to traverse the many-dimensional potential energy surface.

Energy Minimization on Manifolds for Docking Flexible Molecules

In this paper, we extend a recently introduced rigid body minimization algorithm, defined on manifolds, to the problem of minimizing the energy of interacting flexible molecules. The goal is to integrate moving the ligand in six dimensional rotational/translational space with internal rotations around rotatable bonds within the two molecules. We show that adding rotational degrees of freedom to the rigid moves of the ligand results in an overall optimization search space that is a manifold to which our manifold optimization approach can be extended. The effectiveness of the method is shown for three different docking problems of increasing complexity. First, we minimize the energy of fragment-size ligands with a single rotatable bond as part of a protein mapping method developed for the identification of binding hot spots. Second, we consider energy minimization for docking a flexible ligand to a rigid protein receptor, an approach frequently used in existing methods. In the third problem, we account for flexibility in both the ligand and the receptor. Results show that minimization using the manifold optimization algorithm is substantially more efficient than minimization using a traditional all-atom optimization algorithm while producing solutions of comparable quality. In addition to the specific problems considered, the method is general enough to be used in a large class of applications such as docking multidomain proteins with flexible hinges. The code is available under open source license (at http://cluspro.bu.edu/Code/Code_Rigtree.tar) and with minimal effort can be incorporated into any molecular modeling package.

Prediction and clarification of structures of (bio)molecules on surfaces

The design of future materials for biotechnological applications via deposition of molecules on surfaces will require not only exquisite control of the deposition procedure, but of equal importance will be our ability to predict the shapes and stability of individual molecules on various surfaces. Furthermore, one will need to be able to predict the structure patterns generated during the self-organization of whole layers of (bio)molecules on the surface. In this review, we present an overview over the current state of the art regarding the prediction and clarification of structures of biomolecules on surfaces using theoretical and computational methods.

Determining the Energetics of Small β-Sheet Peptides using Adaptive Steered Molecular Dynamics

Mechanically driven unfolding is a useful computational tool for extracting the energetics and stretching pathway of peptides. In this work, two representative β-hairpin peptides, chignolin (PDB: 1UAO ) and trpzip1 (PDB: 1LE0 ), were investigated using an adaptive variant of the original steered molecular dynamics method called adaptive steered molecular dynamics (ASMD). The ASMD method makes it possible to perform energetic calculations on increasingly complex biological systems. Although the two peptides are similar in length and have similar secondary structures, their unfolding energetics are quite different. The hydrogen bonding profile and specific residue pair interaction energies provide insight into the differing stabilities of these peptides and reveal which of the pairs provides the most significant stabilization.

Energy Landscapes and Global Optimization of Self-Assembling Cyclic Peptides

Self-assembled cyclic peptide nanotubes have attracted much attention because of their antimicrobial properties. Here, we present calculations on the formation of cyclic peptide dimers using basin-hopping and discrete path sampling. We present an analysis of the basin-hopping move sets that most efficiently explore the conformations of cyclic peptides. Group rotation moves, in which sections of the ring are rotated as a rigid body, are the most effective for cyclic peptides containing up to 20 residues. For cyclic peptide dimers, we find that a combination of group rotation intramolecular moves and rigid body intermolecular moves performs well. Discrete path sampling calculations on the cyclic peptide dimers show significant differences in the dimerization of hexa- and octapeptides.

Equilibrium molecular thermodynamics from Kirkwood sampling

We present two methods for barrierless equilibrium sampling of molecular systems based on the recently proposed Kirkwood method (J. Chem. Phys.2009, 130, 134102). Kirkwood sampling employs low-order correlations among internal coordinates of a molecule for random (or non-Markovian) sampling of the high dimensional conformational space. This is a geometrical sampling method independent of the potential energy surface. The first method is a variant of biased Monte Carlo, where Kirkwood sampling is used for generating trial Monte Carlo moves. Using this method, equilibrium distributions corresponding to different temperatures and potential energy functions can be generated from a given set of low-order correlations. Since Kirkwood samples are generated independently, this method is ideally suited for massively parallel distributed computing. The second approach is a variant of reservoir replica exchange, where Kirkwood sampling is used to construct a reservoir of conformations, which exchanges conformations with the replicas performing equilibrium sampling corresponding to different thermodynamic states. Coupling with the Kirkwood reservoir enhances sampling by facilitating global jumps in the conformational space. The efficiency of both methods depends on the overlap of the Kirkwood distribution with the target equilibrium distribution. We present proof-of-concept results for a model nine-atom linear molecule and alanine dipeptide.

Exploring energy landscapes: from molecular to mesoscopic systems

We review a comprehensive computational framework to survey the potential energy landscape for systems composed of rigid or partially rigid molecules. Illustrative case studies relevant to a wide range of molecular clusters and soft and condensed matter systems are discussed.

A conformational factorisation approach for estimating the binding free energies of macromolecules

We present a conformational factorization approach. The theory is based on a superposition partition function, decomposed as a sum over contributions from local minima. The factorisation greatly reduces the number of minima that need to be considered, by employing the same local configurations for groups that are sufficiently distant from the binding site. The theory formalises the conditions required to analyse how our definition of the binding site region affects the free energy difference between the apo and holo states. We employ basin-hopping parallel tempering to sample minima that contribute significantly to the partition function, and calculate the binding free energies within the harmonic normal mode approximation. A further significant gain in efficiency is achieved using a recently developed local rigid body framework in both the sampling and the normal mode analysis, which reduces the number of degrees of freedom. We benchmark this approach for human aldose reductase (PDB code 2INE). When varying the size of the rigid region, the free energy difference converges for factorisation of groups at a distance of 14 Å from the binding site, which corresponds to 80% of the protein being locally rigidified.

Exploring Energy Landscapes: Metrics, Pathways, and Normal-Mode Analysis for Rigid-Body Molecules

We present new methodology for exploring the energy landscapes of molecular systems, using angle-axis variables for the rigid-body rotational coordinates. The key ingredient is a distance measure or metric tensor, which is invariant to global translation and rotation. The metric is used to formulate a generalized nudged elastic band method for calculating pathways, and a full prescription for normal-mode analysis is described. The methodology is tested by mapping the potential energy and free energy landscape of the water octamer, described by the TIP4P potential.