High-throughput Binding Affinity Calculations at Extreme Scales

Resistance to chemotherapy and molecularly targeted therapies is a major factor in limiting the effectiveness of cancer treatment. In many cases, resistance can be linked to genetic changes in target proteins, either pre-existing or evolutionarily selected during treatment. Key to overcoming this challenge is an understanding of the molecular determinants of drug binding. Using multi-stage pipelines of molecular simulations we can gain insights into the binding free energy and the residence time of a ligand, which can inform both stratified and personal treatment regimes and drug development. To support the scalable, adaptive and automated calculation of the binding free energy on high-performance computing resources, we introduce the High- throughput Binding Affinity Calculator (HTBAC). HTBAC uses a building block approach in order to attain both workflow flexibility and performance. We demonstrate close to perfect weak scaling to hundreds of concurrent multi-stage binding affinity calculation pipelines. This permits a rapid time-to-solution that is essentially invariant of the calculation protocol, size of candidate ligands and number of ensemble simulations. As such, HTBAC advances the state of the art of binding affinity calculations and protocols

Characterizing early drug resistance-related events using geometric ensembles from HIV protease dynamics:

The use of antiretrovirals (ARVs) has drastically improved the life quality and expectancy of HIV patients since their introduction in health care. Several millions are still afflicted worldwide by HIV and ARV resistance is a constant concern for both healthcare practitioners and patients, as while treatment options are finite, the virus constantly adapts via complex mutation patterns to select for resistant strains under the pressure of drug treatment. The HIV protease is a crucial enzyme for viral maturation and has been a game changing drug target since the first application. Due to similarities in protease inhibitor designs, drug cross-resistance is not uncommon across ARVs of the same class

PLAS-5k: Dataset of Protein-Ligand Affinities from Molecular Dynamics for Machine Learning Applications

Computational methods and recently modern machine learning methods have played a key role in structure-based drug design. Though several benchmarking datasets are available for machine learning applications in virtual screening, accurate prediction of binding affinity for a protein-ligand complex remains a major challenge. New datasets that allow for the development of models for predicting binding affinities better than the state-of-the-art scoring functions are important. For the first time, we have developed a dataset, PLAS-5k comprised of 5000 protein-ligand complexes chosen from PDB database. The dataset consists of binding affinities along with energy components like electrostatic, van der Waals, polar and non-polar solvation energy calculated from molecular dynamics simulations using MMPBSA (Molecular Mechanics Poisson-Boltzmann Surface Area) method. The calculated binding affinities outperformed docking scores and showed a good correlation with the available experimental values. The availability of energy components may enable optimization of desired components during machine learning-based drug design. Further, OnionNet model has been retrained on PLAS-5k dataset and is provided as a baseline for the prediction of binding affinities

Molecular Basis for Drug Resistance in HIV-1 Protease

HIV-1 protease is one of the major antiviral targets in the treatment of patients infected with HIV-1. The nine FDA approved HIV-1 protease inhibitors were developed with extensive use of structure-based drug design, thus the atomic details of how the inhibitors bind are well characterized. From this structural understanding the molecular basis for drug resistance in HIV-1 protease can be elucidated. Selected mutations in response to therapy and diversity between clades in HIV-1 protease have altered the shape of the active site, potentially altered the dynamics and even altered the sequence of the cleavage sites in the Gag polyprotein. All of these interdependent changes act in synergy to confer drug resistance while simultaneously maintaining the fitness of the virus. New strategies, such as incorporation of the substrate envelope constraint to design robust inhibitors that incorporate details of HIV-1 protease’s function and decrease the probability of drug resistance, are necessary to continue to effectively target this key protein in HIV-1 life cycle

Application of machine learning, molecular modelling and structural data mining against antiretroviral drug resistance in HIV-1

Millions are affected with the Human Immunodeficiency Virus (HIV) world wide, even though the death toll is on the decline. Antiretrovirals (ARVs), more specifically protease inhibitors have shown tremendous success since their introduction into therapy since the mid 1990’s by slowing down progression to the Acquired Immune Deficiency Syndrome (AIDS). However, Drug Resistance Mutations (DRMs) are constantly selected for due to viral adaptation, making drugs less effective over time. The current challenge is to manage the infection optimally with a limited set of drugs, with differing associated levels of toxicities in the face of a virus that (1) exists as a quasispecies, (2) may transmit acquired DRMs to drug-naive individuals and (3) that can manifest class-wide resistance due to similarities in design. The presence of latent reservoirs, unawareness of infection status, education and various socio-economic factors make the problem even more complex. Adequate timing and choice of drug prescription together with treatment adherence are very important as drug toxicities, drug failure and sub-optimal treatment regimens leave room for further development of drug resistance. While CD4 cell count and the determination of viral load from patients in resource-limited settings are very helpful to track how well a patient’s immune system is able to keep the virus in check, they can be lengthy in determining whether an ARV is effective. Phenosense assay kits answer this problem using viruses engineered to contain the patient sequences and evaluating their growth in the presence of different ARVs, but this can be expensive and too involved for routine checks. As a cheaper and faster alternative, genotypic assays provide similar information from HIV pol sequences obtained from blood samples, inferring ARV efficacy on the basis of drug resistance mutation patterns. However, these are inherently complex and the various methods of in silico prediction, such as Geno2pheno, REGA and Stanford HIVdb do not always agree in every case, even though this gap decreases as the list of resistance mutations is updated. A major gap in HIV treatment is that the information used for predicting drug resistance is mainly computed from data containing an overwhelming majority of B subtype HIV, when these only comprise about 12% of the worldwide HIV infections. In addition to growing evidence that drug resistance is subtype-related, it is intuitive to hypothesize that as subtyping is a phylogenetic classification, the more divergent a subtype is from the strains used in training prediction models, the less their resistance profiles would correlate. For the aforementioned reasons, we used a multi-faceted approach to attack the virus in multiple ways. This research aimed to (1) improve resistance prediction methods by focusing solely on the available subtype, (2) mine structural information pertaining to resistance in order to find any exploitable weak points and increase knowledge of the mechanistic processes of drug resistance in HIV protease. Finally, (3) we screen for protease inhibitors amongst a database of natural compounds [the South African natural compound database (SANCDB)] to find molecules or molecular properties usable to come up with improved inhibition against the drug target. In this work, structural information was mined using the Anisotropic Network Model, Dynamics Cross-Correlation, Perturbation Response Scanning, residue contact network analysis and the radius of gyration. These methods failed to give any resistance-associated patterns in terms of natural movement, internal correlated motions, residue perturbation response, relational behaviour and global compaction respectively. Applications of drug docking, homology-modelling and energy minimization for generating features suitable for machine-learning were not very promising, and rather suggest that the value of binding energies by themselves from Vina may not be very reliable quantitatively. All these failures lead to a refinement that resulted in a highly sensitive statistically-guided network construction and analysis, which leads to key findings in the early dynamics associated with resistance across all PI drugs. The latter experiment unravelled a conserved lateral expansion motion occurring at the flap elbows, and an associated contraction that drives the base of the dimerization domain towards the catalytic site’s floor in the case of drug resistance. Interestingly, we found that despite the conserved movement, bond angles were degenerate. Alongside, 16 Artificial Neural Network models were optimised for HIV proteases and reverse transcriptase inhibitors, with performances on par with Stanford HIVdb. Finally, we prioritised 9 compounds with potential protease inhibitory activity using virtual screening and molecular dynamics (MD) to additionally suggest a promising modification to one of the compounds. This yielded another molecule inhibiting equally well both opened and closed receptor target conformations, whereby each of the compounds had been selected against an array of multi-drug-resistant receptor variants. While a main hurdle was a lack of non-B subtype data, our findings, especially from the statistically-guided network analysis, may extrapolate to a certain extent to them as the level of conservation was very high within subtype B, despite all the present variations. This network construction method lays down a sensitive approach for analysing a pair of alternate phenotypes for which complex patterns prevail, given a sufficient number of experimental units. During the course of research a weighted contact mapping tool was developed to compare renin-angiotensinogen variants and packaged as part of the MD-TASK tool suite. Finally the functionality, compatibility and performance of the MODE-TASK tool were evaluated and confirmed for both Python2.7.x and Python3.x, for the analysis of normals modes from single protein structures and essential modes from MD trajectories. These techniques and tools collectively add onto the conventional means of MD analysis

Bioinformatics Techniques for Studying Drug Resistance In HIV and Staphylococcus Aureus

The worldwide HIV/AIDS pandemic has been partly controlled and treated by antivirals targeting HIV protease, integrase and reverse transcriptase, however, drug resistance has become a serious problem. HIV-1 drug resistance to protease inhibitors evolves by mutations in the PR gene. The resistance mutations can alter protease catalytic activity, inhibitor binding, and stability. Different machine learning algorithms (restricted boltzmann machines, clustering, etc.) have been shown to be effective machine learning tools for classification of genomic and resistance data. Application of restricted boltzmann machine produced highly accurate and robust classification of HIV protease resistance. They can also be used to compare resistance profiles of different protease inhibitors. HIV drug resistance has also been studied by enzyme kinetics and X-ray crystallography. Triple mutant HIV-1 protease with resistance mutations V32I, I47V and V82I has been used as a model for the active site of HIV-2 protease. The effects of four investigational antiviral inhibitors was measured for Triple mutant. The tested compounds had significantly worse inhibition of triple mutant with Ki values of 17-40 nM compared to 2-10 pM for wild type protease. The crystal structure of triple mutant in complex with GRL01111 was solved and showed few changes in protease interactions with inhibitor. These new inhibitors are not expected to be effective for HIV-2 protease or HIV-1 protease with changes V32I, I47V and V82I. Methicillin-resistant Staphylococcus aureus (MRSA) is an opportunistic pathogen that causes hospital and community-acquired infections. Antibiotic resistance occurs because of newly acquired low-affinity penicillin-binding protein (PBP2a). Transcriptome analysis was performed to determine how MuM (mutated PBP2 gene) responds to spermine and how Mu50 (wild type) responds to spermine and spermine–β-lactam synergy. Exogenous spermine and oxacillin were found to alter some significant gene expression patterns with major biochemical pathways (iron, sigB regulon) in MRSA with mutant PBP2 protein

Molecular dynamics simulations of complex systems including HIV-1 protease

Advances in supercomputer architectures have resulted in a situation where many scienti�fic codes are used on systems whose performance characteristics di�ffer considerably from the platform they were developed and optimised for. This is particularly apparent in the realm of Grid computing, where new technologies such as MPIg allow researchers to connect geographically disparate resources together into virtual parallel machines. Finding ways to exploit these new resources efficiently is necessary both to extract the maximum bene�fit from them, and to provide the enticing possibility of enabling new science. In this thesis, an existing general purpose molecular dynamics code (LAMMPS) is extended to allow it to perform more efficiently in a geographically distributed Grid environment showing considerable performance gains as a result. The technique of replica exchange molecular dynamics is discussed along with its applicability to the Grid model and its bene�fits with respect to increasing sampling of configurational space. The dynamics of two sub-structures of the HIV-1 protease (known as the flaps) are investigated using replica exchange molecular dynamics in LAMMPS showing considerable movement that would have been difficult to investigate by traditional methods. To complement this, a study was carried out investigating the use of computational tools to calculate binding affinity between HIV-1 protease mutants and the drug lopinavir in comparison with results derived experimentally by other research groups. The results demonstrate some promise for computational methods in helping to determine the most eff�ective course of treatment for patients in the future

Computational Studies of the HIV-1 Protease Dimer Interface.

HIV-1 protease (HIVp) is one of four major drug targets to prevent propagation of the infectious HIV virion. Currently, all ten marketed HIVp drugs are inhibitors that target the HIVp active site. However, these drug therapies provide selective pressure resulting in mutations of the protease that escape drug efficacy. Consequently, the development of inhibitors of HIVp that have new modes of action is necessary. The dimer interface is an attractive target due to its highly conserved nature and its importance in forming an active enzyme. Until now, all dissociative inhibitors were created by mimicking residues at the dimer interface, and resulted in several non-drug-like compounds. However, we created several receptor-based pharmacophore models of the HIVp dimer interface, using ensembles of multiple protein structures (MPS). The MPS method was used to map the dimer interface with a series of small-molecule probes – methanol, ethane, and benzene. The maps were translated into pharmacophore models which were used to filter in silico, three-dimensional library of small molecules. The MPS method identified several novel small-molecule inhibitors capable of inhibiting dimerization, with several compounds characterized with less than 50 uM-level affinity. In the clinically relevant multi-drug resistant form of HIVp, these compounds maintained dissociative inhibition with nearly identical inhibition rates. Zhang-Poorman kinetic analysis verified the small molecules inhibit HIVp in a dissociative manner. In addition to creating novel inhibitors, we modeled the protein-ligand interaction of known dissociative inhibitors using Langevin Dynamics. Ten, 10-ns simulations were initiated based on the hypothetical mechanism of ligand binding, but the dynamics simulations showed that the complex was unstable. Although the simulations did not result in a clear mechanism for protein-ligand binding, of the known dimerization inhibitors, they did demonstrate the entropic penalty of the proposed binding mechanism is unfavorable. Finally, we propose to use hydrogen/deuterium exchange (HDX) - mass spectrometry techniques to obtain new structural information and further characterize the molecular recognition between HIVp and dimer inhibitors. HDX can provide the first structural evidence defining the mechanism of HIVp dissociative inhibition by small molecules. HDX could be broadly applicable for a range of active-site and allosteric inhibitors.Ph.D.BiophysicsUniversity of Michigan, Horace H. Rackham School of Graduate Studieshttp://deepblue.lib.umich.edu/bitstream/2027.42/89663/1/jeromeq_1.pd

Binding free energy calculations and molecular dynamics simulations on complexes of viral proteases with their ligands

Ein Ziel der biomolekularen Modellierung ist die Berechnung der Affinität deltaG von Liganden an Proteine, insbesondere Enzyme. Das Spektrum der Methoden, die zu diesem Zweck entwickelt wurden, reicht von theoretisch genauen aber aufwändigen Verfahren zu einfachen, eher qualitativen Verfahren. Während letztere häufig empirische Scoring-Funktionen und eine einzelne Struktur als Eingabe verwenden, wird für kompliziertere Methoden der möglichst vollständige Konformationsraum eines Protein-Ligand-Komplexes benötigt. Dieser wird mit Sampling-Verfahren wie der Molekulardynamik (MD) durchmustert. In dieser Promotionsarbeit sollten Verfahren zur Berechnung von deltaG, insbesondere Varianten der Molecular Mechanics Poisson-Boltzmann Surface Area (MMPBSA) Methode, getestet und nach Möglichkeit weiterentwickelt werden. Desweiteren sollte die Auswirkung bestimmter Resistenzmutationen auf Struktur und Dynamik von Proteinen mit unterschiedlichen Maßen aus MD Simulationen heraus erfasst werden. Der erste Schritt der quantitativen Modellierung mit MD ist die Beschreibung der Moleküle durch die Parametrisierung eines Kraftfelds. Anhand des sulfatierten Tyrosins wurde eine solche molekulare Parametrisierung für ein Nicht-Standard-Molekül durchgeführt. Sodann wurden Varianten der tendenziell weniger aufwändigen MMPBSA-Methode getestet im Hinblick auf ihre Konvergenz und ihre Eignung zur Bestimmung genauer deltaG-Werte oder zumindest verschiedene Enzym-Ligand-Komplexe in eine richtige Rangfolge gemäß ihrer deltaG-Werte zu bringen. Die Varianten unterscheiden sich durch verschiedene Solvatisierungsmodelle und Methoden zur Berechnung der Entropie. Als molekulares Referenzsystem wurden Mutanten der HIV Protease im Komplex mit Wirkstoffen verwendet, da es hierzu experimentelle Daten gibt, mit denen die berechneten Werte verglichen werden können. Am anderen Ende des methodischen Spektrums liegt die aufwändige Thermodynamische Integration (TI). Bei einer guten Kraftfeldparametrisierung sollte TI in der Lage sein, deltaG-Effekte in der Größenordnung weniger kJ/mol quantitativ zu bestimmen. Dies wurde anhand der Mutante L76V der HIVProtease, die für einige Wirkstoffe zu einer Resensitivierung (erhöhte Affinität) führt, getestet. Schließlich sollten MD-Simulationen verwendet werden, um die molekularen Effekte von Mutationen der NS3/4A-Protease des humanen Hepatitis C Virus auf die Bindung von Liganden (Substrat, Inhibitoren) zu verstehen.A major aim of biomolecular modelling is the calculation of binding affinities deltaG of ligands to proteins, especially enzymes. The spectrum of methods that has been developed for this task ranges from theoretically exact but expensive to more simple and qualitative ones. While the latter are often empirical scoring functions using one single structure as an input, the more complex methods require the preferably complete conformational space of a protein-ligand complex which can be sampled using methods such as molecular dynamics (MD). The intention of this thesis was to test and further develop methods for the calculation of deltaG, in particular variants of the molecular mechanics Poisson-Boltzmann surface area (MMPBSA) method. Furthermore, the effects of specific resistance mutations on the structure and dynamics of proteins should be determined using different metrics on MD simulation data. The first step to quantitative modelling using MD is the description of the molecules by parameterizing a forcefield. Such a molecular parameterization was performed for the non-standard amino acid sulpho-tyrosine. Subsequently, variants of the less expensive MMPBSA method were tested with regard to their ability to converge and determine deltaG estimates or at least establish the correct ranking of deltaG values for a set of enzyme-ligand complexes. Different solvation models and procedures to calculate the entropy have been used. As a molecular reference system, mutants of the HIV protease complexed with inhibitors were used. For these systems, experimental data are available to which the calculated values can be compared. At the other end of the methodological spectrum is the more expensive thermodynamic integration (TI). With a proper forcefield parameterization, TI should be able to quantitatively determine deltaG effects in the order of a few kJ/mol. This was tested on the HIV protease mutation L76V which is known to lead to a resensitivation (increased affinity) for some drugs. Eventually, MD simulations were used to understand the molecular effects of mutations of the NS3/4A protease, an enzyme of the human hepatitis C virus, on the binding of ligands (substrate, inhibitors)

Uncertainty quantification in classical molecular dynamics

Molecular dynamics simulation is now a widespread approach for understanding complex systems on the atomistic scale. It finds applications from physics and chemistry to engineering, life and medical science. In the last decade, the approach has begun to advance from being a computer-based means of rationalizing experimental observations to producing apparently credible predictions for a number of real-world applications within industrial sectors such as advanced materials and drug discovery. However, key aspects concerning the reproducibility of the method have not kept pace with the speed of its uptake in the scientific community. Here, we present a discussion of uncertainty quantification for molecular dynamics simulation designed to endow the method with better error estimates that will enable it to be used to report actionable results. The approach adopted is a standard one in the field of uncertainty quantification, namely using ensemble methods, in which a sufficiently large number of replicas are run concurrently, from which reliable statistics can be extracted. Indeed, because molecular dynamics is intrinsically chaotic, the need to use ensemble methods is fundamental and holds regardless of the duration of the simulations performed. We discuss the approach and illustrate it in a range of applications from materials science to ligand-protein binding free energy estimation. This article is part of the theme issue 'Reliability and reproducibility in computational science: implementing verification, validation and uncertainty quantification in silico'