Quantum Mechanical-Cluster Approach to Solve the Bioisosteric Replacement Problem in Drug Design.

Bioisosteres are molecules that differ in substituents but still have very similar shapes. Bioisosteric replacements are ubiquitous in modern drug design, where they are used to alter metabolism, change bioavailability, or modify activity of the lead compound. Prediction of relative affinities of bioisosteres with computational methods is a long-standing task; however, the very shape closeness makes bioisosteric substitutions almost intractable for computational methods, which use standard force fields. Here, we design a quantum mechanical (QM)-cluster approach based on the GFN2-xTB semi-empirical quantum-chemical method and apply it to a set of H → F bioisosteric replacements. The proposed methodology enables advanced prediction of biological activity change upon bioisosteric substitution of -H with -F, with the standard deviation of 0.60 kcal/mol, surpassing the ChemPLP scoring function (0.83 kcal/mol), and making QM-based ΔΔG estimation comparable to ∼0.42 kcal/mol standard deviation of in vitro experiment. The speed of the method and lack of tunable parameters makes it affordable in current drug research.

Novel curcumin derivatives as P-glycoprotein inhibitors: Molecular modeling, synthesis and sensitization of multidrug resistant cells to doxorubicin.

The MDR1/P-glycoprotein (Pgp)/ABCB1 multidrug transporter is being investigated as a druggable target for antitumor therapy for decades. The natural product curcumin is known to provide an efficient scaffold for compounds capable of blocking Pgp mediated efflux and sensitization of multidrug resistant (MDR) cells to the Pgp transported drug doxorubicin (Dox). We performed molecular dynamics simulations and docking of curcumin derivatives into the Pgp model. Based on these calculations, a series of pyrazolocurcumin derivatives with predicted metabolic stability and/or improved binding affinity were proposed for synthesis and evaluation of MDR reversal potency against Dox selected K562/4 subline, a derivative of K562 human chronic myelogenous leukemia cell line. Compounds 16 and 19 which are both dimethylcurcumin pyrazole derivatives bearing an N-p-phenylcarboxylic amide substitution, were the most potent Pgp blockers as determined by intracellular Dox accumulation. Furthermore, at non-toxic submicromolar concentrations 16 and 19 dramatically sensitized K562/4 cells to Dox. Together with good water solubility of 16 and 19, these results indicate that the new pyrazolo derivatives of curcumin are a promising scaffold for development of clinically applicable Pgp antagonists.

The role of human in the loop: lessons from D3R challenge 4.

The rapid development of new machine learning techniques led to significant progress in the area of computer-aided drug design. However, despite the enormous predictive power of new methods, they lack explainability and are often used as black boxes. The most important decisions in drug discovery are still made by human experts who rely on intuitions and simplified representation of the field. We used D3R Grand Challenge 4 to model contributions of human experts during the prediction of the structure of protein-ligand complexes, and prediction of binding affinities for series of ligands in the context of absence or abundance of experimental data. We demonstrated that human decisions have a series of biases: a tendency to focus on easily identifiable protein-ligand interactions such as hydrogen bonds, and neglect for a more distributed and complex electrostatic interactions and solvation effects. While these biases still allow human experts to compete with blind algorithms in some areas, the underutilization of the information leads to significantly worse performance in data-rich tasks such as binding affinity prediction.

PF­114, a novel selective inhibitor of BCR­ABL tyrosine kinase, is a potent inducer of apoptosis in chronic myelogenous leukemia cells.

A t(9;22) chromosomal translocation which forms the chimeric tyrosine kinase breakpoint cluster region (BCR)­Abelson murine leukemia viral oncogene homolog 1 (ABL) is a key mechanism underlying the pathogenesis of chronic myelogenous leukemia (CML). Pharmacological inhibition of BCR­ABL with imatinib (Gleevec) has been reported as an effective targeted therapy; however, mutations (including the kinase domain of ABL) suppress the efficacy of inhibitors. PF­114, a derivative of the third generation BCR­ABL inhibitor ponatinib, demonstrated a high inhibitory activity against wild-type and mutant BCR­ABL forms, such as the clinically important T315I. Furthermore, PF­114 exhibited preferential kinase selectivity, safety, notable pharmacokinetic properties and therapeutic efficacy in a murine model. Investigation into the mechanisms of CML cell death revealed an exceptional potency of PF­114 (at low nanomolar concentrations) for the CML­derived K562 cell line, whereas leukemia cell lines that lack the chimeric tyrosine kinase were markedly more refractory. The molecular ordering of events mechanistically associated with K562 cell death included the dephosphorylation of CrkL adaptor protein followed by inhibition of ERK1/2 and Akt, G1 arrest, a decrease of phosphorylated Bcl­2­associated death promoter, Bcl­2­like protein 11, BH3 interacting­domain death agonist, Bcl­extra large and Bcl­2 family apoptosis regulator, and reduced mitochondrial transmembrane potential. Increased Annexin V reactivity, activation of caspases and poly(ADP­ribose)polymerase cleavage were proposed to lead to internucleosomal DNA fragmentation. Thus, PF­114 may be a potent inducer of apoptosis in CML cells. Nevertheless, activation of STAT3 phosphorylation in response to PF­114 may permit cell rescue; thus, a combination of BCR­ABL and STAT3 inhibitors should be considered for improved therapeutic outcome. Collectively, the targeted killing of BCR­ABL­positive cells, along with other beneficial properties, such as in vivo characteristics, suggests PF­114 as a potential candidate for analysis in clinical trials with CML patients.

Quantifying Possible Routes for SpnF-Catalyzed Formal Diels-Alder Cycloaddition.

The Diels-Alder reaction is a cornerstone of modern organic synthesis. Despite this, it remains essentially inaccessible to biosynthetic approaches. Only a few natural enzymes catalyze even a formal [4 + 2] cycloaddition, and it remains uncertain if any of them proceed via the Diels-Alder mechanism. In this study, we focus on the [4 + 2] cycloaddition step in the biosynthesis of spinosyn A, a reaction catalyzed by SpnF enzyme, one of the most promising "true Diels-Alderase" candidates. The four currently proposed mechanisms (including the Diels-Alder one) for this reaction in water (as a first-order approximation of the enzymatic reaction) are evaluated by an exhaustive quantum mechanical search for possible transition states (728 were found in total). We find that the line between the recently proposed bis-pericyclic [J. Am. Chem. Soc. 2016, 138 (11), 3631] and Diels-Alder routes is blurred, and favorable transition states of both types may coexist. Application of the Curtin-Hammett principle, however, reveals that the bis-pericyclic mechanism accounts for ∼83% of the reaction flow in water, while the classical Diels-Alder mechanism contributes only ∼17%. The current findings provide a route for modeling this reaction inside the SpnF active site and inferring the catalytic architecture of possible Diels-Alderases.

An explicit account of solvation is essential for modeling Suzuki-Miyaura coupling in protic solvents.

We compared explicit and implicit solvation approaches in modeling the free energy profile of the final step of Suzuki-Miyaura coupling. Both approaches produced similar ΔG(≠) in all the studied solvents (benzene, toluene, DMF, ethanol, and water). Solvation free energies of individual reaction components reasonably correlated for explicit and implicit models in aprotic solvents (RMSE = 30-50 kJ mol(-1), R(2) > 0.71). However for ethanol and water the correlation was poor. We attributed this difference to the formation of the PdH-O hydrogen bond with Pd(PPh3)2 which was surprisingly observed in explicit modeling. Further QM calculations of the Pd(PPh3)2-H2O system confirmed the direction (PdH) and stability of this bonding. Therefore we stress the need for considering explicit solvation for modeling Pd-catalyzed reactions in protic solvents.

2,3-Dihydroxy-quinoxaline induces ATPase activity of Herpes Simplex Virus thymidine kinase.

2,3-Dihydroxy-quinoxaline, a small molecule that promotes ATPase catalytic activity of Herpes Simplex Virus thymidine kinase (HSV-TK), was identified by virtual screening. This compound competitively inhibited HSV-TK catalyzed phosphorylation of acyclovir with Ki=250 µM (95% CI: 106-405 µM) and dose-dependently increased the rate of the ATP hydrolysis with KM=112 µM (95% CI: 28-195 µM). The kinetic scheme consistent with this experimental data is proposed.

Molecular modeling of different substrate-binding modes and their role in penicillin acylase catalysis.

Molecular modeling was addressed to understand different substrate-binding modes and clarify the role of two positively charged residues of the penicillin G acylase active site - ßR263 and αR145 - in binding of negatively charged substrates. Although the electrostatic contribution to productive substrate binding was dominated by ßR263 rather than αR145, it was found that productive binding was not the only possible mode of substrate placement in the active site. Two extra binding modes - nonproductive and preproductive - were located by means of molecular docking and dynamics with binding affinities comparable with the productive one. A unique feature of nonproductive and preproductive complexes was that the substrate's acyl group did not penetrate the hydrophobic pocket, but occupied a patch on the protein interface spanning from ßR263 to αR145. Nonproductive and preproductive complexes competed with each other and productive binding mode, giving rise to increased apparent substrate binding. Preproductive complex revealed an ability to switch to a productive one during molecular dynamics simulations, and conformational plasticity of the penicillin G acylase active site was shown to be crucial for that. Nonproductive binding observed at molecular modeling corresponded well with experimentally observed substrate inhibition in penicillin acylase catalysis. By combining estimated free energies of substrate binding in each mode, and accounting for two possible conformations of the penicillin G acylase active site (closed and open) quantitative agreement with experimentally measured K(M) values was achieved. Calculated near-attack conformation frequencies from corresponding molecular dynamics simulations were in a quantitative correlation with experimental k(cat) values and demonstrated adequate application of molecular modeling methods.

Alchemical FEP Calculations of Ligand Conformer Focusing in Explicit Solvent.

Slow rotational degrees of freedom in ligands can make alchemical FEP simulations unreliable due to inadequate sampling. We addressed this problem by introducing a FEP-based protocol of ligand conformer focusing in explicit solvent. Our method involves FEP transformations between conformers using equilibrium dihedral angle as a reaction coordinate and provides the cost of "focusing" on one specific conformational state that binds to a protein. The calculated conformer focusing term made a considerable difference of 5-10 kJ/mol in computed relative binding free energies of studied Syk inhibitors and significantly improved the resulting accuracy of predictions.

Lead Finder docking and virtual screening evaluation with Astex and DUD test sets.

Lead Finder is a molecular docking software. Sampling uses an original implementation of the genetic algorithm that involves a number of additional optimization procedures. Lead Finder's scoring functions employ a set of semi-empiric molecular mechanics functionals that have been parameterized independently for docking, binding energy predictions and rank-ordering for virtual screening. Sampling and scoring both utilize a staged approach, moving from fast but less accurate algorithm versions to computationally more intensive but more accurate versions. Lead Finder includes tools for the preparation of full atom protein and ligand models. In this exercise, Lead Finder achieved 72.9% docking success rate on the Astex test set when the original author-prepared full atom models were used, and 74.1% success rate when the structures were prepared by Lead Finder. The major cause of docking failures were scoring errors resulting from the use of imperfect solvation models. In many cases, docking errors could be corrected by the proper protonation and the use of correct cyclic conformations of ligands. In virtual screening experiments on the DUD test set the early enrichment factor of several tens was achieved on average. However, the area under the ROC curve ("AUC ROC") ranged from 0.70 to 0.74 depending on the screening protocol used, and the separation from the null model was not perfect-0.12-0.15 units of AUC ROC. We assume that effective virtual screening in the whole range of enrichment curve and not just at the early enrichment stages requires more accurate solvation modeling and accounting for the protein backbone flexibility.

Hit clustering can improve virtual fragment screening: CDK2 and PARP1 case studies.

Virtual fragment screening could be a promising alternative to existing experimental screening techniques. However, reliable methods of in silico fragment screening are yet to be established and validated. In order to develop such an approach we first checked how successful the existing molecular docking methods can be in predicting fragment binding affinities and poses. Using our Lead Finder docking software the RMSD of the binding energy prediction was observed to be 1.35 kcal/mol(-1) on a set of 26 experimentally characterized fragment inhibitors, and the RMSD of the predicted binding pose from the experimental one was <1.5 Å. Then, we explored docking of 68 fragments obtained from 39 drug molecules for which co-crystal structures were available from the PDB. It appeared that fragments that participate in oriented non-covalent interactions, such as hydrogen bonds and metal coordination, could be correctly docked in 70-80% of cases suggesting the potential success of rediscovering of corresponding drugs by in silico fragment approach. Based on these findings we've developed a virtual fragment screening technique which involved structural filtration of protein-ligand complexes for specific interactions and subsequent clustering in order to minimize the number of preferable starting fragment candidates. Application of this method led to 2 millimolar-scale fragment PARP1 inhibitors with a new scaffold.

TSAR, a new graph-theoretical approach to computational modeling of protein side-chain flexibility: modeling of ionization properties of proteins.

A new graph-theoretical approach called thermodynamic sampling of amino acid residues (TSAR) has been elaborated to explicitly account for the protein side chain flexibility in modeling conformation-dependent protein properties. In TSAR, a protein is viewed as a graph whose nodes correspond to structurally independent groups and whose edges connect the interacting groups. Each node has its set of states describing conformation and ionization of the group, and each edge is assigned an array of pairwise interaction potentials between the adjacent groups. By treating the obtained graph as a belief-network-a well-established mathematical abstraction-the partition function of each node is found. In the current work we used TSAR to calculate partition functions of the ionized forms of protein residues. A simplified version of a semi-empirical molecular mechanical scoring function, borrowed from our Lead Finder docking software, was used for energy calculations. The accuracy of the resulting model was validated on a set of 486 experimentally determined pK(a) values of protein residues. The average correlation coefficient (R) between calculated and experimental pK(a) values was 0.80, ranging from 0.95 (for Tyr) to 0.61 (for Lys). It appeared that the hydrogen bond interactions and the exhaustiveness of side chain sampling made the most significant contribution to the accuracy of pK(a) calculations.

CSAR scoring challenge reveals the need for new concepts in estimating protein-ligand binding affinity.

The dG prediction accuracy by the Lead Finder docking software on the CSAR test set was characterized by R(2)=0.62 and rmsd=1.93 kcal/mol, and the method of preparation of the full-atom structures of the test set did not significantly affect the resulting accuracy of predictions. The primary factors determining the correlation between the predicted and experimental values were the van der Waals interactions and solvation effects. Those two factors alone accounted for R(2)=0.50. The other factors that affected the accuracy of predictions, listed in the order of decreasing importance, were the change of ligand's internal energy upon binding with protein, the electrostatic interactions, and the hydrogen bonds. It appears that those latter factors contributed to the independence of the prediction results from the method of full-atom structure preparation. Then, we turned our attention to the other factors that could potentially improve the scoring function in order to raise the accuracy of the dG prediction. It turned out that the ligand-centric factors, including Mw, cLogP, PSA, etc. or protein-centric factors, such as the functional class of protein, did not improve the prediction accuracy. Following that, we explored if the weak molecular interactions such as X-H...Ar, X-H...Hal, CO...Hal, C-H...X, stacking and π-cationic interactions (where X is N or O), that are generally of interest to the medicinal chemists despite their lack of proper molecular mechanical parametrization, could improve dG prediction. Our analysis revealed that out of these new interactions only CO...Hal is statistically significant for dG predictions using Lead FInder scoring function. Accounting for the CO...Hal interaction resulted in the reduction of the rmsd from 2.19 to 0.69 kcal/mol for the corresponding structures. The other weak interaction factors were not statistically significant and therefore irrelevant to the accuracy of dG prediction. On the basis of our findings from our participation in the CSAR scoring challenge we conclude that a significant increase of accuracy predictions necessitates breakthrough scoring approaches. We anticipate that the explicit accounting for water molecules, protein flexibility, and a more thermodynamically accurate method of dG calculation rather than single point energy calculation may lead to such breakthroughs.

Improving performance of docking-based virtual screening by structural filtration.

In the current study an innovative method of structural filtration of docked ligand poses is introduced and applied to improve the virtual screening results. The structural filter is defined by a protein-specific set of interactions that are a) structurally conserved in available structures of a particular protein with its bound ligands, and b) that can be viewed as playing the crucial role in protein-ligand binding. The concept was evaluated on a set of 10 diverse proteins, for which the corresponding structural filters were developed and applied to the results of virtual screening obtained with the Lead Finder software. The application of structural filtration resulted in a considerable improvement of the enrichment factor ranging from several folds to hundreds folds depending on the protein target. It appeared that the structural filtration had effectively repaired the deficiencies of the scoring functions that used to overestimate decoy binding, resulting into a considerably lower false positive rate. In addition, the structural filters were also effective in dealing with some deficiencies of the protein structure models that would lead to false negative predictions otherwise. The ability of structural filtration to recover relatively small but specifically bound molecules creates promises for the application of this technology in the fragment-based drug discovery.

Novel antitumor L-arabinose derivative of indolocarbazole with high affinity to DNA.

Novel indolocarbazole derivative 12-(alpha-L-arabinopyranosyl)indolo[2,3-alpha]pyrrolo[3,4-c]carbazole-5,7-dione (AIC) demonstrated high potency (at submicromolar concentrations) against the NCI panel of human tumor cell lines and transplanted tumors in vivo. In search of tentative targets for AIC, we found that the drug formed high affinity intercalative complexes with d(AT)(20), d(GC)(20) and calf thymus DNA (binding constants (1.6x10(6)) M(-1)< or =K(a)< or =(3.3x10(6)) M(-1)). The drug intercalated preferentially into GC pairs of the duplex. Importantly, the concentrations at which AIC formed the intercalative complexes with DNA (C< or =1 microM) were identical to the concentrations that triggered p53-dependent gene reporter transactivation, the replication block, the inhibition of topoisomerase I-mediated DNA relaxation and death of HCT116 human colon carcinoma cells. We conclude that the formation of high affinity intercalative complexes with DNA is an important factor for anticancer efficacy of AIC.

Developing novel approaches to improve binding energy estimation and virtual screening: a PARP case study.

Poly-(ADP-ribose)-polymerase (PARP) is a promising anti-cancer target as it plays a crucial role in the cellular reparation and survival mechanisms. However, the development of a robust and cost effective experimental technique to screen PARP inhibitors is still a scientific challenge owing to the difficulties in quantitative detection of the enzyme activity. In this work we demonstrate that the computational chemistry tools including molecular docking and scoring can perform on par with the experimental studies in assessing binding constants and in the recovery of active compounds in virtual screening. Using the recently introduced Lead Finder software we were able to dock a set of 142 well characterized PARP inhibitors and obtain a good correlation between the calculated and experimentally measured binding energies with the rmsd of 1.67 kcal mol(-1). Additionally, fine-tuning of the energy scaling coefficients within the Lead Finder scoring function has further decreased rmsd to the value of 0.88 kcal mol(-1). Moreover, we were able to reproduce the selectivity of ligand binding between the two isoforms of the enzyme-PARP1 and PARP2-suggesting that the Lead Finder software can be used to design isoform-selective inhibitors of PARP. An impressive enrichment was obtained in the virtual screening experiment, in which the mentioned set of PARP inhibitors was mixed with a commercial library of 300,000 compounds. We also demonstrate that the virtual screening performance can be significantly improved by an additional structural filtration of the docked ligand poses through detection of the crucial hydrogen bonding interactions with the enzyme.

Lead finder: an approach to improve accuracy of protein-ligand docking, binding energy estimation, and virtual screening.

An innovative molecular docking algorithm and three specialized high accuracy scoring functions are introduced in the Lead Finder docking software. Lead Finder's algorithm for ligand docking combines the classical genetic algorithm with various local optimization procedures and resourceful exploitation of the knowledge generated during docking process. Lead Finder's scoring functions are based on a molecular mechanics functional which explicitly accounts for different types of energy contributions scaled with empiric coefficients to produce three scoring functions tailored for (a) accurate binding energy predictions; (b) correct energy-ranking of docked ligand poses; and (c) correct rank-ordering of active and inactive compounds in virtual screening experiments. The predicted values of the free energy of protein-ligand binding were benchmarked against a set of experimentally measured binding energies for 330 diverse protein-ligand complexes yielding rmsd of 1.50 kcal/mol. The accuracy of ligand docking was assessed on a set of 407 structures, which included almost all published test sets of the following programs: FlexX, Glide SP, Glide XP, Gold, LigandFit, MolDock, and Surflex. rmsd of 2 A or less was observed for 80-96% of the structures in the test sets (80.0% on the Glide XP and FlexX test sets, 96.0% on the Surflex and MolDock test sets). The ability of Lead Finder to distinguish between active and inactive compounds during virtual screening experiments was benchmarked against 34 therapeutically relevant protein targets. Impressive enrichment factors were obtained for almost all of the targets with the average area under receiver operator curve being equal to 0.92.