How accurate are continuum solvation models for drug-like molecules?

Synthesis and characterization of gold(I) thiolate derivatives and bimetallic complexes for HIV inhibition

Introduction: The human immunodeficiency virus (HIV) remains a significant global health concern, with a reported high infection rate of 38.4 million cases globally; an estimated 2 million new infections and approximately 700,000 HIV/AIDS-related deaths were reported in 2021. Despite the advent of anti-retroviral therapy (ART), HIV/AIDS persists as a chronic disease. To combat this, several studies focus on developing inhibitors targeting various stages of the HIV infection cycle, including HIV-1 protease. This study aims to synthesize and characterize novel glyco diphenylphosphino metal complexes with potential HIV inhibitory properties. Method: A series of new gold(I) thiolate derivatives and three bimetallic complexes, incorporating amino phosphines and thiocarbohydrate as auxiliary ligands, were synthesized using procedures described by Jiang, et al. (2009) and Coetzee et al. (2007). Structural elucidation and purity assessment of the synthesized compounds (1-11) were conducted using micro-analysis, NMR, and infrared spectrometry. Results and Discussion: Using molecular modeling techniques, three of the metal complexes were identified as potential HIV protease inhibitors, exhibiting strong binding affinity interactions with binding pocket residues. These inhibitors demonstrated an ability to inhibit the flexibility of the flap regions of the HIV protease, similar to the known HIV protease inhibitor, darunavir. This study sheds light on the promising avenues for the development of novel therapeutic agents against HIV/AIDS.

Computed interactions of berenil with restricted foldamers of c-MYC DNA G-quadruplexes

G-quadruplexes (G4s) are secondary four-stranded DNA helical structures made up of guanine-rich nucleic acids that can assemble in the promoter regions of multiple genes under the appropriate conditions. Stabilization of G4 structures by small molecules can regulate transcription in non-telomeric regions, including in proto-oncogenes and promoter regions, contributing to anti-proliferative and anti-tumor activities. Because G4s are detectable in cancer cells but not in normal cells, they make excellent drug discovery targets. Diminazene, DMZ (or berenil), has been shown to be an efficient G-quadruplex binder. Due to the stability of the folding topology, G-quadruplex structures are frequently found in the promotor regions of oncogenes and may play a regulatory role in gene activation. Using molecular docking and molecular dynamics simulations on several different binding poses, we have studied DMZ binding toward multiple G4 topologies of the c-MYC G-quadruplex. DMZ binds preferentially to G4s that have extended loops and flanking bases. This preference arises from its interactions with the loops and the flanking nucleotides, which were not found in the structure lacking extended regions. The binding to the G4s with no extended regions instead occurred mostly through end stacking. All binding sites for DMZ were confirmed by 100 ns molecular dynamics simulations and through binding enthalpies calculated using the MM-PBSA method. The primary driving forces were electrostatic, as the cationic DMZ interacts with the anionic phosphate backbone, and through van der Waals interactions, which primarily contributed in end stacking interactions.Communicated by Ramaswamy H. Sarma.

Probing the Activation Mechanisms of Agonist DPI-287 to Delta-Opioid Receptor and Novel Agonists Using Ensemble-Based Virtual Screening with Molecular Dynamics Simulations

Pain drugs targeting mu-opioid receptors face major addiction problems that have caused an epidemic. The delta-opioid receptor (DOR) has shown to not cause addictive effects when bound to an agonist. While the active conformation of the DOR in complex with agonist DPI-287 has been recently solved, there are still no FDA-approved agonists targeting it, providing the opportunity for structure-based virtual screening. In this study, the conformational plasticity of the DOR was probed using molecular dynamics (MD) simulations, identifying two representative conformations from clustering analysis. The two MD conformations as well as the crystal conformation of DOR were used to screen novel compounds from the ZINC database (17 million compounds), in which 69 drugs were picked as potential compounds based on their docking scores. Notably, 37 out of the 69 compounds were obtained from the simulated conformations. The binding stability of the 69 compounds was further investigated using MD simulations. Based on the MM-GBSA binding energy and the predicted drug properties, eight compounds were chosen as the most favorable, six of which were from the simulated conformations. Using a dynamic network model, the communication between the crystal agonist and the top eight molecules with the receptor was analyzed to confirm if these novel compounds share a similar activation mechanism to the crystal ligand. Encouragingly, docking of these eight compounds to the other two opioid receptors (kappa and mu) suggests their good selectivity toward DOR.

Computational analysis of drug resistance of taxanes bound to human β-tubulin mutant (D26E)

The single-point mutation D26E in human β-tubulin is associated with drug resistance seen with two anti-mitotic taxanes (paclitaxel and docetaxel) when used to treat cancers. The molecular mechanism of this resistance remains elusive. However, docetaxel and a third-generation taxane, cabazitaxel, are thought to overcome this resistance. Here, structural models of both the wildtype (WT) and D26E mutant (MT) human β-tubulin were constructed based on the crystal structure of pig β-tubulin in complex with docetaxel (PDB ID: 1TUB). The three taxanes were docked into the WT and MT β-tubulin, and the resulting complexes were submitted to three independent runs of 200 ns molecular dynamic simulations, which were then averaged. MM/GBSA calculations revealed the binding energy of paclitaxel with WT and MT β-Tubulin to be -101.5 ± 8.4 and -90.4 ± 8.9 kcal/mol, respectively. The binding energy of docetaxel was estimated to be -104.7 ± 7.0 kcal/mol with the WT and -103.8 ± 5.5 kcal/mol with the MT β-tubulin. Interestingly, cabazitaxel was found to have a binding energy of -122.8 ± 10.8 kcal/mol against the WT and -106.2 ± 7.0 kcal/mol against the MT β-tubulin. These results show that paclitaxel and docetaxel bound to the MT less strongly than the WT, suggesting possible drug resistance. Similarly, cabazitaxel displayed a greater binding propensity against WT and MT β-tubulin than the other two taxanes. Furthermore, the dynamic cross-correlation matrices (DCCM) analysis suggests that the single-point mutation D26E induces a subtle dynamical difference in the ligand-binding domain. Overall, the present study revealed how the single-point mutation D26E may reduce the binding affinity of the taxanes, however, the effect of the mutation does not significantly affect the binding of cabazitaxel.

In Silico Investigation of the Molecular Mechanism of PARP1 Inhibition for the Treatment of BRCA-Deficient Cancers

The protein PARP1, which plays a crucial role in DNA repair processes, is an attractive target for cancer therapy, especially for BRCA-deficient cancers. To overcome the acquired drug resistance of PARP1, PARP1 G-quadruplex (G4) identified in the PARP1-promotor region is gaining increasing attention. Aiming to explore the molecular mechanism of PARP1 inhibition with PARP1 G4 and PARP1 as potential targets, a comparative investigation of the binding characteristics of the newly identified G4 stabilizer MTR-106, which showed modest activity against talazoparib-resistant xenograft models and the FDA-approved PARP1 inhibitor (PARPi) talazoparib, were performed through molecular simulations. Combined analyses revealed that, relative to the groove binding of talazoparib, MTR-106 induced the formation of a sandwich framework through stacking with dT1 and the capping G-pair (dG2 and dG14) of PARP1 G4 to present largely enhanced binding affinity. For the binding with PARP1, although both were located in the catalytic pocket of PARP1, MTR-106 formed more extensive interactions with the surrounding PARP1 residues compared to talazoparib, in line with its increased binding strength. Importantly, vdW interaction was recognized as a decisive factor in the bindings with PARP1 G4 and PARP1. Collectively, these findings demonstrated the ascendancy of MTR-106 over talazoparib at the atomic level and revealed that the dual targeting of PARP1 G4 and PARP1 might be pivotal for PARPi that is capable of overcoming acquired drug resistance, providing valuable information for the design and development of novel drugs.

Molecular insights into the selective binding mechanism targeting parallel human telomeric G-quadruplex

Stabilizing human telomere DNA G-quadruplex (G4) proves a promising anti-cancer strategy. Though plenty of G4 stabilizing molecules have been reported, little is known about their selective binding mechanism among various G4s. Recently, a designed monohydrazone derivative (compound 15) was reported to display specific preference in binding and stabilizing parallel human telomeric G4. To reveal the selective binding mechanism, a comparative theoretical investigation was performed on two monohydrazone derivatives (compounds 1 and 15) and three telomeric G4s showing parallel, hybrid-I, and hybrid-II conformations. Two probable binding modes, i.e. the end-stacking binding and the groove binding, were predicted by molecular dockings for each monohydrazone in its binding with the telomeric G4s. Further long-timescale molecular dynamics simulations reveal the conversion from the groove binding to the end-stacking binding for both compounds, indicating the preference of the end-stacking binding mode. Structural analysis together with binding free energy calculations show that the van der Waals interaction plays a leading role in ranking the binding affinity. By forming extensive van der Waals interactions, the parallel G4-15 binding complex shows the highest binding affinity, and the corresponding compound 15 exhibits the strongest stabilizing effect to the telomeric G4. These findings agree well with the experimental observations. Through characterizing the selective binding between monohydrazones and telomeric G4s at the atomic level, the current study provides support to the design of novel selective stabilizers targeting telomeric G4s.

To Probe the Binding Interactions between Two FDA Approved Migraine Drugs (Ubrogepant and Rimegepant) and Calcitonin-Gene Related Peptide Receptor (CGRPR) Using Molecular Dynamics Simulations

Recently, the FDA approved ubrogepant and rimegepant as oral drugs to treat migraines by targeting the calcitonin-gene related peptide receptor (CGRPR). Unfortunately, there is no high-resolution complex structure with these two drugs; thus the detailed interaction between drugs and the receptor remains elusive. This study uses molecular docking and molecular dynamics simulation to model the drug-receptor complex and analyze their binding interactions at a molecular level. The complex crystal structure (3N7R) of the gepant drugs' predecessor, olcegepant, was used for our molecular docking of the two drugs and served as a control system. The three systems, with ubrogepant, rimegepant, and crystal olcegepant, were subject to 3 × 1000 ns molecular dynamics simulations and followed by the simulation interaction diagram (SID), structural clustering, and MM-GBSA binding energy analyses. Our MD data revealed that olcegepant binds most strongly to the CGRPR, followed by ubrogepant and then rimegepant, largely due to changes in hydrophobic and electrostatic interactions. The order of our MM-GBSA binding energies of these three compounds is consistent with their experimental IC₅₀ values. SID analysis revealed the pharmacophore of the gepant class to be the dihydroquinazolinone group derivative. Subtle differences in interaction profile have been noted, including interactions with the W74 and W72 residues. The ubrogepant and rimegepant both contact A70 and M42 of the receptor, while olcegepant does not. The results of this study elucidate the interactions in the binding pocket of CGRP receptor and can assist in further development for orally available antagonists of the CGRP receptor.

Optimized Halogen Atomic Radii for PBSA Calculations Using Off-Center Point Charges

In force-field methods, the usage of off-center point charges, also called extra points (EPs), is a common strategy to tackle the anisotropy of the electrostatic potential of covalently bonded halogens (X), thus allowing the description of halogen bonds (XBs) at the molecular mechanics/molecular dynamics (MM/MD) level. Diverse EP implementations exist in the literature differing on the charge sets and/or the X-EP distances. Poisson-Boltzmann and surface area (PBSA) calculations can be used to obtain solvation free energies (ΔG_solv) of small molecules, often to compute binding free energies (ΔG_bind) at the MM-PBSA level. This method depends, among other parameters, on the empirical assignment of atomic radii (PB radii). Given the multiplicity of off-center point-charge models and the lack of specific PB radii for halogens compatible with such implementations, in this work, we assessed the performance of PBSA calculations for the estimation of ΔG_solv values in water (ΔG_hyd), also conducting an optimization of the halogen PB radii (Cl, Br, and I) for each EP model. We not only expand the usage of EP models in the scope of the general AMBER force field (GAFF) but also provide the first optimized halogen PB radii in the context of the CHARMM general force field (CGenFF), thus contributing to improving the description of halogenated compounds in PBSA calculations.

Analysis of natural compounds against the activity of SARS-CoV-2 NSP15 protein towards an effective treatment against COVID-19: a theoretical and computational biology approach

Coronavirus infectious disease 2019 (COVID-19), a viral infection caused by a novel coronavirus (nCoV), continues to emerge as a serious threat to public health. This pandemic caused by SARS-CoV-2 (severe acute respiratory syndrome-coronavirus-2) has infected globally with 1,550,000 plus deaths to date, representing a high risk to public health. No effective drug or vaccine is available to curb down this deadly virus. The expedition for searching for a potential drug or vaccine against COVID-19 is of massive potential and favour to the community. This study is focused on finding an effective natural compound that can be processed further into a potential inhibitor to check the activity of SARS-CoV-2 with minimal side effects targeting NSP15 protein, which belongs to the EndoU enzyme family. The natural screening suggested two efficient compounds (PubChem ID: 95372568 and 1776037) with dihydroxyphenyl region of the compound, found to be important in the interaction with the viral protein showing promising activity which may act as a potent lead inhibitory molecule against the virus. In combination with virtual screening, modelling, drug likeliness, molecular docking, and 500 ns cumulative molecular dynamics simulations (100 ns for each complex) along with the decomposition analysis to calculate and confirm the stability and fold, we propose 95372568 and 1776037 as novel compounds of natural origin capable of getting developed into potent lead molecules against SARS-CoV-2 target protein NSP15.

An in silico investigation of the binding modes and pathway of APTO-253 on c-KIT G-quadruplex DNA

The stability of c-KIT G-quadruplex DNA via ligands has been a significant concern in the growing field of cancer therapy. Thus, it is very important to understand the mechanism behind the high binding affinity of the small drug molecules on the c-KIT G-quadruplex DNA. In this study, we have investigated the binding mode and pathway of the APTO-253 ligand on the c-KIT G-quadruplex DNA employing a total of 10 μs all atom molecular dynamics simulations and further 8.82 μs simulations via the umbrella sampling method using both OL15 and BSC1 latest force fields for DNA structures. From the cluster structure analysis, mainly three binding pathways i.e., top, bottom and side loop stacking modes are identified. Moreover, RMSD, RMSF and 2D-RMSD values indicate that the c-KIT G-quadruplex DNA and APTO-253 molecules are stable throughout the simulation run. Furthermore, the number of hydrogen bonds in each tetrad and the distance between the two central K⁺ cations confirm that the c-KIT G-quadruplex DNA maintains its conformation in the process of complex formation with the APTO-253 ligand. The binding free energies and the minimum values in the potential of mean forces suggest that the binding processes are energetically favorable. Furthermore, we have found that the bottom stacking mode is the most favorable binding mode among all the three modes for the OL15 force field. However, for the BSC1 force field, both the top and bottom binding modes of the APTO-253 ligand in c-KIT G-quadruplex DNA are comparable to each other. To investigate the driving force for the complex formation, we have noticed that the van der Waals (vdW) and π-π stacking interactions are mainly responsible. Our detailed studies provide useful information for the discovery of novel drugs in the field of stabilization of G-quadruplex DNAs.

In silico approach: biological prediction of nordentatin derivatives as anticancer agent inhibitors in the cAMP pathway

A combination of computational techniques has been carried out to predict the binding of nordentatin derivatives based on pyranocoumarin semi-synthesis with the target protein from the expression of the PDE4B gene. The inhibition of the cAMP pathway is the main target of anti-cancer drugs, which is responsible for uncontrolled cell division in cancer. Modeling was done using a combination of semi-empirical methods and the density functional theory (PM3-DFT/6-31G*/B3LYP) to obtain the optimal structure of a small ligand that could be modeled. Studies on the interaction of the ligands and amino acid residues on protein targets were carried out using a combination of molecular docking and molecular dynamic simulation. Molecular docking based on functional grid scores showed a very good native ligand pose with an RMSD of 0.93 Å in determining the initial coordinates of the ligand-receptor interactions. Furthermore, the amino acid residues responsible for interaction through H-bonds were Tyr103, His104, His177, Met217, and Gln313. The binding free energy (kcal mol-1) results of the candidates were PS-1 (-36.84 ± 0.31), PS-2 (-35.34 ± 0.28), PS-3 (-26.65 ± 0.30), PS-5 (-42.66 ± 0.26), PS-7 (-35.33 ± 0.23), and PS-9 (-32.57 ± 0.20), which are smaller than that of the native ligand Z72 (-24.20 ± 0.19), and thus these have good potential as drugs that can inhibit the cAMP pathway. These results provide theoretical information for the efficient inhibition of the cAMP pathway in the future.

To probe the binding pathway of a selective compound (D089-0563) to c-MYC Pu24 G-quadruplex using free ligand binding simulations and Markov state model analysis

D089-0563 is a highly promising anti-cancer compound that selectively binds the transcription-silencing G-quadruplex element (Pu27) at the promoter region of the human c-MYC oncogene; however, its binding mechanism remains elusive.

Assessing the performance of the MM/PBSA and MM/GBSA methods. 10. Impacts of enhanced sampling and variable dielectric model on protein-protein Interactions

Enhanced sampling has been extensively used to capture the conformational transitions in protein folding, but it attracts much less attention in the studies of protein-protein recognition. In this study, we evaluated the impact of enhanced sampling methods and solute dielectric constants on the overall accuracy of the molecular mechanics/Poisson-Boltzmann surface area (MM/PBSA) and molecular mechanics/generalized Born surface area (MM/GBSA) approaches for the protein-protein binding free energy calculations. Here, two widely used enhanced sampling methods, including aMD and GaMD, and conventional molecular dynamics (cMD) simulations with two AMBER force fields (ff03 and ff14SB) were used to sample the conformations for 21 protein-protein complexes. The MM/PBSA and MM/GBSA calculation results illustrate that the standard MM/GBSA based on the cMD simulations yields the best Pearson correlation (r_p = -0.523) between the predicted binding affinities and the experimental data, which is much higher than that given by MM/PBSA (r_p = -0.212). Two enhanced sampling methods (aMD and GaMD) are indeed more efficient for conformational sampling, but they did not improve the binding affinity predictions for protein-protein systems, suggesting that the aMD or GaMD sampling (at least in short timescale simulations) may not be a good choice for the MM/PBSA and MM/GBSA predictions of protein-protein complexes. The solute dielectric constant of 1.0 is recommended to MM/GBSA, but a higher solute dielectric constant is recommended to MM/PBSA, especially for the systems with higher polarity on the protein-protein binding interfaces. Then, a preliminary assessment of the MM/GBSA calculations based on a variable dielectric generalized Born (VDGB) model was conducted. The results highlight the potential power of VDGB in the free energy predictions for protein-protein systems, but more thorough studies should be done in the future.

Tackling Halogenated Species with PBSA: Effect of Emulating the σ-Hole

To model halogen-bond phenomena using classical force fields, an extra point (EP) of charge is frequently introduced at a given distance from the halogen (X) to emulate the σ-hole. The resulting molecular dynamics (MD) trajectories can be used in subsequent molecular mechanics (MM) combined with Poisson-Boltzmann and surface area calculations (PBSA) to estimate protein-ligand binding free energies (Δ G_bind). While EP addition improves the MM/MD description of halogen-containing systems, its effect on the calculation of solvation free energies (Δ G_solv) using the PBSA approach is yet to be assessed. As the PBSA calculations depend, among other parameters, on the empirical assignment of radii (PB radii), a problematic issue arises, since standard halogen radii are smaller than the typical X···EP distances, thus placing the EP within the solvent dielectric. Herein, we took a common literature EP parametrization scheme, which uses X···EP = R_min and RESP charges in the context of GAFF, and performed a comprehensive study on the performance of PBSA (using three different setups) in the calculation of Δ G_solv values for 142 halogenated compounds (bearing Cl, Br, or I) for which the experimental values are known. By conducting an optimization (minimizing the error against experimental values), we provide a new optimized set of halogen PB radii, for each PBSA setup, that should be used in the context of the aforementioned scenario. A simultaneous optimization of PB radii and X···EP distances shows that a wide range of distance/radius pairs can be used without significant loss of accuracy, therefore laying the basis for expanding this halogen radii optimization strategy to other force fields and EP implementations. As ligand Δ G_solv estimation is an important term in the determination of protein-ligand Δ G_bind, this work is particularly relevant in the framework of structure-based virtual screening and related computer-aided drug design routines.

Binding of BRACO19 to a Telomeric G-Quadruplex DNA Probed by All-Atom Molecular Dynamics Simulations with Explicit Solvent

Although BRACO19 is a potent G-quadruplex binder, its potential for clinical usage is hindered by its low selectivity towards DNA G-quadruplex over duplex. High-resolution structures of BRACO19 in complex with neither single-stranded telomeric DNA G-quadruplexes nor B-DNA duplex are available. In this study, the binding pathway of BRACO19 was probed by 27.5 µs molecular dynamics binding simulations with a free ligand (BRACO19) to a DNA duplex and three different topological folds of the human telomeric DNA G-quadruplex (parallel, anti-parallel and hybrid). The most stable binding modes were identified as end stacking and groove binding for the DNA G-quadruplexes and duplex, respectively. Among the three G-quadruplex topologies, the MM-GBSA binding energy analysis suggested that BRACO19′s binding to the parallel scaffold was most energetically favorable. The two lines of conflicting evidence plus our binding energy data suggest conformation-selection mechanism: the relative population shift of three scaffolds upon BRACO19 binding (i.e., an increase of population of parallel scaffold, a decrease of populations of antiparallel and/or hybrid scaffold). This hypothesis appears to be consistent with the fact that BRACO19 was specifically designed based on the structural requirements of the parallel scaffold and has since proven effective against a variety of cancer cell lines as well as toward a number of scaffolds. In addition, this binding mode is only slightly more favorable than BRACO19s binding to the duplex, explaining the low binding selectivity of BRACO19 to G-quadruplexes over duplex DNA. Our detailed analysis suggests that BRACO19′s groove binding mode may not be stable enough to maintain a prolonged binding event and that the groove binding mode may function as an intermediate state preceding a more energetically favorable end stacking pose; base flipping played an important role in enhancing binding interactions, an integral feature of an induced fit binding mechanism.

Binding of Telomestatin, TMPyP4, BSU6037, and BRACO19 to a Telomeric G-Quadruplex-Duplex Hybrid Probed by All-Atom Molecular Dynamics Simulations with Explicit Solvent

A promising anticancer therapeutic strategy is the stabilization of telomeric G-quadruplexes using G-quadruplex-binding small molecules. Although many G-quadruplex-specific ligands have been developed, their low potency and selectivity to G-quadruplexes over duplex remains unsolved. Recently, a crystal structure of a telomeric 3' quadruplex-duplex hybrid was reported and the quadruplex-duplex interface was suggested to a good target to address the issues. However, there are no high-resolution complex structures reported for G-quadruplex ligands except for a docked BSU6037. In this study, molecular dynamic (MD) binding simulations with a free ligand were used to study binding poses and dynamics of four representative ligands: telomestatin, TMPyP4, BSU6037, and BRACO19. The MD data showed that BSU6037 was able to fully intercalate into the interface whereas TMPyP4 and BRACO19 could only maintain partial intercalation into the interface and telomestatin only binds at the quadruplex and duplex ends. Both linear ligands, BSU6037 and BRACO19, were able to interact with the interface, yet they were not selective over duplex DNA. The DNA geometry, binding modes, and binding pathways were systematically characterized, and the binding energy was calculated and compared for each system. The interaction of the ligands to the interface was by the means of an induced-fit binding mechanism rather than a lock-key mechanism, consisting of the DNA unfolding at the interface to allow entrance of the drug and then the refolding and repacking of the DNA and the ligand to further stabilize the G-quadruplex. On the basis of the findings in this study, modifications were suggested to optimize the interface binding for TMPyp4 and telomestatin.

Elucidating the multi-targeted anti-amyloid activity and enhanced islet amyloid polypeptide binding of β-wrapins

β-wrapins are engineered binding proteins stabilizing the β-hairpin conformations of amyloidogenic proteins islet amyloid polypeptide (IAPP), amyloid-β, and α-synuclein, thus inhibiting their amyloid propensity. Here, we use computational and experimental methods to investigate the molecular recognition of IAPP by β-wrapins. We show that the multi-targeted, IAPP, amyloid-β, and α-synuclein, binding properties of β-wrapins originate mainly from optimized interactions between β-wrapin residues and sets of residues in the three amyloidogenic proteins with similar physicochemical properties. Our results suggest that IAPP is a comparatively promiscuous β-wrapin target, probably due to the low number of charged residues in the IAPP β-hairpin motif. The sub-micromolar affinity of β-wrapin HI18, specifically selected against IAPP, is achieved in part by salt-bridge formation between HI18 residue Glu10 and the IAPP N-terminal residue Lys1, both located in the flexible N-termini of the interacting proteins. Our findings provide insights towards developing novel protein-based single- or multi-targeted therapeutics.

An attempt to incorporate effect of direct interaction between a ligand and explicit water molecules into MM/3D-RISM

Endpoint methods using continuum-solvent models are widely used to estimate protein-ligand affinity. A recently developed method, MM/3D-RISM, estimates the solvation energy using statistical mechanics by 3D-RISM. This method is theoretically expected to accurately describe solvation effects and to also be less dependent on protein-ligand systems. In this study, we examined the performance of MM/3D-RISM for a set of α-thrombin inhibitors with a non-congeneric series of ligands, containing three diverse chemical scaffolds. The standard MM/3D-RISM showed a weak correlation (R²  = 0.191) but correctly estimated affinity for two of the three scaffolds. However, the simplest inhibitor, benzamidine, was not ranked appropriately. From visual inspection of inhibitor-binding modes, an attempt was made to incorporate the direct interaction between a ligand and water molecules into MM/3D-RISM. A model (Model-1) dealing with directly interacting water molecules (Wat) as an independent component of a protein (R)-ligand (L) complex-formation, that is, R + L + Wat → R-L-Wat, showed a better linearity (R²  = 0.422) than that of the standard MM/3D-RISM model and achieved a good ranking of all three scaffolds of α-thrombin inhibitors. Additionally, an attempt was made to model avidin-biotin system with a congeneric series of inhibitors, and results showed that both the standard MM/3D-RISM model (R²  = 0.839) and Model-1 (R²  = 0.695) satisfactorily estimated the affinity.

To probe interaction of morphine and IBNtxA with 7TM and 6TM variants of the human μ-opioid receptor using all-atom molecular dynamics simulations with an explicit membrane

IBNtxA, a morphine derivative, is 10-fold more potent and has a better safety profile than morphine. Animal studies indicate that the analgesic effect of IBNtxA appears to be mediated by the activation of truncated splice variants (6TM) of the Mu opioid receptor (MOR-1) where transmembrane helix 1 (TM1) is removed. Interestingly, morphine is unable to activate 6TM variants. To date, a high resolution structure of 6TM variants is missing, and the interaction of 6TM variants with IBNtxA and morphine remains elusive. In this study we used homology modeling, docking and molecular dynamics (MD) simulations to study a representative 6TM variant (G1) and a full-length 7TM variant of human MOR-1 in complex with IBNtxA and morphine respectively. The structural models of human G1 and 7TM were obtained by homology modeling using the X-ray solved crystal structure of the active mouse 7TM bound to an agonist BU72 (PDB id: ) as the template. Our 6000 ns MD data show that either TM1 truncation (i.e. from 7TM to 6TM) or ligand modification (i.e. from morphine to IBNtxA) alone causes the loss of key morphine-7TM interactions that are well-known to be required for MOR-1 activation. Receptor disruptions are mainly located at TMs 2, 3, 6 and 7 in comparison with the active crystal complex. However, when both perturbations occur in the 6TM-IBNtxA complex, the key ligand-receptor interactions and the receptor conformation are recovered to resemble those in the active 7TM-morphine complex. Our molecular switch analysis further explains well why morphine is not able to activate 6TM variants. The close resemblance between 6TM-IBTtxA and 7TM in complex with PZM21, a G-protein biased 7TM agonist, suggests the possible biased agonism of IBNtxA on G1, which is consistent with its reduced side effects.

Binding modes and pathway of RHPS4 to human telomeric G-quadruplex and duplex DNA probed by all-atom molecular dynamics simulations with explicit solvent

RHPS4, a potent binder to human telomeric DNA G-quadruplex, shows high efficacy in tumor cell growth inhibition. However, it's preferential binding to DNA G-quadruplex over DNA duplex (about 10 fold) remains to be improved toward its clinical application. A high resolution structure of the single-stranded telomeric DNA G-quadruplexes, or B-DNA duplex, in complex with RHPS4 is not available yet, and the binding nature of this ligand to these DNA forms remains to be elusive. In this study, we carried out 40 μs molecular dynamics binding simulations with a free ligand to decipher the binding pathway of RHPS4 to a DNA duplex and three G-quadruplex folders (parallel, antiparallel and hybrid) of the human telomeric DNA sequence. The most stable binding mode identified for the duplex, parallel, antiparallel and hybrid G-quadruplexes is an intercalation, bottom stacking, top intercalation and bottom intercalation mode, respectively. The intercalation mode with similar binding strength to both the duplex and the G-quadruplexes, explains the lack of binding selectivity of RHPS4 to the G-quadruplex form. Therefore, a ligand modification that destabilizes the duplex intercalation mode but stabilizes the G-quadruplex intercalation mode will improve the binding selectivity toward G-quadruplex. The intercalation mode of RHPS4 to both the duplex and the antiparallel and the hybrid G-quadruplex follows a base flipping-insertion mechanism rather than an open-insertion mechanism. The groove binding, the side binding and the intercalation with flipping out of base were observed to be intermediate states before the full intercalation state with paired bases.

Noncovalent Interactions in Specific Recognition Motifs of Protein-DNA Complexes

In view of the importance of protein-DNA interactions in biological processes, we extracted from the Protein Data Bank several one-to-one complexes of amino acids with nucleotides that matched certain geometric and energetic specificity criteria and investigated them using quantum chemistry methods. The CCSD(T)/CBS interaction energies were used as a benchmark to compare the performance of the MP2.5, MP2-F12, DFT-D3, and PM6-D3H4 methods. All methods yielded good agreement with the reference values, with declining accuracy from MP2.5 to PM6-D3H4. Regardless of the site of interaction, the minima found after full optimization in implicit solvent with high dielectric constant were close to the structures experimentally detected in protein-DNA complexes. According to DFT-SAPT analysis, the nature of noncovalent interactions strongly depends on the type of amino acid. The negatively charged sugar-phosphate backbone of DNA heavily influences the strength of interactions and must be included in the computational model, especially in the case of interactions with charged amino acids.

Development of a novel series of non-natural triaryl agonists and antagonists of the Pseudomonas aeruginosa LasR quorum sensing receptor

Bacterial chemical communication, through a process called quorum sensing (QS), plays a central role in infection in numerous bacterial pathogens. Quorum sensing in Pseudomonas aeruginosa employs a series of small molecule receptors including the master QS regulator, LasR. In this study we investigate a non-natural triaryl series of LasR ligands using a combination of structure activity relationship studies and computational modeling. These studies have enabled the identification of key structural requirements for ligand binding and have revealed a new strategy for inducing the therapeutically relevant antagonism of LasR.

dMM-PBSA: A New HADDOCK Scoring Function for Protein-Peptide Docking

Molecular-docking programs coupled with suitable scoring functions are now established and very useful tools enabling computational chemists to rapidly screen large chemical databases and thereby to identify promising candidate compounds for further experimental processing. In a broader scenario, predicting binding affinity is one of the most critical and challenging components of computer-aided structure-based drug design. The development of a molecular docking scoring function which in principle could combine both features, namely ranking putative poses and predicting complex affinity, would be of paramount importance. Here, we systematically investigated the performance of the MM-PBSA approach, using two different Poisson-Boltzmann solvers (APBS and DelPhi), in the currently rising field of protein-peptide interactions (PPIs), identifying the correct binding conformations of 19 different protein-peptide complexes and predicting their binding free energies. First, we scored the decoy structures from HADDOCK calculation via the MM-PBSA approach in order to assess the capability of retrieving near-native poses in the best-scoring clusters and of evaluating the corresponding free energies of binding. MM-PBSA behaves well in finding the poses corresponding to the lowest binding free energy, however the built-in HADDOCK score shows a better performance. In order to improve the MM-PBSA-based scoring function, we dampened the MM-PBSA solvation and coulombic terms by 0.2, as proposed in the HADDOCK score and LIE approaches. The new dampened MM-PBSA (dMM-PBSA) outperforms the original MM-PBSA and ranks the decoys structures as the HADDOCK score does. Second, we found a good correlation between the dMM-PBSA and HADDOCK scores for the near-native clusters of each system and the experimental binding energies, respectively. Therefore, we propose a new scoring function, dMM-PBSA, to be used together with the built-in HADDOCK score in the context of protein-peptide docking simulations.

Mechanism of C-Terminal Fragments of Amyloid β-Protein as Aβ Inhibitors: Do C-Terminal Interactions Play a Key Role in Their Inhibitory Activity?

Targeting the early oligomerization of amyloid β protein (Aβ) is a promising therapeutic strategy for Alzheimer's disease (AD). Recently, certain C-terminal fragments (CTFs) derived from Aβ42 were shown to be potent inhibitors of Aβ-induced toxicity. The shortest peptide studied, Aβ(39-42), has been shown to modulate Aβ oligomerization and inhibit Aβ toxicity. Understanding the mechanism of these CTFs, especially Aβ(39-42), is of significance for future therapeutic development of AD and peptidomimetic-based drug development. Here we used ion mobility spectrometry-mass spectrometry to investigate the interactions between two modified Aβ(39-42) derivatives, VVIA-NH2 and Ac-VVIA, and full-length Aβ42. VVIA-NH2 was previously shown to inhibit Aβ toxicity, whereas Ac-VVIA did not. Our mass spectrometry analysis revealed that VVIA-NH2 binds directly to Aβ42 monomer and small oligomers while Ac-VVIA binds only to Aβ42 monomer. Ion mobility studies showed that VVIA-NH2 modulates Aβ42 oligomerization by not only inhibiting the dodecamer formation but also disaggregating preformed Aβ42 dodecamer. Ac-VVIA also inhibits and removes preformed Aβ42 dodecamer. However, the Aβ42 sample with the addition of Ac-VVIA clogged the nanospray tip easily, indicating that larger aggregates are formed in the solution in the presence of Ac-VVIA. Molecular dynamics simulations suggested that VVIA-NH2 binds specifically to the C-terminal region of Aβ42 while Ac-VVIA binds dispersedly to multiple regions of Aβ42. This work implies that C-terminal interactions and binding to Aβ oligomers are important for C-terminal fragment inhibitors.

The MM/PBSA and MM/GBSA methods to estimate ligand-binding affinities

INTRODUCTION

The molecular mechanics energies combined with the Poisson-Boltzmann or generalized Born and surface area continuum solvation (MM/PBSA and MM/GBSA) methods are popular approaches to estimate the free energy of the binding of small ligands to biological macromolecules. They are typically based on molecular dynamics simulations of the receptor-ligand complex and are therefore intermediate in both accuracy and computational effort between empirical scoring and strict alchemical perturbation methods. They have been applied to a large number of systems with varying success.

AREAS COVERED

The authors review the use of MM/PBSA and MM/GBSA methods to calculate ligand-binding affinities, with an emphasis on calibration, testing and validation, as well as attempts to improve the methods, rather than on specific applications.

EXPERT OPINION

MM/PBSA and MM/GBSA are attractive approaches owing to their modular nature and that they do not require calculations on a training set. They have been used successfully to reproduce and rationalize experimental findings and to improve the results of virtual screening and docking. However, they contain several crude and questionable approximations, for example, the lack of conformational entropy and information about the number and free energy of water molecules in the binding site. Moreover, there are many variants of the method and their performance varies strongly with the tested system. Likewise, most attempts to ameliorate the methods with more accurate approaches, for example, quantum-mechanical calculations, polarizable force fields or improved solvation have deteriorated the results.

g_mmpbsa--a GROMACS tool for high-throughput MM-PBSA calculations

Molecular mechanics Poisson-Boltzmann surface area (MM-PBSA), a method to estimate interaction free energies, has been increasingly used in the study of biomolecular interactions. Recently, this method has also been applied as a scoring function in computational drug design. Here a new tool g_mmpbsa, which implements the MM-PBSA approach using subroutines written in-house or sourced from the GROMACS and APBS packages is described. g_mmpbsa was developed as part of the Open Source Drug Discovery (OSDD) consortium. Its aim is to integrate high-throughput molecular dynamics (MD) simulations with binding energy calculations. The tool provides options to select alternative atomic radii and different nonpolar solvation models including models based on the solvent accessible surface area (SASA), solvent accessible volume (SAV), and a model which contains both repulsive (SASA-SAV) and attractive components (described using a Weeks-Chandler-Andersen like integral method). We showcase the effectiveness of the tool by comparing the calculated interaction energy of 37 structurally diverse HIV-1 protease inhibitor complexes with their experimental binding free energies. The effect of varying several combinations of input parameters such as atomic radii, dielectric constant, grid resolution, solute-solvent dielectric boundary definition, and nonpolar models was investigated. g_mmpbsa can also be used to estimate the energy contribution per residue to the binding energy. It has been used to identify those residues in HIV-1 protease that are most critical for binding a range of inhibitors.

Biological applications of classical electrostatics methods

Continuum electrostatics modeling of solvation based on the Poisson–Boltzmann (PB) equation has gained wide acceptance in biomolecular applications such as energetic analysis and structural visualization. Successful application of the PB solvent models requires careful calibration of the solvation parameters. Extensive testing and validation is also important to ensure accuracy in their applications. Limitation in the continuum modeling of solvation is also a known issue in certain biomolecular applications. Growing interest in membrane systems has further spurred developmental efforts to allow inclusion of membrane in the PB solvent models. Despite their past successes due to careful parameterization, algorithm development and parallel implementation, there is still much to be done to improve their transferability from the small molecular systems upon which they were developed and validated to complex macromolecular systems as advances in technology continue to push forward, providing ever greater computational resources to researchers to study more interesting biological systems of higher complexity.

Structural basis of redox-dependent modulation of galectin-1 dynamics and function

Galectin-1 (Gal-1), a member of a family of multifunctional lectins, plays key roles in diverse biological processes including cell signaling, immunomodulation, neuroprotection and angiogenesis. The presence of an unusual number of six cysteine residues within Gal-1 sequence prompted a detailed analysis of the impact of the redox environment on the functional activity of this lectin. We examined the role of each cysteine residue in the structure and function of Gal-1 using both experimental and computational approaches. Our results show that: (i) only three cysteine residues present in each carbohydrate recognition domain (CRD) (Cys2, Cys16 and Cys88) were important in protein oxidation, (ii) oxidation promoted the formation of the Cys16-Cys88 disulfide bond, as well as multimers through Cys2, (iii) the oxidized protein did not bind to lactose, probably due to poor interactions with Arg48 and Glu71, (iv) in vitro oxidation by air was completely reversible and (v) oxidation by hydrogen peroxide was relatively slow (1.7 ± 0.2 M(-1) s(-1) at pH 7.4 and 25°C). Finally, an analysis of key cysteines in other human galectins is also provided in order to predict their behaviour in response to redox variations. Collectively, our data provide new insights into the structural basis of Gal-1 redox regulation with critical implications in physiology and pathology.

Structural similarities and differences between amyloidogenic and non-amyloidogenic islet amyloid polypeptide (IAPP) sequences and implications for the dual physiological and pathological activities of these peptides

IAPP, a 37 amino-acid peptide hormone belonging to the calcitonin family, is an intrinsically disordered protein that is coexpressed and cosecreted along with insulin by pancreatic islet β-cells in response to meals. IAPP plays a physiological role in glucose regulation; however, in certain species, IAPP can aggregate and this process is linked to β-cell death and Type II Diabetes. Using replica exchange molecular dynamics with extensive sampling (16 replicas per sequence and 600 ns per replica), we investigate the structure of the monomeric state of two species of aggregating peptides (human and cat IAPP) and two species of non-aggregating peptides (pig and rat IAPP). Our simulations reveal that the pig and rat conformations are very similar, and consist of helix-coil and helix-hairpin conformations. The aggregating sequences, on the other hand, populate the same helix-coil and helix-hairpin conformations as the non-aggregating sequence, but, in addition, populate a hairpin structure. Our exhaustive simulations, coupled with available peptide-activity data, leads us to a structure-activity relationship (SAR) in which we propose that the functional role of IAPP is carried out by the helix-coil conformation, a structure common to both aggregating and non-aggregating species. The pathological role of this peptide may have multiple origins, including the interaction of the helical elements with membranes. Nonetheless, our simulations suggest that the hairpin structure, only observed in the aggregating species, might be linked to the pathological role of this peptide, either as a direct precursor to amyloid fibrils, or as part of a cylindrin type of toxic oligomer. We further propose that the helix-hairpin fold is also a possible aggregation prone conformation that would lead normally non-aggregating variants of IAPP to form fibrils under conditions where an external perturbation is applied. The SAR relationship is used to suggest the rational design of therapeutics for treating diabetes.

IAPP, a 37 amino-acid peptide hormone belonging to the calcitonin family, is an intrinsically disordered peptide produced along with insulin by pancreatic islet β-cells in response to meals. In its functional form, IAPP acts as a synergic partner of insulin to reduce blood glucose. IAPP can, however, also play a pathological role, contributing to Type II diabetes (T2D). Knowledge of the structural nature of the physiological and pathological forms of IAPP will facilitate the rational design of novel drugs for therapeutic treatment of T2D. However, because IAPP does not fold to a single structure, but rather co-exists between multiple functional (and toxic) structures, it is extremely challenging for experimental methods to gain detailed structural information. Using a computational approach, we were able to obtain detailed structures of four IAPP variants and propose a novel structural hypothesis for the two opposing roles of this peptide.

Modelling and enhanced molecular dynamics to steer structure-based drug discovery

The ever-increasing gap between the availabilities of the genome sequences and the crystal structures of proteins remains one of the significant challenges to the modern drug discovery efforts. The knowledge of structure-dynamics-functionalities of proteins is important in order to understand several key aspects of structure-based drug discovery, such as drug-protein interactions, drug binding and unbinding mechanisms and protein-protein interactions. This review presents a brief overview on the different state of the art computational approaches that are applied for protein structure modelling and molecular dynamics simulations of biological systems. We give an essence of how different enhanced sampling molecular dynamics approaches, together with regular molecular dynamics methods, assist in steering the structure based drug discovery processes.

Assessing the accuracy and performance of implicit solvent models for drug molecules: conformational ensemble approaches

The accuracy and performance of implicit solvent methods for solvation free energy calculations were assessed on a set of 20 neutral drug molecules. Molecular dynamics (MD) provided ensembles of conformations in water and water-saturated octanol. The solvation free energies were calculated by popular implicit solvent models based on quantum mechanical (QM) electronic densities (COSMO-RS, MST, SMD) as well as on molecular mechanical (MM) point-charge models (GB, PB). The performance of the implicit models was tested by a comparison with experimental water-octanol transfer free energies (ΔG(ow)) by using single- and multiconformation approaches. MD simulations revealed difficulties in a priori estimation of the flexibility features of the solutes from simple structural descriptors, such as the number of rotatable bonds. An increasing accuracy of the calculated ΔG(ow) was observed in the following order: GB1 ~ PB < GB7 ≪ MST < SMD ~ COSMO-RS with a clear distinction identified between MM- and QM-based models, although for the set excluding three largest molecules, the differences among COSMO-RS, MST, and SMD were negligible. It was shown that the single-conformation approach applied to crystal geometries provides a rather accurate estimate of ΔG(ow) for rigid molecules yet fails completely for the flexible ones. The multiconformation approaches improved the performance, but only when the deformation contribution was ignored. It was revealed that for large-scale calculations on small molecules a recent GB model, GB7, provided a reasonable accuracy/speed ratio. In conclusion, the study contributes to the understanding of solvation free energy calculations for physical and medicinal chemistry applications.

Early stage intercalation of doxorubicin to DNA fragments observed in molecular dynamics binding simulations

The intercalation mode between doxorubicin (an anticancer drug) and two 6-base-pair DNA model fragments (d(CGATCG)₂ and d(CGTACG)₂) has been well studied by X-ray crystallography and NMR experimental methods. Yet, the detailed intercalation pathway at molecular level remains elusive. In this study, we conducted molecular dynamics binding simulations of these two systems using AMBER DNA (parmbsc0) and drug (GAFF) force fields starting from the unbound state. We observed outside binding (minor groove binding or end-binding) in all six independent binding simulations (three for each DNA fragment), followed by the complete intercalation of a drug molecule in two simulations (one for each DNA fragment). First, our data directly supported that the minor groove binding is the dominant pre-intercalation step. Second, we observed that the opening and flipping of a local base pair (A3-T10 for d(CGATCG)₂ and C1-G12 for d(CGTACG)₂) in the two intercalation trajectories. This locally cooperative flipping-intercalation mechanism was different from the previously proposed rise-insertion mechanism by which the distance between two neighboring intact base pairs increases to create a space for the drug insertion. Third, our simulations provided the first set of data to support the applicability of the AMBER DNA and drug force fields in drug-DNA atomistic binding simulations. Implications on the kinetics pathway and drug action are also discussed.

Comparison of end-point continuum-solvation methods for the calculation of protein-ligand binding free energies

We have compared the predictions of ligand-binding affinities from several methods based on end-point molecular dynamics simulations and continuum solvation, that is, methods related to MM/PBSA (molecular mechanics combined with Poisson-Boltzmann and surface area solvation). Two continuum-solvation models were considered, viz., the Poisson-Boltzmann (PB) and generalised Born (GB) approaches. The nonelectrostatic energies were also obtained in two different ways, viz., either from the sum of the bonded, van der Waals, nonpolar solvation energies, and entropy terms (as in MM/PBSA), or from the scaled protein-ligand van der Waals interaction energy (as in the linear interaction energy approach, LIE). Three different approaches to calculate electrostatic energies were tested, viz., the sum of electrostatic interaction energies and polar solvation energies, obtained either from a single simulation of the complex or from three independent simulations of the complex, the free protein, and the free ligand, or the linear-response approximation (LRA). Moreover, we investigated the effect of scaling the electrostatic interactions by an effective internal dielectric constant of the protein (ϵ(int) ). All these methods were tested on the binding of seven biotin analogues to avidin and nine 3-amidinobenzyl-1H-indole-2-carboxamide inhibitors to factor Xa. For avidin, the best results were obtained with a combination of the LIE nonelectrostatic energies with the MM+GB electrostatic energies from a single simulation, using ϵ(int) = 4. For fXa, standard MM/GBSA, based on one simulation and using ϵ(int) = 4-10 gave the best result. The optimum internal dielectric constant seems to be slightly higher with PB than with GB solvation.