Binding and selectivity studies of phosphatidylinositol 3-kinase (PI3K) inhibitors

Overexpression of the Phosphatidylinositol 3-kinase (PI3K) proteins have been observed in cancer cells. Targeting the phosphatidylinositol 3-kinase (PI3K) signaling transduction pathway by inhibition of the PI3K substrate recognition sites has been proved to be an effective approach to block cancer progression. Many PI3K inhibitors have been developed. Seven drugs have been approved by the US FDA with a mechanism of targeting the phosphatidylinositol 3-kinase/protein kinase-B/mammalian target of rapamycin (PI3K/AKT/mTOR) signaling pathway. In this study, we used docking tools to investigate selective binding of ligands toward four different subtypes of PI3Ks (PI3Kα, PI3Kβ, PI3Kγ and PI3Kδ). The affinity predicted from both the Glide dock and the Movable-Type (MT)-based free energy calculations agreed well with the experimental data. The validation of our predicted methods with a large dataset of 147 ligands showed very small mean errors. We identified residues that may dictate the subtype-specific binding. Particularly, residues Asp964, Ser806, Lys890 and Thr886 of PI3Kγ might be utilized for PI3Kγ-selective inhibitor design. Residues Val828, Trp760, Glu826 and Tyr813 may be important for PI3Kδ-selective inhibitor binding.

Molecular dynamics study on the behavior and binding mechanism of target protein Transgelin-2 with its agonist TSG12 for anti-asthma drug discovery

Transgelin-2 (TG2) is a novel promising therapeutic target for the treatment of asthma as it plays an important role in relaxing airway smooth muscles and reducing pulmonary resistance in asthma. The compound TSG12 is the only reported TG2 agonist with in vivo anti-asthma activity. However, the dynamic behavior and ligand binding sites of TG2 and its binding mechanism with TSG12 remain unclear. In this study, we performed 12.6 μs molecular dynamics (MD) simulations for apo-TG2 and TG2-TSG12 complex, respectively. The results suggested that the apo-TG2 has 4 most populated conformations, and that its binding of the agonist could expand the conformation distribution space of the protein. The simulations revealed 3 potential binding sites in 3 most populated conformations, one of which is induced by the agonist binding. Free energy decomposition uncovered 8 important residues with contributions stronger than -1 kcal/mol. Computational alanine scanning for the important residues by 100 ns conventional MD simulation for each mutated TG2-TSG12 complexes demonstrated that E27, R49 and F52 are essential residues for the agonist binding. These results should be helpful to understand the dynamic behavior of TG2 and its binding mechanism with the agonist TSG12, which could provide some structural insights into the novel mechanism for anti-asthma drug development.

Receptor-Based Virtual Screening of Large Libraries in a Multi-Level In Silico Approach

Structure-based drug design (SBDD) has become an alternative to high throughput screening (HTS) as it reduces experimental costs and time. It works like a funnel, filtering out compounds that do not show good affinity (or score) toward a particular target, with known 3D structure.Here, we describe a protocol for structure-based drug design using a multi-level in silico approach, combining Molecular Docking, Virtual Screening, Molecular Dynamics Simulations and Free energy calculations to find new lead molecules for experimental testing, predict binding affinities and characterize binding modes.

RestraintMaker: a graph-based approach to select distance restraints in free-energy calculations with dual topology

The calculation of relative binding free energies (RBFE) involves the choice of the end-state/system representation, of a sampling approach, and of a free-energy estimator. System representations are usually termed “single topology” or “dual topology”. As the terminology is often used ambiguously in the literature, a systematic categorization of the system representations is proposed here. In the dual-topology approach, the molecules are simulated as separate molecules. Such an approach is relatively easy to automate for high-throughput RBFE calculations compared to the single-topology approach. Distance restraints are commonly applied to prevent the molecules from drifting apart, thereby improving the sampling efficiency. In this study, we introduce the program RestraintMaker, which relies on a greedy algorithm to find (locally) optimal distance restraints between pairs of atoms based on geometric measures. The algorithm is further extended for multi-state methods such as enveloping distribution sampling (EDS) or multi-site \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda $$\end{document}λ-dynamics. The performance of RestraintMaker is demonstrated for toy models and for the calculation of relative hydration free energies. The Python program can be used in script form or through an interactive GUI within PyMol. The selected distance restraints can be written out in GROMOS or GROMACS file formats. Additionally, the program provides a human-readable JSON format that can easily be parsed and processed further. The code of RestraintMaker is freely available on GitHub https://github.com/rinikerlab/restraintmaker.

Relative free-energy calculations for scaffold hopping-type transformations with an automated RE-EDS sampling procedure

The calculation of relative free-energy differences between different compounds plays an important role in drug design to identify potent binders for a given protein target. Most rigorous methods based on molecular dynamics simulations estimate the free-energy difference between pairs of ligands. Thus, the comparison of multiple ligands requires the construction of a “state graph”, in which the compounds are connected by alchemical transformations. The computational cost can be optimized by reducing the state graph to a minimal set of transformations. However, this may require individual adaptation of the sampling strategy if a transformation process does not converge in a given simulation time. In contrast, path-free methods like replica-exchange enveloping distribution sampling (RE-EDS) allow the sampling of multiple states within a single simulation without the pre-definition of alchemical transition paths. To optimize sampling and convergence, a set of RE-EDS parameters needs to be estimated in a pre-processing step. Here, we present an automated procedure for this step that determines all required parameters, improving the robustness and ease of use of the methodology. To illustrate the performance, the relative binding free energies are calculated for a series of checkpoint kinase 1 inhibitors containing challenging transformations in ring size, opening/closing, and extension, which reflect changes observed in scaffold hopping. The simulation of such transformations with RE-EDS can be conducted with conventional force fields and, in particular, without soft bond-stretching terms.

Mechanistic and thermodynamic characterization of oxathiazolones as potent and selective covalent immunoproteasome inhibitors

The ubiquitin–proteasome system is responsible for the degradation of proteins and plays a critical role in key cellular processes. While the constitutive proteasome (cPS) is expressed in all eukaryotic cells, the immunoproteasome (iPS) is primarily induced during disease processes, and its inhibition is beneficial in the treatment of cancer, autoimmune disorders and neurodegenerative diseases. Oxathiazolones were reported to selectively inhibit iPS over cPS, and the inhibitory activity of several oxathiazolones against iPS was experimentally determined. However, the detailed mechanism of the chemical reaction leading to irreversible iPS inhibition and the key selectivity drivers are unknown, and separate characterization of the noncovalent and covalent inhibition steps is not available for several compounds. Here, we investigate the chemical reaction between oxathiazolones and the Thr1 residue of iPS by quantum mechanics/molecular mechanics (QM/MM) simulations to establish a plausible reaction mechanism and to determine the rate-determining step of covalent complex formation. The modelled binding mode and reaction mechanism are in line with the selective inhibition of iPS versus cPS by oxathiazolones. The k_inact value of several ligands was estimated by constructing the potential of mean force of the rate-determining step by QM/MM simulations coupled with umbrella sampling. The equilibrium constant K_i of the noncovalent complex formation was evaluated by classical force field-based thermodynamic integration. The calculated K_i and k_inact values made it possible to analyse the contribution of the noncovalent and covalent steps to the overall inhibitory activity. Compounds with similar intrinsic reactivities exhibit varying selectivities for iPS versus cPS owing to subtle differences in the binding modes that slightly affect K_i, the noncovalent affinity, and importantly alter k_inact, the covalent reactivity of the bound compounds. A detailed understanding of the inhibitory mechanism of oxathiazolones is useful in designing iPS selective inhibitors with improved drug-like properties.

Sterenin M as a potential inhibitor of SARS-CoV-2 main protease identified from MeFSAT database using molecular docking, molecular dynamics simulation and binding free energy calculation

The disease outbreak of Coronavirus disease-19 (COVID-19), caused by the novel SARS-CoV-2 virus, remains a public health concern. COVID-19 is spreading rapidly with a high mortality rate due to unavailability of effective treatment or vaccine for the disease. The high rate of mutation and recombination in SARS-CoV2 makes it difficult for scientist to develop specific anti-CoV2 drugs and vaccines. SARS-CoV-2-Mpro cleaves the viral polyprotein to produce a variety of non-structural proteins, but in human host it also cleaves the nuclear transcription factor kappa B (NF-κB) essential modulator (NEMO), which suppresses the activation of the NF-κB pathway and weakens the immune response. Since the main protease (Mpro) is required for viral gene expression and replication, it is a promising target for antagonists to treat novel coronavirus disease and discovery of high resolution crystal structure of SARS-CoV-2-Mpro provide an opportunity for in silico identification of its possible inhibitors. In this study we intend to find novel and potential Mpro inhibitors from around 1830 chemically diverse and therapeutically important secondary metabolites available in the MeFSAT database by performing molecular docking against the Mpro structure of SARS-CoV-2 (PDB ID: 6LZE). After ADMET (absorption, distribution, metabolism, excretion, and toxicity) profile and binding energy calculation through MM-GBSA for top five hits, Sterenin M was proposed as a SARS-CoV2-Mpro inhibitor with validation of molecular dynamics (MD) simulation study. Sterenin M seems to have the potential to be a promising ligand against SARS-CoV-2, and thus it requires further validation by in vitro and in vivo studies.

Ensemble-based screening of natural products and FDA-approved drugs identified potent inhibitors of SARS-CoV-2 that work with two distinct mechanisms

The recent outbreak of SARS-CoV-2 is responsible for high morbidity and mortality rate across the globe. This requires an urgent identification of drugs and other interventions to overcome this pandemic. Computational drug repurposing represents an alternative approach to provide a more effective approach in search for COVID-19 drugs. Selected natural product known to have antiviral activities were screened, and based on their hits; a similarity search with FDA approved drugs was performed using computational methods. Obtained drugs from similarity search were assessed for their stability and inhibition against SARS-CoV-2 targets. Diosmin (DB08995) was found to be a promising drug that works with two distinct mechanisms, preventing viral replication and viral fusion into the host cell. Isoquercetin (DB12665) and rutin (DB01698) work by inhibiting viral replication and preventing cell entry, respectively. Our analysis based on molecular dynamics simulation and MM-PBSA binding free energy calculation suggests that diosmin, isoquercetin, rutin and other similar flavone glycosides could serve as SARS-CoV-2 inhibitor, hence an alternative solution to treat COVID-19 upon further clinical validation.

New insights into Hexakis macrocycles as a novel nano-carrier for highly potent anti-cancer treatment: A new challenge in drug delivery

The purpose of this study is to design and evaluate a new tubular assembly structure of Hexakis (m-phenylene ethynylene) (m-PE) macrocycles and to explore its potential application as an innovative drug delivery system. First, we focused on how (m-PE) macrocycles can be self-assembled in both chloroform (CHCl₃) and water solvents for the formation of the assembled nanotube using molecular dynamics (MD) simulations. In contrast to their behavior in water solvent, all ten (m-PE) macrocycles remain aggregated at low concentrations of CHCl₃. We found that these macrocycles carrying chiral side chains and capable of H-bonded self-association, assemble into tubular stacks. Then, the dual delivery strategy for the transport of doxorubicin (DOX) and curcumin (Cur) on the self-assembly system of hexakis (m-PE) nanocarrier is examined using molecular dynamics (MD) simulation and free energy calculation. The obtained results indicated that the binding energy of DOX (- 298.9 kJ/mol) on hexakis (m-PE) in the presence of Cur is higher than free DOX (- 247.7 kJ/mol). Furthermore, in the interaction of the DOX and hexakis m-PE, the contribution of van der Walls (vdW) energy is higher than electrostatic (elec) energy, which can be related to the strong π-π interactions between the drug molecules with the carrier surface. In general, the results indicated that the simultaneous delivery of DOX and Cur through DOX/Cur/hexakis (m-PE) could be a promising vehicle in tumor therapy. Based on the obtained results of the present research, hexakis (m-PE) macrocycle can be used as a drug delivery vehicle for targeted or systemic delivery.

SAMPL6 host-guest binding affinities and binding poses from spherical-coordinates-biased simulations

Host-guest binding is a challenging problem in computer simulation. The prediction of binding affinities between hosts and guests is an important part of the statistical assessment of the modeling of proteins and ligands (SAMPL) challenges. In this work, the volume-based variant of well-tempered metadynamics is employed to calculate the binding affinities of the host-guest systems in the SAMPL6 challenge. By biasing the spherical coordinates describing the relative position of the host and the guest, the initial-configuration-induced bias vanishes and all possible binding poses are explored. The agreement between the predictions and the experimental results and the observation of new binding poses indicate that the volume-based technique serves as a nice candidate for the calculation of binding free energies and the search of the binding poses.

Binding and intracellular transport of 25-hydroxycholesterol by Niemann-Pick C2 protein

Side-chain oxidized cholesterol derivatives, like 25-hydroxycholesterol (25-OH-Chol) are important regulators of cellular cholesterol homeostasis. How transport of oxysterols through the endo-lysosomal pathway contributes to their biological function is not clear. The Niemann-Pick C2 protein (NPC2) is a small lysosomal sterol transfer protein required for export of cholesterol from late endosomes and lysosomes (LE/LYSs). Here, we show that 25-hydroxy-cholestatrienol, (25-OH-CTL), an intrinsically fluorescent analogue of 25-OH-Chol, becomes trapped in LE/LYSs of NPC2-deficient fibroblasts, but can efflux from the cells even in the absence of NPC2 upon removal of the sterol source. Fluorescence recovery after photobleaching (FRAP) of 25-OH-CTL in endo-lysosomes was rapid and extensive and only partially dependent on NPC2 function. Using quenching of NPC2's intrinsic fluorescence, we show that 25-OH-Chol and 25-OH-CTL can bind to NPC2 though with lower affinity compared to cholesterol and its fluorescent analogues, cholestatrienol (CTL) and dehydroergosterol (DHE). This is confirmed by calculations of binding energies which additionally show that 25-OH-CTL can bind in two orientations to NPC2, in stark contrast to cholesterol and its analogues. We conclude that NPC2's affinity for all sterols is energetically favored over their self-aggregation in the lysosomal lumen. Lysosomal export of 25-OH-Chol is not strictly dependent on the NPC2 protein.

Development of sodium glucose co-transporter 2 (SGLT2) inhibitors with novel structure by molecular docking and dynamics simulation

In this study, molecular docking studies were carried out to explore the binding interactions of sodium glucose co-transporter 2 (SGLT2) with its inhibitors. A correlation between the docking scores and the experimental bioactivity was observed (R² = 0.8368, N = 24). The new inhibitors were designed using the 3D quantitative structure activity relationship (3D-QSAR) method, and the activities were predicted by the docking method. In order to understand the structure-activity correlation of compound 1 m (the highest score of docking) and compound 1 t (the lowest score), we carried out a combined molecular dynamics simulation and MM-GBSA method. It was found that, in the system of SGLT2-1 m, the interaction between Gln271 and Val272 exhibited significant effects, which were absent in the SGLT2-1 t system. This study is expected to shed light on the mechanism of action of compound 1 m, leading to development of active drug candidates targeting SGLT2.

Computational insights into the binding of IN17 inhibitors to MELK

The protein kinase MELK is an important kinase in cell signaling and has shown to be a promising anti-cancer target. Recent work has resulted in a novel small molecule scaffold targeting MELK, IN17. However, there has been little structural information or physical understanding of MELK-IN17 interactions. Using Tinker-OpenMM on GPUs, we have performed free energy simulations on MELK binding with IN17 and 11 derivatives. This series of studies provides structural insights into how substitution on IN17 leads to differences in complex structure and binding thermodynamics. In addition, this study serves as an assessment of the current capabilities of the AMOEBA forcefield, accelerated by GPU computing, to serve as a molecular-dynamics-based free energy simulation platform for lead optimization.

Structural and mechanistic insights into modified G-quadruplex thrombin-binding DNA aptamers

Thrombin-binding aptamer (TBA) can fold into a G-quadruplex structure necessary for interacting with thrombin. When one thymidine residue of the TGT loop at position 7 is replaced with unlocked uracil (UNA), d-isothymidine (D-isoT) or l-isothymidine (L-isoT), these modified sequences display different activities. To date, the mechanisms of how D/L-isoT and UNA influence the biological properties of TBA have not been illustrated in the literature. In this paper, we fill this gap by probing the structure variations and binding modes of these modified TBAs via molecular dynamics (MD) simulation and free energy calculation. Comparative structural analyses demonstrated that both D-IsoT and UNA changed the local conformation of TGT loop and formed stronger interactions with the target protein. Particularly, D-IsoT and UNA adopted similar conformation which can well explain their similar biological activities. In addition, the flexibility of the two TT loops were described clearly. In contrast, L-IsoT at position 7 led to an obvious tendency to unfold. Free energy calculation and the analysis of key residues energy contributions eventually provide a clear picture of interactions for further understanding of the structure-activity relationships. Collectively, our findings open the way for a rational design of modified aptamers.

Modeling and Simulation of Oligonucleotide Hybrids: Outlining a Strategy

Molecular dynamics simulations with a state-of-the-art force field provide an atomistic detailed description of the structural and thermodynamic features of biomolecules. Effects of chemical modifications and of the environment such as sequence, solvent, and ionic strength can explicitly be taken into account. Molecular simulation techniques can also provide insight in change in binding affinity, in protonation (pK_a shift) and tautomeric propensity due to changes in the environment or in the molecular system. The quality and reliability of a simulation depend on the quality of the force field and on the reproducibility of the data, and validation depends on the availability of suitable experimental data. Here, we describe the workflow to investigate oligonucleotide hybrids using molecular simulation including hardware and software information.

Selective mechanisms and molecular design of 2,4 Diarylaminopyrimidines as ALK inhibitors

As an attractive therapeutic target for non-small-cell lung cancer (NSCLC), anaplastic lymphoma kinase (ALK) has got increased attention, and the selectivity of ALK inhibitors is an enormous challenge. Recently, 2,4-Diarylaminopyrimidines with high inhibitory activity over InsR/IGF1R were reported as ALK inhibitors, which harboring phosphine oxide moiety. In this work, it is the first time to reveal that the incorporation of dimethylphosphine oxide moiety and the smaller active pocket of ALK is key factor in the selectivity of inhibitor 11q toward ALK over IGF1R/InsR. The results of molecular simulation indicate that the subtle change in the binding pocket of ALK is mainly associated with the flexibility of P-loop and the own residues K1150 and D1270. The replacement of the dimethylphosphine oxide and methylpiperazine of inhibitor 11q would alter the major inhibitory effects of binding and activation. The results further combined 3D-QSAR can not only profile the binding mechanism between the 2,4-Diarylaminopyrimidines inhibitors and ALK, but also supply the useful information for the rational design of a more potential small molecule inhibitor bound to ALK receptor.

A molecular dynamics investigation of CDK8/CycC and ligand binding: conformational flexibility and implication in drug discovery

Abnormal activity of cyclin-dependent kinase 8 (CDK8) along with its partner protein cyclin C (CycC) is a common feature of many diseases including colorectal cancer. Using molecular dynamics (MD) simulations, this study determined the dynamics of the CDK8-CycC system and we obtained detailed breakdowns of binding energy contributions for four type-I and five type-II CDK8 inhibitors. We revealed system motions and conformational changes that will affect ligand binding, confirmed the essentialness of CycC for inclusion in future computational studies, and provide guidance in development of CDK8 binders. We employed unbiased all-atom MD simulations for 500 ns on twelve CDK8-CycC systems, including apoproteins and protein-ligand complexes, then performed principal component analysis (PCA) and measured the RMSF of key regions to identify protein dynamics. Binding pocket volume analysis identified conformational changes that accompany ligand binding. Next, H-bond analysis, residue-wise interaction calculations, and MM/PBSA were performed to characterize protein-ligand interactions and find the binding energy. We discovered that CycC is vital for maintaining a proper conformation of CDK8 to facilitate ligand binding and that the system exhibits motion that should be carefully considered in future computational work. Surprisingly, we found that motion of the activation loop did not affect ligand binding. Type-I and type-II ligand binding is driven by van der Waals interactions, but electrostatic energy and entropic penalties affect type-II binding as well. Binding of both ligand types affects protein flexibility. Based on this we provide suggestions for development of tighter-binding CDK8 inhibitors and offer insight that can aid future computational studies.

Identification of curcumin derivatives as human LMTK3 inhibitors for breast cancer: a docking, dynamics, and MM/PBSA approach

Human lemur tyrosine kinase-3 (LMTK3) is primarily involved in regulation of estrogen receptor-α (ERα) by phosphorylation activity. LMTK3 acts as key biomarker for ERα positive breast cancer and identified as novel drug target for breast cancer. Due to the absence of experimental reports, the computational approach has been followed to screen LMTK3 inhibitors from natural product curcumin derivatives based on rational inhibitor design. The initial virtual screening and re-docking resulted in identification of top three leads with favorable binding energy and strong interactions in critical residues of ATP-binding cavity. ADME prediction confirmed the pharmacological activity of the leads with various properties. The stability and binding affinity of leads were well refined in dynamic system from 25 ns MD simulations. The behavior of protein motion towards closure of ATP-binding cavity was evaluated based on eigenvectors by PCA. In addition, MM/PBSA calculations also confirmed the relative binding free energy of LMTK3-lead complexes in favor of the effective binding. From our study, novel LMTK3 inhibitors tetrahydrocurcumin, curcumin 4,4'-diacetate, and demethoxycurcumin have been proposed with inhibition mechanism. Further experimental evaluation on reported lead candidates might prove its role in breast cancer therapeutics.

Free energy calculations on the stability of the 14-3-3ζ protein

Mutations of cysteine are often introduced to e.g. avoid formation of non-physiological inter-molecular disulfide bridges in in-vitro experiments, or to maintain specificity in labeling experiments. Alanine or serine is typically preferred, which usually do not alter the overall protein stability, when the original cysteine was surface exposed. However, selecting the optimal mutation for cysteines in the hydrophobic core of the protein is more challenging. In this work, the stability of selected Cys mutants of 14-3-3ζ was predicted by free-energy calculations and the obtained data were compared with experimentally determined stabilities. Both the computational predictions as well as the experimental validation point at a significant destabilization of mutants C94A and C94S. This destabilization could be attributed to the formation of hydrophobic cavities and a polar solvation of a hydrophilic side chain. A L12E, M78K double mutant was further studied in terms of its reduced dimerization propensity. In contrast to naïve expectations, this double mutant did not lead to the formation of strong salt bridges, which was rationalized in terms of a preferred solvation of the ionic species. Again, experiments agreed with the calculations by confirming the monomerization of the double mutants. Overall, the simulation data is in good agreement with experiments and offers additional insight into the stability and dimerization of this important family of regulatory proteins.

Multiscale molecular dynamics simulations of rotary motor proteins

Protein functions require specific structures frequently coupled with conformational changes. The scale of the structural dynamics of proteins spans from the atomic to the molecular level. Theoretically, all-atom molecular dynamics (MD) simulation is a powerful tool to investigate protein dynamics because the MD simulation is capable of capturing conformational changes obeying the intrinsically structural features. However, to study long-timescale dynamics, efficient sampling techniques and coarse-grained (CG) approaches coupled with all-atom MD simulations, termed multiscale MD simulations, are required to overcome the timescale limitation in all-atom MD simulations. Here, we review two examples of rotary motor proteins examined using free energy landscape (FEL) analysis and CG-MD simulations. In the FEL analysis, FEL is calculated as a function of reaction coordinates, and the long-timescale dynamics corresponding to conformational changes is described as transitions on the FEL surface. Another approach is the utilization of the CG model, in which the CG parameters are tuned using the fluctuation matching methodology with all-atom MD simulations. The long-timespan dynamics is then elucidated straightforwardly by using CG-MD simulations.

eTOX ALLIES: an automated pipeLine for linear interaction energy-based simulations

BACKGROUND

Computational methods to predict binding affinities of small ligands toward relevant biological (off-)targets are helpful in prioritizing the screening and synthesis of new drug candidates, thereby speeding up the drug discovery process. However, use of ligand-based approaches can lead to erroneous predictions when structural and dynamic features of the target substantially affect ligand binding. Free energy methods for affinity computation can include steric and electrostatic protein-ligand interactions, solvent effects, and thermal fluctuations, but often they are computationally demanding and require a high level of supervision. As a result their application is typically limited to the screening of small sets of compounds by experts in molecular modeling.

RESULTS

We have developed eTOX ALLIES, an open source framework that allows the automated prediction of ligand-binding free energies requiring the ligand structure as only input. eTOX ALLIES is based on the linear interaction energy approach, an efficient end-point free energy method derived from Free Energy Perturbation theory. Upon submission of a ligand or dataset of compounds, the tool performs the multiple steps required for binding free-energy prediction (docking, ligand topology creation, molecular dynamics simulations, data analysis), making use of external open source software where necessary. Moreover, functionalities are also available to enable and assist the creation and calibration of new models. In addition, a web graphical user interface has been developed to allow use of free-energy based models to users that are not an expert in molecular modeling.

CONCLUSIONS

Because of the user-friendliness, efficiency and free-software licensing, eTOX ALLIES represents a novel extension of the toolbox for computational chemists, pharmaceutical scientists and toxicologists, who are interested in fast affinity predictions of small molecules toward biological (off-)targets for which protein flexibility, solvent and binding site interactions directly affect the strength of ligand-protein binding.

CDOCKER and $$\lambda$$-dynamics for prospective prediction in D3R Grand Challenge 2

The opportunity to prospectively predict ligand bound poses and free energies of binding to the Farnesoid X Receptor in the D3R Grand Challenge 2 provided a useful exercise to evaluate CHARMM based docking (CDOCKER) and [Formula: see text]-dynamics methodologies for use in "real-world" applications in computer aided drug design. In addition to measuring their current performance, several recent methodological developments have been analyzed retrospectively to highlight best procedural practices in future applications. For pose prediction with CDOCKER, when the protein structure used for rigid receptor docking was close to the crystallographic holo structure, reliable poses were obtained. Benzimidazoles, with a known holo receptor structure, were successfully docked with an average RMSD of 0.97 [Formula: see text]. Other non-benzimidazole ligands displayed less accuracy largely because the receptor structures we chose for docking were too different from the experimental holo structures. However, retrospective analysis has shown that when these ligands were re-docked into their holo structures, the average RMSD dropped to 1.18 [Formula: see text] for all ligands. When sulfonamides and spiros were docked with the apo structure, which agrees more with their holo structure than the structures we chose, five out of six ligands were correctly docked. These docking results emphasize the need for flexible receptor docking approaches. For [Formula: see text]-dynamics techniques, including multisite [Formula: see text]-dynamics (MS[Formula: see text]D), reasonable agreement with experiment was observed for the 33 ligands investigated; root mean square errors of 2.08 and 1.67 kcal/mol were obtained for free energy sets 1 and 2, respectively. Retrospectively, soft-core potentials, adaptive landscape flattening, and biasing potential replica exchange (BP-REX) algorithms were critical to model large substituent perturbations with sufficient precision and within restrictive timeframes, such as was required with participation in Grand Challenge 2. These developments, their associated benefits, and proposed procedures for their use in future applications are discussed.

Receptor-guided 3D-QSAR studies, molecular dynamics simulation and free energy calculations of Btk kinase inhibitors

BACKGROUND

Bruton tyrosine kinase (Btk) plays an important role in B-cell development, differentiation, and signaling. It is also found be in involved in male immunodeficiency disease such as X-linked agammaglobulinemia (XLA). Btk is considered as a potential therapeutic target for treating autoimmune diseases and hematological malignancies.

RESULTS

In this work, a combined molecular modeling study was performed on a series of thieno [3,2-c] pyridine-4-amine derivatives as Btk inhibitors. Receptor-guided COMFA (q 2 = 0.574, NOC = 3, r 2 = 0.924) and COMSIA (q 2 = 0.646, NOC = 6, r 2 = 0.971) models were generated based on the docked conformation of the most active compound 26. All the developed models were tested for robustness using various validation techniques. Furthermore, a 5-ns molecular dynamics (MD) simulation and binding free energy calculations were carried out to determine the binding modes of the inhibitors and to identify crucial interacting residues. The rationality and stability of molecular docking and 3D-QSAR results were validated by MD simulation. The binding free energies calculated by the MM/PBSA method showed the importance of the van der Waals interaction.

CONCLUSIONS

A good correlation between the MD results, docking studies, and the contour map analysis were observed. The study has identified the key amino acid residues in Btk binding pocket. The results from this study can provide some insights into the development of potent, novel Btk inhibitors.

The importance of protonation and tautomerization in relative binding affinity prediction: a comparison of AMBER TI and Schrödinger FEP

In drug discovery, protonation states and tautomerization are easily overlooked. Through a Merck-Rutgers collaboration, this paper re-examined the initial settings and preparations for the Thermodynamic Integration (TI) calculation in AMBER Free-Energy Workflows, demonstrating the value of careful consideration of ligand protonation and tautomer state. Finally, promising results comparing AMBER TI and Schrödinger FEP+ are shown that should encourage others to explore the value of TI in routine Structure-based Drug Design.

Free energy calculation provides insight into the action mechanism of selective PARP-1 inhibitor

Selective poly (ADP-ribose) polymerase (PARP)-1 inhibitor represents promising therapy against cancers with a good balance between efficacy and safety. Owing to the conserved structure between PARP-1 and PARP-2, most of the clinical and experimental drugs show equivalent inhibition against both targets. Most recently, it's disclosed a highly selective PARP-1 inhibitor (NMS-P118) with promising pharmacokinetic properties. Herein, we combined molecular simulation with free energy calculation to gain insights into the selective mechanism of NMS-P118. Our results suggest the reduction of binding affinity for PARP-2 is attributed to the unfavorable conformational change of protein, which is accompanied by a significant energy penalty. Alanine-scanning mutagenesis study further reveals the important role for a tyrosine residue of donor loop (Tyr889(PARP-1) and Tyr455(PARP-2)) in contributing to the ligand selectivity. Retrospective structural analysis indicates the ligand-induced movement of Tyr455(PARP-2) disrupts the intra-molecule hydrogen bonding network, which partially accounts for the "high-energy" protein conformation in the presence of NMS-P118. Interestingly, such effect isn't observed in other non-selective PARP inhibitors including BMN673 and A861695, which validates the computational prediction. Our work provides energetic insight into the subtle variations in the crystal structures and could facilitate rational design of new selective PARP inhibitor.

A Python tool to set up relative free energy calculations in GROMACS

Free energy calculations based on molecular dynamics (MD) simulations have seen a tremendous growth in the last decade. However, it is still difficult and tedious to set them up in an automated manner, as the majority of the present-day MD simulation packages lack that functionality. Relative free energy calculations are a particular challenge for several reasons, including the problem of finding a common substructure and mapping the transformation to be applied. Here we present a tool, alchemical-setup.py, that automatically generates all the input files needed to perform relative solvation and binding free energy calculations with the MD package GROMACS. When combined with Lead Optimization Mapper (LOMAP; Liu et al. in J Comput Aided Mol Des 27(9):755-770, 2013), recently developed in our group, alchemical-setup.py allows fully automated setup of relative free energy calculations in GROMACS. Taking a graph of the planned calculations and a mapping, both computed by LOMAP, our tool generates the topology and coordinate files needed to perform relative free energy calculations for a given set of molecules, and provides a set of simulation input parameters. The tool was validated by performing relative hydration free energy calculations for a handful of molecules from the SAMPL4 challenge (Mobley et al. in J Comput Aided Mol Des 28(4):135-150, 2014). Good agreement with previously published results and the straightforward way in which free energy calculations can be conducted make alchemical-setup.py a promising tool for automated setup of relative solvation and binding free energy calculations.

Strontium selectivity in sodium nonatitanate Na ₄Ti₉O₂₀·xH₂O

We study the extraction of strontium by sodium nonatitanate powder from nitrate strontium and acetate sodium mixture. Experiments show that adsorption is quantitative. The excess Gibbs free energy has been modeled by various models (ideal, 2D Coulomb, regular solution model) for the solid phase. We find that the free energy of the solid phase is controlled by short-range interactions rather than long-ranged Coulombic forces. The selectivity is the consequence of a competition between the liquid and solid phases: both phases prefer strontium rather than sodium but the solid contribution is predominant.

Correcting for the free energy costs of bond or angle constraints in molecular dynamics simulations

BACKGROUND

Free energy simulations are an important tool in the arsenal of computational biophysics, allowing the calculation of thermodynamic properties of binding or enzymatic reactions. This paper introduces methods to increase the accuracy and precision of free energy calculations by calculating the free energy costs of constraints during post-processing. The primary purpose of employing constraints for these free energy methods is to increase the phase space overlap between ensembles, which is required for accuracy and convergence.

METHODS

The free energy costs of applying or removing constraints are calculated as additional explicit steps in the free energy cycle. The new techniques focus on hard degrees of freedom and use both gradients and Hessian estimation. Enthalpy, vibrational entropy, and Jacobian free energy terms are considered.

RESULTS

We demonstrate the utility of this method with simple classical systems involving harmonic and anharmonic oscillators, four-atomic benchmark systems, an alchemical mutation of ethane to methanol, and free energy simulations between alanine and serine. The errors for the analytical test cases are all below 0.0007kcal/mol, and the accuracy of the free energy results of ethane to methanol is improved from 0.15 to 0.04kcal/mol. For the alanine to serine case, the phase space overlaps of the unconstrained simulations range between 0.15 and 0.9%. The introduction of constraints increases the overlap up to 2.05%. On average, the overlap increases by 94% relative to the unconstrained value and precision is doubled.

CONCLUSIONS

The approach reduces errors arising from constraints by about an order of magnitude. Free energy simulations benefit from the use of constraints through enhanced convergence and higher precision.

GENERAL SIGNIFICANCE

The primary utility of this approach is to calculate free energies for systems with disparate energy surfaces and bonded terms, especially in multi-scale molecular mechanics/quantum mechanics simulations. This article is part of a Special Issue entitled Recent developments of molecular dynamics.

Molecular modeling study on the dynamical structural features of human smoothened receptor and binding mechanism of antagonist LY2940680 by metadynamics simulation and free energy calculation

BACKGROUND

The smoothened (SMO) receptor, one of the Class F G protein coupled receptors (GPCRs), is an essential component of the canonical hedgehog signaling pathway which plays a key role in the regulation of embryonic development in animals. The function of the SMO receptor can be modulated by small-molecule agonists and antagonists, some of which are potential antitumour agents. Understanding the binding mode of an antagonist in the SMO receptor is crucial for the rational design of new antitumour agents.

METHODS

Molecular dynamics (MD) simulation and dynamical network analysis are used to study the dynamical structural features of SMO receptor. Metadynamics simulation and free energy calculation are employed to explore the binding mechanism between the antagonist and SMO receptor.

RESULTS

The MD simulation results and dynamical network analysis show that the conserved KTXXXW motif in helix VIII has strong interaction with helix I. The α-helical extension of transmembrane 6 (TM6) is detected as part of the ligand-binding pocket and dissociation pathway of the antagonist. The metadynamics simulation results illustrate the binding mechanism of the antagonist in the pocket of SMO receptor, and free energy calculation shows the antagonist needs to overcome about 38kcal/mol of energy barrier to leave the binding pocket of SMO receptor.

CONCLUSIONS

The unusually long TM6 plays an important role on the binding behavior of the antagonist in the pocket of SMO receptor.

GENERAL SIGNIFICANCE

The results can not only profile the binding mechanism between the antagonist and Class F GPCRs, but also supply the useful information for the rational design of a more potential small molecule antagonist bound to SMO receptor.

Dual base-flipping of cytosines in a CpG dinucleotide sequence

Simultaneous flipped-out conformation of two neighboring bases on opposite strands of DNAs has been observed in several X-ray structures. It has also been detected for two cytosines on opposite strands in different contexts of CpG sites. In this paper, we study by MD simulations the dual base flipping of the two cytosines in hemi-methylated CpG site. We calculate the potential of mean force of flipping-out the unmethylated cytosine in three model systems. The first is for DNA bound to the regulatory protein UHRF1. In this case, the methyl-cytosine on the complementary strand is flipped-out into the binding pocket of the SRA domain of the protein. The other two systems are for unbound DNAs in which the methyl-cytosine is either intra-helical or extra-helical. We find that when the methyl-cytosine is flipped-out it is easier to flip-out the other (unmethylated) cytosine on the opposite strand by about 14-16kJ/mol. This lower penalty for dual-base flipping is observed for both the bound and unbound states of the DNA. Analyses of the hydrogen bond network and stacking interactions within the CpG site indicate that the lower penalty is due to stabilization of the dual-base flipped-out conformation via interactions involving the orphan guanines. The results presented in this paper suggest that the extra-helical conformation of the methyl-cytosine recognized by UHRF1 can facilitate the base-flipping process of the target cytosine to be methylated by Dnmt1.