Linear interaction energy models for β-secretase (BACE) inhibitors: Role of van der Waals,electrostatic, and continuum-solvation terms

Computing the binding affinity of a protein–ligand complex is one of the most fundamental and difficult tasks in computer-aided drug design. Many approaches for computing binding affinities can be classified as linear interaction energy (LIE) models as they rely on some type of linear fit of computed interaction energies between ligand and protein. We have examined the computed interaction energies of a series of β-secretase (BACE) inhibitors in terms of van der Waals, coulombic, and continuum-solvation contributions to ligand binding. We have also systematically examined the effect of different protonation states of the protein and ligands. We find that the binding affinities are relatively insensitive to the protonation state of the protein when neutral ligands are considered. Inclusion of charged ligands leads to large deviations in the coulomb, solvation, and even van der Waals terms. The latter is due to increased repulsive van der Waals interactions in the complex due to the strong coulomb attraction found between oppositely charged functional groups in the protein and ligand. In general, we find that the best models are obtained when the protein is judiciously charged (e.g. Asp32⁻, Arg235⁺) and the potentially charged ligands are treated as neutral.     

Linear interaction energy models for beta-secretase (BACE) inhibitors: Role of van der Waals, electrostatic, and continuum-solvation terms

Computing the binding affinity of a protein–ligand complex is one of the most fundamental and difficult tasks in computer-aided drug design. Many approaches for computing binding affinities can be classified as linear interaction energy (LIE) models as they rely on some type of linear fit of computed interaction energies between ligand and protein. We have examined the computed interaction energies of a series of β-secretase (BACE) inhibitors in terms of van der Waals, coulombic, and continuum-solvation contributions to ligand binding. We have also systematically examined the effect of different protonation states of the protein and ligands. We find that the binding affinities are relatively insensitive to the protonation state of the protein when neutral ligands are considered. Inclusion of charged ligands leads to large deviations in the coulomb, solvation, and even van der Waals terms. The latter is due to increased repulsive van der Waals interactions in the complex due to the strong coulomb attraction found between oppositely charged functional groups in the protein and ligand. In general, we find that the best models are obtained when the protein is judiciously charged (e.g. Asp32⁻, Arg235⁺) and the potentially charged ligands are treated as neutral.     

Predicting proteinase specificities from free energy calculations

The role of the primary binding residue (P1) in complexes between three different subtilases (subtilisin Carlsberg, thermitase and proteinase K) and their canonical protein inhibitor eglin c have been studied by free energy calculations. Based on the crystal structures of eglin c in complex with subtilisin Carlsberg and thermitase, and a homology model of the eglin c-proteinase K complex, a total of 57 mutants have been constructed and docked into their host proteins. The binding free energy was then calculated using molecular dynamics (MD) simulations combined with the linear interaction energy (LIE) method for all complexes differing only in the nature of the amino acid at the P1 position. LIE calculations for 19 different complexes for each subtilase were thus carried out excluding proline. The effects of substitutions at the P1 position on the binding free energies are found to be very large, and positively charged residues (Arg, Lys and His) are particularly deleterious for all three enzymes. The charged variants of the acidic side chains are found to bind more favorably as compared to their protonated states in all three subtilases. Furthermore, hydrophobic amino acids are accommodated most favorably at the S1-site in all three enzymes. Comparison of the three series of binding free energies shows only minor differences in the 19 computed relative binding free energies among these subtilases. This is further reflected in the correlation coefficient between the 23 relative binding free energies obtained, including the possible protonation states of ionizable side chains, but excluding the P1 Pro, for subtilisin Carlsberg versus thermitase (0.95), subtilisin versus proteinase K (0.94) and thermitase versus proteinase K (0.96).     

Could MM-GBSA be accurate enough for calculation of absolute protein/ligand binding free energies?

Implicit solvation methods such as MM-GBSA, when applied to evaluating protein/ligand binding free energies, are widely believed to be accurate only for the estimation of relative binding free energies for a congeneric series of ligands. In this work, we show that the MM-GBSA flavor of Prime 3.0, VSGB-2.0, with a variable dielectric model and a novel energy function, could be approaching the accuracy required for evaluating absolute binding free energies, albeit, through a linear regression fit. The data-set used for validation includes 106 protein–ligand complexes that were carefully selected to control for variability in the affinity data as well as error in the modeled complexes. Through systematic analysis, we also quantify the degradation in the R² of fit between experimental and calculated values with either greater variability in the affinity data or an increase in error in the modeled protein/ligand complexes. Limitations for its application in drug discovery are discussed along with the identification of areas for future development.     

Monte Carlo simulations of I2- (CO2)16 and I2- (N2O)16 clusters. Minimum energy structures and solvation energy

Structures of the I2- (CO2)16 and I2- (N2O)16 clusters have been studied along the dissociation path of the halide ion. Minimum energy structures and solvation energies have been estimated. The influence of the electrostatic interaction on the I2- charge distribution and the structure of the solvation shell has been discussed. It has been shown that the small dipole moment of the N2O molecule does not have much influence on the cluster properties, and the large quadrupole moments of both solvent clusters give the major contribution into electrostatic interactions. The structure of the clusters is more sensitive to different potential models than to the dipole moment of the N2O molecule.     

Scoring binding affinity of multiple ligands using implicit solvent and a single molecular dynamics trajectory: application to influenza neuraminidase

We explore a perturbative approach to calculation of binding free energy of multiple ligands, based on a single molecular dynamics simulation of a reference ligand-receptor complex and analysis via a hybrid force field/continuum model potential. The methodology is applied to prediction of relative binding free energies of 10 Influenza neuraminidase inhibitors, using Poisson-Boltzmann and generalised Born models of implicit solvent. These single-step MM-PB/SA and MM-GB/SA approaches predict the experimentally most potent ligand as first- or second-ranked according to total binding free energy. Ranking of inhibitors displays only moderate sensitivity to the choice of reference trajectory and ligand partial charge scheme. When ranked according to total electrostatic binding free energy, correlation with experiment improves (r(2) of 0.72); this may be related to underestimated first solvation shell effects by the implicit water models. Therefore, to increase the generality of this single-step approach as part of a potential computational compound optimisation strategy, further development of the treatment of short-range solvent interactions is warranted.     

Identification of non-resistant ROS-1 inhibitors using structure based pharmacophore analysis

Proto-oncogene receptor tyrosine kinase ROS-1 plays a key role in regulating a variety of cancers mainly non-small cell lung cancer (NSCLC). The marketed ROS-1 inhibitors such as Crizotinib suffer from the tribulations of growing resistance due to mutations primarily Gly2032Arg in the ROS-1 protein. To curb the problem of resistance, researchers have developed inhibitors such as Lorlatinib against the mutant protein. The present study was designed to identify inhibitors against wild type (WT) as well as mutant ROS-1 protein that will offer a broader spectrum of activity. Exploring crystal structure of ROS-1 complexed with Lorlatinib, receptor-ligand pharmacophore model was developed using Discovery Studio (DS) software. The developed pharmacophore model consisted of one hydrogen bond acceptor (HBA), one hydrogen bond donor (HBD) and two hydrophobic features (HY), subsequently utilized for virtual screening of commercially available databases and the retrieved hits were further subjected to fitness score and Lipinski’s filter. Thereafter, the retrieved hits were docked in WT and mutated (Gly2032Arg) proteins of ROS-1. Total five molecules were retrieved with good docking scores and good binding interactions within the active site of WT and mutated ROS-1. The binding energies of the ligand-receptor complexes were predicted via calculation of MM-GBSA score. To predict the stability of the ligand receptor complexes with mutant and wild type proteins, molecular dynamic simulation was performed. Thus, these identified hits show good binding affinities with WT and mutant ROS-1 proteins that may be further evaluated for their in-vitro/in-vivo activity.     

Insight into the structural mechanism for PKBα allosteric inhibition by molecular dynamics simulations and free energy calculations

Protein kinase B (PKB/Akt) is an attractive target for the treatment of tumor. Unlike PKB's ATP-competitive inhibitors, its allosteric inhibitors can maintain PKB's inactive state via its binding in a pocket between PH domain and kinase domain, which specifically inhibit PKB by preventing the phosphorylations of Thr308 and Ser473. In the present studies, MD simulations were performed on three allosteric inhibitors with different inhibitory potencies (IC₅₀) to investigate the interaction modes between the inhibitors and PKBα. MM/GB(PB)SA were further applied to calculate the binding free energies of these inhibitors binding to PKBα. The computed binding free energies were consistent with the ranking of their experimental bioactivities. The key residues of PKBα interacting with the allosteric inhibitor were further discussed by analyzing the different interaction modes of these three inhibitors binding to PKBα and by calculating binding free energy contributions of corresponding residues around the binding pocket. The structural requirements were then summarized for the allosteric inhibitor binding to PKBα. A possible structural mechanism of PKBα inhibition induced by the binding of allosteric inhibitor was formulated. The current studies indicate that there should be an optimum balance between the van der Waals and total electrostatic interactions for further designing of PKBα allosteric inhibitors.     

Computational alanine scanning and free energy decomposition for E. coli type I signal peptidase with lipopeptide inhibitor complex

A thorough investigation of different roles of Escherichia coli type I signal peptidase residues binding to lipopeptide inhibitor has been performed by a combination of computational alanine scanning mutagenesis and free energy decomposition methods. PB and GB models are both used to evaluate the binding free energy in computational alanine scanning method and only GB model can be used to decompose the binding free energy on a per-residue basis. The regression analysis between the PB and GB model and also between the computational alanine scanning and free energy decomposition have been reported with a correlation coefficient of 0.96 and 0.83, respectively, which suggest they are both in fair agreement with each other. Moreover, the contribution components from van der Waals, electrostatic interaction, non-polar and polar energy of solvation, have been determined as well as the effects of backbones and side-chains. The results indicate that Lys145 is the most important residue for the binding but also acts as a general base, activating Ser90 to increase its nucleophility, recognizing and stabilizing the binding of lipopeptide inhibitor to the E. coli SPase. The hydroxyl group of Ser88 plays a key role for the binding of the inhibitor. Ser90 contributes more to the intramolecular interaction than to the intermolecular interaction. Tyr143 and Phe84 contribute larger van der Waals interaction energies, indicating that these residues can be important for the selection based on the shape of the inhibitors. The contributions from other several interfacial residues of the E. coli SPase are also analyzed. This study can be a guide for the optimization of lipopeptide inhibitors and future design of new therapeutic agents for the treatment of bacterial infections.     

Computational studies of the resistance patterns of mutant HIV-1 aspartic proteases towards ABT-538 (ritonavir) and design of new derivatives

Kinetic characterization and cross resistance pattern studies of HIV-1 aspartic protase (PR) inhibitors have shown that some mutations cause considerable reduction in inhibition efficiency. We have performed a computational study of the binding of ABT-538 (ritonavir) with wild type (wt) PR and 12 model mutant structures (R8Q, V321, M461, V82A, V82F, V821, I84V, M46I/V82F, M46I/I84V, V32I/I84V, V82F/I84V and V32I/K45I/F53L/A71V/I84V/L89M (6X)) for which inhibition data are available. Our computational studies indicate a significant correlation between computed complexation energies of ABT-538 with the modeled mutant enzyme structures and the corresponding experimental inhibition constants. By evaluating non-bonding interaction energies between the inhibitor and the mutant enzymes, we have carried out a mechanistic analysis to ascertain the reasons underlying the decrease in binding affinities. This analysis indicated that several residues in addition to the mutated residues contribute to the loss of binding. Taking these considerations into account, a number of new derivatives of ABT-538 were designed, so as to increase van der Waal's and hydrogen bonding interactions with selected mutants. A significant improvement in calculated complexation energies towards both mutant and wt PR structures was obtained for several of the redesigned analogues.     

Computational insights into the protonation states of catalytic dyad in BACE1–acyl guanidine based inhibitor complex

Developing small compound based drugs targeting the β-secretase (BACE) enzyme is one of the most promising strategies in treatment of the Alzheimer’s disease. As the enzyme shows the activity based on the acid-base reaction at a very narrow pH range, the protonation state of aspartic acids with the residue number 32 and 228 (Asp32 and Asp228), which forms the active site dyad, along with the protonation state of the ligand (substrate or inhibitor) play very critical role in interactions between the ligand and enzyme. Thus, understanding the nature of the protonation state of both enzyme’s active site dyad and ligand is crucial for drug design in Alzheimer’s disease field. Here we have investigated the protonation state of the Asp32 and Asp228 residues in the presence of a highly potent beta secretase inhibitor, containing acyl guanidine warhead that have recently been devised but not extensively studied. Our Quantum Mechanical, Molecular Dynamics and Docking studies on all the possible protonation states have suggested that the dyad residues are in di-deprotonated states in the presence of protonated inhibitor.     

Molecular dynamics simulations studies and free energy analysis on inhibitors of MDM2–p53 interaction

The oncoprotein MDM2 (murine double minute 2) negatively regulates the activity and stability of tumor suppressor p53. Inactivation of the MDM2–p53 interaction by potent inhibitors offers new possibilities for anticancer therapy. Molecular dynamics (MD) simulations were performed on three inhibitors–MDM2 complexes to investigate the stability and structural transitions. Simulations show that the backbone of MDM2 maintains stable during the whole time. However, slightly structural changes of inhibitors and MDM2 are observed. Furthermore, the molecular mechanics generalized Born surface area (MM-GBSA) approach was introduced to analyze the interactions between inhibitors and MDM2. The results show that binding of inhibitor pDIQ to MDM2 is significantly stronger than that of pMI and pDI to MDM2. The side chains of residues have more contribution than backbone of residues in energy decomposition. The structure–affinity analyses show that L54, I61, M62, Y67, Q72, H73 and V93 produce important interaction energy with inhibitors. The residue W/Y22′ is also very important to the interaction between the inhibitors and MDM2. The three-dimensional structures at different times indicate that the mobility of Y100 influences on the binding of inhibitors to MDM2, and its change has important role in conformations of inhibitors and MDM2.     

Conformational and energy evaluations of novel peptides binding to dengue virus envelope protein

Effective novel peptide inhibitors which targeted the domain III of the dengue envelope (E) protein by blocking dengue virus (DENV) entry into target cells, were identified. The binding affinities of these peptides towards E-protein were evaluated by using a combination of docking and explicit solvent molecular dynamics (MD) simulation methods. The interactions of these complexes were further investigated by using the Molecular Mechanics-Poisson Boltzmann Surface Area (MMPBSA) and Molecular Mechanics Generalized Born Surface Area (MMGBSA) methods. Free energy calculations of the peptides interacting with the E-protein demonstrated that van der Waals (vdW) and electrostatic interactions were the main driving forces stabilizing the complexes. Interestingly, calculated binding free energies showed good agreement with the experimental dissociation constant (K_d) values. Our results also demonstrated that specific residues might play a crucial role in the effective binding interactions. Thus, this study has demonstrated that a combination of docking and molecular dynamics simulations can accelerate the identification process of peptides as potential inhibitors of dengue virus entry into host cells.     

Study of the interaction of Huperzia saururus Lycopodium alkaloids with the acetylcholinesterase enzyme

In the present study, we describe and compare the binding modes of three Lycopodium alkaloids (sauroine, 6-hydroxylycopodine and sauroxine; isolated from Huperzia saururus) and huperzine A with the enzyme acetylcholinesterase. Refinement and rescoring of the docking poses (obtained with different programs) with an all atom force field helped to improve the quality of the protein–ligand complexes. Molecular dynamics simulations were performed to investigate the complexes and the alkaloid's binding modes. The combination of the latter two methodologies indicated that binding in the active site is favored for the active compounds. On the other hand, similar binding energies in both the active and the peripheral sites were obtained for sauroine, thus explaining its experimentally determined lack of activity. MM-GBSA predicted the order of binding energies in agreement with the experimental IC₅₀ values.     

A molecular dynamics study of the BACE1 conformational change from Apo to closed form induced by hydroxyethylamine derived compounds

BACE1 is an aspartyl protease which is a therapeutic target for Alzheimer’s disease (AD) because of its participation in the rate-limiting step in the production of Aβ-peptide, the accumulation of which produces senile plaques and, in turn, the neurodegenerative effects associated with AD. The active site of this protease is composed in part by two aspartic residues (Asp93 and Asp289). Additionally, the catalytic site has been found to be covered by an antiparallel hairpin loop called the flap. The dynamics of this flap are fundamental to the catalytic function of the enzyme. When BACE1 is inactive (Apo), the flap adopts an open conformation, allowing a substrate or inhibitor to access the active site. Subsequent interaction with the ligand induces flap closure and the stabilization of the macromolecular complex. Further, the protonation state of the aspartic dyad is affected by the chemical nature of the species entering the active site, so that appropriate selection of protonation states for the ligand and the catalytic residues will permit the elucidation of the inhibitory pathway for BACE1. In the present study, comparative analysis of different combinations of protonation states for the BACE1-hydroxyethylamine (HEA) system is reported. HEAs are potent inhibitors of BACE1 with favorable pharmacological and kinetic properties, as well as oral bioavailability. The results of Molecular Dynamics (MD) simulations and population density calculations using 8 different parameters demonstrate that the LnAsp289 configuration (HEA with a neutral amine and the Asp289 residue protonated) is the only one which permits the expected conformational change in BACE1, from apo to closed form, after flap closure. Additionally, differences in their capacities to establish and maintain interactions with residues such as Asp93, Gly95, Thr133, Asp289, Gly291, and Asn294 during this step allow differentiation among the inhibitory activities of the HEAs. The results and methodology here reported will serve to elucidate the inhibitory pathway of other families of compounds that act as BACE1 inhibitors, as well as the design of better leader compounds for the treatment of AD.     

Design and computational support for the binding stability of a new CCR5/CXCR4 dual tropic inhibitor: Computational design of a CCR5/CXCR4 drug

The human immunodeficiency virus (HIV) infects healthy human cells by binding to the glycoprotein cluster of differentiation 4 receptors on the surface of helper T-cells, along with either of two chemokine receptors, CC chemokine receptor type 5 (CCR5) or C-X-C chemokine receptor type 4 (CXCR4). Recently, a pyrazolo-piperdine ligand was synthesized and the corresponding biological data showed good binding to both chemokine receptors, effectively blocking HIV-1 entry. Here, we exhaustively assess the atomistic binding interactions of this compound with both CCR5 and CXCR4, and we find that binding is driven by π-stacking interactions between aromatic rings on the ligand and receptor residues, as well as electrostatic interactions involving the protonated piperidine nitrogen. However, these favorable binding interactions were partially offset by unfavorable desolvation of active site glutamates and aspartates, prompting our proposal of a new, synthetically-accessible derivative designed to increase the electrostatic interactions without compromising the π-stacking features.