Enhancing backbone sampling in Monte Carlo simulations using Internal Coordinates Normal Mode Analysis

Normal mode methods are becoming a popular alternative to sample the conformational landscape of proteins. In this study, we describe the implementation of an internal coordinate normal mode analysis method and its application in exploring protein flexibility by using the Monte Carlo method PELE. This new method alternates two different stages, a perturbation of the backbone through the application of torsional normal modes, and a resampling of the side chains. We have evaluated the new approach using two test systems, ubiquitin and c-Src kinase, and the differences to the original ANM method are assessed by comparing both results to reference molecular dynamics simulations. The results suggest that the sampled phase space in the internal coordinate approach is closer to the molecular dynamics phase space than the one coming from a Cartesian coordinate anisotropic network model. In addition, the new method shows a great speedup (∼∼5-7x), making it a good candidate for future normal mode implementations in Monte Carlo methods.The authors thank D. E. Shaw Research lab. for providing the kinase MD coordinates and Dr. López Blanco for sharing the code developed in his thesis and for providing useful comments. This work was supported by the CTQ-48287-R projects of the Spanish Ministry of Economy and Competitiveness (MINECO) and the grant SEV-2011-00067 of Severo Ochoa Program, awarded by the Spanish Government.Peer ReviewedPostprint (author's final draft

Ligand binding mechanism in steroid receptors; from conserved plasticity to differential evolutionary constraints

Steroid receptor drugs have been available for more than half a century, but details 24 of the ligand binding mechanism has remained elusive. We solved X-ray structures of 25 the glucocorticoid and mineralocorticoid receptors to identify a conserved plasticity at 26 helix 6-7 region that extend the ligand binding pocket towards the receptor surface. 27 Since none of the endogenous ligands exploit this region, we hypothesized that it 28 constitutes an integral part of the binding event. Extensive all atom unbiased ligand 29 exit and entrance simulations corroborate a ligand binding pathway that gives the 30 observed structural plasticity a key functional role. Kinetic measurements reveal that 31 the receptor residence time correlate with structural rearrangements observed in both 32 structures and simulations. Ultimately, our findings reveal why nature has conserved 33 the capacity to open up this region and highlight how differences in the details of the 34 ligand entry process result in differential evolutionary constraints across the steroid 35 receptors.This study was supported by The European Research Council (2009-Adg25027-535 PELE) to V.G and by the SEV-2011-00067 grant of the Severo Ochoa Program. We 536 would like to acknowledge our AstraZeneca colleagues J. Hartleib, R.Unwin and 537 R.Knöll for helpful discussions. We also thank N. Blomberg (ELIXIR) and R. Neutze 538 (University of Gothenburg) for careful reading of the manuscript.Peer ReviewedPostprint (author's final draft

Neue Tabu-Search Algorithmen zur Untersuchung von Energielandschaften molekularer Systeme

The visualization of energy functions is based on the possibility of separating different degrees of freedom. The most important one is the Born-Oppenheimer-approximation, which separates nucleus and electron movements. This allows the illustration of the potential energy as a function of the nuclei coordinates. Minima of the surface correspond to stable points like isomers or conformers. They are important for predicting the stability or thermodynamical of a system. Stationary points of first order correspond to transition points. They describe phase transitions, chemical reaction, or conformational changes. Furthermore, the partition function connects the potential hypersurface to the free energy of the system. The aim of the present work is the development and application of new approaches for the efficient exploration of multidimensional hypersurfaces. Initially, the Conformational Analysis and Search Tool (CAST) program was developed to create a basis for the new methods and algorithms. The development of CAST in object oriented C++ included, among other things, the implementation of a force field, different interfaces to external programs, analysis tools, and optimization libraries. Descriptions of an energy landscape require knowledge about the most stable minima. The Gradient Only Tabu Search (GOTS) has been shown to be very efficient in the optimization of mathematical test functions. Therefore, GOTS was taken as a starting point. Tabu-Search is based on the steepest descent - modest ascent strategy. The steepest descent is used for finding local minima, while the modest ascent is taken for leaving a minimum quickly. Furthermore, Tabu-Search is combined with an adaptive memory design to avoid cycling or returning. The highly accurate exploration of the phase space by Tabu-Search is often too expensive for complex optimization problems. Therefore, an algorithm for diversification of the search is required. After exploration of the proximity of the search space, the algorithm would guide the search to new and hopefully promising parts of the phase space. First application of GOTS to conformational search revealed weaknesses in the diversification search and the modest ascent part. On the one hand, the original methodology for diversification is insufficiently diverse. The algorithm is considerably improved by combining the more local GOTS with the wider searching Basin Hopping (BH) approach. The second weak point is a too inaccurate and inefficient modest ascent strategy. Analysis of common transition state search algorithms lead to the adaption of the Dimer-method to the Tabu-Search approach. The Dimer-method only requires the first derivatives for locating the closest transition state. For conformational search, dihedral angles are usually the most flexible degrees of freedom. Therefore, only those are used in the Dimer-method for leaving a local minimum. Furthermore, the exact localization of the reaction pathway and the transition state is not necessary as the local minimum position should only be departed as fast as possible. This allows for larger step sizes during the Dimer-search. In the following optimization step, all coordinates are relaxed to remove possible strains in the system. The new Tabu-Search method with Dimer-search delivers more and improved minima. Furthermore, the approach is faster for larger systems. For a system with approximately 1200 atoms, an acceleration of 40 was measured. The new approach was compared to Molecular Dynamics with optimization (MD), Simulated Annealing (SA), and BH with the help of conformational search problems of bio-organic systems. In all cases, a better performance was found. A comparison to the Monte Carlo Multiple Minima/Low Mode Sampling (MCMM/LM) method proved the outstanding performance of the new Tabu-Search approach. The solvation of the chignolin protein further revealed the possibility of uncovering discrepancies between the employed theoretical model and the experimental starting structure. Ligand optimization for improvement of x-ray structures was one further new application field. Besides the global optimization, the search for transition states and reaction pathways is also of paramount importance. These points describe different transitions of stable states. Therefore, a new approach for the exploration of such cases was developed. The new approach is based on a global minimization of a hyperplane being perpendicular to the reaction coordinate. Minima of this reduced phase space belong to traces of transition states between reactant and product states on the unchanged hypersurface. Optimization to the closest transition state using the Dimer-method delivers paths lying between the initial and the final state. An iterative approach finally yields complex reaction pathways with many intermediate local minima. The PathOpt algorithm was tested by means of rearrangements of argon clusters showing very promising results.Die visuelle Darstellung von Energiefunktionen basiert auf der Möglichkeit, verschiedene Freiheitsgrade zu separieren. Die wichtigste Näherung ist dabei die Born-Oppenheimer-Näherung. Sie erlaubt damit die Darstellung der potentiellen Energie als Funktion der Kernkoordinaten. Die daraus entstehende mehrdimensionale Hyperfläche entspricht der Summenformel eines beliebigen Systems. Minima der Fläche entsprechen stabilen Punkten wie Isomeren oder Konformeren. Diese sind wichtig für Aussagen über die Stabilität oder die Thermodynamik eines Systems. Stationäre Punkte erster Ordnung entsprechen Übergangsstrukturen und beschreiben Phasenübergänge, chemische Reaktionen aber auch Konformationsänderungen. Über die Zustandssumme ist die Hyperfläche zudem mit der freien Energie verknüpft. Das Ziel dieser Arbeit ist die Entwicklung und Anwendung neuer Methoden zur effizienten Untersuchung mehrdimensionaler Hyperflächen. Dabei wurde zunächst das Conformational Analysis and Search Tool (CAST)-Programm entwickelt. Die Entwicklung des CAST-Programms in objektorientiertem C++ beinhaltete unter anderem die Implementierung eines Kraftfeldes, verschiedene Schnittstellen zu externen Programmen, Analysealgorithmen und verschiedene Optimierungsmodule. Um Aussagen über eine Energielandschaft treffen zu können, müssen zuerst die stabilsten Minima gefunden werden. Der Gradient Only Tabu Search (GOTS) hat sich als sehr effizient in der Optimierung von mathematischen Funktionen erwiesen. Daher wurde GOTS als Startpunkt verwendet. Tabu-Search basiert auf dem steepest descent – modest ascent Prinzip. Zum Finden neuer Minima wird der steilste Abstieg (steepest descent) verwendet, ein Minimum wird auf dem Weg des geringsten Anstiegs (modest ascent) wieder verlassen. Tabu-Search ist zudem mit einem lernfähigen Speicherdesign kombiniert, wodurch ein Zurück- und im Kreis laufen vermieden wird. Der Phasenraum wird von Tabu-Search sehr genau untersucht, was für komplexere Probleme zu aufwendig wird. Daher bedarf es eines Diversifizierungsschritts, welcher nach Absuchen eines Teils des Phasenraums, die Suche in neue vielversprechende Bereiche bringt. Erste Anwendungen auf Konformationssuchen zeigten, dass GOTS Schwächen im Diversifizierungsschritt und der modest ascent Strategie besitzt. Zum einen ist die ursprünglich verwendete Methodik für die Diversifizierung zu wenig divers. Eine Kombination des mehr lokalen GOTS mit der weiträumiger suchenden Basin Hopping (BH) Methode brachte eine erhebliche Verbesserung. Der zweite Schwachpunkt besteht aus einer zu ungenauen und ineffizienten modest ascent Methode. Daher wurde die Dimer-Methode für Tabu-Search adaptiert. Diese benötigt lediglich die erste Ableitung, um zum Übergangszustand erster Ordnung zu konvergieren. Dabei werden in der Dimer-Methode nur Diederwinkel variiert. Zudem muss der Reaktionspfad und der Übergangszustand nicht exakt getroffen werden, da das Minimum nur möglichst schnell verlassen werden soll. Dies erlaubt größere Schrittweiten in der Dimer-Suche. Im nachfolgenden Optimierungsschritt werden alle Koordinaten relaxiert. Die neue Tabu-Search-Methode mit Dimer-Suche liefert mehr und deutlich verbesserte Minima. Zudem ist sie für größere Systeme deutlich schneller. Für ein System mit circa 1200 Atomen wurde eine Beschleunigung um den Faktor 40 erzielt. Die neue Methode wurde am Beispiel der Konformationssuche von bio-organischen Systemen mit Molekulardynamik mit Optimierung (MD), Simulated Annealing (SA) und BH verglichen, wobei sich in allen Fällen eine bessere Effizienz zeigte. Ein Vergleich zur Monte Carlo Multiple Minima/Low Mode Sampling Methode anhand der Optimierung von peptidischen Ligand-Rezeptor-Komplexen belegte ebenfalls die hervorragende Effizienz des neuen Ansatzes. Die Solvatisierung des Chignolin-Proteins mit Tabu-Search deckte die Möglichkeit auf, Differenzen zwischen der verwendeten theoretischen Methode und der experimentellen Startstruktur aufzudecken. Als weiterer neuer Anwendungsbereich wurde die Optimierung von Ligand-Enzym-Komplexen zur Verbesserung von Röntgenstrukturen untersucht. Neben der globalen Optimierung ist auch die Suche nach Übergangszuständen und Reaktionspfaden von größter Wichtigkeit. Diese beschreiben verschiedene Übergänge zwischen stabilen Zuständen. Daher wurde ein neuer Ansatz zur Untersuchung dieser Fragestellungen entwickelt. Dieser basiert auf einer globalen Minimierung einer Hyperfläche, welche senkrecht zum Reaktionspfad steht. Die Minima des reduzierten Phasenraums gehören zu Spuren zu Übergangszuständen zwischen dem Edukt und dem Produkt. Durch Optimierung dieser Punkte mittels der Dimer-Methode werden also Pfade gefunden, die zwischen Anfangs- und Endpunkt liegen. Ein iteratives Vorgehen liefert letztendlich komplexe Reaktionspfade. PathOpt wurde an Umlagerungen von Argon-Clustern evaluiert, welche sehr vielversprechende Ergebnisse lieferten

Molecular Rift: Virtual Reality for Drug Designers

Recent advances in interaction design have created new ways to use computers. One example is the ability to create enhanced 3D environments that simulate physical presence in the real world a virtual reality. This is relevant to drug discovery since molecular models are frequently used to obtain deeper understandings of, say, ligand protein complexes. We have developed a tool (Molecular Rift), which creates a virtual reality environment steered with hand movements. Oculus Rift, a head-mounted display, is used to create the virtual settings. The program is controlled by gesture-recognition, using the gaming sensor MS Kinect v2, eliminating the need for standard input devices. The Open Babel toolkit was integrated to provide access to powerful cheminformatics functions. Molecular Rift was developed with a focus on usability, including iterative test-group evaluations. We conclude with reflections on virtual reality's future capabilities in chemistry and education. Molecular Rift is open source and can be downloaded from GitHub

Enhancing backbone sampling in Monte Carlo simulations using Internal Coordinates Normal Mode Analysis

Normal mode methods are becoming a popular alternative to sample the conformational landscape of proteins. In this study, we describe the implementation of an internal coordinate normal mode analysis method and its application in exploring protein flexibility by using the Monte Carlo method PELE. This new method alternates two different stages, a perturbation of the backbone through the application of torsional normal modes, and a resampling of the side chains. We have evaluated the new approach using two test systems, ubiquitin and c-Src kinase, and the differences to the original ANM method are assessed by comparing both results to reference molecular dynamics simulations. The results suggest that the sampled phase space in the internal coordinate approach is closer to the molecular dynamics phase space than the one coming from a Cartesian coordinate anisotropic network model. In addition, the new method shows a great speedup (∼∼5-7x), making it a good candidate for future normal mode implementations in Monte Carlo methods.The authors thank D. E. Shaw Research lab. for providing the kinase MD coordinates and Dr. López Blanco for sharing the code developed in his thesis and for providing useful comments. This work was supported by the CTQ-48287-R projects of the Spanish Ministry of Economy and Competitiveness (MINECO) and the grant SEV-2011-00067 of Severo Ochoa Program, awarded by the Spanish Government.Peer Reviewe

A New Tabu-Search-Based Algorithm for Solvation of Proteins

The proper description of explicit water shells is of enormous importance for all-atom calculations. We propose a new approach for the setup of water shells around proteins based on Tabu-Search global optimization and compare its efficiency with standard molecular dynamics protocols using the chignolin protein as a test case. Both algorithms generate reasonable water shells, but the new approach provides solvated systems with an increased water–enzyme interaction and offers further advantages. It enables a stepwise buildup of the solvent shell, so that the more important inner part can be prepared more carefully. It also allows the generation of solute structures which can be biased either toward the (experimental) starting structure or the underlying theoretical model, i.e., the employed force field

Impact of applicability domains to generative artificial intelligence

Molecular generative artificial intelligence is drawing significant attention in the drug design community, with several experimentally validated proofs of concepts already published. Nevertheless, generative models are known for sometimes generating unrealistic, unsynthesizable or unstable structures. This calls for methods to constrain those algorithms to generate structures in reasonable portions of the chemical space. While the concept of applicability domains (AD) for predictive models is well studied, its counterpart for generative models is not yet defined. In this work, we examine empirically various possibilities and propose applicability domains suited for generative models. Using both public and internal datasets, we use state-of-the-art generative methods to generate novel structures that are predicted actives by a corresponding QSAR model, while constraining the generative model to stay within a given applicability domain. Our work looks at several applicability domain definitions, combining various criteria, such as structural similarity to the training set, similarity of physico-chemical properties, unwanted substructures, and Quantitative Estimate of Drug- Likeness (QED). We assess both from a qualitative and quantitative point of view the structures generated, and find that the applicability domain definitions have a strong influence on the chemical beauty of generated molecules. An extensive analysis of our results allows us to identify applicability domain definitions that are best suited for generating drug-like molecules with generative models. We anticipate that this work will help foster the adoption of generative models in an industrial context

Binding Mode and Induced Fit Predictions for Prospective Computational Drug Design

Computer-aided drug design plays an important role in medicinal chemistry to obtain insights into molecular mechanisms and to prioritize design strategies. Although significant improvement has been made in structure based design, it still remains a key challenge to accurately model and predict induced fit mechanisms. Most of the current available techniques either do not provide sufficient protein conformational sampling or are too computationally demanding to fit an industrial setting. The current study presents a systematic and exhaustive investigation of predicting binding modes for a range of systems using PELE (Protein Energy Landscape Exploration), an efficient and fast protein–ligand sampling algorithm. The systems analyzed (cytochrome P, kinase, protease, and nuclear hormone receptor) exhibit different complexities of ligand induced fit mechanisms and protein dynamics. The results are compared with results from classical molecular dynamics simulations and (induced fit) docking. This study shows that ligand induced side chain rearrangements and smaller to medium backbone movements are captured well in PELE. Large secondary structure rearrangements, however, remain challenging for all employed techniques. Relevant binding modes (ligand heavy atom RMSD < 1.0 Å) can be obtained by the PELE method within a few hours of simulation, positioning PELE as a tool applicable for rapid drug design cycles