Consistent Prediction of Mutation Effect on Drug Binding in HIV-1 Protease Using Alchemical Calculations

Machine learning-guided design of potent darunavir analogs targeting HIV-1 proteases: A computational approach for antiretroviral drug discovery

In the pursuit of novel antiretroviral therapies for human immunodeficiency virus type-1 (HIV-1) proteases (PRs), recent improvements in drug discovery have embraced machine learning (ML) techniques to guide the design process. This study employs ensemble learning models to identify crucial substructures as significant features for drug development. Using molecular docking techniques, a collection of 160 darunavir (DRV) analogs was designed based on these key substructures and subsequently screened using molecular docking techniques. Chemical structures with high fitness scores were selected, combined, and one-dimensional (1D) screening based on beyond Lipinski's rule of five (bRo5) and ADME (absorption, distribution, metabolism, and excretion) prediction implemented in the Combined Analog generator Tool (CAT) program. A total of 473 screened analogs were subjected to docking analysis through convolutional neural networks scoring function against both the wild-type (WT) and 12 major mutated PRs. DRV analogs with negative changes in binding free energy (ΔΔ G bind ) compared to DRV could be categorized into four attractive groups based on their interactions with the majority of vital PRs. The analysis of interaction profiles revealed that potent designed analogs, targeting both WT and mutant PRs, exhibited interactions with common key amino acid residues. This observation further confirms that the ML model-guided approach effectively identified the substructures that play a crucial role in potent analogs. It is expected to function as a powerful computational tool, offering valuable guidance in the identification of chemical substructures for synthesis and subsequent experimental testing.

FMO-guided design of darunavir analogs as HIV-1 protease inhibitors

The prevalence of HIV-1 infection continues to pose a significant global public health issue, highlighting the need for antiretroviral drugs that target viral proteins to reduce viral replication. One such target is HIV-1 protease (PR), responsible for cleaving viral polyproteins, leading to the maturation of viral proteins. While darunavir (DRV) is a potent HIV-1 PR inhibitor, drug resistance can arise due to mutations in HIV-1 PR. To address this issue, we developed a novel approach using the fragment molecular orbital (FMO) method and structure-based drug design to create DRV analogs. Using combinatorial programming, we generated novel analogs freely accessible via an on-the-cloud mode implemented in Google Colab, Combined Analog generator Tool (CAT). The designed analogs underwent cascade screening through molecular docking with HIV-1 PR wild-type and major mutations at the active site. Molecular dynamics (MD) simulations confirmed the assess ligand binding and susceptibility of screened designed analogs. Our findings indicate that the three designed analogs guided by FMO, 19-0-14-3, 19-8-10-0, and 19-8-14-3, are superior to DRV and have the potential to serve as efficient PR inhibitors. These findings demonstrate the effectiveness of our approach and its potential to be used in further studies for developing new antiretroviral drugs.

Modelling the transmission dynamics and optimal control strategies for HIV infection in China

In order to end the AIDS epidemic by 2030 that was put forward by the Joint United Nations Programme on HIV/AIDS in 2014, China needs to take more effective measures to achieve the three 90% goals (90-90-90). We establish a compartmental model to study the dynamics of HIV transmission with control strategies. The analytical results show the existence and stability of the disease-free equilibrium and endemic equilibrium. An optimal control model is constructed to evaluate the impacts of control measures. The simulation results show that the optimal control strategy proposed in this work can eradicate AIDS by 2030. The cost-effectiveness analysis indicates that the cost of the control strategy that combines screening for latent individuals and enhancing education for unaware infected individuals is the lowest. Our findings can provide guidance for public health authorities on effective mitigation strategies to achieve the goals proposed by the United Nations Program on HIV/AIDS.

Accurately Predicting Protein pKa Values Using Nonequilibrium Alchemy

The stability, solubility, and function of a protein depend on both its net charge and the protonation states of its individual residues. pKa is a measure of the tendency for a given residue to (de)protonate at a specific pH. Although pKa values can be resolved experimentally, theory and computation provide a compelling alternative. To this end, we assess the applicability of a nonequilibrium (NEQ) alchemical free energy method to the problem of pKa prediction. On a data set of 144 residues that span 13 proteins, we report an average unsigned error of 0.77 ± 0.09, 0.69 ± 0.09, and 0.52 ± 0.04 pK for aspartate, glutamate, and lysine, respectively. This is comparable to current state-of-the-art predictors and the accuracy recently reached using free energy perturbation methods (e.g., FEP+). Moreover, we demonstrate that our open-source, pmx-based approach can accurately resolve the pKa values of coupled residues and observe a substantial performance disparity associated with the lysine partial charges in Amber14SB/Amber99SB*-ILDN, for which an underused fix already exists.

Identifying and Overcoming the Sampling Challenges in Relative Binding Free Energy Calculations of a Model Protein:Protein Complex

Relative alchemical binding free energy calculations are routinely used in drug discovery projects to optimize the affinity of small molecules for their drug targets. Alchemical methods can also be used to estimate the impact of amino acid mutations on protein:protein binding affinities, but these calculations can involve sampling challenges due to the complex networks of protein and water interactions frequently present in protein:protein interfaces. We investigate these challenges by extending a graphics processing unit (GPU)-accelerated open-source relative free energy calculation package (Perses) to predict the impact of amino acid mutations on protein:protein binding. Using the well-characterized model system barnase:barstar, we describe analyses for identifying and characterizing sampling problems in protein:protein relative free energy calculations. We find that mutations with sampling problems often involve charge-changes, and inadequate sampling can be attributed to slow degrees of freedom that are mutation-specific. We also explore the accuracy and efficiency of current state-of-the-art approaches─alchemical replica exchange and alchemical replica exchange with solute tempering─for overcoming relevant sampling problems. By employing sufficiently long simulations, we achieve accurate predictions (RMSE 1.61, 95% CI: [1.12, 2.11] kcal/mol), with 86% of estimates within 1 kcal/mol of the experimentally determined relative binding free energies and 100% of predictions correctly classifying the sign of the changes in binding free energies. Ultimately, we provide a model workflow for applying protein mutation free energy calculations to protein:protein complexes, and importantly, catalog the sampling challenges associated with these types of alchemical transformations. Our free open-source package (Perses) is based on OpenMM and is available at https://github.com/choderalab/perses.

Interactions between selectivity filter and pore helix control filter gating in the MthK channel

K+ channel activity can be limited by C-type inactivation, which is likely initiated in part by dissociation of K+ ions from the selectivity filter and modulated by the side chains that surround it. While crystallographic and computational studies have linked inactivation to a "collapsed" selectivity filter conformation in the KcsA channel, the structural basis for selectivity filter gating in other K+ channels is less clear. Here, we combined electrophysiological recordings with molecular dynamics simulations, to study selectivity filter gating in the model potassium channel MthK and its V55E mutant (analogous to KcsA E71) in the pore-helix. We found that MthK V55E has a lower open probability than the WT channel, due to decreased stability of the open state, as well as a lower unitary conductance. Simulations account for both of these variables on the atomistic scale, showing that ion permeation in V55E is altered by two distinct orientations of the E55 side chain. In the "vertical" orientation, in which E55 forms a hydrogen bond with D64 (as in KcsA WT channels), the filter displays reduced conductance compared to MthK WT. In contrast, in the "horizontal" orientation, K+ conductance is closer to that of MthK WT; although selectivity filter stability is lowered, resulting in more frequent inactivation. Surprisingly, inactivation in MthK WT and V55E is associated with a widening of the selectivity filter, unlike what is observed for KcsA and reminisces recent structures of inactivated channels, suggesting a conserved inactivation pathway across the potassium channel family.

Comprehending the Structure, Dynamics, and Mechanism of Action of Drug-Resistant HIV Protease

Since the emergence of the Human Immunodeficiency Virus (HIV) in the 1980s, strategies to combat HIV-AIDS are continuously evolving. Among the many tested targets to tackle this virus, its protease enzyme (PR) was proven to be an attractive option that brought about numerous research publications and ten FDA-approved drugs to inhibit the PR activity. However, the drug-induced mutations in the enzyme made these small molecule inhibitors ineffective with prolonged usage. The research on HIV PR, therefore, remains a thrust area even today. Through this review, we reiterate the importance of understanding the various structural and functional components of HIV PR in redesigning the structure-based small molecule inhibitors. We also discuss at length the currently available FDA-approved drugs and how these drug molecules induced mutations in the enzyme structure. We then recapitulate the reported mechanisms on how these drug-resistant variants remain sufficiently active to cleave the natural substrates. We end with the future scope covering the recently proposed strategies that show promise to deal with the mutations.

Identifying and overcoming the sampling challenges in relative binding free energy calculations of a model protein:protein complex

Relative alchemical binding free energy calculations are routinely used in drug discovery projects to optimize the affinity of small molecules for their drug targets. Alchemical methods can also be used to estimate the impact of amino acid mutations on protein:protein binding affinities, but these calculations can involve sampling challenges due to the complex networks of protein and water interactions frequently present in protein:protein interfaces. We investigate these challenges by extending a GPU-accelerated opensource relative free energy calculation package (Perses) to predict the impact of amino acid mutations on protein:protein binding. Using the well-characterized model system barnase:barstar, we describe analyses for identifying and characterizing sampling problems in protein:protein relative free energy calculations. We find that mutations with sampling problems often involve charge-changes, and inadequate sampling can be attributed to slow degrees of freedom that are mutation-specific. We also explore the accuracy and efficiency of current state-of-the-art approaches-alchemical replica exchange and alchemical replica exchange with solute tempering-for overcoming relevant sampling problems. By employing sufficiently long simulations, we achieve accurate predictions (RMSE 1.61, 95% CI: [1.12, 2.11] kcal/mol), with 86% of estimates within 1 kcal/mol of the experimentally-determined relative binding free energies and 100% of predictions correctly classifying the sign of the changes in binding free energies. Ultimately, we provide a model workflow for applying protein mutation free energy calculations to protein:protein complexes, and importantly, catalog the sampling challenges associated with these types of alchemical transformations. Our free open-source package (Perses) is based on OpenMM and available at https://github.com/choderalab/perses .

Accurate calculation of absolute free energy of binding for SHP2 allosteric inhibitors using free energy perturbation

Accurate prediction of binding affinity is a primary objective in structure-based drug discovery. A free energy perturbation (FEP) method based on molecular dynamics simulation shows great promise for protein-ligand binding affinity predictions. However, accurate calculation of binding affinity for allosteric inhibitors remains unknown and elusive, which hampers the discovery of allosteric inhibitors. Allosteric inhibitors exhibit several significant advantages over orthosteric inhibitors including higher specificity and lower side effects. Allosteric inhibitors against SHP2 are thought to be beneficial not only for diseases related to metabolism, but also for cancer, which make SHP2 a potential drug target. However, high structural sensitivity makes structural optimization of SHP2 allosteric inhibitors face challenges. Herein, we calculated the absolute binding free energy of SHP2 allosteric inhibitors using the FEP method by employing different λ-windows/simulation time sampling strategies. A simulation run with 32 λ-windows/64 ps sampling strategy delivered an excellent correlation (r = 0.96) and an unprecedented low mean absolute error of 0.5 kcal mol^-1 between predicted binding free energies and experimental ones, outperforming the MM/PBSA method. Our study demonstrates the possibility to accurately calculate the absolute binding free energy of allosteric inhibitors using FEP, which offers exciting prospects for the discovery of more effective allosteric inhibitors.

Hybrid QM/MM Free-Energy Evaluation of Drug-Resistant Mutational Effect on the Binding of an Inhibitor Indinavir to HIV-1 Protease

A human immunodeficiency virus-1 (HIV-1) protease is a homodimeric aspartic protease essential for the replication of HIV. The HIV-1 protease is a target protein in drug discovery for antiretroviral therapy, and various inhibitor molecules of transition state analogues have been developed. However, serious drug-resistant mutants have emerged. For understanding the molecular mechanism of the drug resistance, an accurate examination of the impacts of the mutations on ligand binding and enzymatic activity is necessary. Here, we present a molecular simulation study on the ligand binding of indinavir, a potent transition state analogue inhibitor, to the wild-type protein and a V82T/I84V drug-resistant mutant of the HIV-1 protease. We employed a hybrid ab initio quantum mechanical/molecular mechanical (QM/MM) free-energy optimization technique which combines a highly accurate QM description of the ligand molecule and its interaction with statistically ample conformational sampling of the MM protein environment by long-time molecular dynamics simulations. Through the free-energy calculations of protonation states of catalytic groups at the binding pocket and of the ligand-binding affinity changes upon the mutations, we successfully reproduced the experimentally observed significant reduction of the binding affinity upon the drug-resistant mutations and elucidated the underlying molecular mechanism. The present study opens the way for understanding the molecular mechanism of drug resistance through the direct quantitative comparison of ligand binding and enzymatic reaction with the same accuracy.

Accurate absolute free energies for ligand-protein binding based on non-equilibrium approaches

The accurate calculation of the binding free energy for arbitrary ligand-protein pairs is a considerable challenge in computer-aided drug discovery. Recently, it has been demonstrated that current state-of-the-art molecular dynamics (MD) based methods are capable of making highly accurate predictions. Conventional MD-based approaches rely on the first principles of statistical mechanics and assume equilibrium sampling of the phase space. In the current work we demonstrate that accurate absolute binding free energies (ABFE) can also be obtained via theoretically rigorous non-equilibrium approaches. Our investigation of ligands binding to bromodomains and T4 lysozyme reveals that both equilibrium and non-equilibrium approaches converge to the same results. The non-equilibrium approach achieves the same level of accuracy and convergence as an equilibrium free energy perturbation (FEP) method enhanced by Hamiltonian replica exchange. We also compare uni- and bi-directional non-equilibrium approaches and demonstrate that considering the work distributions from both forward and reverse directions provides substantial accuracy gains. In summary, non-equilibrium ABFE calculations are shown to yield reliable and well-converged estimates of protein-ligand binding affinity.

A shift of dynamic equilibrium between the KIT active and inactive states causes drug resistance

Tyrosine phosphorylation, a highly regulated post-translational modification, is carried out by the enzyme tyrosine kinase (TK). TKs are important mediators in signaling cascades, facilitating diverse biological processes in response to stimuli. TKs may acquire mutations leading to malignancy and are viable targets for anti-cancer drugs. Mast/stem cell growth factor receptor KIT is a TK involved in cell differentiation, whose dysregulation leads to various types of cancer, including gastrointestinal stromal tumors, leukemia, and melanoma. KIT can be targeted by a range of inhibitors that predominantly bind to the inactive state of the enzyme. A mutation Y823D in the activation loop of KIT is known to be responsible for the loss of sensitivity to some drugs in metastatic tumors. We used all-atom molecular dynamics simulations to study the impact of Y823D on the KIT conformation and dynamics and compared it to the effect of phosphorylation of Y823. We simulated in total 6.4 μs of wild-type, mutant and phosphorylated KIT in the active- and inactive-state conformations. We found that Y823D affects the protein dynamics differently: in the active state, the mutation increases the protein stability, whereas in the inactive state it induces local destabilization, thus shifting the dynamic equilibrium towards the active state, altering the communication between distant regulatory regions. The observed dynamics of the Y823D mutant is similar to the dynamics of KIT phosphorylated at position Y823, thus we hypothesize that this mutation mimics a constitutively active kinase, which is not responsive to inhibitors that bind its inactive conformation.

Non-active site mutants of HIV-1 protease influence resistance and sensitisation towards protease inhibitors

Background

HIV-1 can develop resistance to antiretroviral drugs, mainly through mutations within the target regions of the drugs. In HIV-1 protease, a majority of resistance-associated mutations that develop in response to therapy with protease inhibitors are found in the protease’s active site that serves also as a binding pocket for the protease inhibitors, thus directly impacting the protease-inhibitor interactions. Some resistance-associated mutations, however, are found in more distant regions, and the exact mechanisms how these mutations affect protease-inhibitor interactions are unclear. Furthermore, some of these mutations, e.g. N88S and L76V, do not only induce resistance to the currently administered drugs, but contrarily induce sensitivity towards other drugs. In this study, mutations N88S and L76V, along with three other resistance-associated mutations, M46I, I50L, and I84V, are analysed by means of molecular dynamics simulations to investigate their role in complexes of the protease with different inhibitors and in different background sequence contexts.

Results

Using these simulations for alchemical calculations to estimate the effects of mutations M46I, I50L, I84V, N88S, and L76V on binding free energies shows they are in general in line with the mutations’ effect on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$IC_{50}$$\end{document}IC50 values. For the primary mutation L76V, however, the presence of a background mutation M46I in our analysis influences whether the unfavourable effect of L76V on inhibitor binding is sufficient to outweigh the accompanying reduction in catalytic activity of the protease. Finally, we show that L76V and N88S changes the hydrogen bond stability of these residues with residues D30/K45 and D30/T31/T74, respectively.

Conclusions

We demonstrate that estimating the effect of both binding pocket and distant mutations on inhibitor binding free energy using alchemical calculations can reproduce their effect on the experimentally measured \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$IC_{50}$$\end{document}IC50 values. We show that distant site mutations L76V and N88S affect the hydrogen bond network in the protease’s active site, which offers an explanation for the indirect effect of these mutations on inhibitor binding. This work thus provides valuable insights on interplay between primary and background mutations and mechanisms how they affect inhibitor binding.

Effect of mutations on binding of ligands to guanine riboswitch probed by free energy perturbation and molecular dynamics simulations

Riboswitches can regulate gene expression by direct and specific interactions with ligands and have recently attracted interest as potential drug targets for antibacterial. In this work, molecular dynamics (MD) simulations, free energy perturbation (FEP) and molecular mechanics generalized Born surface area (MM-GBSA) methods were integrated to probe the effect of mutations on the binding of ligands to guanine riboswitch (GR). The results not only show that binding free energies predicted by FEP and MM-GBSA obtain an excellent correlation, but also indicate that mutations involved in the current study can strengthen the binding affinity of ligands GR. Residue-based free energy decomposition was applied to compute ligand-nucleotide interactions and the results suggest that mutations highly affect interactions of ligands with key nucleotides U22, U51 and C74. Dynamics analyses based on MD trajectories indicate that mutations not only regulate the structural flexibility but also change the internal motion modes of GR, especially for the structures J12, J23 and J31, which implies that the aptamer domain activity of GR is extremely plastic and thus readily tunable by nucleotide mutations. This study is expected to provide useful molecular basis and dynamics information for the understanding of the function of GR and possibility as potential drug targets for antibacterial.

Large scale relative protein ligand binding affinities using non-equilibrium alchemy

Ligand binding affinity calculations based on molecular dynamics (MD) simulations and non-physical (alchemical) thermodynamic cycles have shown great promise for structure-based drug design. However, their broad uptake and impact is held back by the notoriously complex setup of the calculations. Only a few tools other than the free energy perturbation approach by Schrödinger Inc. (referred to as FEP+) currently enable end-to-end application. Here, we present for the first time an approach based on the open-source software pmx that allows to easily set up and run alchemical calculations for diverse sets of small molecules using the GROMACS MD engine. The method relies on theoretically rigorous non-equilibrium thermodynamic integration (TI) foundations, and its flexibility allows calculations with multiple force fields. In this study, results from the Amber and Charmm force fields were combined to yield a consensus outcome performing on par with the commercial FEP+ approach. A large dataset of 482 perturbations from 13 different protein-ligand datasets led to an average unsigned error (AUE) of 3.64 ± 0.14 kJ mol-1, equivalent to Schrödinger's FEP+ AUE of 3.66 ± 0.14 kJ mol-1. For the first time, a setup is presented for overall high precision and high accuracy relative protein-ligand alchemical free energy calculations based on open-source software.

Accurate Estimation of Ligand Binding Affinity Changes upon Protein Mutation

The design of proteins with novel ligand-binding functions holds great potential for application in biomedicine and biotechnology. However, our ability to engineer ligand-binding proteins is still limited, and current approaches rely primarily on experimentation. Computation could reduce the cost of the development process and would allow rigorous testing of our understanding of the principles governing molecular recognition. While computational methods have proven successful in the early stages of the discovery process, optimization approaches that can quantitatively predict ligand affinity changes upon protein mutation are still lacking. Here, we assess the ability of free energy calculations based on first-principles statistical mechanics, as well as the latest Rosetta protocols, to quantitatively predict such affinity changes on a challenging set of 134 mutations. After evaluating different protocols with computational efficiency in mind, we investigate the performance of different force fields. We show that both the free energy calculations and Rosetta are able to quantitatively predict changes in ligand binding affinity upon protein mutations, yet the best predictions are the result of combining the estimates of both methods. These closely match the experimentally determined ΔΔG values, with a root-mean-square error of 1.2 kcal/mol for the full benchmark set and of 0.8 kcal/mol for a subset of protein systems providing the most reproducible results. The currently achievable accuracy offers the prospect of being able to employ computation for the optimization of ligand-binding proteins as well as the prediction of drug resistance.

Exploring molecular mechanism of allosteric inhibitor to relieve drug resistance of multiple mutations in HIV-1 protease by enhanced conformational sampling

Recently, allosteric regulations of HIV-1 protease (PR) are suggested as a promising approach to relieve drug resistance of mutations toward inhibitors targeting the active site of PR. Replica-exchange molecular dynamics (REMD) simulations and normal mode analysis (NMA) are integrated to enhance conformational sampling of PR. Molecular mechanics generalized Born surface area (MM-GBSA) method was applied to calculate binding free energies of three inhibitors APV, DRV, and NIT to the wild-type (WT) and multidrug resistance (MDR) PRs. The results suggest that binding free energies of APV and DRV are decreased in the MDR PR relative to the WT PR, suggesting drug resistance of mutations on these two inhibitors. However, the binding ability of the allosteric inhibitor NIT is not impaired in the MDR PR. In addition, internal dynamics analysis based on REMD simulations proves that mutations hardly produce obvious effect on the conformation of the MDR PR in comparison to the WT PR. Scanning of hydrophobic contacts and hydrogen bond contacts of inhibitors with residues of PRs on the concatenated trajectories of REMD demonstrates that mutations change the symmetric interaction networks of APV and DRV with PR, but do not generate obvious influence on the asymmetric interaction network of NIT with PR. In summary, allosteric inhibitor NIT can adapt the MDR PR better than those inhibitors toward the active site of PR, thus allosteric inhibitors of PR may be a possible channel to overcome drug resistance of PR.