1
|
Hattori S, Zhu Q. Revisiting Aspirin Polymorphic Stability Using a Machine Learning Potential. ACS OMEGA 2024; 9:36589-36599. [PMID: 39220495 PMCID: PMC11360032 DOI: 10.1021/acsomega.4c04782] [Citation(s) in RCA: 0] [Impact Index Per Article: 0] [Reference Citation Analysis] [Abstract] [Track Full Text] [Download PDF] [Figures] [Subscribe] [Scholar Register] [Received: 05/21/2024] [Revised: 07/26/2024] [Accepted: 07/30/2024] [Indexed: 09/04/2024]
Abstract
In this study, we present a systematic computational investigation to analyze the long-debated free energy stability of two well-known aspirin polymorphs, denoted as Form I and Form II. Specifically, we developed a strategy to collect training configurations covering diverse interatomic interactions between representative functional groups in aspirin crystals. Utilizing a state-of-the-art neural network interatomic potential (NNIP) model, we trained an accurate machine learning potential to simulate aspirin crystal dynamics under finite temperature conditions with ∼0.46 kJ/mol/molecule accuracy. Employing the trained NNIP model, we performed thermodynamic integration to assess the free energy difference between aspirins I and II, accounting for the anharmonic effects in a large supercell consisting of 512 molecules. For the first time, our results convincingly demonstrated that Form I is more stable than Form II at 300 K, ranging from 0.74 to 1.83 kJ/mol/molecule, aligning with experimental observations. Unlike the majority of previous simulations based on (quasi)harmonic approximations in a small super cell, which often found degenerate energies between aspirins I and II, our findings underscore the importance of anharmonic effects in determining polymorphic stability ranking. Furthermore, we proposed the use of the rotational degrees of freedom of methyl and ester/phenyl groups in aspirin crystals as characteristic motions to highlight rotational entropic contribution that favors the stability of Form I. From the structural perspective, we also found that the subtle free energy difference can be used to explain the distinct thermal expansion responses as observed in both experimental and simulation data. Beyond the aspirin polymorphism, we anticipate that such entropy-driven stabilization can be broadly applicable to many other organic systems, suggesting that our approach holds great promise for stability studies in small-molecule drug design.
Collapse
Affiliation(s)
- Shinnosuke Hattori
- Advanced
Research Laboratory, Research Platform, Sony Group Corporation, 4−14−1 Asahi-cho, Atsugi-shi 243−0014, Japan
| | - Qiang Zhu
- Department
of Mechanical Engineering and Engineering Science, University of North Carolina at Charlotte, Charlotte, North Carolina 28223, United States
| |
Collapse
|
2
|
Tan S, Gong X, Liu H, Yao X. Identification of novel LRRK2 inhibitors by structure-based virtual screening and alchemical free energy calculation. Phys Chem Chem Phys 2024; 26:19775-19786. [PMID: 38984923 DOI: 10.1039/d4cp01762e] [Citation(s) in RCA: 0] [Impact Index Per Article: 0] [Reference Citation Analysis] [Abstract] [MESH Headings] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 07/11/2024]
Abstract
The Leucine-rich repeat kinase 2 (LRRK2) target has been identified as a promising drug target for Parkinson's disease (PD) treatment. This study focuses on optimizing the activity of LRRK2 inhibitors using alchemical relative binding free energy (RBFE) calculations. Initially, we assessed various free energy calculation methods across different LRRK2 kinase inhibitor scaffolds. The results indicate that alchemical free energy calculations are promising for prospective predictions on LRRK2 inhibitors, especially for the aminopyrimidine scaffold with an RMSE of 1.15 kcal mol-1 and Rp of 0.83. Following this, we optimized a potent LRRK2 kinase inhibitor identified from previous virtual screenings, featuring a novel scaffold. Guided by RBFE predictions using alchemical methods, this optimization led to the discovery of compound LY2023-001. This compound, with a [1,2,4]triazolo[5,6-b]indole scaffold, exhibited enhanced inhibitory activity against G2019S LRRK2 (IC50 = 12.9 nM). Molecular dynamics (MD) simulations revealed that LY2023-001 formed stable hydrogen bonds with Glu1948, and Ala1950 in the G2019S LRRK2 protein. Additionally, its phenyl substituents engage in strong electrostatic interactions with Lys1906 and van der Waals interactions with Leu1885, Phe1890, Val1893, Ile1933, Met1947, Leu1949, Leu2001, Ala2016, and Asp2017. Our findings underscore the potential of computational methods in the successful optimization of small molecules, offering important insights for the development of novel LRRK2 inhibitors.
Collapse
Affiliation(s)
- Shuoyan Tan
- College of Chemistry and Chemical Engineering, Lanzhou University, Lanzhou 730000, China
| | - Xiaoqing Gong
- Faculty of Applied Sciences, Macao Polytechnic University, Macao SAR, 999078, China.
| | - Huanxiang Liu
- Faculty of Applied Sciences, Macao Polytechnic University, Macao SAR, 999078, China.
| | - Xiaojun Yao
- Faculty of Applied Sciences, Macao Polytechnic University, Macao SAR, 999078, China.
| |
Collapse
|
3
|
Lagardère L, Maurin L, Adjoua O, El Hage K, Monmarché P, Piquemal JP, Hénin J. Lambda-ABF: Simplified, Portable, Accurate, and Cost-Effective Alchemical Free-Energy Computation. J Chem Theory Comput 2024; 20:4481-4498. [PMID: 38805379 DOI: 10.1021/acs.jctc.3c01249] [Citation(s) in RCA: 0] [Impact Index Per Article: 0] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 05/30/2024]
Abstract
We introduce the lambda-Adaptive Biasing Force (lambda-ABF) method for the computation of alchemical free-energy differences. We propose a software implementation and showcase it on biomolecular systems. The method arises from coupling multiple-walker adaptive biasing force with λ-dynamics. The sampling of the alchemical variable is continuous and converges toward a uniform distribution, making manual optimization of the λ schedule unnecessary. Contrary to most other approaches, alchemical free-energy estimates are obtained immediately without any postprocessing. Free diffusion of λ improves orthogonal relaxation compared to fixed-λ thermodynamic integration or free-energy perturbation. Furthermore, multiple walkers provide generic orthogonal space coverage with minimal user input and negligible computational overhead. We show that our high-performance implementations coupling the Colvars library with NAMD and Tinker-HP can address real-world cases including ligand-receptor binding with both fixed-charge and polarizable models, with a demonstrably richer sampling than fixed-λ methods. The implementation is fully open-source, publicly available, and readily usable by practitioners of current alchemical methods. Thanks to the portable Colvars library, lambda-ABF presents a unified user interface regardless of the back-end (NAMD, Tinker-HP, or any software to be interfaced in the future), sparing users the effort of learning multiple interfaces. Finally, the Colvars Dashboard extension of the visual molecular dynamics (VMD) software provides an interactive monitoring and diagnostic tool for lambda-ABF simulations.
Collapse
Affiliation(s)
- Louis Lagardère
- Sorbonne Université, Laboratoire de Chimie Théorique, UMR 7616 CNRS, Paris 75005, France
- Sorbonne Université, Institut Parisien de Chimie Physique et Théorique, FR2622 CNRS, 75005 Paris, France
- Qubit Pharmaceuticals, 29 rue du Faubourg Saint Jacques, 75014 Paris, France
| | - Lise Maurin
- Sorbonne Université, Laboratoire de Chimie Théorique, UMR 7616 CNRS, Paris 75005, France
- Sorbonne Université, Laboratoire Jacques-Louis Lions, UMR 7589 CNRS, 75005 Paris, France
| | - Olivier Adjoua
- Sorbonne Université, Laboratoire de Chimie Théorique, UMR 7616 CNRS, Paris 75005, France
| | - Krystel El Hage
- Qubit Pharmaceuticals, 29 rue du Faubourg Saint Jacques, 75014 Paris, France
| | - Pierre Monmarché
- Sorbonne Université, Laboratoire de Chimie Théorique, UMR 7616 CNRS, Paris 75005, France
- Sorbonne Université, Laboratoire Jacques-Louis Lions, UMR 7589 CNRS, 75005 Paris, France
| | - Jean-Philip Piquemal
- Sorbonne Université, Laboratoire de Chimie Théorique, UMR 7616 CNRS, Paris 75005, France
- Qubit Pharmaceuticals, 29 rue du Faubourg Saint Jacques, 75014 Paris, France
| | - Jérôme Hénin
- Laboratoire de Biochimie Théorique, Université Paris Cité, CNRS, UPR 9080, 75005 Paris, France
| |
Collapse
|
4
|
Wang X, Nayak S, Wilson RE, Soderholm L, Servis MJ. Solvent effects on extractant conformational energetics in liquid-liquid extraction: a simulation study of molecular solvents and ionic liquids. Phys Chem Chem Phys 2024; 26:2877-2886. [PMID: 38048065 DOI: 10.1039/d3cp04680j] [Citation(s) in RCA: 0] [Impact Index Per Article: 0] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 12/05/2023]
Abstract
Extractant design in liquid-liquid extraction (LLE) is a research frontier of metal ion separations that typically focuses on the direct extractant-metal interactions. However, a more detailed understanding of energetic drivers of separations beyond primary metal coordination is often lacking, including the role of solvent in the extractant phase. In this work, we propose a new mechanism for enhancing metal-complexant energetics with nanostructured solvents. Using molecular dynamics simulations with umbrella sampling, we find that the organic solvent can reshape the energetics of the extractant's intramolecular conformational landscape. We calculate free energy profiles of different conformations of a representative bidentate extractant, n-octyl(phenyl)-N,N-diisobutyl carbamoyl methyl phosphinoxide (CMPO), in four different solvents: dodecane, tributyl phosphate (TBP), and dry and wet ionic liquid (IL) 1-ethyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imide ([EMIM][Tf2N]). By promoting reorganization of the extractant molecule into its binding conformation, our findings reveal how particular solvents can ameliorate this unfavorable step of the metal separation process. In particular, the charge alternating nanodomains formed in ILs substantially reduce the free energy penalty associated with extractant reorganization. Importantly, using alchemical free energy calculations, we find that this stabilization persists even when we explicitly include the extracted cation. These findings provide insight into the energetic drivers of metal ion separations and potentially suggest a new approach to designing effective separations using a molecular-level understanding of solvent effects.
Collapse
Affiliation(s)
- Xiaoyu Wang
- Chemical Sciences and Engineering Division, Argonne National Laboratory, 9700 S Cass Ave, Lemont, IL 60439, USA.
| | - Srikanth Nayak
- Chemical Sciences and Engineering Division, Argonne National Laboratory, 9700 S Cass Ave, Lemont, IL 60439, USA.
| | - Richard E Wilson
- Chemical Sciences and Engineering Division, Argonne National Laboratory, 9700 S Cass Ave, Lemont, IL 60439, USA.
| | - L Soderholm
- Chemical Sciences and Engineering Division, Argonne National Laboratory, 9700 S Cass Ave, Lemont, IL 60439, USA.
| | - Michael J Servis
- Chemical Sciences and Engineering Division, Argonne National Laboratory, 9700 S Cass Ave, Lemont, IL 60439, USA.
| |
Collapse
|
5
|
York DM. Modern Alchemical Free Energy Methods for Drug Discovery Explained. ACS PHYSICAL CHEMISTRY AU 2023; 3:478-491. [PMID: 38034038 PMCID: PMC10683484 DOI: 10.1021/acsphyschemau.3c00033] [Citation(s) in RCA: 1] [Impact Index Per Article: 1.0] [Reference Citation Analysis] [Abstract] [Key Words] [Grants] [Track Full Text] [Download PDF] [Figures] [Subscribe] [Scholar Register] [Received: 07/17/2023] [Revised: 09/12/2023] [Accepted: 09/13/2023] [Indexed: 12/02/2023]
Abstract
This Perspective provides a contextual explanation of the current state-of-the-art alchemical free energy methods and their role in drug discovery as well as highlights select emerging technologies. The narrative attempts to answer basic questions about what goes on "under the hood" in free energy simulations and provide general guidelines for how to run simulations and analyze the results. It is the hope that this work will provide a valuable introduction to students and scientists in the field.
Collapse
Affiliation(s)
- Darrin M. York
- Laboratory for Biomolecular
Simulation Research, Institute for Quantitative Biomedicine, and Department
of Chemistry and Chemical Biology, Rutgers
University, Piscataway, New Jersey 08854, United States
| |
Collapse
|
6
|
Gnyawali SC, Denune JA, Hockman B, Kristjánsdóttir JV, Ragnarsdóttir MS, Timsina LR, Ghatak S, Lechler K, Sen CK, Roy S. Moisture mitigation using a vented liner and a vented socket system for individuals with transfemoral amputation. Sci Rep 2023; 13:16557. [PMID: 37783779 PMCID: PMC10545693 DOI: 10.1038/s41598-023-43572-2] [Citation(s) in RCA: 0] [Impact Index Per Article: 0] [Reference Citation Analysis] [Abstract] [Key Words] [MESH Headings] [Grants] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Received: 05/04/2023] [Accepted: 09/26/2023] [Indexed: 10/04/2023] Open
Abstract
Sweating and heat buildup at the skin-liner interface is a major challenge for persons with limb loss. Liners made of heat-non-conducting materials may cause sweating of the residual limb and may result in liners slipping off the skin surface especially on a warm day or during high activity, causing skin breakdown and affecting limb health. To address this, we evaluated the efficacy of the vented liner-socket system (VS, Össur) compared to Seal-In silicone liner and non-vented socket (nVS, Össur) in reducing relative humidity (RH) during increased sweat. Nine individuals with limb loss using nVS were randomized to VS or nVS and asked for activity in a 20-min treadmill walk. RH was significantly attenuated (p = 0.0002) and perceived sweating, as reported by prosthesis users, improved (p = 0.028) with VS, patient-reported comprehensive lower limb amputee socket survey (CLASS) outcomes to determine the suspension, stability, and comfort were not significantly different between VS and nVS. There are limited rigorous scientific studies that clearly provide evidence-based guidelines to the prosthetist in the selection of liners from numerous available options. The present study is innovative in clearly establishing objective measures for assessing humidity and temperatures at the skin-liner interface while performing activity. As shown by the measured data and perceived sweat scores provided by the subjects based on their daily experience, this study provided clear evidence establishing relative humidity at the skin-liner interface is reduced with the use of a vented liner-socket system when compared to a similar non-vented system.
Collapse
Affiliation(s)
- Surya C Gnyawali
- McGowan Institute for Regenerative Medicine (MIRM), Department of Surgery, University of Pittsburgh School of Medicine, Pittsburgh, PA, USA.
| | - Jeffrey A Denune
- Indiana Centre for Regenerative Medicine and Engineering (ICRME), Indiana University Health Comprehensive Wound Centre, Indiana University School of Medicine, Indianapolis, IN, USA
| | - Bryce Hockman
- Indiana Centre for Regenerative Medicine and Engineering (ICRME), Indiana University Health Comprehensive Wound Centre, Indiana University School of Medicine, Indianapolis, IN, USA
| | | | | | - Lava R Timsina
- Indiana Centre for Regenerative Medicine and Engineering (ICRME), Indiana University Health Comprehensive Wound Centre, Indiana University School of Medicine, Indianapolis, IN, USA
| | - Subhadip Ghatak
- McGowan Institute for Regenerative Medicine (MIRM), Department of Surgery, University of Pittsburgh School of Medicine, Pittsburgh, PA, USA
- Indiana Centre for Regenerative Medicine and Engineering (ICRME), Indiana University Health Comprehensive Wound Centre, Indiana University School of Medicine, Indianapolis, IN, USA
| | - Knut Lechler
- Össur Ehf., R&D, Medical Office, Reykjavik, Iceland
| | - Chandan K Sen
- McGowan Institute for Regenerative Medicine (MIRM), Department of Surgery, University of Pittsburgh School of Medicine, Pittsburgh, PA, USA
- Indiana Centre for Regenerative Medicine and Engineering (ICRME), Indiana University Health Comprehensive Wound Centre, Indiana University School of Medicine, Indianapolis, IN, USA
| | - Sashwati Roy
- McGowan Institute for Regenerative Medicine (MIRM), Department of Surgery, University of Pittsburgh School of Medicine, Pittsburgh, PA, USA.
- Indiana Centre for Regenerative Medicine and Engineering (ICRME), Indiana University Health Comprehensive Wound Centre, Indiana University School of Medicine, Indianapolis, IN, USA.
| |
Collapse
|
7
|
Lundborg M, Lidmar J, Hess B. On the Path to Optimal Alchemistry. Protein J 2023; 42:477-489. [PMID: 37651042 PMCID: PMC10480267 DOI: 10.1007/s10930-023-10137-1] [Citation(s) in RCA: 0] [Impact Index Per Article: 0] [Reference Citation Analysis] [Abstract] [Key Words] [Grants] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Accepted: 07/04/2023] [Indexed: 09/01/2023]
Abstract
Alchemical free energy calculations have become a standard and widely used tool, in particular for calculating and comparing binding affinities of drugs. Although methods to compute such free energies have improved significantly over the last decades, the choice of path between the end states of interest is usually still the same as two decades ago. We will show that there is a fundamentally arbitrary, implicit choice of parametrization of this path. To address this, the notion of the length of a path or a metric is required. A metric recently introduced in the context of the accelerated weight histogram method also proves to be very useful here. We demonstrate that this metric can not only improve the efficiency of sampling along a given path, but that it can also be used to improve the actual choice of path. For a set of relevant use cases, the combination of these improvements can increase the efficiency of alchemical free energy calculations by up to a factor 16.
Collapse
Affiliation(s)
| | - Jack Lidmar
- Department of Physics, KTH Royal Institute of Technology, 10691, Stockholm, Sweden
| | - Berk Hess
- Department of Applied Physics, KTH Royal Institute of Technology, 10691, Stockholm, Science for Life Laboratory, Solna, Sweden.
| |
Collapse
|
8
|
Liao J, Shu Z, Gao J, Wu M, Chen C. SurfPB: A GPU-Accelerated Electrostatic Calculation and Visualization Tool for Biomolecules. J Chem Inf Model 2023; 63:4490-4496. [PMID: 37500509 DOI: 10.1021/acs.jcim.3c00745] [Citation(s) in RCA: 0] [Impact Index Per Article: 0] [Reference Citation Analysis] [Abstract] [MESH Headings] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 07/29/2023]
Abstract
In this work, we present SurfPB as a useful tool for the study of biomolecules. It can do many typical calculations, including the molecular surface, electrostatic potential, solvation free energy, entropy, and binding free energy. Among all of the calculations, the entropy calculation is the most time-consuming one. In SurfPB, the calculation can be performed in a vacuum or implicit solvent and accelerated on GPU. The Poisson-Boltzmann equation solver is accelerated on GPU as well. Moreover, we developed a graphical user interface for SurfPB. It allows users to input the parameters and complete the whole calculation in a visual way. The calculated electrostatic potentials are shown on the molecular surface in a three-dimensional scene.
Collapse
Affiliation(s)
- Jun Liao
- Biomolecular Physics and Modeling Group, School of Physics Huazhong University of Science and Technology, Wuhan 430074, Hubei, China
| | - Zirui Shu
- Biomolecular Physics and Modeling Group, School of Physics Huazhong University of Science and Technology, Wuhan 430074, Hubei, China
| | - Junyong Gao
- Biomolecular Physics and Modeling Group, School of Physics Huazhong University of Science and Technology, Wuhan 430074, Hubei, China
| | - Mincong Wu
- Biomolecular Physics and Modeling Group, School of Physics Huazhong University of Science and Technology, Wuhan 430074, Hubei, China
| | - Changjun Chen
- Biomolecular Physics and Modeling Group, School of Physics Huazhong University of Science and Technology, Wuhan 430074, Hubei, China
| |
Collapse
|
9
|
Bassani D, Moro S. Past, Present, and Future Perspectives on Computer-Aided Drug Design Methodologies. Molecules 2023; 28:molecules28093906. [PMID: 37175316 PMCID: PMC10180087 DOI: 10.3390/molecules28093906] [Citation(s) in RCA: 25] [Impact Index Per Article: 25.0] [Reference Citation Analysis] [Abstract] [Key Words] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Received: 03/22/2023] [Revised: 04/28/2023] [Accepted: 05/02/2023] [Indexed: 05/15/2023] Open
Abstract
The application of computational approaches in drug discovery has been consolidated in the last decades. These families of techniques are usually grouped under the common name of "computer-aided drug design" (CADD), and they now constitute one of the pillars in the pharmaceutical discovery pipelines in many academic and industrial environments. Their implementation has been demonstrated to tremendously improve the speed of the early discovery steps, allowing for the proficient and rational choice of proper compounds for a desired therapeutic need among the extreme vastness of the drug-like chemical space. Moreover, the application of CADD approaches allows the rationalization of biochemical and interactive processes of pharmaceutical interest at the molecular level. Because of this, computational tools are now extensively used also in the field of rational 3D design and optimization of chemical entities starting from the structural information of the targets, which can be experimentally resolved or can also be obtained with other computer-based techniques. In this work, we revised the state-of-the-art computer-aided drug design methods, focusing on their application in different scenarios of pharmaceutical and biological interest, not only highlighting their great potential and their benefits, but also discussing their actual limitations and eventual weaknesses. This work can be considered a brief overview of computational methods for drug discovery.
Collapse
Affiliation(s)
- Davide Bassani
- Pharmaceutical Research & Early Development, Roche Innovation Center Basel, F. Hoffmann-La Roche Ltd., 4070 Basel, Switzerland
- Molecular Modeling Section (MMS), Department of Pharmaceutical and Pharmacological Sciences, University of Padova, Via Marzolo 5, 35131 Padova, Italy
| | - Stefano Moro
- Molecular Modeling Section (MMS), Department of Pharmaceutical and Pharmacological Sciences, University of Padova, Via Marzolo 5, 35131 Padova, Italy
| |
Collapse
|
10
|
Amsler J, Plessow PN, Studt F, Bučko T. Anharmonic Correction to Free Energy Barriers from DFT-Based Molecular Dynamics Using Constrained Thermodynamic Integration. J Chem Theory Comput 2023; 19:2455-2468. [PMID: 37043693 DOI: 10.1021/acs.jctc.3c00169] [Citation(s) in RCA: 2] [Impact Index Per Article: 2.0] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 04/14/2023]
Abstract
For the calculation of anharmonic contributions to free energy barriers, constrained thermodynamic λ-path integration (λ-TI) from a harmonic reference force field to density functional theory is presented as an alternative to the established Blue Moon ensemble method (ξ-TI), in which free energy gradients along the reaction coordinate ξ are integrated. With good agreement in all cases, the λ-TI method is benchmarked against the ξ-TI method for several reactions, including the internal CH3 group rotation in ethane, a nucleophilic substitution of CH3Cl, a retro-Diels-Alder reaction, and a proton transfer in zeolite H-SSZ-13. An advantage of λ-TI is that one can use virtually any reference state to compute anharmonic contributions to reaction free energies or free energy barriers. This is particularly relevant for catalysis, where it is now possible to compute anharmonic corrections to the free energy of a transition state relative to any reference, for example, the most stable state of the active site and the reactants in the gas phase. This is in contrast to ξ-TI, where free energy barriers can only be computed relative to an initial state with all reactants coadsorbed. Finally, the Bennett acceptance ratio method combined with λ-TI is demonstrated to reduce the number of required integration grid points with tolerable accuracy, favoring thus λ-TI over ξ-TI in terms of computational efficiency.
Collapse
Affiliation(s)
- Jonas Amsler
- Institute of Catalysis Research and Technology, Karlsruhe Institute of Technology, Hermann-von-Helmholtz-Platz 1, 76344 Eggenstein-Leopoldshafen, Germany
| | - Philipp N Plessow
- Institute of Catalysis Research and Technology, Karlsruhe Institute of Technology, Hermann-von-Helmholtz-Platz 1, 76344 Eggenstein-Leopoldshafen, Germany
| | - Felix Studt
- Institute of Catalysis Research and Technology, Karlsruhe Institute of Technology, Hermann-von-Helmholtz-Platz 1, 76344 Eggenstein-Leopoldshafen, Germany
- Institute for Chemical Technology and Polymer Chemistry, Karlsruhe Institute of Technology, Kaiserstr. 12, 76131 Karlsruhe, Germany
| | - Tomáš Bučko
- Department of Physical and Theoretical Chemistry, Faculty of Natural Sciences, Comenius University in Bratislava, Ilkovičova 6, SK-84215 Bratislava, Slovakia
- Institute of Inorganic Chemistry, Slovak Academy of Sciences, Dúbravská cesta 9, SK-84236 Bratislava, Slovakia
| |
Collapse
|
11
|
Liu X, Zheng L, Qin C, Zhang JZH, Sun Z. Comprehensive evaluation of end-point free energy techniques in carboxylated-pillar[6]arene host-guest binding: I. Standard procedure. J Comput Aided Mol Des 2022; 36:735-752. [PMID: 36136209 DOI: 10.1007/s10822-022-00475-0] [Citation(s) in RCA: 0] [Impact Index Per Article: 0] [Reference Citation Analysis] [Abstract] [Key Words] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Received: 07/06/2022] [Accepted: 09/06/2022] [Indexed: 10/14/2022]
Abstract
Despite the massive application of end-point free energy methods in protein-ligand and protein-protein interactions, computational understandings about their performance in relatively simple and prototypical host-guest systems are limited. In this work, we present a comprehensive benchmark calculation with standard end-point free energy techniques in a recent host-guest dataset containing 13 host-guest pairs involving the carboxylated-pillar[6]arene host. We first assess the charge schemes for solutes by comparing the charge-produced electrostatics with many ab initio references, in order to obtain a preliminary albeit detailed view of the charge quality. Then, we focus on four modelling details of end-point free energy calculations, including the docking procedure for the generation of initial condition, the charge scheme for host and guest molecules, the water model used in explicit-solvent sampling, and the end-point methods for free energy estimation. The binding thermodynamics obtained with different modelling schemes are compared with experimental references, and some practical guidelines on maximizing the performance of end-point methods in practical host-guest systems are summarized. Further, we compare our simulation outcome with predictions in the grand challenge and discuss further developments to improve the prediction quality of end-point free energy methods. Overall, unlike the widely acknowledged applicability in protein-ligand binding, the standard end-point calculations cannot produce useful outcomes in host-guest binding and thus are not recommended unless alterations are performed.
Collapse
Affiliation(s)
- Xiao Liu
- School of Mathematics, Physics and Statistics, Shanghai University of Engineering Science, Shanghai, 201620, China.
| | - Lei Zheng
- NYU-ECNU Center for Computational Chemistry at NYU Shanghai, Shanghai, 200062, China
| | - Chu Qin
- School of Mathematics, Physics and Statistics, Shanghai University of Engineering Science, Shanghai, 201620, China
| | - John Z H Zhang
- NYU-ECNU Center for Computational Chemistry at NYU Shanghai, Shanghai, 200062, China.,School of Chemistry and Molecular Engineering, East China Normal University, Shanghai, 200062, China.,Shenzhen Institute of Advanced Technology, Chinese Academy of Sciences, Shenzhen, Guangdong, China.,Department of Chemistry, New York University, New York, NY, 10003, USA
| | - Zhaoxi Sun
- College of Chemistry and Molecular Engineering, Peking University, Beijing, 100871, China.
| |
Collapse
|
12
|
Akkus E, Tayfuroglu O, Yildiz M, Kocak A. Accurate Binding Free Energy Method from End-State MD Simulations. J Chem Inf Model 2022; 62:4095-4106. [PMID: 35972783 PMCID: PMC9472276 DOI: 10.1021/acs.jcim.2c00601] [Citation(s) in RCA: 1] [Impact Index Per Article: 0.5] [Reference Citation Analysis] [Abstract] [Track Full Text] [Download PDF] [Figures] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 01/11/2023]
Abstract
![]()
Herein, we introduce a new strategy to estimate binding
free energies
using end-state molecular dynamics simulation trajectories. The method
is adopted from linear interaction energy (LIE) and ANI-2x neural
network potentials (machine learning) for the atomic simulation environment
(ASE). It predicts the single-point interaction energies between ligand–protein
and ligand–solvent pairs at the accuracy of the wb97x/6-31G*
level for the conformational space that is sampled by molecular dynamics
(MD) simulations. Our results on 54 protein–ligand complexes
show that the method can be accurate and have a correlation of R = 0.87–0.88 to the experimental binding free energies,
outperforming current end-state methods with reduced computational
cost. The method also allows us to compare BFEs of ligands with different
scaffolds. The code is available free of charge (documentation and
test files) at https://github.com/otayfuroglu/deepQM.
Collapse
Affiliation(s)
- Ebru Akkus
- Department of Bioengineering, Gebze Technical University, 41400 Gebze, Kocaeli, Turkey
| | - Omer Tayfuroglu
- Department of Chemistry, Gebze Technical University, 41400 Gebze, Kocaeli, Turkey
| | - Muslum Yildiz
- Department of Molecular Biology and Genetics, Gebze Technical University, 41400 Gebze, Kocaeli, Turkey
| | - Abdulkadir Kocak
- Department of Chemistry, Gebze Technical University, 41400 Gebze, Kocaeli, Turkey
| |
Collapse
|
13
|
Waibl F, Kraml J, Fernández-Quintero ML, Loeffler JR, Liedl KR. Explicit solvation thermodynamics in ionic solution: extending grid inhomogeneous solvation theory to solvation free energy of salt-water mixtures. J Comput Aided Mol Des 2022; 36:101-116. [PMID: 35031880 PMCID: PMC8907097 DOI: 10.1007/s10822-021-00429-y] [Citation(s) in RCA: 6] [Impact Index Per Article: 3.0] [Reference Citation Analysis] [Abstract] [Key Words] [MESH Headings] [Grants] [Track Full Text] [Download PDF] [Figures] [Journal Information] [Subscribe] [Scholar Register] [Received: 09/14/2021] [Accepted: 10/28/2021] [Indexed: 12/03/2022]
Abstract
Hydration thermodynamics play a fundamental role in fields ranging from the pharmaceutical industry to environmental research. Numerous methods exist to predict solvation thermodynamics of compounds ranging from small molecules to large biomolecules. Arguably the most precise methods are those based on molecular dynamics (MD) simulations in explicit solvent. One theory that has seen increased use is inhomogeneous solvation theory (IST). However, while many applications require accurate description of salt-water mixtures, no implementation of IST is currently able to estimate solvation properties involving more than one solvent species. Here, we present an extension to grid inhomogeneous solvation theory (GIST) that can take salt contributions into account. At the example of carbazole in 1 M NaCl solution, we compute the solvation energy as well as first and second order entropies. While the effect of the first order ion entropy is small, both the water-water and water-ion entropies contribute strongly. We show that the water-ion entropies are efficiently approximated using the Kirkwood superposition approximation. However, this approach cannot be applied to the water-water entropy. Furthermore, we test the quantitative validity of our method by computing salting-out coefficients and comparing them to experimental data. We find a good correlation to experimental salting-out constants, while the absolute values are overpredicted due to the approximate second order entropy. Since ions are frequently used in MD, either to neutralize the system or as a part of the investigated process, our method greatly extends the applicability of GIST. The use-cases range from biopharmaceuticals, where many assays require high salt concentrations, to environmental research, where solubility in sea water is important to model the fate of organic substances.
Collapse
Affiliation(s)
- Franz Waibl
- Department of General, Inorganic, and Theoretical Chemistry, University of Innsbruck, Innrain 80/82, 6020, Innsbruck, Austria
| | - Johannes Kraml
- Department of General, Inorganic, and Theoretical Chemistry, University of Innsbruck, Innrain 80/82, 6020, Innsbruck, Austria
| | - Monica L Fernández-Quintero
- Department of General, Inorganic, and Theoretical Chemistry, University of Innsbruck, Innrain 80/82, 6020, Innsbruck, Austria
| | - Johannes R Loeffler
- Department of General, Inorganic, and Theoretical Chemistry, University of Innsbruck, Innrain 80/82, 6020, Innsbruck, Austria
| | - Klaus R Liedl
- Department of General, Inorganic, and Theoretical Chemistry, University of Innsbruck, Innrain 80/82, 6020, Innsbruck, Austria.
| |
Collapse
|
14
|
Sun Z, Kalhor P, Xu Y, Liu J. Extensive numerical tests of leapfrog integrator in middle thermostat scheme in molecular simulations. CHINESE J CHEM PHYS 2021. [DOI: 10.1063/1674-0068/cjcp2111242] [Citation(s) in RCA: 2] [Impact Index Per Article: 0.7] [Reference Citation Analysis] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 12/19/2022]
Affiliation(s)
- Zhaoxi Sun
- Beijing National Laboratory for Molecular Sciences, College of Chemistry and Molecular Engineering, Institute of Theoretical and Computational Chemistry, Peking University, Beijing 100871, China
| | - Payam Kalhor
- Beijing National Laboratory for Molecular Sciences, College of Chemistry and Molecular Engineering, Institute of Theoretical and Computational Chemistry, Peking University, Beijing 100871, China
| | - Yang Xu
- Beijing National Laboratory for Molecular Sciences, College of Chemistry and Molecular Engineering, Institute of Theoretical and Computational Chemistry, Peking University, Beijing 100871, China
| | - Jian Liu
- Beijing National Laboratory for Molecular Sciences, College of Chemistry and Molecular Engineering, Institute of Theoretical and Computational Chemistry, Peking University, Beijing 100871, China
| |
Collapse
|
15
|
Petrov D. Perturbation Free-Energy Toolkit: An Automated Alchemical Topology Builder. J Chem Inf Model 2021; 61:4382-4390. [PMID: 34415755 PMCID: PMC8479811 DOI: 10.1021/acs.jcim.1c00428] [Citation(s) in RCA: 10] [Impact Index Per Article: 3.3] [Reference Citation Analysis] [Abstract] [MESH Headings] [Grants] [Track Full Text] [Download PDF] [Figures] [Journal Information] [Subscribe] [Scholar Register] [Received: 04/15/2021] [Indexed: 11/30/2022]
Abstract
Free-energy calculations play an important role in the application of computational chemistry to a range of fields, including protein biochemistry, rational drug design, or materials science. Importantly, the free-energy difference is directly related to experimentally measurable quantities such as partition and adsorption coefficients, water activity, and binding affinities. Among several techniques aimed at predicting free-energy differences, perturbation approaches, involving the alchemical transformation of one molecule into another through intermediate states, stand out as rigorous methods based on statistical mechanics. However, despite the importance of free-energy calculations, the applicability of the perturbation approaches is still largely impeded by a number of challenges, including the definition of the perturbation path, i.e., alchemical changes leading to the transformation of one molecule to the other. To address this, an automatic perturbation topology builder based on a graph-matching algorithm is developed, which can identify the maximum common substructure (MCS) of two or multiple molecules and provide the perturbation topologies suitable for free-energy calculations using the GROMOS and the GROMACS simulation packages. Various MCS search options are presented leading to alternative definitions of the perturbation pathway. Moreover, perturbation topologies generated using the default multistate MCS search are used to calculate the changes in free energy between lysine and its two post-translational modifications, 3-methyllysine and acetyllysine. The pairwise free-energy calculations performed on this test system led to a cycle closure of 0.5 ± 0.3 and 0.2 ± 0.2 kJ mol-1, with GROMOS and GROMACS simulation packages, respectively. The same relative free energies between the three states are obtained by employing the enveloping distribution sampling (EDS) approach when compared to the pairwise perturbations. Importantly, this toolkit is made available online as an open-source Python package (https://github.com/drazen-petrov/SMArt).
Collapse
Affiliation(s)
- Drazen Petrov
- Department of Material Sciences
and Process Engineering, Institute of Molecular Modeling and Simulation, University of Natural Resources and Life Sciences
Vienna, Muthgasse 18, A-1190 Vienna, Austria
| |
Collapse
|
16
|
Sun Z, Liu Z. BAR‐Based Multi‐Dimensional Nonequilibrium Pulling for Indirect Construction of QM/MM Free Energy Landscapes: Varying the QM Region. ADVANCED THEORY AND SIMULATIONS 2021. [DOI: 10.1002/adts.202100185] [Citation(s) in RCA: 8] [Impact Index Per Article: 2.7] [Reference Citation Analysis] [Track Full Text] [Subscribe] [Scholar Register] [Indexed: 12/26/2022]
Affiliation(s)
- Zhaoxi Sun
- Beijing National Laboratory for Molecular Sciences, College of Chemistry and Molecular Engineering Peking University Beijing 100871 China
| | - Zhirong Liu
- Beijing National Laboratory for Molecular Sciences, College of Chemistry and Molecular Engineering Peking University Beijing 100871 China
| |
Collapse
|
17
|
King E, Aitchison E, Li H, Luo R. Recent Developments in Free Energy Calculations for Drug Discovery. Front Mol Biosci 2021; 8:712085. [PMID: 34458321 PMCID: PMC8387144 DOI: 10.3389/fmolb.2021.712085] [Citation(s) in RCA: 49] [Impact Index Per Article: 16.3] [Reference Citation Analysis] [Abstract] [Key Words] [Grants] [Track Full Text] [Download PDF] [Figures] [Journal Information] [Subscribe] [Scholar Register] [Received: 05/19/2021] [Accepted: 07/27/2021] [Indexed: 01/11/2023] Open
Abstract
The grand challenge in structure-based drug design is achieving accurate prediction of binding free energies. Molecular dynamics (MD) simulations enable modeling of conformational changes critical to the binding process, leading to calculation of thermodynamic quantities involved in estimation of binding affinities. With recent advancements in computing capability and predictive accuracy, MD based virtual screening has progressed from the domain of theoretical attempts to real application in drug development. Approaches including the Molecular Mechanics Poisson Boltzmann Surface Area (MM-PBSA), Linear Interaction Energy (LIE), and alchemical methods have been broadly applied to model molecular recognition for drug discovery and lead optimization. Here we review the varied methodology of these approaches, developments enhancing simulation efficiency and reliability, remaining challenges hindering predictive performance, and applications to problems in the fields of medicine and biochemistry.
Collapse
Affiliation(s)
- Edward King
- Department of Molecular Biology and Biochemistry, University of California, Irvine, CA, United States
| | - Erick Aitchison
- Department of Molecular Biology and Biochemistry, University of California, Irvine, CA, United States
| | - Han Li
- Department of Chemical and Biomolecular Engineering, University of California, Irvine, CA, United States
| | - Ray Luo
- Department of Molecular Biology and Biochemistry, University of California, Irvine, CA, United States
- Department of Chemical and Biomolecular Engineering, University of California, Irvine, CA, United States
- Department of Materials Science and Engineering, University of California, Irvine, CA, United States
- Department of Biomedical Engineering, University of California, Irvine, CA, United States
| |
Collapse
|
18
|
Kashefolgheta S, Wang S, Acree WE, Hünenberger PH. Evaluation of nine condensed-phase force fields of the GROMOS, CHARMM, OPLS, AMBER, and OpenFF families against experimental cross-solvation free energies. Phys Chem Chem Phys 2021; 23:13055-13074. [PMID: 34105547 PMCID: PMC8207520 DOI: 10.1039/d1cp00215e] [Citation(s) in RCA: 5] [Impact Index Per Article: 1.7] [Reference Citation Analysis] [Abstract] [Grants] [Track Full Text] [Figures] [Journal Information] [Subscribe] [Scholar Register] [Received: 01/15/2021] [Accepted: 04/28/2021] [Indexed: 12/02/2022]
Abstract
Experimental solvation free energies are nowadays commonly included as target properties in the validation of condensed-phase force fields, sometimes even in their calibration. In a previous article [Kashefolgheta et al., J. Chem. Theory. Comput., 2020, 16, 7556-7580], we showed how the involved comparison between experimental and simulation results could be made more systematic by considering a full matrix of cross-solvation free energies . For a set of N molecules that are all in the liquid state under ambient conditions, such a matrix encompasses N×N entries for considering each of the N molecules either as solute (A) or as solvent (B). In the quoted study, a cross-solvation matrix of 25 × 25 experimental value was introduced, considering 25 small molecules representative for alkanes, chloroalkanes, ethers, ketones, esters, alcohols, amines, and amides. This experimental data was used to compare the relative accuracies of four popular condensed-phase force fields, namely GROMOS-2016H66, OPLS-AA, AMBER-GAFF, and CHARMM-CGenFF. In the present work, the comparison is extended to five additional force fields, namely GROMOS-54A7, GROMOS-ATB, OPLS-LBCC, AMBER-GAFF2, and OpenFF. Considering these nine force fields, the correlation coefficients between experimental values and simulation results range from 0.76 to 0.88, the root-mean-square errors (RMSEs) from 2.9 to 4.8 kJ mol-1, and average errors (AVEEs) from -1.5 to +1.0 kJ mol-1. In terms of RMSEs, GROMOS-2016H66 and OPLS-AA present the best accuracy (2.9 kJ mol-1), followed by OPLS-LBCC, AMBER-GAFF2, AMBER-GAFF, and OpenFF (3.3 to 3.6 kJ mol-1), and then by GROMOS-54A7, CHARM-CGenFF, and GROMOS-ATB (4.0 to 4.8 kJ mol-1). These differences are statistically significant but not very pronounced, and are distributed rather heterogeneously over the set of compounds within the different force fields.
Collapse
Affiliation(s)
- Sadra Kashefolgheta
- Laboratorium für Physikalische Chemie, ETH Zürich, ETH-Hönggerberg, HCICH-8093 ZürichSwitzerland+41 44 632 55 03
| | - Shuzhe Wang
- Laboratorium für Physikalische Chemie, ETH Zürich, ETH-Hönggerberg, HCICH-8093 ZürichSwitzerland+41 44 632 55 03
| | - William E. Acree
- Department of Chemistry, University of North Texas1155 Union Circle Drive #305070DentonTexas 76203USA
| | - Philippe H. Hünenberger
- Laboratorium für Physikalische Chemie, ETH Zürich, ETH-Hönggerberg, HCICH-8093 ZürichSwitzerland+41 44 632 55 03
| |
Collapse
|
19
|
Huai Z, Yang H, Sun Z. Binding thermodynamics and interaction patterns of human purine nucleoside phosphorylase-inhibitor complexes from extensive free energy calculations. J Comput Aided Mol Des 2021; 35:643-656. [PMID: 33759016 DOI: 10.1007/s10822-021-00382-w] [Citation(s) in RCA: 6] [Impact Index Per Article: 2.0] [Reference Citation Analysis] [Abstract] [Key Words] [Journal Information] [Subscribe] [Scholar Register] [Received: 10/12/2020] [Accepted: 03/13/2021] [Indexed: 11/29/2022]
Abstract
Human purine nucleoside phosphorylase (hPNP) plays a significant role in the catabolism of deoxyguanosine. The trimeric protein is an important target in the treatment of T-cell cancers and autoimmune disorders. Experimental studies on the inhibition of the hPNP observe that the first ligand bound to one of three subunits effectively inhibits the protein, while the binding of more ligands to the subsequent sites shows negative cooperativities. In this work, we performed extensive end-point and alchemical free energy calculations to determine the binding thermodynamics of the trimeric protein-ligand system. 13 Immucillin inhibitors with experimental results are under calculation. Two widely accepted charge schemes for small molecules including AM1-BCC and RESP are adopted for ligands. The results of RESP are in better agreement with the experimental reference. Further investigations of the interaction networks in the protein-ligand complexes reveal that several residues play significant roles in stabilizing the complex structure. The most commonly observed ones include PHE200, GLU201, MET219, and ASN243. The conformations of the protein in different protein-ligand complexes are observed to be similar. We expect these insights to aid the development of potent drugs targeting hPNP.
Collapse
Affiliation(s)
- Zhe Huai
- State Key Laboratory of Precision Spectroscopy, School of Physics and Electronic Science, East China Normal University, Shanghai, 200062, China
| | - Huaiyu Yang
- College of Engineering, Hebei Normal University, Shijiazhuang, 050024, China
| | - Zhaoxi Sun
- State Key Laboratory of Precision Spectroscopy, School of Physics and Electronic Science, East China Normal University, Shanghai, 200062, China.
| |
Collapse
|
20
|
Bieniek M, Bhati AP, Wan S, Coveney PV. TIES 20: Relative Binding Free Energy with a Flexible Superimposition Algorithm and Partial Ring Morphing. J Chem Theory Comput 2021; 17:1250-1265. [PMID: 33486956 PMCID: PMC7876800 DOI: 10.1021/acs.jctc.0c01179] [Citation(s) in RCA: 14] [Impact Index Per Article: 4.7] [Reference Citation Analysis] [Abstract] [Grants] [Track Full Text] [Download PDF] [Figures] [Journal Information] [Subscribe] [Scholar Register] [Received: 11/10/2020] [Indexed: 12/14/2022]
Abstract
The TIES (Thermodynamic Integration with Enhanced Sampling) protocol is a formally exact alchemical approach in computational chemistry to the calculation of relative binding free energies. The validity of TIES relies on the correctness of matching atoms across compared pairs of ligands, laying the foundation for the transformation along an alchemical pathway. We implement a flexible topology superimposition algorithm which uses an exhaustive joint-traversal for computing the largest common component(s). The algorithm is employed to enable matching and morphing of partial rings in the TIES protocol along with a validation study using 55 transformations and five different proteins from our previous work. We find that TIES 20 with the RESP charge system, using the new superimposition algorithm, reproduces the previous results with mean unsigned error of 0.75 kcal/mol with respect to the experimental data. Enabling the morphing of partial rings decreases the size of the alchemical region in the dual-topology transformations resulting in a significant improvement in the prediction precision. We find that increasing the ensemble size from 5 to 20 replicas per λ window only has a minimal impact on the accuracy. However, the non-normal nature of the relative free energy distributions underscores the importance of ensemble simulation. We further compare the results with the AM1-BCC charge system and show that it improves agreement with the experimental data by slightly over 10%. This improvement is partly due to AM1-BCC affecting only the charges of the atoms local to the mutation, which translates to even fewer morphed atoms, consequently reducing issues with sampling and therefore ensemble averaging. TIES 20, in conjunction with the enablement of ring morphing, reduces the size of the alchemical region and significantly improves the precision of the predicted free energies.
Collapse
Affiliation(s)
- Mateusz
K. Bieniek
- Centre for Computational Science, Department
of Chemistry, University College London, 20 Gordon Street, London WC1H 0AJ, United Kingdom
| | - Agastya P. Bhati
- Centre for Computational Science, Department
of Chemistry, University College London, 20 Gordon Street, London WC1H 0AJ, United Kingdom
| | - Shunzhou Wan
- Centre for Computational Science, Department
of Chemistry, University College London, 20 Gordon Street, London WC1H 0AJ, United Kingdom
| | - Peter V. Coveney
- Centre for Computational Science, Department
of Chemistry, University College London, 20 Gordon Street, London WC1H 0AJ, United Kingdom
| |
Collapse
|
21
|
Amsler J, Plessow PN, Studt F, Bučko T. Anharmonic Correction to Adsorption Free Energy from DFT-Based MD Using Thermodynamic Integration. J Chem Theory Comput 2021; 17:1155-1169. [PMID: 33482059 DOI: 10.1021/acs.jctc.0c01022] [Citation(s) in RCA: 16] [Impact Index Per Article: 5.3] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 11/29/2022]
Abstract
Adsorption processes are often governed by weak interactions for which the estimation of entropy contributions by means of the harmonic approximation is prone to be inaccurate. Thermodynamic integration (TI) from the harmonic to the fully interacting system (λ-path integration) can be used to compute anharmonic corrections. Here, we combine TI with (curvilinear) internal coordinates in periodic systems to make the formalism available in computational studies. Our implementation of ab initio molecular dynamics in VASP is independent of the reaction path and can be thus applied to study adsorption processes relative to the gas phase and does hence provide a useful tool for computational catalysis. We discuss the application of the approach on three model systems for which exact semianalytical solutions exist and illustrate and quantify the importance of anharmonic vibrations, hindered rotations, and hindered translations (dissociation). Eventually, we apply the method to study the adsorption of small adsorbates in a zeolite (H-SSZ-13).
Collapse
Affiliation(s)
- Jonas Amsler
- Institute of Catalysis Research and Technology, Karlsruhe Institute of Technology, Hermann-von-Helmholtz-Platz 1, 76344 Eggenstein-Leopoldshafen, Germany
| | - Philipp N Plessow
- Institute of Catalysis Research and Technology, Karlsruhe Institute of Technology, Hermann-von-Helmholtz-Platz 1, 76344 Eggenstein-Leopoldshafen, Germany
| | - Felix Studt
- Institute of Catalysis Research and Technology, Karlsruhe Institute of Technology, Hermann-von-Helmholtz-Platz 1, 76344 Eggenstein-Leopoldshafen, Germany.,Institute for Chemical Technology and Polymer Chemistry, Karlsruhe Institute of Technology, Kaiserstr. 12, 76131 Karlsruhe, Germany
| | - Tomáš Bučko
- Department of Physical and Theoretical Chemistry, Faculty of Natural Sciences, Comenius University in Bratislava, Ilkovičova 6, SK-84215 Bratislava, Slovakia.,Institute of Inorganic Chemistry, Slovak Academy of Sciences, Dúbravská cesta 9, SK-84236 Bratislava, Slovakia
| |
Collapse
|
22
|
Huai Z, Shen Z, Sun Z. Binding Thermodynamics and Interaction Patterns of Inhibitor-Major Urinary Protein-I Binding from Extensive Free-Energy Calculations: Benchmarking AMBER Force Fields. J Chem Inf Model 2020; 61:284-297. [PMID: 33307679 DOI: 10.1021/acs.jcim.0c01217] [Citation(s) in RCA: 15] [Impact Index Per Article: 3.8] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 12/23/2022]
Abstract
Mouse major urinary protein (MUP) plays a key role in the pheromone communication system. The one-end-closed β-barrel of MUP-I forms a small, deep, and hydrophobic central cavity, which could accommodate structurally diverse ligands. Previous computational studies employed old protein force fields and short simulation times to determine the binding thermodynamics or investigated only a small number of structurally similar ligands, which resulted in sampled regions far from the experimental structure, nonconverged sampling outcomes, and limited understanding of the possible interaction patterns that the cavity could produce. In this work, extensive end-point and alchemical free-energy calculations with advanced protein force fields were performed to determine the binding thermodynamics of a series of MUP-inhibitor systems and investigate the inter- and intramolecular interaction patterns. Three series of inhibitors with a total of 14 ligands were simulated. We independently simulated the MUP-inhibitor complexes under two advanced AMBER force fields. Our benchmark test showed that the advanced AMBER force fields including AMBER19SB and AMBER14SB provided better descriptions of the system, and the backbone root-mean-square deviation (RMSD) was significantly lowered compared with previous computational studies with old protein force fields. Surprisingly, although the latest AMBER force field AMBER19SB provided better descriptions of various observables, it neither improved the binding thermodynamics nor lowered the backbone RMSD compared with the previously proposed and widely used AMBER14SB. The older but widely used AMBER14SB actually achieved better performance in the prediction of binding affinities from the alchemical and end-point free-energy calculations. We further analyzed the protein-ligand interaction networks to identify important residues stabilizing the bound structure. Six residues including PHE38, LEU40, PHE90, ALA103, LEU105, and TYR120 were found to contribute the most significant part of protein-ligand interactions, and 10 residues were found to provide favorable interactions stabilizing the bound state. The two AMBER force fields gave extremely similar interaction networks, and the secondary structures also showed similar behavior. Thus, the intra- and intermolecular interaction networks described with the two AMBER force fields are similar. Therefore, AMBER14SB could still be the default option in free-energy calculations to achieve highly accurate binding thermodynamics and interaction patterns.
Collapse
Affiliation(s)
- Zhe Huai
- State Key Laboratory of Precision Spectroscopy, School of Physics and Electronic Science, East China Normal University, Shanghai 200062, China
| | - Zhaoxi Shen
- Institute of Applied Physics and Materials Engineering, University of Macau, Avenida da Universidade, Taipa, Macau 999078, China
| | - Zhaoxi Sun
- State Key Laboratory of Precision Spectroscopy, School of Physics and Electronic Science, East China Normal University, Shanghai 200062, China
| |
Collapse
|
23
|
Lee TS, Allen BK, Giese TJ, Guo Z, Li P, Lin C, McGee TD, Pearlman DA, Radak BK, Tao Y, Tsai HC, Xu H, Sherman W, York DM. Alchemical Binding Free Energy Calculations in AMBER20: Advances and Best Practices for Drug Discovery. J Chem Inf Model 2020; 60:5595-5623. [PMID: 32936637 PMCID: PMC7686026 DOI: 10.1021/acs.jcim.0c00613] [Citation(s) in RCA: 175] [Impact Index Per Article: 43.8] [Reference Citation Analysis] [Abstract] [MESH Headings] [Grants] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 01/03/2023]
Abstract
Predicting protein-ligand binding affinities and the associated thermodynamics of biomolecular recognition is a primary objective of structure-based drug design. Alchemical free energy simulations offer a highly accurate and computationally efficient route to achieving this goal. While the AMBER molecular dynamics package has successfully been used for alchemical free energy simulations in academic research groups for decades, widespread impact in industrial drug discovery settings has been minimal because of the previous limitations within the AMBER alchemical code, coupled with challenges in system setup and postprocessing workflows. Through a close academia-industry collaboration we have addressed many of the previous limitations with an aim to improve accuracy, efficiency, and robustness of alchemical binding free energy simulations in industrial drug discovery applications. Here, we highlight some of the recent advances in AMBER20 with a focus on alchemical binding free energy (BFE) calculations, which are less computationally intensive than alternative binding free energy methods where full binding/unbinding paths are explored. In addition to scientific and technical advances in AMBER20, we also describe the essential practical aspects associated with running relative alchemical BFE calculations, along with recommendations for best practices, highlighting the importance not only of the alchemical simulation code but also the auxiliary functionalities and expertise required to obtain accurate and reliable results. This work is intended to provide a contemporary overview of the scientific, technical, and practical issues associated with running relative BFE simulations in AMBER20, with a focus on real-world drug discovery applications.
Collapse
Affiliation(s)
- Tai-Sung Lee
- Rutgers, the State University of New Jersey, Laboratory for Biomolecular Simulation Research, and Department of Chemistry and Chemical Biology, United States
| | - Bryce K. Allen
- Silicon Therapeutics, Boston, Massachusetts 02210, United States
| | - Timothy J. Giese
- Rutgers, the State University of New Jersey, Laboratory for Biomolecular Simulation Research, and Department of Chemistry and Chemical Biology, United States
| | - Zhenyu Guo
- Silicon Therapeutics, Boston, Massachusetts 02210, United States
| | - Pengfei Li
- Silicon Therapeutics, Boston, Massachusetts 02210, United States
| | - Charles Lin
- Silicon Therapeutics, Boston, Massachusetts 02210, United States
| | - T. Dwight McGee
- Silicon Therapeutics, Boston, Massachusetts 02210, United States
| | - David A. Pearlman
- QSimulate Incorporated, Cambridge, Massachusetts 02139, United States
| | - Brian K. Radak
- Silicon Therapeutics, Boston, Massachusetts 02210, United States
| | - Yujun Tao
- Rutgers, the State University of New Jersey, Laboratory for Biomolecular Simulation Research, and Department of Chemistry and Chemical Biology, United States
| | - Hsu-Chun Tsai
- Rutgers, the State University of New Jersey, Laboratory for Biomolecular Simulation Research, and Department of Chemistry and Chemical Biology, United States
| | - Huafeng Xu
- Silicon Therapeutics, Boston, Massachusetts 02210, United States
| | - Woody Sherman
- Silicon Therapeutics, Boston, Massachusetts 02210, United States
| | - Darrin M. York
- Rutgers, the State University of New Jersey, Laboratory for Biomolecular Simulation Research, and Department of Chemistry and Chemical Biology, United States
| |
Collapse
|
24
|
Monroe JI, Hatch HW, Mahynski NA, Shell MS, Shen VK. Extrapolation and interpolation strategies for efficiently estimating structural observables as a function of temperature and density. J Chem Phys 2020; 153:144101. [PMID: 33086808 DOI: 10.1063/5.0014282] [Citation(s) in RCA: 9] [Impact Index Per Article: 2.3] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 12/28/2022] Open
Abstract
Thermodynamic extrapolation has previously been used to predict arbitrary structural observables in molecular simulations at temperatures (or relative chemical potentials in open-system mixtures) different from those at which the simulation was performed. This greatly reduces the computational cost in mapping out phase and structural transitions. In this work, we explore the limitations and accuracy of thermodynamic extrapolation applied to water, where qualitative shifts from anomalous to simple-fluid-like behavior are manifested through shifts in the liquid structure that occur as a function of both temperature and density. We present formulas for extrapolating in volume for canonical ensembles and demonstrate that linear extrapolations of water's structural properties are only accurate over a limited density range. On the other hand, linear extrapolation in temperature can be accurate across the entire liquid state. We contrast these extrapolations with classical perturbation theory techniques, which are more conservative and slowly converging. Indeed, we show that such behavior is expected by demonstrating exact relationships between extrapolation of free energies and well-known techniques to predict free energy differences. An ideal gas in an external field is also studied to more clearly explain these results for a toy system with fully analytical solutions. We also present a recursive interpolation strategy for predicting arbitrary structural properties of molecular fluids over a predefined range of state conditions, demonstrating its success in mapping qualitative shifts in water structure with density.
Collapse
Affiliation(s)
- Jacob I Monroe
- National Institute of Standards and Technology, Gaithersburg, Maryland 20899, USA
| | - Harold W Hatch
- National Institute of Standards and Technology, Gaithersburg, Maryland 20899, USA
| | - Nathan A Mahynski
- National Institute of Standards and Technology, Gaithersburg, Maryland 20899, USA
| | - M Scott Shell
- University of California - Santa Barbara, Santa Barbara, California 93106, USA
| | - Vincent K Shen
- National Institute of Standards and Technology, Gaithersburg, Maryland 20899, USA
| |
Collapse
|
25
|
SAMPL7 TrimerTrip host-guest binding affinities from extensive alchemical and end-point free energy calculations. J Comput Aided Mol Des 2020; 35:117-129. [PMID: 33037549 DOI: 10.1007/s10822-020-00351-9] [Citation(s) in RCA: 14] [Impact Index Per Article: 3.5] [Reference Citation Analysis] [Abstract] [Key Words] [Journal Information] [Subscribe] [Scholar Register] [Received: 07/29/2020] [Accepted: 09/30/2020] [Indexed: 12/12/2022]
Abstract
The prediction of host-guest binding affinities with computational modelling is still a challenging task. In the 7th statistical assessment of the modeling of proteins and ligands (SAMPL) challenge, a new host named TrimerTrip was synthesized and the thermodynamic parameters of 16 structurally diverse guests binding to the host were characterized. In the TrimerTrip-guest challenge, only structures of the host and the guests are provided, which indicates that the predictions of both the binding poses and the binding affinities are under assessment. In this work, starting from the binding poses obtained from our previous enhanced sampling simulations in the configurational space, we perform extensive alchemical and end-point free energy calculations to calculate the host-guest binding affinities retrospectively. The alchemical predictions with two widely accepted charge schemes (i.e. AM1-BCC and RESP) are in good agreement with the experimental reference, while the end-point estimates perform poorly in reproducing the experimental binding affinities. Aside from the absolute value of the binding affinity, the rank of binding free energies is also crucial in drug design. Surprisingly, the end-point MM/PBSA method seems very powerful in reproducing the experimental rank of binding affinities. Although the length of our simulations is long and the intermediate spacing is dense, the convergence behavior is not very good, which may arise from the flexibility of the host molecule. Enhanced sampling techniques in the configurational space may be required to obtain fully converged sampling. Further, as the length of sampling in alchemical free energy calculations already achieves several hundred ns, performing direct simulations of the binding/unbinding event in the physical space could be more useful and insightful. More details about the binding pathway and mechanism could be obtained in this way. The nonequilibrium method could also be a nice choice if one insists to use the alchemical method, as the intermediate sampling is avoided to some extent.
Collapse
|
26
|
König G, Glaser N, Schroeder B, Kubincová A, Hünenberger PH, Riniker S. An Alternative to Conventional λ-Intermediate States in Alchemical Free Energy Calculations: λ-Enveloping Distribution Sampling. J Chem Inf Model 2020; 60:5407-5423. [PMID: 32794763 DOI: 10.1021/acs.jcim.0c00520] [Citation(s) in RCA: 18] [Impact Index Per Article: 4.5] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 01/01/2023]
Abstract
Alchemical free energy calculations typically rely on intermediate states to bridge between the relevant phase spaces of the two end states. These intermediate states are usually created by mixing the energies or parameters of the end states according to a coupling parameter λ. The choice of the procedure has a strong impact on the efficiency of the calculation, as it affects both the encountered energy barriers and the phase space overlap between the states. The present work builds on the connection between the minimum variance pathway (MVP) and enveloping distribution sampling (EDS). It is shown that both methods can be regarded as special cases of a common scheme referred to as λ-EDS, which can also reproduce the behavior of conventional λ-intermediate states. A particularly attractive feature of λ-EDS is its ability to emulate the use of soft core potentials (SCP) while avoiding the associated computational overhead when applying efficient free energy estimators such as the multistate Bennett's acceptance ratio (MBAR). The method is illustrated for both relative and absolute free energy calculations considering five benchmark systems. The first two systems (charge inversion and cavity creation in a dipolar solvent) demonstrate the use of λ-EDS as an alternative coupling scheme in the context of thermodynamic integration (TI). The three other systems (change of bond length, change of dihedral angles, and cavity creation in water) investigate the efficiency and optimal choice of parameters in the context of free energy perturbation (FEP) and Bennett's acceptance ratio (BAR). It is shown that λ-EDS allows larger steps along the alchemical pathway than conventional intermediate states.
Collapse
Affiliation(s)
- Gerhard König
- Laboratory of Physical Chemistry, Department of Chemistry and Applied Biosciences, ETH Zürich, Vladimir-Prelog-Weg 2, 8093 Zürich, Switzerland
| | - Nina Glaser
- Laboratory of Physical Chemistry, Department of Chemistry and Applied Biosciences, ETH Zürich, Vladimir-Prelog-Weg 2, 8093 Zürich, Switzerland
| | - Benjamin Schroeder
- Laboratory of Physical Chemistry, Department of Chemistry and Applied Biosciences, ETH Zürich, Vladimir-Prelog-Weg 2, 8093 Zürich, Switzerland
| | - Alžbeta Kubincová
- Laboratory of Physical Chemistry, Department of Chemistry and Applied Biosciences, ETH Zürich, Vladimir-Prelog-Weg 2, 8093 Zürich, Switzerland
| | - Philippe H Hünenberger
- Laboratory of Physical Chemistry, Department of Chemistry and Applied Biosciences, ETH Zürich, Vladimir-Prelog-Weg 2, 8093 Zürich, Switzerland
| | - Sereina Riniker
- Laboratory of Physical Chemistry, Department of Chemistry and Applied Biosciences, ETH Zürich, Vladimir-Prelog-Weg 2, 8093 Zürich, Switzerland
| |
Collapse
|
27
|
Sun Z. SAMPL7 TrimerTrip host-guest binding poses and binding affinities from spherical-coordinates-biased simulations. J Comput Aided Mol Des 2020; 35:105-115. [PMID: 32776199 DOI: 10.1007/s10822-020-00335-9] [Citation(s) in RCA: 14] [Impact Index Per Article: 3.5] [Reference Citation Analysis] [Abstract] [Key Words] [Journal Information] [Subscribe] [Scholar Register] [Received: 04/04/2020] [Accepted: 08/04/2020] [Indexed: 12/21/2022]
Abstract
Host-guest binding remains a major challenge in modern computational modelling. The newest 7th statistical assessment of the modeling of proteins and ligands (SAMPL) challenge contains a new series of host-guest systems. The TrimerTrip host binds to 16 structurally diverse guests. Previously, we have successfully employed the spherical coordinates as the collective variables coupled with the enhanced sampling technique metadynamics to enhance the sampling of the binding/unbinding event, search for possible binding poses and calculate the binding affinities in all three host-guest binding cases of the 6th SAMPL challenge. In this work, we report a retrospective study on the TrimerTrip host-guest systems by employing the same protocol to investigate the TrimerTrip host in the SAMPL7 challenge. As no binding pose is provided by the SAMPL7 host, our simulations initiate from randomly selected configurations and are proceeded long enough to obtain converged free energy estimates and search for possible binding poses. The calculated binding affinities are in good agreement with the experimental reference, and the obtained binding poses serve as a nice starting point for end-point or alchemical free energy calculations. Note that as the work is performed after the close of the SAMPL7 challenge, we do not participate in the challenge and the results are not formally submitted to the SAMPL7 challenge.
Collapse
Affiliation(s)
- Zhaoxi Sun
- State Key Laboratory of Precision Spectroscopy, School of Chemistry and Molecular Engineering, East China Normal University, Shanghai, 200062, China.
| |
Collapse
|
28
|
Li Y, Nam K. Repulsive Soft-Core Potentials for Efficient Alchemical Free Energy Calculations. J Chem Theory Comput 2020; 16:4776-4789. [PMID: 32559374 DOI: 10.1021/acs.jctc.0c00163] [Citation(s) in RCA: 9] [Impact Index Per Article: 2.3] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 01/07/2023]
Abstract
In alchemical free energy (FE) simulations, annihilation and creation of atoms are generally achieved with the soft-core potential that shifts the interparticle separations. While this soft-core potential eliminates the numerical instability occurring near the two end states of the transformation, it makes the hybrid Hamiltonian vary nonlinearly with respect to the parameter λ, which interpolates between the Hamiltonians representing the two end states. This complicates FE estimation by Bennett acceptance ratio (BAR), free energy perturbation (FEP), and thermodynamic integration (TI) methods, thus reducing their calculation efficiency. In this work, we develop a new type of repulsive soft-core potential, called Gaussian soft-core (GSC) potential, with two parameters controlling its maximum and width. The main advantage of this potential is the linearity of the hybrid Hamiltonian with respect to λ, thus permitting the direct application of BAR, FEP, TI, and other variant FE methods. The accuracy and efficiency of the GSC potential are demonstrated by comparing the free energies of annihilation determined for 13 small molecules and an alchemical mutation of a protein side chain. In addition, in combination with a TI integrand (∂H/∂λ) estimation strategy, we show that GSC can considerably reduce the number of λ simulations compared to the commonly used separation-shifted soft-core potential.
Collapse
Affiliation(s)
- Yaozong Li
- Department of Chemistry, Umeå University, SE-901 87 Umeå, Sweden.,Department of Biochemistry, University of Zurich, Winterthurerstrasse 190, CH-8057 Zurich, Switzerland
| | - Kwangho Nam
- Department of Chemistry, Umeå University, SE-901 87 Umeå, Sweden.,Department of Chemistry and Biochemistry, University of Texas at Arlington, Arlington, Texas 76019-0065, United States
| |
Collapse
|
29
|
Zanetti-Polzi L, Daidone I, Amadei A. Fully Atomistic Multiscale Approach for p Ka Prediction. J Phys Chem B 2020; 124:4712-4722. [PMID: 32427481 DOI: 10.1021/acs.jpcb.0c01752] [Citation(s) in RCA: 9] [Impact Index Per Article: 2.3] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 01/29/2023]
Abstract
The ionization state of titratable amino acids strongly affects proteins structure and functioning in a large number of biological processes. It is therefore essential to be able to characterize the pKa of ionizable groups inside proteins and to understand its microscopic determinants in order to gain insights into many functional properties of proteins. A big effort has been devoted to the development of theoretical approaches for the prediction of deprotonation free energies, yet the accurate theoretical/computational calculation of pKa values is recognized as a current challenge. A methodology based on a hybrid quantum/classical approach is here proposed for the computation of deprotonation free energies. The method is applied to calculate the pKa of formic acid, methylammonium, and methanethiol, providing results in good agreement with the corresponding experimental estimates. The pKa is also calculated for aspartic acid and lysine as single residues in solution and for three aspartic/glutamic acids inside a well-characterized protein: hen egg white lysozyme. While for small molecules the method is able to deal with multiple protonation states of all titratable groups, this becomes computationally very expensive for proteins. The calculated pKa values for the single amino acids (except for the zwitterionic aspartic acid) and inside the protein display a systematic shift with respect to the experimental values that suggests that the fine balance between hydrophobic and polar interactions might be not accurately reproduced by the usual classical force-fields, thus affecting the computation of deprotonation free energies. The calculated pKa shifts inside the protein are in good agreement with the corresponding experimental ones (within 1 pKa unit), well reproducing the pKa changes due to the protein environment even in the case of large pKa shifts.
Collapse
Affiliation(s)
| | - Isabella Daidone
- Department of Physical and Chemical Sciences, University of L'Aquila, Via Vetoio, I-67010 L'Aquila, Italy
| | - Andrea Amadei
- Department of Chemical and Technological Sciences, University of Rome "Tor Vergata", Via della Ricerca Scientifica, I-00185 Rome, Italy
| |
Collapse
|
30
|
Rondina GG, Böhm MC, Müller-Plathe F. Predicting the Mobility Increase of Coarse-Grained Polymer Models from Excess Entropy Differences. J Chem Theory Comput 2020; 16:1431-1447. [DOI: 10.1021/acs.jctc.9b01088] [Citation(s) in RCA: 13] [Impact Index Per Article: 3.3] [Reference Citation Analysis] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 11/30/2022]
Affiliation(s)
- Gustavo G. Rondina
- Eduard-Zintl-Institut für Anorganische und Physikalische Chemie, Technische Universität Darmstadt, Alarich-Weiss-Straße 8, 64287 Darmstadt, Germany
| | - Michael C. Böhm
- Eduard-Zintl-Institut für Anorganische und Physikalische Chemie, Technische Universität Darmstadt, Alarich-Weiss-Straße 8, 64287 Darmstadt, Germany
| | - Florian Müller-Plathe
- Eduard-Zintl-Institut für Anorganische und Physikalische Chemie, Technische Universität Darmstadt, Alarich-Weiss-Straße 8, 64287 Darmstadt, Germany
| |
Collapse
|
31
|
Hahn DF, König G, Hünenberger PH. Overcoming Orthogonal Barriers in Alchemical Free Energy Calculations: On the Relative Merits of λ-Variations, λ-Extrapolations, and Biasing. J Chem Theory Comput 2020; 16:1630-1645. [DOI: 10.1021/acs.jctc.9b00853] [Citation(s) in RCA: 11] [Impact Index Per Article: 2.8] [Reference Citation Analysis] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 01/05/2023]
Affiliation(s)
- David F. Hahn
- Laboratory of Physical Chemistry, Department of Chemistry and Applied Biosciences, ETH Zürich, Vladimir-Prelog-Weg 2, 8093 Zürich, Switzerland
| | - Gerhard König
- Laboratory of Physical Chemistry, Department of Chemistry and Applied Biosciences, ETH Zürich, Vladimir-Prelog-Weg 2, 8093 Zürich, Switzerland
| | - Philippe H. Hünenberger
- Laboratory of Physical Chemistry, Department of Chemistry and Applied Biosciences, ETH Zürich, Vladimir-Prelog-Weg 2, 8093 Zürich, Switzerland
| |
Collapse
|
32
|
Sun Z, Wang X, Zhang JZH. Theoretical understanding of the thermodynamics and interactions in transcriptional regulator TtgR-ligand binding. Phys Chem Chem Phys 2019; 22:1511-1524. [PMID: 31872826 DOI: 10.1039/c9cp05980f] [Citation(s) in RCA: 10] [Impact Index Per Article: 2.0] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 02/06/2023]
Abstract
The transcriptional regulator TtgR belongs to the TetR family of transcriptional repressors. It depresses the transcription of the TtgABC operon and itself and thus regulates the extrusion of noxious chemicals with efflux pumps in bacterial cells. As the ligand-binding domain of TtgR is rather flexible, it can bind with a number of structurally diverse ligands, such as antibiotics, flavonoids and aromatic solvents. In the current work, we perform equilibrium and nonequilibrium alchemical free energy simulation to predict the binding affinities of a series of ligands targeting the TtgR protein and an agreement between the theoretical prediction and the experimental result is observed. End-point methods MM/PBSA and MM/GBSA are also employed for comparison. We further study the interaction maps and contacts between the protein and the ligand and identify important interactions in the protein-ligand binding cases. The dynamics fluctuation and secondary structures are also investigated. The current work sheds light on atomic and thermodynamic understanding of the TtgR-ligand interactions.
Collapse
Affiliation(s)
- Zhaoxi Sun
- Computational Biomedicine (IAS-5/INM-9), Forschungszentrum Jülich, Jülich 52425, Germany. and State Key Laboratory of Precision Spectroscopy, School of Chemistry and Molecular Engineering, East China Normal University, Shanghai 200062, China
| | - Xiaohui Wang
- State Key Laboratory of Precision Spectroscopy, School of Chemistry and Molecular Engineering, East China Normal University, Shanghai 200062, China and Institute of Computational Science, Università della Svizzera italiana (USI), Via Giuseppe Buffi 13, CH-6900, Lugano, Ticino, Switzerland
| | - John Z H Zhang
- State Key Laboratory of Precision Spectroscopy, School of Chemistry and Molecular Engineering, East China Normal University, Shanghai 200062, China and NYU-ECNU Center for Computational Chemistry at NYU Shanghai, Shanghai 200062, China and Department of Chemistry, New York University, NY, NY 10003, USA.
| |
Collapse
|
33
|
Applicability of a thermodynamic cycle approach for a force field parametrization targeting non-aqueous solvation free energies. J Comput Aided Mol Des 2019; 34:71-82. [PMID: 31781991 DOI: 10.1007/s10822-019-00261-5] [Citation(s) in RCA: 0] [Impact Index Per Article: 0] [Reference Citation Analysis] [Abstract] [Key Words] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Received: 06/26/2019] [Accepted: 11/21/2019] [Indexed: 10/25/2022]
Abstract
Accurate solvation free energy ΔGsolv predictions require well parametrized force fields. In order to refit Lennard-Jones (LJ) parameters for improved ΔGsolv predictions for a variety of compound classes and chemical environments, a large number of ΔGsolv calculations is required. As the calculation of ΔGsolv is computational expensive, there is need for efficient but precise calculation methods. In this work, we focus on the computation of non-aqueous solvation free energies. We compare ΔGsolv results from highly precise reference simulations for transferring a solute from the vacuum into a condensed phase and results obtained from a thermodynamic cycle implementation. As test systems, we alter LJ parameters ε and σ of widely used GAFF atom types ca, cl, n1, oh and os in various solute/solvent combinations. We examine the degree of configurational space overlap and find an impact by hydrogen bonds and the solvent accessible surface area. We conclude that the application of a thermodynamic cycle for the parametrization of force fields targeting ΔGsolv is useful if the adaptation of LJ parameters is limited to atom types in the solute or if only ε is changed.
Collapse
|
34
|
Sun Z. BAR-based multi-dimensional nonequilibrium pulling for indirect construction of QM/MM free energy landscapes: from semi-empirical to ab initio. Phys Chem Chem Phys 2019; 21:21942-21959. [PMID: 31552953 DOI: 10.1039/c9cp04113c] [Citation(s) in RCA: 12] [Impact Index Per Article: 2.4] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 01/23/2023]
Abstract
The indirect method for the construction of quantum mechanics (QM)/molecular mechanics (MM) free energy landscapes provides a cheaper alternative for free energy simulations at the QM level. The indirect method features a direct calculation of the free energy profile with a computationally efficient but less accurate Hamiltonian (i.e. low-level Hamiltonian) and a low-level-to-high-level correction. In the thermodynamic cycle, the direct low-level calculation along the physically meaningful reaction coordinate is corrected via the alchemical method, which is often achieved with perturbation-based techniques. In our previous work, a multi-dimensional nonequilibrium pulling framework is proposed for the indirect construction of QM/MM free energy landscapes. Previously, we focused on obtaining semi-empirical QM (SQM) results indirectly from direct MM simulations and MM to SQM corrections. In this work, we apply this method to obtain results under ab initio QM Hamiltonians by combining direct SQM results and SQM to QM corrections. A series of SQM and QM Hamiltonians are benchmarked. It is observed that PM6 achieves the best performance among the low-level Hamiltonians. Therefore, we recommend using PM6 as the low-level theory in the indirect free energy simulation. Considering its higher similarity to the high-level Hamiltonians, PM6 corrected with the bond charge correction could be more accurate than the existing AM1-BCC model. Another central result in the current work is a basic protocol of choosing the strength of restraints and an appropriate time step in nonequilibrium free energy simulation at the stiff spring limit. We provide theoretical derivations to emphasize the importance of using a sufficiently large force constant and choosing an appropriate time step. It is worth noting that a general rule of thumb for choosing the time step, according to our derivation, is that a time step of 1 fs or smaller should be used, as long as the stiff spring approximation is employed, even in simulations with constraints on bonds involving hydrogen atoms.
Collapse
Affiliation(s)
- Zhaoxi Sun
- State Key Laboratory of Precision Spectroscopy, School of Chemistry and Molecular Engineering, East China Normal University, Shanghai 200062, China
| |
Collapse
|
35
|
Gurina DL, Antipova ML, Petrenko VE. Hydroxycinnamic acids in supercritical carbon dioxide. The dependence of cosolvent-induced solubility enhancement on the selective solvation. J Supercrit Fluids 2019. [DOI: 10.1016/j.supflu.2019.04.014] [Citation(s) in RCA: 2] [Impact Index Per Article: 0.4] [Reference Citation Analysis] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 10/27/2022]
|
36
|
Abstract
In practical free energy estimation, the bias is often neglected once it has been shown to vanish in the large-sample limit. Yet finite-sample bias always exists and ought to be considered in any rigorous study. This work develops a metric for bias in a broad class of free energy "bridge estimators" (e.g., Bennett's method). The framework complements existing variance estimation methods and provides a means for comparing systematic and statistical errors. Examples show that, contrary to what is often assumed, the bias can be quite substantial when the sample size is modest.
Collapse
Affiliation(s)
- Brian K Radak
- Theoretical and Computational Biophysics Group, NIH Center for Macromolecular Modeling and Bioinformatics, Beckman Institute for Advanced Science and Technology, University of Illinois at Urbana-Champaign, Urbana, Illinois, 61801-2325, USA
| |
Collapse
|
37
|
Hudson PS, Woodcock HL, Boresch S. Use of Interaction Energies in QM/MM Free Energy Simulations. J Chem Theory Comput 2019; 15:4632-4645. [PMID: 31142113 DOI: 10.1021/acs.jctc.9b00084] [Citation(s) in RCA: 21] [Impact Index Per Article: 4.2] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 12/14/2022]
Abstract
The use of the most accurate (i.e., QM or QM/MM) levels of theory for free energy simulations (FES) is typically not possible. Primarily, this is because the computational cost associated with the extensive configurational sampling needed for converging FES is prohibitive. To ensure the feasibility of QM-based FES, the "indirect" approach is generally taken, necessitating a free energy calculation between the MM and QM/MM potential energy surfaces. Ideally, this step is performed with standard free energy perturbation (Zwanzig's equation) as it only requires simulations be carried out at the low level of theory; however, work from several groups over the past few years has conclusively shown that Zwanzig's equation is ill-suited to this task. As such, many approximations have arisen to mitigate difficulties with Zwanzig's equation. One particularly popular notion is that the convergence of Zwanzig's equation can be improved by using interaction energy differences instead of total energy differences. Although problematic numerical fluctuations (a major problem when using Zwanzig's equation) are indeed reduced, our results and analysis demonstrate that this "interaction energy approximation" (IEA) is theoretically incorrect, and the implicit approximation invoked is spurious at best. Herein, we demonstrate this via solvation free energy calculations using IEA from two different low levels of theory to the same target high level. Results from this proof-of-concept consistently yield the wrong results, deviating by ∼1.5 kcal/mol from the rigorously obtained value.
Collapse
Affiliation(s)
- Phillip S Hudson
- Department of Chemistry , University of South Florida , 4202 East Fowler Avenue, CHE205 , Tampa , Florida 33620-5250 , United States.,Laboratory of Computational Biology , National Institutes of Health, National Heart, Lung and Blood Institute , 12 South Drive, Rm 3053 , Bethesda , Maryland 20892-5690 , United States
| | - H Lee Woodcock
- Department of Chemistry , University of South Florida , 4202 East Fowler Avenue, CHE205 , Tampa , Florida 33620-5250 , United States
| | - Stefan Boresch
- Faculty of Chemistry, Department of Computational Biological Chemistry , University of Vienna , Währingerstraße 17 , Vienna A-1090 , Austria
| |
Collapse
|
38
|
Jaganade T, Chattopadhyay A, Pazhayam NM, Priyakumar UD. Energetic, Structural and Dynamic Properties of Nucleobase-Urea Interactions that Aid in Urea Assisted RNA Unfolding. Sci Rep 2019; 9:8805. [PMID: 31217494 PMCID: PMC6584539 DOI: 10.1038/s41598-019-45010-8] [Citation(s) in RCA: 6] [Impact Index Per Article: 1.2] [Reference Citation Analysis] [Abstract] [Key Words] [MESH Headings] [Grants] [Track Full Text] [Download PDF] [Figures] [Journal Information] [Subscribe] [Scholar Register] [Received: 01/11/2019] [Accepted: 05/28/2019] [Indexed: 01/21/2023] Open
Abstract
Understanding the structure-function relationships of RNA has become increasingly important given the realization of its functional role in various cellular processes. Chemical denaturation of RNA by urea has been shown to be beneficial in investigating RNA stability and folding. Elucidation of the mechanism of unfolding of RNA by urea is important for understanding the folding pathways. In addition to studying denaturation of RNA in aqueous urea, it is important to understand the nature and strength of interactions of the building blocks of RNA. In this study, a systematic examination of the structural features and energetic factors involving interactions between nucleobases and urea is presented. Results from molecular dynamics (MD) simulations on each of the five DNA/RNA bases in water and eight different concentrations of aqueous urea, and free energy calculations using the thermodynamic integration method are presented. The interaction energies between all the nucleobases with the solvent environment and the transfer free energies become more favorable with respect to increase in the concentration of urea. Preferential interactions of urea versus water molecules with all model systems determined using Kirkwood-Buff integrals and two-domain models indicate preference of urea by nucleobases in comparison to water. The modes of interaction between urea and the nucleobases were analyzed in detail. In addition to the previously identified hydrogen bonding and stacking interactions between urea and nucleobases that stabilize the unfolded states of RNA in aqueous solution, NH-π interactions are proposed to be important. Dynamic properties of each of these three modes of interactions have been presented. The study provides fundamental insights into the nature of interaction of urea molecules with nucleobases and how it disrupts nucleic acids.
Collapse
Affiliation(s)
- Tanashree Jaganade
- Center for Computational Natural Sciences and Bioinformatics, International Institute of Information Technology, Hyderabad, 500032, India
| | - Aditya Chattopadhyay
- Center for Computational Natural Sciences and Bioinformatics, International Institute of Information Technology, Hyderabad, 500032, India
| | - Nila M Pazhayam
- Center for Computational Natural Sciences and Bioinformatics, International Institute of Information Technology, Hyderabad, 500032, India
| | - U Deva Priyakumar
- Center for Computational Natural Sciences and Bioinformatics, International Institute of Information Technology, Hyderabad, 500032, India.
| |
Collapse
|
39
|
Sun Z, Wang X, Zhao Q, Zhu T. Understanding Aldose Reductase-Inhibitors interactions with free energy simulation. J Mol Graph Model 2019; 91:10-21. [PMID: 31128525 DOI: 10.1016/j.jmgm.2019.05.011] [Citation(s) in RCA: 9] [Impact Index Per Article: 1.8] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Received: 04/24/2019] [Revised: 05/13/2019] [Accepted: 05/14/2019] [Indexed: 12/15/2022]
Abstract
Aldose Reductase (AR) reduces a variety of substrates, such as aldehydes, aldoses and corticosteroids. It is the first and rate-limiting enzyme of the polyol pathway and is an important target enzyme for diabetes-associated complications, including retinopathy, neuropathy, and nephropathy. Inhibitors targeting this enzyme are structurally different and some of them have side effects. In existing publications, computational techniques are applied to investigate the binding affinities of existing inhibitors and predicting the affinities of newly designed ligands. However, these calculations only employ coarse and approximated methods such as docking and MM/PBSA. Brute force simulations are employed to study the dynamics of the system but no converged statistics is obtained. As a result, these computations provide results not consistent with experimental values and large discrepancies exist. In the current work, we employ the enhanced sampling technique of alchemical free energy simulation to calculate the binding affinities of several ligands targeting AR. The statistical error is defined with care and the correlation in the time-series data is fully considered. The statistically optimal estimators are used to extract the free energy estimates and the predicted results are in agreement with the experimental values. Less computationally demanding end-point free energy methods are also performed to compare their efficiency with the alchemical methods. As is expected, the end-point methods are of less accuracy and reliability compared with the alchemical free energy methods. The decomposition of the free energy difference in each alchemical transformation into the enthalpic and entropic components gives further insights on the thermodynamics. The enthalpy-entropy compensation is observed in this case. As the structural data obtained from experiments are only snapshots and more details are needed to understand the dynamics of the protein-ligand system, the conformational ensemble is analyzed. We identify important residues involved in the protein-ligand binding case and short-lived interactions formed due to fluctuations in the conformational ensemble. The current work shed light on the atomic detailed understanding of the dynamics of AR-inhibitors interactions.
Collapse
Affiliation(s)
- Zhaoxi Sun
- State Key Laboratory of Precision Spectroscopy, School of Chemistry and Molecular Engineering, East China Normal University, Shanghai, 200062, China; Computational Biomedicine (IAS-5/INM-9), Forschungszentrum Jülich, Jülich, 52425, Germany.
| | - Xiaohui Wang
- State Key Laboratory of Precision Spectroscopy, School of Chemistry and Molecular Engineering, East China Normal University, Shanghai, 200062, China; Institute of Computational Science, Università della Svizzera italiana (USI), Via Giuseppe Buffi 13, CH-6900, Lugano, Ticino, Switzerland
| | - Qianqian Zhao
- Computational Biomedicine (IAS-5/INM-9), Forschungszentrum Jülich, Jülich, 52425, Germany; College of Chemistry, Fuzhou University, Fuzhou, 350002, China
| | - Tong Zhu
- State Key Laboratory of Precision Spectroscopy, School of Chemistry and Molecular Engineering, East China Normal University, Shanghai, 200062, China
| |
Collapse
|
40
|
Taddese T, Kitabata M, Okazaki S. All-atom molecular dynamics study on the non-solvent induced phase separation: Thermodynamics of adding water to poly(vinylidene fluoride)/N-methyl-2-pyrrolidone solution. J Chem Phys 2019; 150:184505. [PMID: 31091903 DOI: 10.1063/1.5094088] [Citation(s) in RCA: 1] [Impact Index Per Article: 0.2] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 02/04/2023] Open
Abstract
The change in the thermodynamics when adding water in poly(vinylidene fluoride) (PVDF)/N-methyl-2-pyrrolidone (NMP) solution is studied from all atom molecular dynamics (MD) simulations. This is done by estimating the free energy of mixing of PVDF/NMP solution with increasing volume fraction of water (ϕw) using an appropriately chosen thermodynamic cycle and the Bennett acceptance ratio method. The MD calculations predict the thermodynamic phase separation point of water/NMP/PVDF to be at ϕw = 0.08, in close agreement with the experimental cloud point measurement (ϕw = 0.05). Examining the enthalpic and entropic components of the free energy of mixing reveals that at low concentrations of water, the enthalpy term has the most significant contribution to the miscibility of the ternary system, whereas at higher concentrations of water, the entropy term dominates. Finally, the free energy of mixing was compared with the Flory-Huggins (FH) free energy of mixing by computing the concentration-dependent interaction parameters from MD simulations. The FH model inadequately predicted the miscibility of the PVDF solution, mainly due to its negligence of the excess entropy of mixing.
Collapse
Affiliation(s)
- Tseden Taddese
- Department of Materials Chemistry, Graduate School of Engineering, Nagoya University, Furo-Cho, Chikusa-Ku, Nagoya, Aichi 464-8603, Japan
| | - Masahiro Kitabata
- Department of Materials Chemistry, Graduate School of Engineering, Nagoya University, Furo-Cho, Chikusa-Ku, Nagoya, Aichi 464-8603, Japan
| | - Susumu Okazaki
- Department of Materials Chemistry, Graduate School of Engineering, Nagoya University, Furo-Cho, Chikusa-Ku, Nagoya, Aichi 464-8603, Japan
| |
Collapse
|
41
|
Wang X, Sun Z. Understanding PIM-1 kinase inhibitor interactions with free energy simulation. Phys Chem Chem Phys 2019; 21:7544-7558. [PMID: 30895980 DOI: 10.1039/c9cp00070d] [Citation(s) in RCA: 14] [Impact Index Per Article: 2.8] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 12/19/2022]
Abstract
The proviral integration site of the Moloney leukemia virus (PIM) family includes three homologous members. PIM-1 kinase is an important target in effective therapeutic interventions of lymphomas, prostate cancer and leukemia. In the current work, we performed free energy calculations to calculate the binding affinities of several inhibitors targeting this protein. The alchemical method with integration and perturbation-based estimators and the end-point methods were compared. The computational results indicated that the alchemical method can accurately predict the binding affinities, while the end-point methods give relatively unreliable predictions. Decomposing the free energy difference into enthalpic and entropic components with MBAR reweighting enabled us to investigate the detailed thermodynamic parameters with which the entropy-enthalpy compensation in this protein-ligand binding case is identified. We then studied the conformational ensemble, and the important protein-ligand interactions were identified. The current work sheds light on the understanding of the PIM-1-kinase-inhibitor interactions at the atomic level and will be useful in the further development of potential drugs.
Collapse
Affiliation(s)
- Xiaohui Wang
- State Key Laboratory of Precision Spectroscopy, School of Chemistry and Molecular Engineering, East China Normal University, Shanghai 200062, China
| | | |
Collapse
|
42
|
Hahn DF, Hünenberger PH. Alchemical Free-Energy Calculations by Multiple-Replica λ-Dynamics: The Conveyor Belt Thermodynamic Integration Scheme. J Chem Theory Comput 2019; 15:2392-2419. [PMID: 30821973 DOI: 10.1021/acs.jctc.8b00782] [Citation(s) in RCA: 13] [Impact Index Per Article: 2.6] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 02/06/2023]
Abstract
A new method is proposed to calculate alchemical free-energy differences based on molecular dynamics (MD) simulations, called the conveyor belt thermodynamic integration (CBTI) scheme. As in thermodynamic integration (TI), K replicas of the system are simulated at different 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 K/2 equispaced replicas moving forward and K/2 equispaced replicas 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 simultaneously, 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 or λ-local elevation umbrella sampling (λ-LEUS) 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 barostating schemes, and bypasses the need to carefully preselect a λ-ladder and a swapping-attempt frequency. Compared to λ-LEUS, it eliminates (or drastically reduces) the dead time associated with the preoptimization of a biasing potential. The goal of this article is to provide the mathematical/physical formulation of the proposed CBTI scheme, along with an initial application of the method to the calculation of the hydration free energy of methanol.
Collapse
Affiliation(s)
- David F Hahn
- Laboratory of Physical Chemistry, Department of Chemistry and Applied Biosciences , ETH Zürich , Vladimir-Prelog-Weg 2 , 8093 Zürich , Switzerland
| | - Philippe H Hünenberger
- Laboratory of Physical Chemistry, Department of Chemistry and Applied Biosciences , ETH Zürich , Vladimir-Prelog-Weg 2 , 8093 Zürich , Switzerland
| |
Collapse
|
43
|
Wang X, He Q, Sun Z. BAR-based multi-dimensional nonequilibrium pulling for indirect construction of a QM/MM free energy landscape. Phys Chem Chem Phys 2019; 21:6672-6688. [PMID: 30855611 DOI: 10.1039/c8cp07012a] [Citation(s) in RCA: 17] [Impact Index Per Article: 3.4] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 02/06/2023]
Abstract
Construction of free energy landscapes at the quantum mechanics (QM) level is computationally demanding. As shown in previous studies, by employing an indirect scheme (i.e. constructing a thermodynamic cycle connecting QM states via an alchemical pathway), simulations are converged with much less computational burden. The indirect scheme makes QM/molecular mechanics (MM) free energy simulation orders of magnitude faster than the direct QM/MM schemes. However, the indirect QM/MM simulations were mostly equilibrium sampling based and the nonequilibrium methods were merely exploited in one-dimensional alchemical QM/MM end-state correction at two end states. In this work, we represent a multi-dimensional nonequilibrium pulling scheme for indirect QM/MM free energy simulations, where the whole free energy simulation is performed only with nonequilibrium methods. The collective variable (CV) space we explore is a combination of one alchemical CV and one physically meaningful CV. The current nonequilibrium indirect QM/MM simulation method can be seen as the generalization of equilibrium perturbation based indirect QM/MM methods. The test systems include one backbone dihedral case and one distance case. The two cases are significantly different in size, enabling us to investigate the dependence of the speedup of the indirect scheme on the size of the system. It is shown that the speedup becomes larger when the size of the system becomes larger, which is consistent with the scaling behavior of QM Hamiltonians.
Collapse
Affiliation(s)
- Xiaohui Wang
- State Key Laboratory of Precision Spectroscopy, School of Chemistry and Molecular Engineering, East China Normal University, Shanghai 200062, China
| | | | | |
Collapse
|
44
|
Wang X, Tu X, Deng B, Zhang JZH, Sun Z. BAR-based optimum adaptive steered MD for configurational sampling. J Comput Chem 2019; 40:1270-1289. [PMID: 30762879 DOI: 10.1002/jcc.25784] [Citation(s) in RCA: 18] [Impact Index Per Article: 3.6] [Reference Citation Analysis] [Abstract] [Key Words] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Received: 06/26/2018] [Revised: 11/05/2018] [Accepted: 01/06/2019] [Indexed: 11/08/2022]
Abstract
The equilibrium and nonequilibrium adaptive alchemical free energy simulation methods optimum Bennett's acceptance ratio and optimum crooks' equation (OCE), based on the statistically optimal bidirectional reweighting estimator named Bennett's Acceptance Ratio or Crooks' equation, perform initial sampling in the staging alchemical transformation and then determine the importance rank of different states via the time-derivative of the variance. The method is proven to give speedups compared with the equal time rule. In the current work, we extend the time derivative of variance guided adaptive sampling method to the configurational space, falling in the term of steered MD (SMD). The SMD approach biasing physically meaningful collective variable (CV) such as one dihedral or one distance to pulling the system from one conformational state to another. By minimizing the variance of the free energy differences along the pathway in an optimized way, a new type of adaptive SMD (ASMD) is introduced. As exhibits in the alchemical case, this adaptive sampling method outperforms the traditional equal-time SMD in nonequilibrium stratification. Also, the method gives much more efficient calculation of potential of mean force than the selection criterion-based ASMD scheme, which is proven to be more efficient than traditional SMD. The OCE workflow is periodicity-of-CV dependent while ASMD is not. The performance is demonstrated in a dihedral flipping case and two distance pulling cases, accounting for periodic and nonperiodic CVs, respectively. © 2019 Wiley Periodicals, Inc.
Collapse
Affiliation(s)
- Xiaohui Wang
- State Key Laboratory of Precision Spectroscopy, East China Normal University, Shanghai 200062, China.,Institute of Computational Science, Università della Svizzera italiana (USI), CH-6900, Lugano, Ticino, Switzerland
| | - Xingzhao Tu
- Institute of Organic Chemistry, University of Leipzig, Leipzig 04103, Germany
| | - Boming Deng
- Laboratory of Oil Analysis, Beijing Hangfengkewei Equipment Technology Co., Ltd., Beijing 100141, China
| | - John Z H Zhang
- State Key Laboratory of Precision Spectroscopy, East China Normal University, Shanghai 200062, China.,NYU-ECNU Center for Computational Chemistry at NYU Shanghai, Shanghai 200062, China.,Department of Chemistry, New York University, New York, New York, 10003
| | - Zhaoxi Sun
- State Key Laboratory of Precision Spectroscopy, East China Normal University, Shanghai 200062, China.,Computational Biomedicine (IAS-5/INM-9), Forschungszentrum Juelich, Jülich 52425, Germany
| |
Collapse
|
45
|
König G, Pickard FC, Huang J, Thiel W, MacKerell AD, Brooks BR, York DM. A Comparison of QM/MM Simulations with and without the Drude Oscillator Model Based on Hydration Free Energies of Simple Solutes. Molecules 2018; 23:E2695. [PMID: 30347691 PMCID: PMC6222909 DOI: 10.3390/molecules23102695] [Citation(s) in RCA: 25] [Impact Index Per Article: 4.2] [Reference Citation Analysis] [Abstract] [Key Words] [MESH Headings] [Grants] [Track Full Text] [Download PDF] [Figures] [Journal Information] [Subscribe] [Scholar Register] [Received: 09/08/2018] [Revised: 10/15/2018] [Accepted: 10/16/2018] [Indexed: 12/01/2022] Open
Abstract
Maintaining a proper balance between specific intermolecular interactions and non-specific solvent interactions is of critical importance in molecular simulations, especially when predicting binding affinities or reaction rates in the condensed phase. The most rigorous metric for characterizing solvent affinity are solvation free energies, which correspond to a transfer from the gas phase into solution. Due to the drastic change of the electrostatic environment during this process, it is also a stringent test of polarization response in the model. Here, we employ both the CHARMM fixed charge and polarizable force fields to predict hydration free energies of twelve simple solutes. The resulting classical ensembles are then reweighted to obtain QM/MM hydration free energies using a variety of QM methods, including MP2, Hartree⁻Fock, density functional methods (BLYP, B3LYP, M06-2X) and semi-empirical methods (OM2 and AM1 ). Our simulations test the compatibility of quantum-mechanical methods with molecular-mechanical water models and solute Lennard⁻Jones parameters. In all cases, the resulting QM/MM hydration free energies were inferior to purely classical results, with the QM/MM Drude force field predictions being only marginally better than the QM/MM fixed charge results. In addition, the QM/MM results for different quantum methods are highly divergent, with almost inverted trends for polarizable and fixed charge water models. While this does not necessarily imply deficiencies in the QM models themselves, it underscores the need to develop consistent and balanced QM/MM interactions. Both the QM and the MM component of a QM/MM simulation have to match, in order to avoid artifacts due to biased solute⁻solvent interactions. Finally, we discuss strategies to improve the convergence and efficiency of multi-scale free energy simulations by automatically adapting the molecular-mechanics force field to the target quantum method.
Collapse
Affiliation(s)
- Gerhard König
- Laboratory of Computational Biology, National Heart, Lung, and Blood Institute, National Institutes of Health, Bethesda, MD 20892, USA.
- Laboratory for Biomolecular Simulation Research, Institute for Quantitative Biomedicine, Department of Chemistry and Chemical Biology, Rutgers University, Piscataway, NJ 08854, USA.
- Max-Planck-Institut für Kohlenforschung, 45470 Mülheim an der Ruhr, Germany.
| | - Frank C Pickard
- Laboratory of Computational Biology, National Heart, Lung, and Blood Institute, National Institutes of Health, Bethesda, MD 20892, USA.
| | - Jing Huang
- Laboratory of Computational Biology, National Heart, Lung, and Blood Institute, National Institutes of Health, Bethesda, MD 20892, USA.
- Department of Pharmaceutical Science, School of Pharmacy, University of Maryland, 20 Penn Street, Baltimore, MD 21201, USA.
- School of Life Sciences, Westlake University, 18 Shilongshan Street, Hangzhou 310024, China.
| | - Walter Thiel
- Max-Planck-Institut für Kohlenforschung, 45470 Mülheim an der Ruhr, Germany.
| | - Alexander D MacKerell
- Department of Pharmaceutical Science, School of Pharmacy, University of Maryland, 20 Penn Street, Baltimore, MD 21201, USA.
| | - Bernard R Brooks
- Laboratory of Computational Biology, National Heart, Lung, and Blood Institute, National Institutes of Health, Bethesda, MD 20892, USA.
| | - Darrin M York
- Laboratory for Biomolecular Simulation Research, Institute for Quantitative Biomedicine, Department of Chemistry and Chemical Biology, Rutgers University, Piscataway, NJ 08854, USA.
| |
Collapse
|
46
|
König G, Brooks BR, Thiel W, York DM. On the convergence of multi-scale free energy simulations. MOLECULAR SIMULATION 2018; 44:1062-1081. [PMID: 30581251 DOI: 10.1080/08927022.2018.1475741] [Citation(s) in RCA: 29] [Impact Index Per Article: 4.8] [Reference Citation Analysis] [Abstract] [Key Words] [Track Full Text] [Subscribe] [Scholar Register] [Indexed: 10/16/2022]
Abstract
In this work we employ simple model systems to evaluate the relative performance of two of the most important free energy methods: The Zwanzig equation (also known as "Free energy perturbation") and Bennett's acceptance ratio method (BAR). Although our examples should be transferable to other kinds of free energy simulations, we focus on applications of multi-scale free energy simulations. Such calculations are especially complex, since they connect two different levels of theory with very different requirements in terms of speed, accuracy, sampling and parallelizability. We try to reconcile all those different factors by developing some simple criteria to guide the early stages of the development of a free energy protocol. This is accomplished by quantifying how many λ intermediate steps and how many potential energy evaluations are necessary in order to reach a certain level of convergence.
Collapse
Affiliation(s)
- Gerhard König
- Max-Planck-Institut für Kohlenforschung, 45470 Mülheim an der Ruhr, Germany, EU.,Laboratory for Biomolecular Simulation Research, Center for Integrative Proteomics Research, and Department of Chemistry and Chemical Biology, Rutgers University, Piscataway, New Jersey 08854, USA.,Laboratory of Computational Biology, National Heart Lung and Blood Institute, National Institutes of Health, Bethesda, Maryland 20892, USA
| | - Bernard R Brooks
- Laboratory of Computational Biology, National Heart Lung and Blood Institute, National Institutes of Health, Bethesda, Maryland 20892, USA
| | - Walter Thiel
- Max-Planck-Institut für Kohlenforschung, 45470 Mülheim an der Ruhr, Germany, EU
| | - Darrin M York
- Laboratory for Biomolecular Simulation Research, Center for Integrative Proteomics Research, and Department of Chemistry and Chemical Biology, Rutgers University, Piscataway, New Jersey 08854, USA
| |
Collapse
|
47
|
Giese TJ, York DM. A GPU-Accelerated Parameter Interpolation Thermodynamic Integration Free Energy Method. J Chem Theory Comput 2018; 14:1564-1582. [PMID: 29357243 PMCID: PMC5849537 DOI: 10.1021/acs.jctc.7b01175] [Citation(s) in RCA: 38] [Impact Index Per Article: 6.3] [Reference Citation Analysis] [Abstract] [Grants] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 12/26/2022]
Abstract
There has been a resurgence of interest in free energy methods motivated by the performance enhancements offered by molecular dynamics (MD) software written for specialized hardware, such as graphics processing units (GPUs). In this work, we exploit the properties of a parameter-interpolated thermodynamic integration (PI-TI) method to connect states by their molecular mechanical (MM) parameter values. This pathway is shown to be better behaved for Mg2+ → Ca2+ transformations than traditional linear alchemical pathways (with and without soft-core potentials). The PI-TI method has the practical advantage that no modification of the MD code is required to propagate the dynamics, and unlike with linear alchemical mixing, only one electrostatic evaluation is needed (e.g., single call to particle-mesh Ewald) leading to better performance. In the case of AMBER, this enables all the performance benefits of GPU-acceleration to be realized, in addition to unlocking the full spectrum of features available within the MD software, such as Hamiltonian replica exchange (HREM). The TI derivative evaluation can be accomplished efficiently in a post-processing step by reanalyzing the statistically independent trajectory frames in parallel for high throughput. We also show how one can evaluate the particle mesh Ewald contribution to the TI derivative evaluation without needing to perform two reciprocal space calculations. We apply the PI-TI method with HREM on GPUs in AMBER to predict p Ka values in double stranded RNA molecules and make comparison with experiments. Convergence to under 0.25 units for these systems required 100 ns or more of sampling per window and coupling of windows with HREM. We find that MM charges derived from ab initio QM/MM fragment calculations improve the agreement between calculation and experimental results.
Collapse
Affiliation(s)
- Timothy J. Giese
- Laboratory for Biomolecular Simulation Research, Center for Integrative Proteomics Research, and Department of Chemistry and Chemical Biology, Rutgers University, Piscataway, New Jersey 08854-8087, United States
| | - Darrin M. York
- Laboratory for Biomolecular Simulation Research, Center for Integrative Proteomics Research, and Department of Chemistry and Chemical Biology, Rutgers University, Piscataway, New Jersey 08854-8087, United States
| |
Collapse
|
48
|
Manjari SR, Banavali NK. Structural Articulation of Biochemical Reactions Using Restrained Geometries and Topology Switching. J Chem Inf Model 2018; 58:453-463. [PMID: 29357231 DOI: 10.1021/acs.jcim.7b00699] [Citation(s) in RCA: 2] [Impact Index Per Article: 0.3] [Reference Citation Analysis] [Abstract] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 11/29/2022]
Abstract
A strategy named "restrained geometries and topology switching" (RGATS) is presented to obtain detailed trajectories for complex biochemical reactions using molecular mechanics (MM) methods. It enables prediction of realistic dynamical pathways for chemical reactions, especially for accurately characterizing the structural adjustments of highly complex environments to any proximal biochemical reaction. It can be used to generate reactive conformations, model stepwise or concerted reactions in complex environments, and probe the influence of changes in the environment. Its ability to take reactively nonoptimal conformations and generate favorable starting conformations for a biochemical reaction is illustrated for a proton transfer between two model compounds. Its ability to study concerted reactions in explicit solvent is illustrated using proton transfers between an ammonium ion and two conserved histidines in an ammonia transporter channel embedded in a lipid membrane. Its ability to characterize the changes induced by subtle differences in the active site environment is illustrated using nucleotide addition by a DNA polymerase in the presence of two versus three Mg2+ ions. RGATS can be employed within any MM program and requires no additional software implementation. This allows the full assortment of computational methods implemented in all available MM programs to be used to tackle virtually any question about biochemical reactions that is answerable without using a quantum mechanical (QM) model. It can also be applied to generate reasonable starting structures for more detailed and expensive QM or QM/MM methods. In particular, this strategy enables rapid prediction of reactant, intermediary, or product state structures in any macromolecular context, with the only requirement being that the structure in any one of these states is either known or can be accurately modeled.
Collapse
Affiliation(s)
- Swati R Manjari
- Laboratory of Computational and Structural Biology, Division of Genetics, NYS Department of Health , CMS 2008, Biggs Laboratory, Wadsworth Center, Empire State Plaza, Albany, New York 12201-0509, United States
| | - Nilesh K Banavali
- Laboratory of Computational and Structural Biology, Division of Genetics, NYS Department of Health , CMS 2008, Biggs Laboratory, Wadsworth Center, Empire State Plaza, Albany, New York 12201-0509, United States.,Department of Biomedical Sciences, School of Public Health, State University of New York at Albany , Albany, New York 12222, United States
| |
Collapse
|
49
|
Wang X, Tu X, Zhang JZH, Sun Z. BAR-based optimum adaptive sampling regime for variance minimization in alchemical transformation: the nonequilibrium stratification. Phys Chem Chem Phys 2018; 20:2009-2021. [PMID: 29299568 DOI: 10.1039/c7cp07573a] [Citation(s) in RCA: 26] [Impact Index Per Article: 4.3] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 11/21/2022]
Abstract
Following the previously proposed equilibrate-state sampling based adaptive sampling regime Optimum Bennett Acceptance Ratio (OBAR), we introduce its nonequilibrium extension, the Optimum Crooks' Equation (OCE) in the current work. The efficiency of the NonEquilibrium Work (NEW) stratification is improved by adaptively manipulating the significance of each nonequilibrium realization followed by importance sampling. As is exhibited in the equilibrium case, the nonequilibrium extension outperforms the simple equal time rule used in nonequilibrium stratification in the sense of minimizing the total variance of the free energy estimate. The speedup of this non-equal time rule is more than 1-fold. The Time Derivative of total Variance (TDV) proposed for the OBAR criterion is extended to determine the importance of each nonequilibrium transformation, which is linearly dependent on the variance. The TDV in the nonequilibrium case gives a totally different importance rank from the standard errors of the free energy differences and OBAR TDV due to the duration of nonequilibrium pulling being added into the OCE equation. The performance of the OCE workflow is demonstrated in the solvation of several small molecules with a series of lambda increments and relaxation times between successive perturbations. To the best of our knowledge, such a nonequilibrium adaptive sampling regime in alchemical transformation is unprecedented.
Collapse
Affiliation(s)
- Xiaohui Wang
- State Key Laboratory of Precision Spectroscopy, School of Chemistry and Molecular Engineering, East China Normal University, Shanghai 200062, China.
| | | | | | | |
Collapse
|
50
|
Mecklenfeld A, Raabe G. Comparison of RESP and IPolQ-Mod Partial Charges for Solvation Free Energy Calculations of Various Solute/Solvent Pairs. J Chem Theory Comput 2017; 13:6266-6274. [PMID: 29125770 DOI: 10.1021/acs.jctc.7b00692] [Citation(s) in RCA: 15] [Impact Index Per Article: 2.1] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 11/30/2022]
Abstract
The calculation of solvation free energies ΔGsolv by molecular simulations is of great interest as they are linked to other physical properties such as relative solubility, partition coefficient, and activity coefficient. However, shortcomings in molecular models can lead to ΔGsolv deviations from experimental data. Various studies have demonstrated the impact of partial charges on free energy results. Consequently, calculation methods for partial charges aimed at more accurate ΔGsolv predictions are the subject of various studies in the literature. Here we compare two methods to derive partial charges for the general AMBER force field (GAFF), i.e. the default RESP as well as the physically motivated IPolQ-Mod method that implicitly accounts for polarization costs. We study 29 solutes which include characteristic functional groups of drug-like molecules, and 12 diverse solvents were examined. In total, we consider 107 solute/solvent pairs including two water models TIP3P and TIP4P/2005. Comparison with experimental results yields better agreement for TIP3P, regardless of the partial charge scheme. The overall performance of GAFF/RESP and GAFF/IPolQ-Mod is similar, though specific shortcomings in the description of ΔGsolv for both RESP and IPolQ-Mod can be identified. However, the high correlation between free energies obtained with GAFF/RESP and GAFF/IPolQ-Mod demonstrates the compatibility between the modified charges and remaining GAFF parameters.
Collapse
Affiliation(s)
- Andreas Mecklenfeld
- Institut für Thermodynamik, Technische Universität Braunschweig , Hans-Sommer-Strasse 5, 38106 Braunschweig, Germany.,Center of Pharmaceutical Engineering, Technische Universität Braunschweig , Franz-Liszt-Strasse 35a, 38106 Braunschweig, Germany
| | - Gabriele Raabe
- Institut für Thermodynamik, Technische Universität Braunschweig , Hans-Sommer-Strasse 5, 38106 Braunschweig, Germany
| |
Collapse
|