Computational structure‐based drug design: Predicting target flexibility

The role of molecular modeling in drug design has experienced a significant revamp in the last decade. The increase in computational resources and molecular models, along with software developments, is finally introducing a competitive advantage in early phases of drug discovery. Medium and small companies with strong focus on computational chemistry are being created, some of them having introduced important leads in drug design pipelines. An important source for this success is the extraordinary development of faster and more efficient techniques for describing flexibility in three‐dimensional structural molecular modeling. At different levels, from docking techniques to atomistic molecular dynamics, conformational sampling between receptor and drug results in improved predictions, such as screening enrichment, discovery of transient cavities, etc. In this review article we perform an extensive analysis of these modeling techniques, dividing them into high and low throughput, and emphasizing in their application to drug design studies. We finalize the review with a section describing our Monte Carlo method, PELE, recently highlighted as an outstanding advance in an international blind competition and industrial benchmarks.We acknowledge the BSC-CRG-IRB Joint Research Program in Computational Biology. This work was supported by a grant from the Spanish Government CTQ2016-79138-R.J.I. acknowledges support from SVP-2014-068797, awarded by the Spanish Government.Peer ReviewedPostprint (author's final draft

Protein kinases: Structure modeling, inhibition, and protein-protein interactions

Human protein kinases belong to a large and diverse enzyme family that contains more than 500 members. Deregulation of protein kinases is associated with many disorders, and this is why protein kinases are attractive targets for drug discovery. Due to the high conservation of the ATP binding pocket among this family, designing specific and/or selective inhibitors against certain member(s) is challenging. Several studies have been conducted on protein kinases to validate them as suitable drug targets. Although there are numerous target-validated protein kinases, the efforts to develop small molecule inhibitors have so far led to only a limited number of therapeutic agents and drug candidates. In our studies, we tried to understand the basic structural features of protein kinases using available computational tools. There are wide structural variations between different states of the same protein kinase that affect the enzyme specificity and inhibition. Many protein kinases do not yet have an available X-ray crystal structure and have not yet been validated to be drug targets. For these reasons, we developed a new homology modeling approach to facilitate modeling non-crystallized protein kinases and protein kinase states. Our homology modeling approach was able to model proteins having long amino acid sequences and multiple protein domains with reliable model quality and a manageable amount of computational time. Then, we checked the applicability of different docking algorithms (the routinely used computational methodology in virtual screening) in protein kinase studies. After performing the basic study of kinase structure modeling, we focused our research on cyclin dependent kinase 2 (CDK2) and glycogen synthase kinase-3β (GSK-3β). We prepared a non-redundant database from 303 available CDK2 PDB structures. We removed all structural anomalies and proceeded to use the CDK2 database in studying CDK2 structure in its different states, upon ATP, ligand and cyclin binding. We clustered the database based on our findings, and the CDK2 clusters were used to generate protein ligand interaction fingerprints (PLIF). We generated a PLIF-based pharmacophore model which is highly selective for CDK2 ligands. A virtual screening workflow was developed making use of the PLIF-based pharmacophore model, ligand fitting into the CDK2 active site and selective CDK2 shape scoring. We studied the structural basis for selective inhibition of CDK2 and GSK-3β. We compared the amino acid sequence, the 3D features, the binding pockets, contact maps, structural geometry, and Sphoxel maps. From this study we found 1) the ligand structural features that are required for the selective inhibition of CDK2 and GSK-3β, and 2) the amino acid residues which are essential for ligand binding and selective inhibition. We used the findings of this study to design a virtual screening workflow to search for selective inhibitors for CDK2 and GSK-3β. Because protein–protein interactions are essential in the function of protein kinases, and in particular of CDK2, we used protein–protein docking knowledge and binding energy calculations to examine CDK2 and cyclin binding. We applied this study to the voltage dependent calcium channel 1 (VDAC1) binding to Bax. We were able to provide important data relevant to future experimental researchers such as on the possibility of Bax to cross biological membranes and the most relevant amino acid residues in VDAC1 that interact with Bax

Identifying Ligand Binding Conformations of the β2-Adrenergic Receptor by Using Its Agonists as Computational Probes

Recently available G-protein coupled receptor (GPCR) structures and biophysical studies suggest that the difference between the effects of various agonists and antagonists cannot be explained by single structures alone, but rather that the conformational ensembles of the proteins need to be considered. Here we use an elastic network model-guided molecular dynamics simulation protocol to generate an ensemble of conformers of a prototypical GPCR, β2-adrenergic receptor (β2AR). The resulting conformers are clustered into groups based on the conformations of the ligand binding site, and distinct conformers from each group are assessed for their binding to known agonists of β2AR. We show that the select ligands bind preferentially to different predicted conformers of β2AR, and identify a role of β2AR extracellular region as an allosteric binding site for larger drugs such as salmeterol. Thus, drugs and ligands can be used as "computational probes" to systematically identify protein conformers with likely biological significance. © 2012 Isin et al

Insights into the Development of Chemotherapeutics Targeting PFKFB Enzymes

The PFKFB enzymes control the primary checkpoint in the glycolytic pathway and are implicated in a multitude of diseases: from cancer, to schizophrenia, to diabetes, and heart disease. The inducible isoform, PFKFB3, is known to be associated with the upregulation of glycolysis in many cancers. The first study within this work investigates the potential for using tier-based approaches of virtual screening to target small molecule kinases, with PFKFB3 serving as a case study. For this investigation, bioactive compounds for PFKFB3 were identified from a compound library of 1364 compounds via high-throughput screening, with bioactive compounds being further characterized as either competitive or non-competitive for F6P. Using the F6P-competitive compounds, several structure based docking programs were assessed individually and in conjunction with a pharmacophore screening. The results showed that the tiered virtual screening approach, using pharmacophore screening in addition to structure-based docking, improved enrichments rates in 80% of cases, reduced CPU costs up to 7-fold, and lessened variability among different structure-based docking methods. The second study investigates the structural and kinetic characteristics of citrate inhibition on the heart PFKFB isoenzyme, PFKFB2. High levels of citrate, an intermediate of the TCA cycle, signify an abundance of biosynthetic precursors and that additional glucose need not be degraded for this purpose. Previous studies have noted that citrate acts as an important negative feed-back mechanism to limit glycolytic activity by inhibiting PFKFB enzymes, yet the structural and mechanistic details of citrate’s inhibition had not been determined. To study the molecular basis for citrate inhibition, the three-dimensional structures of the human and bovine PFKFB2 orthologues were solved, each in complex with citrate. For both cases, citrate primarily occupied the binding site of Fructose-6-phosphate (F6P), competitively blocking F6P from binding. Additionally, a carboxy arm of citrate extended into the γ-phosphate binding site of ATP, sterically and electrostatically blocking the catalytic binding mode for ATP. In the human orthologue, which utilized AMPPNP as an ATP analogue, conformational changes were observed in the 2-kinase domain as well as the binding mode for AMPPNP. This study gives new insights as to how the citrate-mediate negative feedback loop influences glycolytic flux through PFKFB enzymes

Insights from Free-Energy Calculations: Protein Conformational Equilibrium, Driving Forces, and Ligand-Binding Modes

AbstractAccurate free-energy calculations provide mechanistic insights into molecular recognition and conformational equilibrium. In this work, we performed free-energy calculations to study the thermodynamic properties of different states of molecular systems in their equilibrium basin, and obtained accurate absolute binding free-energy calculations for protein-ligand binding using a newly developed M2 algorithm. We used a range of Asp-Phe-Gly (DFG)-in/out p38α mitogen-activated protein kinase inhibitors as our test cases. We also focused on the flexible DFG motif, which is closely connected to kinase activation and inhibitor binding. Our calculations explain the coexistence of DFG-in and DFG-out states of the loop and reveal different components (e.g., configurational entropy and enthalpy) that stabilize the apo p38α conformations. To study novel ligand-binding modes and the key driving forces behind them, we computed the absolute binding free energies of 30 p38α inhibitors, including analogs with unavailable experimental structures. The calculations revealed multiple stable, complex conformations and changes in p38α and inhibitor conformations, as well as balance in several energetic terms and configurational entropy loss. The results provide relevant physics that can aid in designing inhibitors and understanding protein conformational equilibrium. Our approach is fast for use with proteins that contain flexible regions for structure-based drug design

Development of a normal mode-based geometric simulation approach for investigating the intrinsic mobility of proteins

Specific functions of biological systems often require conformational transitions of macromolecules. Thus, being able to describe and predict conformational changes of biological macromolecules is not only important for understanding their impact on biological function, but will also have implications for the modelling of (macro)molecular complex formation and in structure-based drug design approaches. The “conformational selection model” provides the foundation for computational investigations of conformational fluctuations of the unbound protein state. These fluctuations may reveal conformational states adopted by the bound proteins. The aim of this work is to incorporate directional information in a geometry-based approach, in order to sample biologically relevant conformational space extensively. Interestingly, coarse-grained normal mode (CGNM) approaches, e.g., the elastic network model (ENM) and rigid cluster normal mode analysis (RCNMA), have emerged recently and provide directions of intrinsic motions in terms of harmonic modes (also called normal modes). In my previous work and in other studies it has been shown that conformational changes upon ligand binding occur along a few low-energy modes of unbound proteins and can be efficiently calculated by CGNM approaches. In order to explore the validity and the applicability of CGNM approaches, a large-scale comparison of essential dynamics (ED) modes from molecular dynamics (MD) simulations and normal modes from CGNM was performed over a dataset of 335 proteins. Despite high coarse-graining, low frequency normal modes from CGNM correlate very well with ED modes in terms of directions of motions (average maximal overlap is 0.65) and relative amplitudes of motions (average maximal overlap is 0.73). In order to exploit the potential of CGNM approaches, I have developed a three-step approach for efficient exploration of intrinsic motions of proteins. The first two steps are based on recent developments in rigidity and elastic network theory. Initially, static properties of the protein are determined by decomposing the protein into rigid clusters using the graph-theoretical approach FIRST at an all-atom representation of the protein. In a second step, dynamic properties of the molecule are revealed by the rotations-translations of blocks approach (RTB) using an elastic network model representation of the coarse-grained protein. In the final step, the recently introduced idea of constrained geometric simulations of diffusive motions in proteins is extended for efficient sampling of conformational space. Here, the low-energy (frequency) normal modes provided by the RCNMA approach are used to guide the backbone motions. The NMSim approach was validated on hen egg white lysozyme by comparing it to previously mentioned simulation methods in terms of residue fluctuations, conformational space explorations, essential dynamics, sampling of side-chain rotamers, and structural quality. Residue fluctuations in NMSim generated ensemble is found to be in good agreement with MD fluctuations with a correlation coefficient of around 0.79. A comparison of different geometry-based simulation approaches shows that FRODA is restricted in sampling the backbone conformational space. CONCOORD is restricted in sampling the side-chain conformational space. NMSim sufficiently samples both the backbone and the side-chain conformations taking experimental structures and conformations from the state of the art MD simulation as reference. The NMSim approach is also applied to a dataset of proteins where conformational changes have been observed experimentally, either in domain or functionally important loop regions. The NMSim simulations starting from the unbound structures are able to reach conformations similar to ligand bound conformations (RMSD 0.7) between the RMS fluctuations derived from NMSim generated structures and two experimental structures are observed. Furthermore, intrinsic fluctuations in NMSim simulation correlate with the region of loop conformational changes observed upon ligand binding in 2 out of 3 cases. The NMSim generated pathway of conformational change from the unbound structure to the ligand bound structure of adenylate kinase is validated by a comparison to experimental structures reflecting different states of the pathway as proposed by previous studies. Interestingly, the generated pathway confirms that the LID domain closure precedes the closing of the NMPbind domain, even if no target conformation is provided in NMSim. Hence, the results in this study show that, incorporating directional information in the geometry-based approach NMSim improves the sampling of biologically relevant conformational space and provides a computationally efficient alternative to state of the art MD simulations.Konformationsänderungen von Proteinen sind häufig eine grundlegende Voraussetzung für deren biologische Funktion. Die genaue Charakterisierung und Vorhersage dieser Konformationsänderungen ist für das Verständnis ihres Einflusses auf die Funktion erforderlich. Eines der dafür am häufigsten verwendeten und genauesten computergestützten Verfahren ist die Molekulardynamik-Simulationen (MD Simulationen). Diese sind jedoch nach wie vor sehr rechenintensiv und durchmustern den Konformationsraum nur in begrenztem Maße. Daher wurden Anstrengungen unternommen, alternative geometriebasierte Methoden (wie etwa CONCOORD oder FRODA) zu entwickeln, die auf einer reduzierten Darstellung von Proteinen beruhen. Das Ziel dieser Arbeit ist es, Richtungsinformationen in einen geometriebasierten Ansatz zu integrieren, und so den biologisch relevanten Konformationsraum erschöpfend zu durchmustern. Diese Idee führte kürzlich zur Entwicklung von „coarse-grained normal mode“ (CGNM) Methoden, wie zum Beispiel dem „elastic network model“ (ENM) und der von mir in vorangegangenen Arbeiten entwickelte „rigid cluster normal mode analysis“ (RCNMA). Beide Methoden liefern die gewünschte Richtungsinformation der intrinsischen Bewegungen eines Proteins in Form von harmonischen Moden (auch Normalmoden). Um die Aussagekraft, Robustheit und breite Anwendbarkeit solcher CGNM Verfahren zu untersuchen, wurde im Rahmen dieser Dissertation ein umfangreicher Vergleich zwischen „essential dynamics“ (ED) Moden aus MD Simulationen und Normalmoden aus CGNM Berechnungen durchgeführt. Der zugrundeliegende Datensatz enthielt 335 Proteine. Obwohl die CGNM Verfahren eine stark vereinfachte Darstellung für Proteine verwenden, korrelieren die niederfrequenten Moden dieser Verfahren bezüglich ihrer Bewegungs-Richtung (durchschnittliche maximale Überschneidung: 0,65) und -Amplitude (durchschnittliche maximale Überschneidung: 0,73) sehr gut mit ED Moden. Im Durchschnitt beschreibt das erste Viertel der Normalmoden 85 % des Raumes, der durch die ersten fünf ED Moden aufgespannt wird. Um die Leistungsfähigkeit von CGNM Verfahren genauer zu bestimmen, wurde im Rahmen der vorliegenden Studie eine dreistufige Methode zur Untersuchung der intrinsischen Dynamik von Proteinen entwickelt. Die ersten beiden Stufen basieren auf neusten Entwicklungen in der Rigiditäts-Theorie und der Beschreibung von elastischen Netzwerken. Diese sind im RCNMA Ansatz verwirklich und ermöglichen die Bestimmung der Normalmoden. Im letzten Schritt werden die Bewegungen des Proteinrückgrates entlang der mittels RCNMA erzeugten niederenergetischen Normalmoden ausgerichtet. Die Seitenkettenkonformrationen werden dabei durch Diffusionsbewegungen hin zu energetisch günstigen Rotameren erzeugt. Dies ist ein iterativer Prozess, bestehend aus mehreren kleineren Schritten, in denen jeweils intermediäre Konformationen erzeugt werden. Zur Validierung des NMSim Ansatzes wurde dieser mit den anderen zuvor genannten Simulationsmethoden am Beispiel von Lysozym verglichen. Die Fluktuationen der Aminosäurereste aus dem mit NMSim erzeugten Ensemble stimmen mit berechneten Fluktuationen aus der MD Simulation gut überein (Korrelationskoeffizient R = 0,79). Ein Vergleich der unterschiedlichen geometriebasierten Simulationsansätze zeigt, dass bei FRODA die Durchmusterung des Konformationsraumes des Proteinrückrates unzureichend ist. Bei CONCOORD ist hingegen die Durchmusterung des Konformationsraumes der Seitenketten unzureichend. NMSim hingegen durchmustert sowohl den Konformationsraum des Proteinrückrates als auch den der Seitenketten angemessen, wenn man die experimentell und mittels MD Simulationen erzeugten Konformationen als Referenz verwendet. Der NMSim Ansatz wurde ebenfalls auf einen Datensatz von Proteinen angewendet, für die Konformationsänderungen in Domänen oder in funktionell wichtigen Schleifenregionen experimentell beobacht wurden. In Übereinstimmung mit dem Konformations-Selektions-Modell ist der NMSim Ansatz bei vier von fünf Proteinen, die eine Domänenbewegung aufweisen, in der Lage, ausgehend von der ungebundenen Struktur neue Konformationen zu erzeugen, die der ligandgebundenen Konformation entsprechen (RMSD 0,7) zwischen der RMS Fluktuation der durch NMSim erzeugten Konformationen und jeweils zwei experimentellen Strukturen erreicht. Hingegen korrelieren die intrinischen Fluktuationen der NMSim Simulation in zwei von drei Fällen mit dem Bereich der ligandinduzierten Konformationsänderung in den Schleifen. Der mit NMSim generierte Pfad für die Konformationsänderungen von der ungebundenen Struktur zur ligandgebundenen Struktur der Adenylat-Kinase wurde durch den Vergleich zu experimentellen Strukturen validiert, die verschiedene Zustände des Pfades widerspiegeln. Die unterschiedlichen Kristallstrukturen, die entlang der Konformationsänderungen von der ungebundenen zur ligandgebundenen Struktur liegen, werden auf dem von NMSim erzeugten Pfad durchmustert. Interessanterweise bestätigt der generierte Pfad, dass die Schließbewegung der LID Domäne derjenigen der NMPbind Domäne vorangeht, sogar wenn keine Zielkonformation für die NMSim Simulation verwendet wurde

Computational studies on protein-ligand docking

This thesis describes the development and refinement of a number of techniques for molecular docking and ligand database screening, as well as the application of these techniques to predict the structures of several protein-ligand complexes and to discover novel ligands of an important receptor protein. Global energy optimisation by Monte-Carlo minimisation in internal co-ordinates was used to predict bound conformations of eight protein-ligand complexes. Experimental X-ray crystallography structures became available after the predictions were made. Comparison with the X-ray structures showed that the docking procedure placed 30 to 70% of the ligand molecule correctly within 1.5A from the native structure. The discrimination potential for identification of high-affinity ligands was derived and optimised using a large set of available protein-ligand complex structures. A fast boundary-element solvation electrostatic calculation algorithm was implemented to evaluate the solvation component of the discrimination potential. An accelerated docking procedure utilising pre-calculated grid potentials was developed and tested. For 23 receptors and 63 ligands extracted from X-ray structures, the docking and discrimination protocol was capable of correct identification of the majority of native receptor-ligand couples. 51 complexes with known structures were predicted. 35 predictions were within 3A from the native structure, giving correct overall positioning of the ligand, and 26 were within 2A, reproducing a detailed picture of the receptor-ligand interaction. Docking and ligand discrimination potential evaluation was applied to screen the database of more than 150000 commercially available compounds for binding to the fibroblast growth factor receptor tyrosine kinase, the protein implicated in several pathological cell growth aberrations. As expected, a number of compounds selected by the screening protocol turned out to be known inhibitors of the tyrosine kinases. 49 putative novel ligands identified by the screening protocol were experimentally tested and five compounds have shown inhibition of phosphorylation activity of the kinase. These compounds can be used as leads for further drug development

Rational Design of Small-Molecule Inhibitors of Protein-Protein Interactions: Application to the Oncogenic c-Myc/Max Interaction

Protein-protein interactions (PPIs) constitute an emerging class of targets for pharmaceutical intervention pursued by both industry and academia. Despite their fundamental role in many biological processes and diseases such as cancer, PPIs are still largely underrepresented in today's drug discovery. This dissertation describes novel computational approaches developed to facilitate the discovery/design of small-molecule inhibitors of PPIs, using the oncogenic c-Myc/Max interaction as a case study.First, we critically review current approaches and limitations to the discovery of small-molecule inhibitors of PPIs and we provide examples from the literature.Second, we examine the role of protein flexibility in molecular recognition and binding, and we review recent advances in the application of Elastic Network Models (ENMs) to modeling the global conformational changes of proteins observed upon ligand binding. The agreement between predicted soft modes of motions and structural changes experimentally observed upon ligand binding supports the view that ligand binding is facilitated, if not enabled, by the intrinsic (pre-existing) motions thermally accessible to the protein in the unliganded form.Third, we develop a new method for generating models of the bioactive conformations of molecules in the absence of protein structure, by identifying a set of conformations (from different molecules) that are most mutually similar in terms of both their shape and chemical features. We show how to solve the problem using an Integer Linear Programming formulation of the maximum-edge weight clique problem. In addition, we present the application of the method to known c-Myc/Max inhibitors.Fourth, we propose an innovative methodology for molecular mimicry design. We show how the structure of the c-Myc/Max complex was exploited to designing compounds that mimic the binding interactions that Max makes with the leucine zipper domain of c-Myc.In summary, the approaches described in this dissertation constitute important contributions to the fields of computational biology and computer-aided drug discovery, which combine biophysical insights and computational methods to expedite the discovery of novel inhibitors of PPIs

Simulating Molecular Mechanisms Of the Mdm2-Mediated Regulatory Interactions: a Conformational Selection Model Of the Mdm2 Lid Dynamics

Diversity and complexity of MDM2 mechanisms govern its principal function as the cellular antagonist of the p53 tumor suppressor. Structural and biophysical studies have demonstrated that MDM2 binding could be regulated by the dynamics of a pseudo-substrate lid motif. However, these experiments and subsequent computational studies have produced conflicting mechanistic models of MDM2 function and dynamics. We propose a unifying conformational selection model that can reconcile experimental findings and reveal a fundamental role of the lid as a dynamic regulator of MDM2-mediated binding. In this work, structure, dynamics and energetics of apo-MDM2 are studied as a function of posttranslational modifications and length of the lid. We found that the dynamic equilibrium between closed and semi-closed lid forms may be a fundamental characteristic of MDM2 regulatory interactions, which can be modulated by phosphorylation, phosphomimetic mutation as well as by the lid size. Our results revealed that these factors may regulate p53-MDM2 binding by fine-tuning the thermodynamic equilibrium between preexisting conformational states of apo-MDM2. In agreement with NMR studies, the effect of phosphorylation on MDM2 interactions was more pronounced with the truncated lid variant that favored the thermodynamically dominant closed form. The phosphomimetic mutation S17D may alter the lid dynamics by shifting the thermodynamic equilibrium towards the ensemble of semi-closed conformations. The dominant semi-closed lid form and weakened dependence on the phosphorylation seen in simulations with the complete lid can provide a rationale for binding of small p53-based mimetics and inhibitors without a direct competition with the lid dynamics. The results suggested that a conformational selection model of preexisting MDM2 states may provide a robust theoretical framework for understanding MDM2 dynamics. Probing biological functions and mechanisms of MDM2 regulation would require further integration of computational and experimental studies and may help to guide drug design of novel anti-cancer therapeutics