Keynote Talk
FEP for Drug Design Including Computation of Absolute Free Energies of Binding
William L. Jorgensen | Department of Chemistry, Yale University, New Haven, USABoth Monte Carlo statistical mechanics and molecular dynamics are being used to compute relative and absolute free energies of binding for protein-ligand complexes. Applications have included computing absolute free energies of binding for L99A lysozyme with benzene and seven analogs, and a drug-like molecule MIF180 with human macrophage migration inhibitory factor (MIF). It is noted that FEP calculations for molecular creations are much more efficient than for annihilations, and strikingly as few as 10 charge and Lennard-Jones windows are needed for creation of MIF180, which has 22 non-hydrogen atoms. Additional computational issues are considered including performance of OPLS-AA/M, CHARMM, and AMBER force fields, atomic charge models, water equilibration, and convergence. Significant advantages in computation speed are reported with MC over MD for similar extents of configurational sampling.
References
[1] Computer-Aided Discovery of Anti-HIV Agents. Jorgensen, W. L. Bioorg. Med. Chem. 2016, 24, 4768-4788.
[2] Enhanced Monte Carlo Methods for Modeling Proteins Including Computation of Absolute Free Energies of Binding. Cabeza de Vaca, I.; Qian, Y.; Vilseck, J. Z.; Tirado-Rives, J.; Jorgensen, W. L. J. Chem. Theory Comput. 2018, 14, 3279-3288.
[3] Robust FEP Protocols for Creating Molecules in Solution. Cabeza de Vaca, I.; Zarzuela, R.; Tirado-Rives, J.; Jorgensen, W. L. J. Chem. Theory Comput. 2019, 14, 0000-0000.
Invited Talks
Adaptive free energy calculations and free-energy guided design
John D. Chodera | Memorial Sloan Kettering Cancer Center, New York City, USAA major challenge in scaling up alchemical free energy methods for jointly optimizing affinity, selectivity, and other physical properties is how we can optimally utilize finite computational resources to most rapidly address the questions of interest. Here, we consider algorithms that will enable us to scale up to make adaptive, optimal use of distributed computing resources, utilizing both relative and absolute free energy calculations.
Improving the accuracy of binding free energy calculations via enhanced binding mode sampling and better force fields
David L. Mobley | University of California, Irvine, USAPredicting binding affinities of new compounds is challenging for a number of reasons. One challenge we encounter particularly often is the difficulty of knowing the binding mode in advance. While new compounds in a congeneric series often may share a binding mode, this is not always the case. Additionally, we seek computational tools which can work across compound series, or even across protein targets (e.g. to predict selectivity), so we need to be able to predict binding modes in addition to computing binding free energies of specific poses. I give an update on our work applying nonequilibrium candidate Monte Carlo (NCMC) to the ligand binding mode sampling problem in the context of our BLUES package, and discuss how this can impact binding free energy calculations. Additionally, I will briefly highlight some of our work on improving small molecule force fields in connection with the Open Force Field Initiative.
Using Long-Timescale Molecular Dynamics Simulations to Benchmark Enhanced Sampling Methods
Albert Pan | D.E. Shaw Research, New York City, USAAll-atom molecular dynamics (MD) simulation is a valuable technique for providing detailed information about the dynamics of biomolecules, but its computational expense can often be prohibitive. This limitation has motivated the development of “enhanced sampling” simulation methods - purely algorithmic changes to conventional MD that aim to accelerate the sampling of configurational states. Although many such methods are available, relatively few systematic studies of their performance have been done. Quantitative claims about the performance of enhanced sampling simulations are typically limited to (i) comparisons with conventional MD simulations of small, model systems, which may present qualitatively different sampling challenges than complex biological systems, or (ii) comparisons with experimental data, which may be complicated by discrepancies between the actual experimental conditions and the modeling of those experimental conditions in the enhanced sampling simulation, including errors in the physical model, or force field, used in simulation. An effective alternative approach to quantifying performance of enhanced sampling methods would be to compare the results they obtain on complex biological systems directly to those obtained by conventional MD simulations using the same force field. Here, we use long-timescale conventional MD simulations to assess the performance of certain commonly used enhanced sampling methods in accelerating the sampling of protein conformational changes, protein folding, and ligand binding.
Reproducible Free Energy Workflows with BioSimSpace
Antonia Mey | University of Edinburgh, UKUsing biomolecular simulations at an atomistic level to estimate free energies of binding between small ligands and a protein have become a popular approach in a computer-aided drug discovery pipelines over the past few years. While the general reliability of affinity predictions is improving, choices made during simulation setup and analysis can still have a strong influence on the results. These choices can be as simple as what molecular simulation engine to use, e.g. Gromacs, Amber, SOMD or more complex choices such as deciding on forcefield parameters for the ligands of interest. Another aspect is the choice of perturbation paths taken in the case of relative free energy calculations using alchemical free energy perturbations as an underlying simulation method. These choices make sharing automated workflows difficult and rely heavily on an expert computational chemist setting up simulations. BioSimSpace was designed with the goal of increasing shareability and reproducibility of biomolecular simulation workflows. BioSimSpace is an interoperable software framework building on existing simulation (e.g. Gromacs, Amber, SOMD) and analysis tools (mdtraj, mdanalysis) providing a high-level Python API to generate interoperable workflow components. The workflow components may be assembled into workflows, either by creating command line pipelines or using a workflow tool such as Knime. New workflows can be built up from a library of workflow components (or nodes) provided with BioSimSpace, or written from scratch in Python. Nodes are executed with best practice recommendations defined by node developers, though these may be overridden according to user preference. I will demonstrate how BioSimSpace has been used to prototype an automated alchemical free energy workflow on a D3R grand challenge dataset. In particular, I will highlight how the workflow is not only readily shareable, but also adjustable to be re-run with a different set of ligand parameters, perturbation map, or using a different underlying simulation engine. In this way, I will share insights gained from comparing different simulation strategies for the same data set.
References
[1] biosimspace.org (2018) Hedges, Mey, Michel, Woods.
[2] Abraham, M. et al. (2015) SoftwareX, 1, 19-25
[3] R. Salomon-Ferrer, D.A. Case, R.C. Walker. (2013) "An overview of the Amber biomolecular simulation package." WIREs Comput. Mol. Sci. 3, 198-210
[4] Sire molecular simulations framework Woods, C. et al. (2016)
[5] McGibbon, R. T. et al. (2015) Biophysical Journal 109, 1528-1532
[6] Gowers, R. J. et al. (2016) Proceedings of the 15th Python in Science Conference, 98-105
Performance of Alchemical Relative Free Energy Calculations With and Without Replica-Exchange
Peter V. Coveney | University College London, UKWe apply a systematic protocol for the evaluation of free energy calculations with and without replica-exchange. The method is based on ensemble averaging to generate accurate assessments of the uncertainties in the predictions. We thus compare FEP+ with TIES, the latter with and without so-called “enhanced sampling” based on replica-exchange protocols which is intrinsic to FEP+. Standard TIES performs best for a reference set of targets and compounds—previously used for comparison purposes between TIES and Amber pmemdGTI—wherein no benefits are forthcoming from replica-exchange methods. TIES incorporating the enhanced sampling technique known as REST2 (replica exchange with solute tempering), TIES-with-REST, captures large conformational changes and generates correct free energy differences, which are not observed using standard TIES. Evaluation of FEP+ and TIES-with-REST reveals a significant systematic underestimation of free energy differences in FEP+ compared to TIES, in turn raising a number of questions pertaining to the accuracy of predictions with the REST2 technique not hitherto widely discussed. Finally, a further source of error, intrinsic to the chaotic nature of molecular dynamics trajectories, will be briefly discussed.
In silico Lead Optimization via Rapid Exploration of Synthetically Tractable Chemical Space with Combined Free Energy and Machine Learning Approaches
Thomas Steinbrecher | Schrödinger GmbH, Mannheim, GermanyImproving or maintaining the potency of lead compounds, while simultaneously optimizing multiple other properties required for safety and biological efficacy, is a primary objective of lead optimization in small molecule drug discovery. A series of computational techniques have recently become available to aid in these tasks. We will discuss how predictive free energy calculations, synthesis-aware compound enumeration and next-generation machine learning techniques can rapidly accelerate the identification, under realistic conditions, of chemical matter that fulfills appropriate property constraints. Specifically, we will recount how the Pathfinder enumeration workflow, FEP+ binding free energy calculations and DeepChem machine learning models have been combined to rapidly generate and profile a library of novel compounds, based on commercially available building blocks, to explore both R-group and core-hopping modifications of an initial CDK2 inhibitor. The entire workflow provides a library of high-affinity, novel CDK2 inhibitors at a high enough throughput to positively impact pharmaceutical industry project timelines. Surprisingly, only a small fraction of the suggested cores had existing matches in BindingDB, despite featuring typical kinase hinge binding features and being achievable by straightforward synthesis. This indicates that exhaustive mining of easily accessible chemical space is typically not done in drug discovery today, even in a highly explored target class like kinases. Overall, it appears that the combination of the complementary computational techniques of binding free energy calculations, machine learning and compound enumeration offers huge synergistic effects. We will present evidence that a true de novo computational design solution benefiting from these advances is now feasible.
The Role of Water in Mediating Biomolecular Binding: From Water Locations to Their Impact on Binding Affinity
Jonathan W. Essex | University of Southampton, UKWater plays an intimate role in protein-ligand binding, not only through solvation/desolvation effects, but more subtly through the formation of direct interactions between the protein and ligand in the binding site. The targeting of bound water molecules for displacement as part of ligand optimization is a long invoked paradigm based around the release of configurational entropy, but there are many examples where displacing water leads to a loss in ligand binding affinity. Quantitatively accurate approaches to address this problem are arguable inadequate – water displacement and ligand interactions are intimately related and difficult to disentangle both experimentally and, hitherto, computationally. We have a long-standing interest in developing and using Grand Canonical Monte Carlo (GCMC) simulation approaches to explore water binding in protein-ligand systems. Through GCMC we are able to locate water molecules with good accuracy when compared against crystal structures. More significantly, the simulations clearly demonstrate the important role of water cooperativity; the mutual stabilization of water molecules means that individual water molecules cannot always be considered in isolation, but rather as part of a network. GCMC allows water binding sites and network binding free energies to be simultaneously calculated. In addition, by combining GCMC with alchemical perturbations of the ligand, networks of bound water molecules are able to adapt and maintain equilibrium with bulk water as the perturbation proceeds. Furthermore, the ability to extract active-site hydration free energies allows the deconvolution of protein-ligand binding free energies into separate protein- and water-mediated components, thereby providing rich, additional detail to the structure-activity relationship (SAR). In this presentation, our underlying methodology GCMC methodology will be described, together with examples of its application to water placement, binding free energy calculations, and protein-ligand affinity prediction.
Free Energy Calculation Using a Reference State: Replica-Exchange Enveloping Distribution Sampling (RE-EDS)
Gerhard König | ETH Zürich, SwitzerlandFree energy calculations require sufficient phase space overlap between the end states in order to converge. So far, this is achieved by dividing the free energy calculation into multiple substeps by introducing intermediate states along a user-defined “reaction coordinate” lambda. The potential energy function of such states is usually generated by linearly mixing the potential energies of the end states, or using soft-core potentials. Multiple improvements have been proposed for this scheme, but a particularly promising approach is the use of replica-exchange enveloping distribution sampling (RE-EDS). RE-EDS removes the need for conventional intermediate states (and associated problems such as the van der Waals endpoint problem) by simulating a single reference state that is able to sample the relevant configurations of all end states at the same time. Instead of investing computer time to create phase-space overlap, transitions between end states are achieved in RE-EDS by coupling reference states with different degree of smoothing of energy barriers in a replica-exchange scheme, thus enhancing sampling in an efficient manner. The choice of smoothed reference states for replica exchange can be performed in a highly automatized fashion. We will illustrate this with several examples, starting from simple model systems and proceeding to biomedical applications.
References
[1] C. Christ, W. F. van Gunsteren, J. Chem. Theory Comput. 5, 276 (2009).
[2] D. Sidler, A. Schwaninger, S. Riniker, J. Chem. Phys. 145, 154114 (2016).
[3] D. Sidler, M. Cristòfol-Clough, S. Riniker, J. Chem. Theory Comput. 13, 3020 (2017)
Applications and Advances in the One-Step Perturbation Approach
Chris Oostenbrink | BOKU, Vienna, AustriaPhysically relevant free-energy differences can be calculated from simulations of unphysical reference states. Over the years, we have successfully applied the one-step perturbation method by the judicious design of reference states. These could involve soft atoms with a mixed character [1], or molecules with chiral centers that could dynamically adjust their chemical configuration [2]. The diversity of end states in terms of their molecular dipoles or charge state was limited, due to a fixed polarization of the environment due to the reference state. A combination with a linear response or a third-power fitting approach alleviates this problem, still with a limited number of additional simulations [3]. Here, we apply this approach to the saturation mutagenesis of a small peptide and tackle the problem of conformational sampling of the sidechains by using appropriate rotamer libraries [4]. An alternative elegant solution to ensure overlap between the reference state and the end states is obtained by using enveloping distribution sampling (EDS) [5]. Here, the reference state is constructed from the end states directly, potentially maximizing the overlap in sampling [6]. We have recently proposed an alternative formulation of the EDS Hamiltonian, which maintains energy minima and lowers energy barriers by accelerating the EDS potential [7]. The accelerated EDS approach was applied to various model systems and realistic examples.
References
[1] A.B. Nørholm, P. Francotte, E. Goffin, I. Botez, L. Danober, P. Lestage, B. Pirotte, J.S. Kastrup, L. Olsen, and C. Oostenbrink. Thermodynamic characterization of new positive allosteric modulators binding to the glutamate receptor A2 ligand-binding domain: combining experimental and computational methods unravels differences in driving forces. J. Chem. Inf. Model 54 (2014) 3404−3416.
[2] M. M. H. Graf, Z. Lin, U. Bren, D. Haltrich, W. F. van Gunsteren, and C. Oostenbrink. Pyranose Dehydrogenase Ligand Promiscuity: A Generalized Approach to Simulate Monosaccharide Solvation, Binding, and Product Formation. PLOS Comput. Biol. 10 (2014) e1003995.
[3] A. de Ruiter and C. Oostenbrink. Efficient and accurate free energy calculations on trypsin inhibitors. J. Chem. Theory Comput. 8 (2012) 3686–3695.
[4] Z. Jandova, D. Fast, M. Setz, M. Pechlaner, C. Oostenbrink. Saturation mutagenesis by efficient free-energy calculation. J. Chem. Theory Comput. 14 (2018) 894–904.
[5] C.D. Christ, and W.F. van Gunsteren, Enveloping distribution sampling: a method to calculate free energy differences from a single simulation. J. Chem. Phys. (2007) 184110.
[6] M. Maurer, N. Hansen and C. Oostenbrink. Comparison of Free-Energy Methods Using a Tripeptide-Water Model System. J. Comput. Chem. 39 (2018) 2226–2242.
[7] J.W. Perthold and C. Oostenbrink. Accelerated enveloping distribution sampling: enabling sampling of multiple end-states while preserving local energy minima. J. Phys. Chem. B, 122 (2018) 5030–5037.
Predicting Antimicrobial Resistance: The Role of Computational Modelling in Translating Genetics into Clinical Microbiology
Philip W. Fowler | University of Oxford, UKThe growth of antimicrobial resistance (AMR) is widely recognised as a threat to much of modern medicine. To determine the suitable antibiotic to treat a serious infection, standard clinical practice is to culture a sample and then examine how well the bacterium grows in the presence of a number of different antibiotics. In the last five years, there has been a global push to replace culture-based clinical microbiology with whole genome sequencing where a sample of the pathogen is directly sequenced and its genome examined for mutations known to confer resistance to specific antibiotics. Whilst laudable, one key weakness of this approach is it is inferential so cannot return a result when a rare or novel genetic variant is encountered. In this talk I will describe the progress we have made in deploying computational methods, mainly alchemical free energy methods [1] but also machine learning [2], to predict de novo the effect of genetic variants on the action of specific antibiotics on S. aureus and M. tuberculosis.
References
[1] Fowler PW, Cole K, Gordon NC, Kearns AM, Llewelyn MJ, et al. Robust Prediction of Resistance to Trimethoprim in Staphylococcus aureus. Cell. Chem. Biol. 2018, 25:339.
[2] Carter JJ, Walker TM, Walker AS, Whitfield MG, Peto TEA, Crook DW, Fowler PW. Prediction of pyrazinamide resistance in Mycobacterium tuberculosis using structure-based machine learning approaches. bioRxiv 2019, DOI: 10.1101/518142.
ABFE, Water and Chemists: A Useful Combination?
Philip C. Biggin | University of Oxford, UKWe have previously shown how absolute binding free energies can be accurately computed both retrospectively and prospectively for compounds that target bromodomain proteins. Since then we have been interested to learn just how useful, this approach would be in an ongoing medicinal chemistry program. Bromodomains are acetylated-lysine ‘reader’ domains and are often found within larger macromolecules that are involved in gene transcription and degradation. We are interested in TRIM33, which is an E3-ubiquitin ligase, containing a C-terminal bromomdomain (BRD) and known to play a role in the repair of single-strand breaks (SSBs) in DNA. TRIM33’s implication in certain cancers makes it an interesting therapeutic target, an aspect that can be further explored by developing chemical probes for the BRD. Several compounds that bound to the BRD of TRIM33 were identified through a high throughput screen, one of which has been selected for optimisation. As there is no crystal structure of the TRIM-ligand complex, computational methods were employed to enable structure-based drug optimisation. Molecular dynamics and absolute binding free energies were utilized to increase confidence in ranking binding modes. Following our previous work, the energies of binding-site water molecules were probed with Grand Canonical Monte Carlo (GCMC) calculations. Combining all aspects has led to an improved binding model that explains all the available SAR and is now being used to yield rational designs for future compounds. In this presentation, I will discuss some of the challenges in taking this forward and where some of the challenges still lie.
From Seven Million to Four
Clara Christ | Bayer AG, Berlin, GermanyFree energy calculations are now routinely used in drug discovery to focus synthesis efforts onto the most promising molecules. In this talk we will show how large chemical spaces of virtual molecules can be filtered down with a cascade of methods using large scale binding affinity estimation as the final selection criterion in the cascade. We will point out strengths and weaknesses of the approach and give an outlook on future developments.
Free Energy Perturbation Calculations: Going Beyond SAR Analysis?
Zara Sands | UCB BioPharma, Braine-l'Alleud, BelgiumIn this talk we will describe instances where we have applied FEP technology to live pharmaceutical projects. We will flag some of the challenges and successes we have encountered, as well as areas we believe need to be addressed if FEP+ is to be used routinely for prospective drug design.
Prediction of Activity Cliffs Using FEP+ and Gromacs FEP
Herman van Vlijmen | Janssen Pharmaceutica, Antwerp, BelgiumActivity cliffs (ACs) are an important type of structure−activity relationship in medicinal chemistry where small structural changes result in unexpectedly large differences in biological activity. They are often difficult to predict with ligand-based approaches, as analogous effects are often not present in an available training data set. If no major structural changes in the protein occur between the binding of compounds that have an activity cliff, free energy based methods should in theory be able to predict the change in binding affinity. In this study we collected 148 compound pairs across 18 different target proteins from public and internal programs, with clear activity cliffs between the compounds. We performed FEP calculations with FEP+ and Gromacs, and showed that in the vast majority of cases the direction of affinity change was correctly predicted, with errors of about 1.4 kcal/mol (~1 log unit of affinity). Results from the FEP+ and Gromacs calculations were strongly correlated, with equal performance on one internal Janssen data set, and better performance of FEP+ on another one.
Fragment Binding Pose Predictions Using Unbiased Simulations and Markov-State-Models
Daniel Seeliger | Boehringer Ingelheim, Biberach, GermanyPredicting the co-structure of small molecule ligands and their respective target proteins has been a long-standing problem in drug discovery. For weak binding compounds typically identified in fragment-based screening (FBS) campaigns, determination of the correct binding site and correct binding mode is usually done experimentally via X-ray crystallography. For many targets of pharmaceutical interest, however, establishing an X-ray system which allows for sufficient throughput to support a drug discovery project is not possible. In this case exploration of fragment hits becomes a very laborious and consequently slow process with the generation of protein/ligand co-crystal structures as the bottleneck of the entire process. In this work we introduce a computational method which is able to reliably predict binding sites and binding modes of fragment-like small molecules using solely the structure of the apo-protein and the ligand’s chemical structure as input information. The method is based on molecular dynamics simulations and Markov-state models and can be run as a fully automated protocol requiring minimal human intervention. We describe the application of the method to a representative subset of different target classes and fragments from historical fragment-based screening efforts at Boehringer Ingelheim and discuss its potential integration into the overall fragment-based drug discovery workflow.
Contributed Talks
Probing different aspects of stability for the Pin1-WW domain
Niels Hansen | University of Stuttgart, GermanyWe present an in-silico mutagenesis study for the Pin1-WW domain based on alchemical free-energy molecular dynamics simulations. The class of WW domains has been studied intensively over the years, both in experiment and simulation, mainly due to their complex thermodynamic character and their small size. The considered data set of nearly 100 point mutations comprises side-chain [1] and amide-to-ester backbone mutations [2]. We show that subtle local variations in the protein starting structure may lead to substantial deviations in the calculated free-energy changes. For such cases, the usage of multiple starting structures in combination with Hamiltonian replica exchange was found to be essential to obtain converged free-energy estimates [3]. However, even then surprising deviations with some experimental data points occur despite an overall good agreement. Therefore, further investigation regarding the comparability between experiment and computation and the validity of the regularly made two-state assumption in the analysis of experimental data is reported.
References
[1] K. Dave, M. Jäger, H. Nguyen, J. W. Kelly, M. Gruebele, J. Mol. Biol. 428, 2016, 1617-1636.
[2] S. Deechongkit, P. Dawson, J. Kelly, J. Am. Chem. Soc. 126, 2004, 16762-16771.
[3] D. Markthaler, H. Kraus, N. Hansen, J. Chem. Inf. Model. 58, 2018, 2305-2318.
Integrated design of novel biologics
Philip Kim | University of Toronto, CanadaI will present our integrated technology platform to develop novel biologics as a use case for thermodynamic integration approaches. Briefly, we integrate computational library design with modern in-cell selection strategies to uncover novel therapeutics. Having obtained promising leads with this methodology, we optimize these using a variety of approaches, including thermodynamic integration. To obtain higher affinity compounds, we make use of non-canonical amino acids as well as macrocyclization techniques.
Free energy calculations in active drug discovery projects: use cases and challenges
Christina Schindler | Merck KGaA, Darmstadt, GermanyFree energy calculations have become a powerful addition to the computational chemist’s toolbox to support structure-based drug design in hit-to-lead and lead optimization stages of drug discovery projects. Here, we present a large prospective data set from using FEP+ in active drug discovery projects at Merck KGaA and compare this performance to results obtained on a new, challenging benchmark of pharmaceutically relevant targets. We discuss opportunities and challenges and highlight use cases and conditions that can maximize the impact of the method in projects.
Enhancing relative protein-ligand binding free energy calculations with grand canonical Monte Carlo simulations of water molecules
Gregory A. Ross | Schrödinger Inc., New York City, USABuried water molecules are frequently a sampling bottleneck in molecular dynamics simulations and protein-ligand free energy calculations. Grand canonical Monte Carlo (GCMC) can greatly enhance the sampling of water molecules, particularly those in buried pockets. Indeed, recent studies have reported very promising results when using GCMC in protein-ligand free energy calculations. We report the development of a GPU-accelerated implementation of GCMC within the FEP+ relative free energy calculation workflow. Our implementation speeds up the sampling efficiency of water by orders of magnitude at a negligible increase in the computational cost. With transformations that displace water in buried pockets, FEP+ predictions that use GCMC are consistent with calculations performed using the double decoupling method. When applied to a large dataset of protein-ligand binding, using FEP+ with GCMC can significantly improve accuracy and the rate of convergence, and removes the dependency of the prediction on the starting solvation structure. The improved sampling efficiency afforded by GCMC also highlights systems that require greater care in their set-up, as well as areas for force-field improvement.
QligFEP: an automated workflow for small molecule free energy calculations in Q
Willem Jespers | Uppsala University, SwedenThe process of ligand binding to a biological target can be represented as the equilibrium between the relevant solvated and bound states of the ligand, which is the basis of structure-based, rigorous methods such as the estimation of relative binding affinities by free energy perturbation (FEP). Despite the growing capacity of computing power and the development of more accurate force fields, a high throughput application of FEP is currently hampered due to the need, in the current schemes, of an expert user definition of the “alchemical” transformations between molecules in the series explored. Here, we present QligFEP, a solution to this problem using an automated workflow for FEP calculations based on a dual topology approach. In this scheme, the starting poses of each of the two ligands, for which the relative affinity is to be calculated, are explicitly present in the MD simulations associated with the (dual topology) FEP transformation, making the connection pathway univocal. We show that this generalized method can be applied to accurately estimate solvation free energies for amino acid sidechain mimics, as well as the binding affinity shifts due to the chemical changes typical of lead optimization processes. This is illustrated in a number of protein systems extracted from other FEP studies in the literature: inhibitors of CDK2 kinase and a series of A2A adenosine G-protein coupled receptor antagonists, where the results obtained with QligFEP are in excellent agreement with experimental data. In addition, our protocol allows for scaffold hopping perturbations to identify the binding affinities between different core scaffolds, which we illustrate with a series of Chk1 kinase inhibitors. QligFEP is implemented in the open-source MD package Q, and works with the most common family of force fields: OPLS, CHARMM and AMBER.
Using Free Energy Perturbation Calculations to Predict Relative Binding Affinities in Drug Design
Zoe Cournia | Biomedical Research Foundation, Academy of Athens, GreeceComputer-aided drug design has become an integral part of drug discovery and is nowadays extensively used in the lead identification and lead optimization phases. The process of delivering optimized leads of higher affinity than the parent compound and providing coherent structure–activity relationships that can efficiently guide synthetic efforts can now be both time and cost-efficient using reliable computational methods to calculate protein–ligand binding affinities. A rigorous method for calculating relative protein–ligand binding affinities is the free energy perturbation (FEP) framework, which is currently the most accurate qualitative link between experimental and computational studies [1], when using sufficiently small mutations and adequate sampling. In this talk, we present our methodology for predicting relative binding affinities of ligands bound to proteins. Several case studies are presented on optimizing M2TM inhibitors [2,3], Arp2/3 allosteric ligands [4], and FXR ligands within the D3R Grand Challenge 2 [5]. Based on our results, applying FEP calculations has a high success rate if appropriately used, thus saving valuable resources in a lead optimization project. Practical considerations and critical aspects of running relative binding free energy calculations are presented. Comparison of FEP results between different software programs is also performed. Although FEP/MD calculations are possible with some free MD engines such as NAMD, no automated tool has been developed to streamline the process, making the calculations tedious and unfeasible for a large number of molecules. Here, we present a web-based tool, FEPprepare, which automates the set-up procedure for performing NAMD/FEP simulations [6].
Acknowledgements
This research has been co‐financed by the European Union and Greek national funds through the Operational Program Competitiveness, Entrepreneurship and Innovation, under the call RESEARCH – CREATE – INNOVATE (project code:T1EDK-03532).
References
[1] Cournia Z*, Allen B, Sherman W, J. Chem. Inf. Model. 2017, 57(12), 2911-2937.
[2] Gkeka P, Eleftheratos S, Kolocouris A, Cournia Z*. J Chem Theor Comp. 2013, 9(2), 1272–1281.
[3] Ioannidis H, Drakopoulos A, Tzitzoglaki C, Homeyer N, Kolarov F, Gkeka P, Freudenberger K, Liolios C, Gauglitz G, Cournia Z*, Gohlke H*, Kolocouris A*. J Chem Inf Model. 2016, 56(5), 862-876.
[4] Baggett AW*, Cournia Z*, Han MS, Patargias G, Glass AC, Liu SY, Nolen BJ. ChemMedChem. 2012, 7(7):1286-94.
[5] Athanasiou C, Vasilakaki S, Dellis D, Cournia Z*. J Comput Aid Mol Des, 2018, 32(1):21-44.
[6] http://feprepare.vi-seem.eu.
The SAMPL6 SAMPLing challenge: Assessing the reliability and efficiency of binding free energy calculations
Andrea Rizzi | Weill Cornell University, New York City, USAWe present the conceptualization and results for the first SAMPLing challenge from the SAMPL series focusing on the assessment of convergence properties and reproducibility of binding free energy methodologies. We provided parameter files and multiple initial geometries for two octa-acid (OA) and one cucurbit[8]uril (CB8) host-guest systems, for which it is computationally feasible to obtain converged binding affinity estimates in a matter of hours or a few days. Participants submitted binding free energy predictions as a function of the computational effort for six different alchemical- and physical-pathway (e.g. molecular dynamics and potential of mean force) methodologies based on GROMACS, AMBER, and OpenMM implementations. For the two small OA binders, the free energy estimates computed with alchemical and potential of mean force approaches show relatively similar variance and bias as a function of the number of energy/force evaluations, with the attach-pull-release (APR) and GROMACS expanded ensemble methodologies performing particularly well. The differences between the methods widen when analyzing the CB8-quinine system, where both the guest size and correlation times are greater. For this system, coupled topologies non-equilibrium switching (CP-NS) obtained the overall highest efficiency followed by Hamiltonian replica exchange (HREX). Among the conclusions emerging from the data, we found that CP-NS convergence can be enhanced by increasing the length of the non-equilibrium protocol, that HREX, while displaying very small variance, can incur into substantial bias that depends on the initial population of the replicas, and that the Berendsen barostat introduces non-negligible artifacts in expanded ensemble simulations. Surprisingly, the results suggest that specifying the force field parameters and charges is insufficient to ensure reproducibility to better than ~0.5 kcal/mol. Further work will be required to identify the exact source of these discrepancies.
Large scale pmx/Gromacs based non-equilibrium free energy calculations
Vytautas Gapsys | Max Planck Institute for Biophysical Chemistry, Göttingen, GermanyIn the talk I will present a method to generating hybrid structures/topologies for amino acid, nucleic acid and ligand modifications. Automation of the process allows for a large scale thermodynamic scanning over various biomolecular complexes employing diverse simulation setups and force field variants. I will highlight several applications of the non-equilibrium free energy calculations in estimating relative free energy differences for protein thermostability, drug resistance, protein-DNA interactions and ligand binding affinities. Comparing the obtained results to experimental measurements benchmarks the accuracy achievable with the contemporary molecular mechanics force fields. An attractive opportunity to combine estimates from different force fields into a consensus result emerges, yielding further improvement in accuracy.
Posters
Developing a Robust Method for Automated Assessment of Binding Affinity via Free Energy Perturbation
Maximilian Kuhn[a,b], Paolo Tosco[a], Antonia Mey[b], Mark Mackey[a], Julien Michel[b] | Cresset[a] and University of Edinburgh[b], UKThere is a rising interest towards FEP (Free Energy Perturbation) calculations in the drug discovery community.[1,2] FEP calculations are performed to predict the relative binding affinity changes (ΔΔG) within a congeneric ligand series. Such calculations consist of non-physical (“alchemical”) transformations, in which a molecule (A) is gradually converted into a structurally related molecule (B) through a number of discrete steps, the so-called λ windows. The ligand simulated in each window can be thought of as an alchemical (i.e., hybrid) molecule consisting of a 1-λ fraction of A and a λ fraction of B. The free energy difference between the end states of the transformation can be assessed by a variety of methods, such as the Multistate Bennett Acceptance Ratio (MBAR) and corresponds to the binding affinity difference between the two molecules. Extending this approach to a network of congeneric molecules leads to assessing the binding affinities of the corresponding ligand series. Recently, alchemical free energy calculations have been applied to predict ligand binding affinities of large data sets, yielding accuracies as close as 0.8 to 1.1 kcal/mol compared to experimental values[2]. However, apart from intrinsic errors like inaccuracies of the force field as well as insufficient phase space sampling, the lack of automation is still a major obstacle to routine application of these methods. For instance, the initial generation of the perturbation network requires comparison of the ligand structures and subsequent connection of the molecules in the network graph based on their structural similarity. This step is especially problematic when the input structures are relatively heterogenous, as this requires the insertion or deletion of a large number of atoms during the simulations. Such perturbations are prone to errors, and usually need a larger number of λ windows to complete successfully, thus increasing calculation time. Therefore, to increase reliability of the results and avoid wasting computing time on poorly designed perturbations, early recognition and subsequent modification of problematic alchemical networks are necessary. A number of open-source applications are available to assist in various parts of the FEP workflow. These include various applications implemented in the Sire3 molecular simulation package for network setup and results evaluation, LOMAP[4] for network generation, AmberTools[5] for preparing the topology files and SOMD[6] for running the simulations. However, installing and using these relatively complex tools requires expertise. In order to make FEP more accessible, a fully automated workflow for performing FEP calculations is under development for inclusion in version 3.0 of Cresset’s drug design software Flare. The workflow can be accessed through a graphical user interface (GUI) or a Python application programming interface (API). The original open-source code base was enhanced to automatically create additional molecules serving as intermediates for ligand pairs which are too dissimilar to be directly mutated into one another. To increase the workflow’s reliability, frequent cycle closures are ensured during the generation of the perturbation network, and energy convergence checks during postprocessing are enforced. Furthermore, an effort is underway to deduce heuristics to infer the optimal number of λ windows for each perturbation, in order to avoid unnecessary calculations and minimise the probability of insufficient sampling. A variety of datasets were processed using predefined default parameters, including the FEP+ dataset[2]. Results obtained with our method were broadly comparable to published reports, considering the overall reduced simulation time. During our validation of the method we found that the ideal settings for a given set of ligands and their target protein is difficult to predict in advance. Therefore, depending on the project and computational resources available, further optimisation of the results may be desirable once knowledge on the system under observation has been gathered through preliminary simulations. Power users of this FEP implementation will have full control on simulation parameters via Flare’s GUI and Python API.
References
[1] Cournia, Z.; Allen, B.; Sherman, W. J. Chem. Inf. Model. 2017, 57 (12), 2911-2937.
[2] Wang, L.; Wu, Y.; Deng, Y. et al. J. Am. Chem. Soc. 2015, 137 (7), 2695-2703.
[3] Woods, C.; Mey, A.S.J.S.; Calabrò, G.; Bosisio, S.; Michel, J. Sire molecular simulations framework (2016). http://siremol.org. Accessed 28 January 2019.
[4] Liu, S.; Wu, Y.; Lin, T. et al. J. Comput.-Aided Mol. Des. 2013, 27 (9), 755-770
[5] Case, D.A.; Ben-Shalom, I.Y.; Brozell, S.R. et al. AMBER 2018, University of California, San Francisco.
[6] Mey, A.S.J.S.; Jiménez, J.J. Michel, J. J. Comput.-Aided Mol. Des. 2018, 32, 199-210.
Analysis of Thermodynamic Coupling via Alchemical Free Energy Calculations
Martin Werner | Max Planck Institute for Biophysical Chemistry, Göttingen, GermanyInteractions between combined amino acid mutations play a key role in protein engineering affecting properties such as cooperativity and allostery. The effect of a double mutation on properties like thermostability or binding affinity can be significantly different from the sum of effects of the separate single mutations. Such a situation is revealed by a nonadditivity of the corresponding free energies and indicates the mutations involved to be thermodynamically coupled. While this type of epistatic interaction is expected for residues in close spatial contact, additivity of the effects of both single amino acid mutations is widely observed for distant pairs. Indeed, even mutations on opposite sides of a protein can exhibit significant nonadditivity in their free energies indicating a correlation of both mutations persisting over large distances. Molecular dynamics simulations with alchemical amino acid mutations grant access to the free energies and nonadditivities of double mutants via different pathways. The comparison of the resulting thermodynamic couplings to experimental free energies affirm an accurate and robust approach. Thus, alchemical free energy calculations provide a powerful tool for the analysis of correlated amino acid mutations and their interaction over large distances. In fact, the access to accurate mutation free energies and nonadditivities potentially enables monitoring of the involved interaction pathway throughout the entire protein.
Unveiling the molecular mechanism underlying a bleeding disorder by using free energy calculations
Gloria Angelica Sandoval-Pérez | Universidad de los Andes, Bogotá, ColombiaThe von Willebrand Disease (vWD) is the most common hereditary bleeding disorder. It is caused by the decrease of plasma levels or defects of the extracellular protein von Willebrand Factor (vWF). These perturbations impair the normal hemostatic function of adhering platelets at sites of vascular injury of the vWF. In particular, naturally occurring mutations increase or decrease the binding of platelets, resulting in a bleeding condition known as von Willebrand Disease type 2. Understanding the molecular mechanism leading to this pathology, remains largely unknown and it is restricted to only few, out of many, identified mutations. Here we tackle this issue by a computational systematic study of 18 vWF mutants located at the A1 domain, the domain where the platelets bind. Molecular dynamics (MD) simulations are being carried out to determine alterations in structure, dynamics, and mechanical stress of the vWF A1 domain, imposed by this set of mutants, both when the domain is in isolation and in complex with the GPIb-alpha platelet receptor. MD-based free energy calculations are performed to elucidate the thermostability of the vWF A1-GPIb-alpha complex upon mutation. The MD simulations are complemented with experimental biophysical data thereby elucidate the molecular mechanism underlying the vWD-type 2. In particular, this project aims to reveal whether the mutants obey universal principles or instead each mutant distorts platelet binding following its own unique mechanism.
Calcium sensor proteins: Can we model ion binding free energies as a function of metastable Markov states?
Tim Hempel | Freie Universität Berlin, GermanyThe Synaptotagmin-1 (Syt-1) C2a domain is a calcium sensor in the neurotransmitter release machinery [1]. We investigate the interplay of conformational dynamics and binding free energies of the Syt-1 C2a domain. We aim at modeling the calcium binding free energies as a function of the protein metastable state. On the one hand, this is achieved by identifying metastable states and modeling the dynamics between them with hidden Markov models [2, 3]. On the other hand, Alchemical free energy perturbation is used to assess calcium ion binding free energies [4] in each of those metastable states. We plan on computing the equilibrium binding free energies by Boltzmann-weighting the energies obtained for each metastable state. Our analysis is concluded by thoroughly discussing reproducibility, quantify deviations and identify possible sources of errors. Our findings can help to interpret Syt-1's function as a molecular switch, a contribution to elucidating synaptic exocitosis mechanisms on the atomic scale.
References
[1] Südhof. C. Neuron 80, no. 3 (2013)
[2] Prinz, Wu, Sarich, Keller, Senne, Held, Chodera, Schütte, Noé. J. Chem. Phys. 134, no. 17 (2011).
[3] Noé, Wu, Prinz, Plattner. J. Chem. Phys. 139, no. 18 (2013).
[4] Hansen and Gunsteren. J. Chem. Theory Comput. 10, no. 7 (2014).
Repulsive soft-core potentials for efficient alchemical free energy calculations
Yaozong Li | University of Zurich, SwitzerlandIn alchemical free energy (FE) simulations, the separation-shifted soft-core potential is typically used to prevent the end-point instability, which is caused by physical overlap of a transformed solute with surrounding solvents. While this soft-core potential eliminates the end-point issue, the resulting hybrid Hamiltonian becomes non-linear with respect to the parameter λ that interpolates linearly between the two end-state Hamiltonians. This non-linearity complicates the FE estimation by the Bennett’s acceptance ratio (BAR), free energy perturbation (FEP) and thermodynamic integration (TI) methods, and ultimately lowers their efficiency. In this work, we propose a new soft-core potential, called Gaussian soft-core (GSC) potential, which retains the linearity of the hybrid Hamiltonian with respect to λ, and thereby allows straightforward determination of FE values using BAR, FEP, TI and their variants. In addition, the developed GSC potential can be combined with a strategy that estimates the TI integrand (∂H/∂λ) based on FEP formulae to substantially decrease the number of λ simulations in comparison with the commonly used alchemical simulation approaches. Finally, the efficiency and accuracy of the GSC potential are demonstrated by comparing the free energies of annihilation of 13 small molecules and alchemical mutation of a protein side chain.
Towards using absolute binding free energy calculations in the early stages of ligand development
Joe Bluck | University of Oxford, UKAbsolute binding free energy (ABFE) calculations have been shown to accurately rank ligands for a given protein, as well as correctly determine the lowest energy binding poses from docking studies.[1,2] These predictions are extremely applicable to the early stages of a ligand optimisation project, where a lack of co-crystal structures and poor experimental data make structure-based optimisation difficult. Weak binders, often the starting point of a ligand optimisation project, present several challenges for performing ABFE calculations. A lack of ligand-protein shape complementarity and ligand flexibility means a greater degree of sampling is needed for calculations to converge. However, given the recent advances in GPU acceleration in MD packages it is now feasible to perform the longer simulation lengths required with more repeats. The work presented in this poster looks at how we are trying to refine an ABFE calculation protocol in parallel with an active ligand optimisation project. Particularly focusing on the challenges encountered with flexible groups and their impact on convergence.
References
[1] Aldeghi, M., Bodkin, M. J., Knapp, S. & Biggin, P. C. A statistical analysis on the performance of MMPBSA versus absolute binding free energy calculations: bromodomains as a case study. J. Chem. Inf. Model. 57, 2203−2221 (2017).
[2] Aldeghi, M., Heifetz, A., Bodkin, M. J., Knapp, S. & Biggin, P. C. Predictions of ligand selectivity from absolute binding free energy calculations. J. Am. Chem. Soc. 139, 946−957 (2016).
Mutation Free Energies of Human Interferon-Gamma Analogues
Elena Lilkova | Institute of Information and Communication Technologies at the Bulgarian Academy, Sofia, BulgariaHuman interferon-gamma (hIFNg) is a crucial immunomodulating cytokine, which biological effects may range from proliferation to apoptosis. Its aberrant expression is associated with multiple autoimmune diseases. This makes it a target of rational protein design aiming at developing biologically inactive analogues that have preserved affinity for the hIFNg cell-surface receptor. We will report preliminary results on calculating mutation free energies of two hIFNg analogues -- one with a K88Q point mutation and a double K88Q,T27Y mutant. To this end we employ slow growth free energy perturbation calculations using the PMX single topology approach. The calculated mutation free energies will be used to estimate via a thermodynamic cycle the differences in the binding free energies of the native and mutant forms of hIFNg towards its extracellular receptor hIFNgR1. The results of our modeling approach will be compared to thermodynamic parameters, that were obtained experimentally using isothermal titration calorimetry. The goal of our study is to provide an accurate and robust simulation protocol for the in silico evaluation of potential candidates for biologicals that can be used to compete with the native endogenous hIFNg in order to inhibit its activity.
Acknowledgements
This work was supported in part under Grant DN-11/20/2017 of the Bulgarian Science Fund. Computational resources were provided by the HPC Cluster at the Faculty of Physics of Sofia University "St. Kliment Ohridski".
Insights into Chemistry through the Computation of Free Energy Hot-Spots
Johannes Dietschreit | LMU Munich, GermanyPredicting free energies is one of the key challenges in modern quantum chemistry. For small, unimolecular systems this is usually done via a frequency analysis of the molecule using quantum mechanics (QM) calculations. For large systems, one focuses on free energy differences computed from sampled energies applying, e.g., Bennett’s acceptance ratio method. The interpretation of free energy results is in most cases not straightforward, because it is not possible to distinguish contributions from different atoms or residues and, therefore, to understand the underlying effects (e.g., bond weakening, sterical clashes, new non-covalent interactions) causing the free energy to change. Based on an approach originally introduced by Berens et al. (1983) to estimate quantum corrections to thermodynamic properties, we present a method that calculates the vibrational part of the free energy from the vibrational density of states function and locates changes in the potential energy surface. Those are the so-called hot-spots, which we interpreted as the locations causing the overall free energy change.
Estimation of absolute and relative binding free energies with PELE
Joan Francesc Gilabert | Barcelona Supercomputing Center, SpainThe estimation of protein-ligand binding free energies is one of the most important problems in computational biophysics. Despite its importance for fields such as in silico drug discovery, it remains an open problem mainly due to two factors: i) the large amount of sampling needed and ii) the limited accuracy of the models used. The difficulty in obtaining the estimates have forced researchers to adopt several different strategies, which can be classified into absolute or relative approaches. While the former estimate the binding free energy of each system independently, the latter give as a result differences in binding free energy between different systems. These different methods have complementary strengths and weaknesses, which is why in our group we are exploring the use of our Monte Carlo code, Protein Energy Landscape Exploration (PELE), for both relative and absolute strategies. We obtain absolute binding free energies by combining enhanced sampling, standard PELE simulations and Markov State Models. On the other hand, relative binding free energies of different ligands are estimated using Free Energy Perturbation methods. We apply local sampling to compute the energetic difference when following an unphysical pathway to transform one ligand into another smoothly. Through the combination of the different approaches to binding free energy estimation with PELE, we hope to produce a versatile and innovative technique that will help accelerate the drug discovery process.
How Phosphorylation Affects the Binding of C-terminal Peptides to PDZ Domains
Nicolas Künzel | University of Saarland, Saarbrücken, GermanyThis project aims to understand how phosphorylation affects protein-peptide interactions. For this molecular dynamics (MD) simulations and alchemical simulations are used to study the binding of the EQVSAV peptide (C-terminal of RAPGEF6) to the PDZ2 domain of the human PTP1E (PTPN13, FAP-1, PTPL1) protein as well as the PDZ1 domain of the human MAGI1 protein (AIP-3, BAP-1). The results are compared to experiment in order to mechanistically explain the binding preference for the unphosphorylated form of this peptide.
14-3-3 Protein Peptide-binding Pathways from Distance Field Replica Exchange
Gabor Nagy | Max Planck Institute for Biophysical Chemistry, Göttingen, GermanyThe distance field approach provides a grid-based collective coordinate that allows reversible binding to a designated binding site, while avoiding clashes with the binding partner. Combining distance field with Hamiltonian replica exchange molecular dynamics allowed us to identify the most likely peptide binding pathway in the 14-3-3Z protein, and determine free energy profiles and conformational changes of the system along the pathway.
Skin permeability modeling using the Martini force field
Christian Wennberg | ERCO Pharma AB, Stockholm, SwedenAccurate skin permeation modeling requires a good understanding of the partitioning of the topically applied components between their vehicle and the stratum corneum barrier structure. One powerful tool when studying this partitioning is free energy calculations in molecular dynamics (MD). However, due to the possible synergistic effects between multiple components in topical formulations during their partitioning, all-atom MD simulations might be too expensive to accurately model this event. An alternative solution is to use coarse-grained MD simulations in order to perform a preliminary screening of possible partition- and permeability-enhancing effects of various topical formulations. Here we present preliminary results with the new beta-release of the Martini 3.0 forcefield for the permeability and partitioning of a selected group of permeability enhancers into our stratum corneum barrier model.
Funnel metadynamics on the Kelch domain of Keap1
Cecilia Chavez Garcia | The University of Western Ontario, London, CanadaFunnel metadynamics (FM) is a method that allows a ligand to enhance the sampling of the target binding sites and its solvated states. In this method, a funnel-shaped potential is applied to the system, reducing the space to explore in the unbound state. During the simulation, the system visits the bound and unbound states several times, allowing an accurate estimation of the binding free-energy surface within a reasonable simulation time. In this work, we performed funnel metadynamics on a DLGex peptide and the Kelch domain of Keap1. The peptide contains the "DLG" motif necessary in Neh2 for Kelch binding. The Neh2 domain is an important transcription factor responsible for the cell's defense against oxidative stress. The initial structure for the peptide was taken from a 1 μs trajectory with the Amber99SB*-ILDNP force field, and the Kelch domain was retrieved from the Protein Data Bank.
Estimation of free energy surfaces for substrate translocation – Metadynamics versus umbrella sampling
Vinaya Kumar Golla, Manas Joshi, Jigneshkumar Prajapati, and Ulrich Kleinekathöfer | Jacobs University Bremen, GermanyThe permeation of antibiotics across Gram-negative bacterial outer membranes is essential to reach the site of action of the antibacterial drug. The advancement of molecular dynamics methods enables us to understand these translocation processes firmly at an atomic level. To this extent, the calculation of free energy surfaces and the underlying antibiotic translocation mechanisms along the outer membrane channels are highly significant. As examples, the translocations of the antibiotics fosfomycin, ciprofloxacin and enrofloxacin have been computed recently [1-3]. Herein, we investigated a variety of substrates and their free energy surfaces along the Pseudomonas aeruginosa outer membrane channel OprO using the well advanced metadynamics and umbrella sampling free energy methods. The free energy calculations have been performed to illustrate the difference of computed free energies by virtue of increasing the substrate complexity along with the involved asymmetry of the outer membrane channels during the substrate translocation process. Thus, the resulting comparative analysis helps us to choose the most appropriate method for future calculations of similar case studies.
References
[1] V. K. Golla et al., Biophys. J. 116, 258–269 (2019)
[2] J. D. Prajapati et al., JPCB 122, 1417 (2018)
[3] J. D. Prajapati et al., JCTC 13, 4553 (2017)
Combining Force Field-based Hybrid Schemes for Grand-canonical Simulations of Low-resolution Models of GPCRs
Ksenia Korshunova | FZ Jülich, GermanyAll-atom molecular dynamics models for predicting poses and binding affinities of ligands for GPCRs are detailed but computationally expensive. Coarse-grained molecular dynamics is faster and reflects the dynamics of the chains behavior well, but doesn't account for interactions on atomistic level, which are important for the ligand binding processes. As an alternative, Molecular Mechanics/Coarse-Grained (MM/CG) method combined with Hamiltonian Adaptive Resolution Scheme (H-AdResS), which includes both all-atom and coarse-grained resolutions for different parts of the model, was adopted. In this method, only the ligand-binding site is treated on atomistic level, whereas the side chains far from the active site are coarse-grained. Areas of the protein with different resolutions are mechanically coupled. Droplet of atomistic water surrounding extracellular part of the protein is coupled to a reservoir of coarse-grained water so that energy and particles can be exchanged between the regions with different resolution. This ensures the correctness of the grand-canonical ensemble in the all-atom region and thus allows free energy calculation of the ligand binding. This poster presents the current state of the project concerning developing a statistically correct and efficient description for the atomistic region and lastly implementing free energy calculations in the MM/CG/H-AdResS setup.
A free energy convergence study for proteins and molecules in solutions through Monte Carlo methods
Israel Cabeza de Vaca Lopez | Yale University, New Haven, USAStatistical mechanics approaches are the most accurate methods to determine free energies. The accuracy of these predictions relies on two main parts: energy estimation and proper sampling. In the general, molecular mechanics force fields are able to provide good free energy estimations with an accurate configurational sampling. A common technique to estimate free energies with molecular dynamics (MD) or Monte Carlo (MC) is the free energy perturbation method (FEP). In general, FEP and other similar methods such as thermodynamic integration (TI) perform an alchemical transformation from the initial to the final state in stages (stratification or multiple lambda windows) to increase the convergence of the calculations. The main sampling problem is to determine the minimum amount of lambda windows to ensure enough configurational space overlap to converge the free energy independently of the transformation path. A common practice to determine the accuracy of the free energy estimations is to perform the transformation using the same path in both ways. A free energy difference between both directions indicates insufficient sampling or lack of configurational space overlap between the initial and final state. In this poster, we present a free energy convergence analysis between both directions respect to the number of lambda windows and sampling configurations using FEP and MC methods for solvents. We will analyze the hydration free energy of a test set composed by small solutes ranging for argon to larger molecules with rotatable bonds using different creation and annihilation strategies. As is a common practice, decoupling of the perturbation in electrostatic and Lennard-jones interactions are performed showing interesting behavior differences. Remarkably, in all cases, the molecular creation direction is more efficient than the annihilation direction. In particular, the Lennard-jones creation direction requires more sampling to converge than annihilation direction but is less dependent on the number of lambda windows. On the contrary, annihilation direction requires always more lambda windows and the results may converge to wrong values in many cases. Moreover, we will present a MC method applied to determine absolute binding free energies for protein-ligand systems.
Identifying determinants of RNA affinity and specificity in nuclear proteins regulating organelle gene expression
Charles Robert | CNRS Laboratoire de Biochimie Théorique, Paris, FrancePreliminary MD studies aimed at validating existing features and identifying new determinants of affinity and specificity of ssRNA binding in alpha-solenoid proteins.
Lipid-Protein Interactions in Potassium Channel Permeation and Gating
Ruo-Xu Gu | Max Planck Institute for Biophysical Chemistry, Göttingen, GermanyPotassium channels are widely distributed, highly selective ion channels gated by diverse physiological signals. There has been great interest in studies of their structure and function, not only due to their fundamental role in electric activity of neurons and muscle contraction but also because they are clinical targets for treating cardiovascular disease and neurological dysfunction. Investigating how bilayer properties affect ion conductance of potassium channels has biological relevance, because the spatially and temporally heterogeneous lipid distribution (e.g., the putative lipid raft) in cell membranes may regulate protein structure and function. In fact, a variety of experiments have suggested effects of lipids on potassium channel conductance, although the mechanism remains unclear. We use molecular dynamics simulations to simulate ion current across potassium channels in the presence of an electric field, which mimic electrophysiological experiments. Specifically, we are interested in the current through channels as a function of membrane thickness and charges. Lipids with different tails and headgroups are used to modify these membrane properties, and MthK are used as examples of potassium channels. We found lower ion conductance in bilayers with smaller bilayer thickness, due to bending of transmembrane helices which in turn affect the opening of the selectivity filter. In addition, we also found that membrane surface charges accelerate ion permeation in the “opening state” of the selectivity filter.
A metadynamics protocol to predict the correct ligand binding pose
Lucia Fusani | GlaxoSmithKline and University of Strathclyde, UKThe prediction of the correct protein-ligand binding pose or poses is important in structure based drug design and crucial for the evaluation of protein-ligand binding affinity. Most three-dimensional ligand-protein structures are obtained from X-ray crystallography experiments which result in a single static model of an ensemble of conformations. In cases where crystal structures cannot be obtained, computational tools have proven to be a powerful alternative, but they are still not accurate enough to replace experimental techniques. The Binding Pose Metadynamics (BPMD)[1] tool allows the study of ligand stability in a full atomistic details in a computational efficient manner averaging over 10 × 10 ns metadynamics runs with the root-mean square deviation of the ligand heavy atoms as the collective variable. The basic principle of BPMD is that ligand poses which are unstable (average RMSD > 2 Å) under the bias of the metadynamics simulation are likely to be infrequently occupied in the energy landscape and make minimal contributions to the protein-ligand binding affinity. Here the robustness of the method is studied using test cases selected from the Twilight database[2] and by rescoring Induced-Fit Docking (IFD) poses. Results show that bad poses are clearly discriminated from the good ones. In particular, in our investigations of crystal structures, those compounds whose binding pose in not supported by the electron density are separated from those with well-defined electron density.
References
[1] A. J. Clark, P. Tiwary, K. Borrelli, S. Feng, E. B. Miller, R. Abel, R. A. Friesner and B. J. Berne, J Chem Theory Comput, 2016, 12, 2990-2998.
[2] C. X. Weichenberger, E. Pozharski and B. Rupp, Acta Crystallographica Section F, 2013, 69, 195-200.
An automated, efficient, and scalable framework for the benchmarking of molecular force fields, and estimation of physical properties from molecular simulation
Simon Boothroyd | MSKCC, New York City, USADeveloping an accurate force field representation of molecules is key to realizing the full potential of molecular simulations, both as a powerful route to gaining fundamental insight into a broad spectrum of chemical and biological phenomena, and for predicting physicochemical and mechanical properties of substances. The Open Force Field Consortium is an industry-funded open science effort to this end, developing open source tools for rapidly generating new, high-quality small molecule force fields. An integral aspect of this is the parameterization and assessment of force fields against high-quality, condensed-phase physical property data (curated from open data sources such the NIST ThermoML Archive) alongside quantum chemical data. The quantity of data for pure and mixed systems available in the ThermoML Archive alone however would require an onerous amount of human and compute resources to estimate manually, especially when estimations must be made for numerous sets of force field parameters. Here we present an entirely automated, highly scalable framework for evaluating physical properties and their gradients in terms of force field parameters. It is written as a modular and extensible Python framework, which employs a multi-fidelity approach that allows for rapid estimation of properties from cached simulation data, and a pluggable API for estimation of new properties. We demonstrate the utility of the framework by benchmarking the SMIRNOFF99Frosst force field against a set of densities and static dielectric constants, utilizing both direct simulation and MBAR reweighting approaches. In future, the framework will be extended to the rapid estimation of more involved properties, including host-guest and protein-ligand binding affinities, and solvation and hydration free energies.
The Open Force Field Consortium - Open Force Fields in Industrial Pharmaceutical R&D
Daniel Kuhn | Merck, Darmstadt, GermanyThe Open Force Field Consortium is an academic-industry consortium which aims to develop open, modern and extensible biomolecular force fields. The development is based on an open source model where force fields, software and tools are made freely available. This approach should improve predictive design by more accurate force fields, along with establishing the necessary infrastructure to make these force fields easier to build and use. In this poster we present an overview of the Open Force Field Consortium and will show how to use the newly developed Open Force Field Toolkit and highlight use cases in an industrial R&D setting.
Alchemical Free Energy Calculations of Aminoadamantanes Bound to the Closed State of Influenza A/M2TM
Dimitrios Stamatis | University of Athens, GreeceAdamantane derivatives, such as amantadine and rimantadine, have been reported to block the transmembrane domain (TM) of the M2 protein of influenza A virus (A/M2) but their clinical use has been discontinued due to evolved resistance in humans. Although experiments and simulations have provided adequate information about the binding interaction of amantadine or rimantadine to the M2 protein, methods for predicting binding affinities of whole series of M2 inhibitors have so far been scarcely applied. Here we show that alchemical free energy calculations of ligand binding using the free energy perturbation (FEP) method and the thermodynamic integration (TI) are valuable for determining the relative binding affinity of aminoadamantane derivatives against M2TM WT in 1,2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC) lipids to mimic the membrane environment. The binding affinities of 20 compounds were measured by isothermal titration calorimetry (ITC) against the A/M2TM WT tetramer in its closed form at pH 8 and used as experimental probes covering a binding affinity range of only ~ 2 kcal mol-1. However, a fair correlation was found which demonstrates that binding free energy calculations can be used to predict relative binding affinities of aminoadamantane derivatives towards M2TM with good accuracy. Such methods could assist in the development of novel potent inhibitors that overcome A/M2 resistance.
Future plans for the Statistical Assessment of the Modeling of Proteins and Ligands project
Paul Czodrowski[a], Danielle Bergazin[b], David Mobley[b] | [a]TU Dortmund University, Germany; [b]University of California, Irvine, USAThe Statistical Assessment of the Modeling of Proteins and Ligands (SAMPL) project is a set of community-wide blind challenges aimed to advance computational techniques as standard predictive tools in rational drug design. Here, we present an overview of future plans for the SAMPL challenges, which will focus on physical property prediction, including logP and logD values, pKa prediction, host–guest binding, and other properties, as well as broadening to include a protein-ligand component. Previous SAMPL challenges have focused on a broad range of biologically relevant systems with different sizes and levels of complexities including proteins, host–guest complexes, and drug-like small molecules to test the latest modeling methods and force fields. Experimental results, such as binding affinities and hydration free energies, are withheld from participants until after the prediction submission deadline, so that the true predictive power of methods can be revealed. The most recent SAMPL6 challenge included challenges based on aqueous host-guest binding data (binding free energies and optionally binding enthalpies) for three different host molecules; and on pKa prediction, for a set of fragment-like molecules. Our current physical property challenge, SAMPL6 Part II, focuses on a small octanol-water log P prediction set.
Exploring the complex free-energy landscape of hGBP1 GTPase and its mutants utilizing ab initio QM/MM metadynamics simulations
Ravi Tripathi, Rachel Glaves, Jan Noetzel, Dominik Marx | Ruhr-Universität Bochum, GermanyMetadynamics has emerged as a powerful method for sampling free energy pathways in complex processes[1, 2]. Here I will present the versatility of metadynamics technique in exploring the unknown reaction pathways of GTP/GDP hydrolysis in a GTPase, human guanine nucleotide binding protein-1 (hGBP1). The significance of accelerated ab-initio QM/MM metadynam- ics simulations in scrutinizing the molecular level details of the GTP and GDP hydrolysis in hGBP1 will be discussed[3]. Furthermore, the computed minimum free energy path, among various possible path, will also be demonstrated. Last but not least I will show the efficacy of this approach in reveling multiple reaction pathways in mutant hGBP1 which feature very similar energy barriers and thus reaction rates[4]. The molecular level information obtained by revealing such complex reaction mechanisms can be utilized in rational drug design.
References
[1] A. Laio and M. Parrinello, Proc. Natl. Acad. Sci. U.S.A., 2002, 99, 12562–12566.
[2] M. Iannuzzi, A. Laio, and M. Parrinello, Phys. Rev. Lett., 2003, 90, 238302–1–4.
[3] R. Tripathi, R. Glaves, and D. Marx, Chem. Sci., 2017, 8, 371–380.
[4] R. Tripathi, J. Noetzel, and D. Marx, Phys. Chem. Chem. Phys., 2019, 21, 859–867.
Automated Free Energy Calculation for Drug Design: Accelerated Enveloping Distribution Sampling
Jan Walther Perthold | Institute of Molecular Modeling and Simulation, University of Natural Resources, Vienna, AustriaFree-energy perturbation (FEP) methods are commonly used in drug design to calculate relative binding free energies between different ligands. Alchemical ligand transformations are usually performed in multiple steps which need to be chosen carefully to ensure sufficient phase-space overlap between neighboring states. With one-step or single-step FEP techniques, a single reference state is designed that samples phase-space not only representative of a full transformation, but ideally resembles multiple ligand end-states and hence allows for efficient multi-state perturbations [1]. Enveloping distribution sampling (EDS) is one example for such a method in which the reference-state is created by mathematical combination of the different ligand end-states based on solid statistical mechanics [2]. We have recently proposed a novel approach to EDS which enables efficient barrier-crossing between the different end-states, termed accelerated EDS (A-EDS) [3]. In this work, we further simplify the parametrization of the A-EDS reference-state and demonstrate the automated calculation of multiple free-energy differences between different ligands from a single simulation in three different well-described drug design model systems.
References
[1] Oostenbrink, C., Free energy calculations from one-step perturbations. Methods Mol Biol 2012, 819, 487-99.
[2] Christ, C. D.; van Gunsteren, W. F., Enveloping distribution sampling: a method to calculate free energy differences from a single simulation. J Chem Phys 2007, 126 (18), 184110.
[3] Perthold, J. W.; Oostenbrink, C., Accelerated Enveloping Distribution Sampling: Enabling Sampling of Multiple End States while Preserving Local Energy Minima. J Phys Chem B 2018, 122 (19), 5030-5037.
Determining Free Energy Differences Through Variational Morphing
Martin Reinhardt | Max Planck Institute for Biophysical Chemistry, Göttingen, GermanyThe Bennett Acceptance Ratio (BAR) method is widely used to estimate the free energy difference between two given states of a many-body system, such as a macromolecule, where sampling in these states is commonly conducted using Molecular Dynamics (MD) or Monte Carlo simulations. As insufficient phase space overlap of such states can lead to large biases in the results, the calculation is commonly split up into several smaller steps using intermediate states defined along a so-called morphing path. To date, such paths are primarily set as linear interpolations of the start and the end states' Hamiltonians. Here, using BAR, we develop an optimization scheme to find the path with the highest accuracy from a more general class of paths that also includes non-linear interpolations. For a broad range of conceptual one-dimensional systems, we show that the free energy sampling error can be decreased by 25 % to 35 % percent compared to linear interpolations, with a theoretical upper limit of 50 %. In addition, the solution can be altered to account for correlations arising from using the same sample in one state to evaluate differences to more than one other state, as commonly done. Furthermore, we show that the optimization is still valid for a high number of intermediate states with very low sampling for each. We conclude by discussing how the optimized pathways can be applied to MD simulations.
Relaxation-Augmented Free Energy Perturbation
Ying-Chih Chiang | University of Southampton, UKWe re-examine the underlying theory of the widely-employed free energy perturbation (FEP) method, which is also known as the Zwanzig equation. FEP frequently fails to converge, and this poor performance reflects the problem of insufficient sampling inherent in the theory. Using a harmonic oscillator model, we prove analytically that the deviation between the exact solution and the falsely-converged FEP result is related to a relaxation energy that the system would release, if it were allowed to relax after being perturbed. Upon introducing this nonequilibrium relaxation energy as a correction term to the standard FEP, we obtain the exact free energy difference between two states, even for cases where the standard approach with Bennett’s acceptance ratio fails. Our new approach - the relaxation augmented free energy perturbation (RAFEP) - hence sheds a new light on the insufficient sampling problem in free energy calculations.
Relative Principal Components Analysis
Mazen Ahmad | Max-Planck Institute for Informatics, Saarbrücken, GermanyA new method termed “Relative Principal Components analysis” (RPCA) is introduced that extracts optimal relevant principal components to describe the change between two data samples representing two macroscopic states. The method is widely applicable in data-driven science. Calculating the components is based on a physical framework which introduces the objective function (the Kullback-Leibler divergence) appropriate for quantifying the change of the macroscopic state effected by the changes in the microscopic features. To demonstrate the applicability of RPCA, we analyze the thermodynamically relevant conformational changes of the protein HIV-1 protease upon binding to different drug molecules. In this case, the RPCA method provides a sound thermodynamic foundation for analyzing the binding process and thus characterizing both the collective and the locally relevant conformational changes. Moreover, the relevant collective conformational changes can be reconstructed from the informative latent variables to exhibit both the enhanced and the restricted conformational fluctuations upon ligand association.
Calculating permeability through skin
Magnus Lundborg | ERCO Pharma, Stockholm, SwedenThe main protective barrier of the human body is the lipid matrix in the stratum corneum in skin. This barrier prevents excessive water loss and also foreign substances from entering the body. In turn, this means that distribution of pharmaceuticals through skin is challenging. The transdermal delivery route is otherwise very appealing since it avoids, e.g., first-pass metabolism. We present calculations of the permeability coefficient, based on the potential of mean force and the local diffusion coefficient, through our cryo-EM validated model of the stratum corneum lipid barrier. The calculations are based on nonequilibrium forward-reverse pulling molecular dynamics simulations. We also use the free energy of transfer from the solvent into the lipid barrier system to calibrate the potential of mean force. The solvent can be water, to enable comparison to common in vitro/ex vivo experiments, or more complex formulations.
Proton-control of transitions in an amino-acid transporter
Zhiyi Wu | Univeristy of Oxford, UKAmino acid transport into the cell is often coupled to the proton gradient, as found in the solute carrier (SLC) 7 family of cationic amino acid transporters (CATs). Although no structure of a human CAT exists, the structure of a bacterial homolog, GkApcT, has been recently solved in an inward occluded state {Jungnickel, 2018 #8899} and allows an opportunity to examine how protons are coupled to amino acid transport. Our working hypothesis is that release of the amino acid substrate is facilitated by deprotonation of a key glutamate residue, located at the bottom of the binding pocket and which acts as the intracellular gate, allowing the protein to transition from an inward-occluded to an inward-open conformation. To explore this hypothesis, we have used molecular dynamics and umbrella sampling and demonstrate that the transition from inward-occluded to inward-open is more energetically favourable when this key glutamate is deprotonated. We also used absolute free energy calculations to show that the binding of amino acids in the binding pocket is more unfavourable in the inward-open state compared than in the inward-occluded state, consistent with the view that the amino acid would more likely be released to the intracellular solution from the inward-open state rather than from the inward-occluded state.
Net Charge Changes in the Calculation of Protein-Ligand Binding Free Energies via Classical Atomistic MD Simulation
Christoph Öhlknecht | University of Natural Resources and Life Sciences, Vienna, AustriaThe calculation of binding free energies of a charged ligand to a target protein is an often encountered problem in molecular dynamics (MD) studies. Absolute free energies as well as their relative changes (compared to a ligand with a different net charge) are both affected by artifacts caused by the representation of electrostatic interactions in a non-macroscopic simulation box using periodic boundary conditions.[1] In particular, electrostatic interactions are calculated either using a lattice-summation scheme or a cutoff-truncation scheme with Barker-Watts reaction-field correction rather than summing over individual potentials. The application of correction terms had already been proposed but so far only been applied for simple systems, including polyatomic molecules binding to a functionalized C60-buckyball [2] or the binding of a ligand to a cytochrome c peroxidase where the MD simulation was performed with fixed solute atom positions.[3] Here, we present two studies of complex biological model systems: the absolute binding free energy of acetylsalicylic acid coupling to a phospholipase A2 [4] and relative binding free energies of Netropsin and Distamycin binding to DNA [5] were succesfully derived by applying the correction terms to the alchemically computed free energy estimates. The corrections split up into individual terms for the direct ligand-protein interactions, for the P-summation over the individual solvent molecules and for the wrong solvent polarization around the charged ligand. For both systems, the corrected estimates were found to be in very good agreement with experiments.
References
[1] Kastenholz MA, Hünenberger PH (2006) J. Chem. Phys. 2006, 124:124106.
[2] Reif MM, Oostenbrink C (2013) J. Comput. Chem. 2014, 3: 227-243.
[3] Rocklin GJ, Mobley DL, Dill KA, Hünenberger PH (2013) J. Chem. Phys. 2014, 139:184103.
[4] de Ruiter A, Oostenbrink C (2013) J. Chem. Theory Comput. 2013, 9:883-892.
[5] Jožica Dolenc, Chris Oostenbrink, Jože Koller, Wilfred F. van Gunsteren (2005) Nucleic Acids Research 2005, 33(2): 725–733
Minimum free energy pathways of ciprofloxacin and enrofloxacin across a bacterial pore
Jigneshkumar Dahyabhai Prajapati | Jacobs University Bremen, GermanyThe outer membrane of the Gram-negative bacterium E. coli forms a potent physical barrier to antimicrobial agents. The porin OmpC in the outer membrane opens a potential passage for antibiotics to reach the periplasmic space of the cell [1]. However, there is a significant lack of knowledge regarding the permeation routes taken by different classes of antibiotic molecules through this pore. With the availability of metadynamics simulations technique and a string method, we developed a protocol to characterize the minimum free energy pathway and the rate-limiting molecular interactions underlying the transport of ciprofloxacin and enrofloxacin across the OmpC porin [1,2]. Moreover, we were able to reproduce these permeation paths in qualitative manner using our in-house implementation of a temperature accelerated Brownian dynamics simulation technique [3,4]. Overall, we have performed ~40 µs and ~600 µs of MD and BD simulations, respectively, to obtain an atomistic scale insight. In general, the procedure has the potential to build a structure-function relationship between different porins and antibiotics molecules and thereby to aid the development of future antibiotics.
References
[1] J. D. Prajapati et al., JCTC 13, 4553 (2017)
[2] J. D. Prajapati et al., JPCB 122, 1417 (2018)
[3] C. J. F. Solano et al., JCTC 12, 2401 (2016)
[4] C. J. F. Solano et al., JCTC 14, 6701 (2018)
Exploration Of The Activation Mechanism Of Small GTPase RhoA
Ennys Gheyouche | University of Nantes, FranceThe Ras superfamily is one of the largest protein families involved in cell signal transduction, with 166 members. This superfamily shares a common activation mechanism where a GTP-bound is involved in signal transduction and a GDP-bound ras member is found in the unactivated form. Once the GTP is dephosphorylated in GDP, ras members have to be reactivated through a nucleotide exchange mechanism catalyzed by a ras-specific Guanine Exchange Factor (GEF). In our study, we focus our attention on the detailed mechanism of exchange of one member of the Ras family, rhoA, and its specific partner p115-GEF (also known as Arhgef-1). These proteins play a role in the vascular smooth muscles cells in response to angiotensin-II stimuli. The protein-protein interaction involves two region switches on rhoA, the nucleotide binding pocket and a magnesium. In order to understand the activation mechanism at an atomic level, an in silico approach using molecular dynamics (MD) on all-atom simulations is used to investigate the dynamic of the interaction between Rhoa and Arhgef1 and when does the nucleotide exchange take part. Therefore this molecular approach will provide us with informations on the key elements involved in this switch process and decipher this challenging mechanism.
Machine Learning of Hydration Free Energies
Clemens Rauer | Max Planck Institute for Polymer Research, Mainz, GermanyHydration free energies are an important molecular property which can provide an insight into the thermodynamic state of the respective system. However, in order to build a free energy surface for a large molecular compound set, many high level simulations are necessary, making the whole process computationally very expensive. We overcome this problem by building a kernel ridge regression model on the calculated hydration free energies. Using this model, we are able to predict free energies for a large set of small druglike molecules[1]. Moreover, by applying a delta-machine learning[2] model, we can use a "cheap" low level method to predict free energies and learn the correction to a higher level method. Then, we can predict high level free energies for significantly larger compound sets than was used in the training of the original model. We show that by using only limited high level data, highly accurate free energies can be calculated using this approach.
References
[1] Mobley, D.L. & Gilson, M.K. Annu. Rev. Biophys. 2017, 46:531-58
[2] Ramakrishnan et al. J. Chem. Theory Comput. 2015, 11, 2087-2096
Optimization of drug-membrane selectivity from free-energy calculations
Bernadette Mohr | Max Planck Institute for Polymer Research, Mainz, GermanyCurrent virtual screening approaches reduce the computational cost by using simplified representations and approximations, statistical mechanical effects are excluded. We investigate one way to include the effects determining the selectivity of a molecule to a specific target into virtual screening by using physics-based models. Coarse-grained simulations are introduced as a preliminary step to allow the effective screening of a large number of molecules. As a test case, the physical and chemical properties of the fluorescent dye 10-N-nonyl acridine orange are modified to increase its binding affinity to Cardiolipin. In eukaryotes, Cardiolipin is almost exclusively found in the inner mitochondrial membrane. Mitochondria are the location of important metabolic pathways and are linked to diseases and apoptosis. Therefore, targeting a lipid specific to mitochondria is of great interest. In the coarse-grained representation, the bead types at selected positions of the acridine orange molecule are changed and the difference in binding affinity between the original and modified structures are determined by free energy calculations. The present results from the coarse-grained simulations indicate that hydrogen bonding has the desired effect. As a next step, the information lost in the coarse-graining process is to be reintroduced to the most promising subset of the screened molecules by repeating the free energy calculations in atomistic detail.
Standard binding free energies of SIM SUMO complexes
Alexander Kötter | Universität Münster, GermanyThe interaction between the small ubiquitin related modifier (SUMO) and different sumo interacting motifs (SIMs) has been the subject of a number of studies in the past years. We focus on the SIM regions of PIAS and RNF4. NMR studies on diSUMO and peptides containing two SIMs of RNF4 revealed the existence of an antiparallel binding mode for the SIMs [1]. This contradicts the notion, that these SIMs only bind parallel to SUMO [2]. Using atomistic molecular dynamics simulations we show that for this kind of SIMs the parallel binding mode is strongly favored. These simulations are very challenging due to the high flexibility of the SUMO binding peptides and require optimized free energy methods [3]. However, an antiparallel binding mode exists as well and also in this mode the SIM extends the beta sheets of SUMO, which is typical for SIM SUMO complexes. Our simulations furthermore point towards the importance of the neighborhood of considered SIMs for high affinity binding.
References
[1] Xu et al. Nat. Comm. 2014, 5
[2] Song et al. J. Biol. Chem. 2005, 280
[3] Gumbart et al. JCTC 2012, 9
Molecular Dynamics Study of Pressure-induced Effects on the Self-cleavage Catalysis Reaction of Hairpin Ribozyme
Narendra Kumar | Ruhr-University Bochum, GermanyThe discovery of ribozymes, i.e. RNAs with an ability to catalyze biological reactions, gave strong support to "RNA world hypothesis" suggesting RNA as a key molecular precursor of life. The early life might have evolved under extreme conditions such as high pressures found in the deep sea. In a joint effort with experiment, we study the self-cleavage reaction, which is key to ribozyme function, under high pressure conditions. Using large-scale replica-exchange molecular dynamics in conjunction with a state-of-the-art RNA force field we sample the free energy landscapes of various structures along the reaction pathway at ambient conditions and at a pressure of 10 kbar. The simulations disclose that active configurations in the so-called activated precursor state, being close to the transition state of the self-cleavage step and prepared for in-line attack, becomes energetically preferred upon compression. Moreover, the free energy barrier from inactive toward active configurations becomes very low. These findings not only suggest that the self-cleavage step is favored both thermodynamically and kinetically, but are qualitatively consistent with the experimental results from our partner group. At molecular level, our simulations disclose that favorable H-bonding interactions between the reaction site and specific neighboring nucleobases are enhanced under pressure [1]. These results suggest that, given favorable environmental conditions, ribozymes could also accelerate the overall cleavage process, hence, act as piezophilic biocatalysts in addition to their hitherto known properties.
References
[1] C Schuabb*, N Kumar*, S Pataraia*, D Marx, R Winter, Nat. Commun. 2017, 8, 14661. *equal contributions
Dissipation-corrected targeted molecular dynamics for the calculation of free energies and friction factors from non-equilibrium simulations of protein-ligand unbinding
Steffen Wolf | University of Freiburg, GermanyWe present dissipation-corrected targeted MD (dcTMD), [1] a new method for the calculation of free energy profiles and friction factors from constant velocity constraint simulations. Performing a cumulant expansion of the Jarzynski identity [2] and comparing it with a Markovian Langevin equation, we derive an expression for the dissipative work and underlying friction factors, which allow on-the-fly correction of the non-equilibrium work to yield the underlying free energy. Besides friction factors giving valuable insight into protein-ligand and ligand-solvent interaction such as the formation and rearrangement of hydration shells, the combination of both free energy and friction factors allow post-MD Langevin dynamics coarse graining, and thus potentially the calculation of absolute equilibrium process kinetics. Having established the technique with the example of a NaCl-ion pair dissociation in water, we present the extension of the method to the unbinding of several inhibitors of the N-terminal domain of Hsp90 and of agonists/antagonists of the beta(2) adrenergic receptor. In these “real life” examples, the application of dcTMD requires the detection of pathways taken by the ligand through the protein, and sorting of trajectories according to such pathways. For this purpose, we developed two sorting methods, with the first being based on machine learning – non-equilibrium PCA [3,4] and the second one on RMSD distance matrices evaluated by the NeighborNet algorithm. [5] We discuss both benefits and limits / pitfalls of these approaches. [6] Last, the method can be used as a “quick-and-dirty” technique based on “deep” non-equilibrium simulations for the scoring of ligands into fast and slowly unbinding compounds. [7]
References
[1] S. Wolf and G. Stock, J. Chem. Theory Comput. 14, 6175 (2018).
[2] C. Jarzynski, Phys. Rev. Lett. 78, 2690 (1997).
[3] S. Brandt, F. Sittel, M. Ernst, and G. Stock, J Phys Chem Lett 2144 (2018).
[4] M. Post, S. Wolf, and G. Stock, under rewiev (2019).
[5] D. Bryant and V. Moulton, Mol Biol Evol 21, 255 (2004).
[6] S. Bray and S. Wolf, in prep. (2019)
[7] S. Wolf, M. Amaral et al., in revision (2019).
Structure-Based Ligand Design by Targeting an Ordered Water: Interact, Displace, or Replace?
Pierre Matricon | Uppsala University, SwedenFormation of a protein-ligand complex results in the displacement of water molecules from the binding site, which is one of the major driving forces of association. Crystal structures have also revealed that hydrogen bonds to ordered waters that bridge protein-ligand interactions can be essential for molecular recognition and contribute to affinity. Thus, the possibility to identify and exploit binding site waters in ligand optimization would be extremely valuable in the design of drug candidates. A highly ordered water molecule that locks the bound ligand conformation was identified in an adenosine-bound crystal structure of the A2A adenosine receptor. In this work, alchemical free energy calculations were used to evaluate the effect of altering hydration sites around a reference ligand with different substitutions. Inhomogeneous Fluid Solvation Theory (IFST) was also used to unravel the contribution of different hydration sites to relative binding thermodynamics.
Building fitness landscapes for antimicrobial resistance from free energy calculations
Martin Carballo-Pacheco | University of Edinburgh, GermanyMicrobes are evolving resistance to drugs rendering them inefficient. De novo evolution of resistance is caused by mutations in genes that lead to drug-resistant phenotypes. Evolution of resistance can be understood as trajectories in a fitness landscape, which map a phenotypic effect (resistance to drug) onto the highly multidimensional genotypic space. Features from these landscapes determine how evolution ensues and which phenotypes are accessible in an specific environment. However, measuring fitness landscapes experimentally is challenging due to their high dimensionality. Here, we propose a computational model to predict the binding of antimicrobial drugs to their target protein and build the fitness landscape for antimicrobial resistance by taking into account the atomistic details of both the target and the drug.
Development of binding free energy computational models for GPCRs using the MM-PBSA method
E.Vrontaki, P.Lagarias, E.Tzortzini, M. Stampelou, A.Kolocouris | University of Athens, GreeceG-Protein Coupled Receptors (GPCRs) are implicated in many diseases and are targeted by 34% of FDA approved drugs. Important members in the GPCRs family A are adenosine receptors, adrenergic receptors, muscarinic acetylcholine receptors, angiotensin II receptors, dopamine receptors etc. Adenosine receptors are important regulators of many physiological and pathophysiological functions and therefore important targets for drug development. Pharmaceutical companies and academic laboratories perform intense efforts to identify selective antagonists for each adenosine receptor subtype (AR) as potential clinical candidates for "soft" treatment of diseases. A2AAR antagonists can be useful for treating cancer, CNS disorders; A1AR antagonists can provide antiasthamatic and CNS agents; A3AR antagonists are promising for therapeutic applications in asthma, glaucoma and A2BAR antagonists for diabetes, asthma. Aim of our project was the development of computational models that could foresight the potency of a molecule based on computer-aided free energy calculations between a set of molecules with different structure covering a range of three pKi unites before any would be attempted. For adenosine receptors sub-types A1, A2A and A3 representative ligands have been selected from literature and using MD simulations and MM-PBSA calculations the binding free energies were calculated and a promising correlation R2 > 0.8 was found with experimenta pKi values. The models can be a useful predictive tool in the drug discovery against ARs, helping understanding deeper ligand-receptor interactions and reducing the cost for drug discovery process.
Exploiting Enhanced Sampling Methods to Design Novel Allosteric Inhibitors of the Aromatase Enzyme
Angelo Spinello | SISSA, Trieste, ItalyHuman Aromatase (HA) is a cytochrome P450 of key importance as therapeutic target against breast cancer. In fact, one of the most currently used therapeutic strategy relies on HA inhibition. HA catalyzes the highly specific conversion of androgens to estrogens in a three-step reaction using molecular oxygen and electrons. Surprisingly, our study of the last reaction step, using QM/MM metadynamics (MTD), revealed a novel mechanism in which HA performs its catalytic function similarly to all other cytochome P450s. One allosteric pocket lies exactly between Asp309 and Arg192, two residues which proved to be crucial for the formation of the reactive species. Moreover, we have exploited the information obtained in order to develop a novel generation of inhibitors, targeting aromatase allosteric pockets, via a sophisticate computational protocol based on in silico screening of different level of complexity across commercial databases of molecules, followed by cumulative multi-μs molecular dynamics and free energy simulations. Finally, we have employed MTD simulations to provide a refinement of the inhibitors binding pose and furnish information about their dissociation kinetics, which is nowadays increasingly being recognized as a key tract in drug efficacy optimization study.
PIP2 Modulation of KATP Channel and Disease
Tanadet Pipatpolkai, R.A. Corey, F.M. Ashcroft, P.J. Stansfeld | University of Oxford, UKThe pancreatic ATP-regulated potassium (KATP) channel is required for insulin secretion during fed state to maintain blood glucose level. The KATP channel is an octameric protein made of four Kir6.2 channel subunits and four SUR1 regulatory subunits. The Kir6.2 channel is activated by phosphatidyl inositol-4,5-bisphosphate (PIP2), a low concentration anionic lipid. Here, we have used free energy calculations to investigate the contribution of the PIP2 molecules to channel affinity. Using coarse-grained molecular dynamics simulation (CG-MD), we were able to show that the PIP2 binding free energy is mainly a contribution from the glycerol and the inositol headgroup. With an extension to this approach, we were also able to show the variation of binding free energy of other inositide lipids (PIP3, PIP2 and PI4P) to Kir6.2. Lastly, we were also able to illustrate the relationship between neonatal diabetes and congenital hyperinsulinism mutations on PIP2 binding to Kir6.2. This provides us with new insights into how defects in lipid binding may cause human disease.
Multiple-replica lambda-dynamics for the Calculation of Alchemical Free-Energies: The Conveyor Belt Thermodynamic Integration (CBTI) scheme
David F. Hahn, Philippe H. Hünenberger | ETH Zurich, SwitzerlandA new method is proposed to calculate alchemical free-energy differences based on molec- ular dynamics (MD) simulations, called the conveyor belt thermodynamic integration (CBTI) scheme. As in thermodynamic integration (TI), K replicas of the system are simulated at dif- ferent values of the alchemical coupling parameter λ. The number K is taken to be even and the replicas are equally spaced on a forward-turn-backward-turn path, akin to a conveyor belt (CB) between the two physical end-states. And as in λ-dynamics (λD), the λ-values associated with the individual systems evolve in time along the simulation. However, they do so in a concerted fashion, determined by the evolution of a single dynamical variable Λ of period 2π controlling the advance of the entire CB. Thus, a change of Λ is always associated with half of the replicas moving forward and the other half moving backward along λ. As a result, the effective free-energy profile of the replica system along Λ is periodic of period 2πK−1 and the magnitude of its variations decreases rapidly upon increasing K, at least as K−1 in the limit of large K. When a sufficient number of replicas is used, these variations become small, which enables a complete and quasi-homogeneous coverage of the λ-range by the replica system, without application of any biasing potential. If desired, a memory-based biasing potential can still be added to further homogenize the sampling, the preoptimization of which is computationally inexpensive. The final free-energy profile along λ is calculated similarly to TI, by binning of the Hamiltonian λ-derivative as a function of λ considering all replicas jointly, followed by quadrature integration. The associated quadrature error can be kept very low owing to the continuous and quasi-homogeneous λ-sampling. The CBTI scheme can be viewed as a continuous/deterministic/dynamical analog of the Hamiltonian replica-exchange/permutation (HRE/HRP) schemes, or as a correlated multiple-replica analog of the λD schemes. Compared to TI, it shares the advantage of the latter schemes in terms of enhanced orthogonal sampling, i.e. the availability of variable-λ paths to circumvent conformational barriers present at specific λ-values. Compared to HRE/HRP, it permits a deterministic and continuous sampling of the λ-range, is expected to be less sensitive to possible artifacts of the thermo- and barostatting schemes, and bypasses the need to carefully preselect a λ-ladder and a swapping-attempt fre- quency. Compared to λ-D, it facilitates barrier crossing along λ. CBTI is applied to the calculation of the hydration free energy of methanol, to the mutation free energy of guanosine triphosphate (GTP) to 8-Br-GTP and to the side chain mutation of a tripeptide.
Exploring Strategies for Labeling Viruses with Gold Nanoclusters through Non-equilibrium Molecular Dynamics Simulations
Emmi Pohjolainen, Gerrit Groenhof, Hannu Häkkinen | University of Jyväskylä, FinlandBiocompatible gold nanoclusters can be utilized as contrast agents in virus imaging. The labeling of viruses can be achieved noncovalently but site-specifically by linking the cluster to the hydrophobic pocket of a virus via a lipid-like pocket factor. We have estimated the binding affinities of three different pocket factors of echovirus 1 (EV1) in molecular dynamics simulations combined with non-equilibrium free-energy calculations [2]. We have also studied the effects on binding affinities with a pocket factor linked to the Au102pMBA44 nanocluster in different protonation states. Although the absolute binding affinities are over-estimated for all the systems, the trend is in agreement with recent experiments [1]. Our results suggest that the natural pocket factor (palmitic acid) can be replaced by molecules pleconaril (drug) and its derivative Kirtan1 that have higher estimated binding affinities. Our results also suggest that including the gold nanocluster does not decrease the affinity of the pocket factor to the virus, but the affinity is sensitive to the protonation state of the nanocluster, i.e., to pH conditions. The methodology introduced in this work helps in the design of optimal strategies for gold–virus bioconjugation for virus detection and manipulation.
References
[1] Martikainen, M., Salorinne, K., Lahtinen, T., Malola, S., Permi, P., Häkkinen, H., and Marjomäki, V. Nanoscale 2015 7, 17457−17467.
[2] Pohjolainen, E., Malola, S., Groenhof, G., and Häkkinen, H. Simulations. Bioconjugate Chemistry 2017 28 (9), 2327-2339.
Ligand and sidechain mutations by FEP simulations: two sides of the same coin
Willem Jespers | Uppsala University, SwedenThe characterization of ligand binding affinities against a set of mutant proteins (Site-directed mutagenesis, SDM), combined with structure-affinity relationships (SAR) of a ligand series, can provide insight in key receptor-ligand interactions. These data can be mapped onto a structure of the protein, resulting in a detailed understanding of the molecular interactions responsible for affinity, selectivity and receptor function. Based on this philosophy, we have developed a pipeline of computational approaches to integrate the available experimental information, including the emerging landscape of GPCR structures, with structure based drug design methods, in particular free energy perturbation (FEP) simulations. Our method is based on molecular dynamics (MD) sampling of the protein-ligand binding site using spherical boundary conditions in our open source MD software Q. The workflow includes protocols for protein side chain mutations (QresFEP) as well as ligand mutations (QligFEP) written in Python. Here I will describe how we use the synergy of these methods in the design of new chemical moieties targeting the G Protein-Coupled Receptor (GPCR) family, with a particular emphasis on adenosine receptors. These examples include: 1) FEP-guided modelling, synthesis and pharmacological characterization of a series of A2A antagonists to predict the binding mode, later confirmed by x-ray crystallography; 2) The conformational selectivity of agonists and antagonists for the active and inactive states of the receptor and 3) The optimized design of different series of specific chemical modulators of the A2B receptor.
Our methodology is freely available under http://open.gpcr-modsim.org/ and https://github.com/qusers/Q6
References
[1] a) Keranen, H.; Aqvist, J.; Gutierrez-de-Teran, H. Chem Comm 2015, 51:3522; b) Boukharta, L.; Gutierrez-de-Teran, H.; Aqvist, J. PlOS Comp. Biol 2014, 10:e1003585
[2] Jespers, W et al. J. Chem Inf 2019, 11:26
[3] Jespers, W et al. Trends Pharmacol. Sci. 2018, 39:75-89
[4] Jespers, W et al. Molecules 2017, 22 (11), 1945
[5] Jespers, W et al. J. Med. Chem., 2017, 60 (17), pp 7502–7511
Free Energy Calculations as a Tool for Lead Optimization in Drug Discovery
Dimitris Ntekoumes, Zoe Cournia | Biomedical Research Foundation of the Academy of Athens, GreeceOne of the most demanding tasks in the lead optimization phase of the drug design process is to predict, among a series of lead candidates, which ones will bind more strongly to the therapeutic target. One of the most promising approaches to calculate the binding free energy is alchemical Free Energy Perturbation (FEP) in the framework of absolute or relative free energy of binding, due to its rigorous statistical mechanics background. However, the high computational costs, the limited force field accuracy, and the technical challenges to setup, run and analyse free energy simulations, have hampered the success and throughput of free energy simulations in drug discovery. In order to resolve these issues during the past few decades, the computational chemistry community has made a significant effort to accurately estimate the binding affinities through validation studies, blind challenge predictions and prospective applications. In this direction, we focused on two different cases, in order to further contribute to this endeavour. First, we study a series of 12 analogs of CK-666, a micromolar Arp2/3 inhibitor, and compare the difference in the free energy of binding estimates for the analog series between two popular molecular simulation packages, NAMD and GROMACS, and against experimental data obtained by our collaborators. Our data reveal that the predicted binding affinities vs experimental results, for both packages, lie within the error of the method (1 – 2 kcal/mol) that current studies suggest. Moreover, GROMACS simulations reveal improved correlation than NAMD, with the mean unsigned error being 1.20 and 2.53 kcal/mol, respectively. Second, we participate in the SAMPL6 challenge, which is focused on evaluating the convergence and reproducibility (across codes) of absolute free energy of binding predictions. To this end, we are running absolute free energy calculations for three host – guest systems using NAMD to assess the convergence of the calculation, and to also compare the binding affinity estimates against experimental data and other simulation platforms. Our predictions are in agreement with the reference calculations performed with YANK, with the deviation being ~1 kcal/mol, albeit both packages deviate significantly from the experimental results. Following the calculations, we will investigate further the underlying reasons for these deviations with an eye on force field parameters, convergence of the FEP calculations, and more realistic treatment of the environment of the protein/ligand systems.