Prediction of ligand binding mode among multiple cross-docking poses by molecular dynamics simulations.

We propose a method to identify the correct binding mode of a ligand with a protein among multiple predicted docking poses. Our method consists of two steps. First, five independent MD simulations with different initial velocities are performed for each docking pose, in order to evaluate its stability. If the root-mean-square deviations (RMSDs) of heavy atoms from the docking pose are larger than a given threshold (2.0 Å) in all five parallel runs, the pose is filtered out and discarded. Then, we perform accurate all-atom binding free energy calculations for the residual poses only. The pose with the lowest binding free energy is identified as the correct pose. As a test case, we applied our method to a previously built cross-docking test set, which included 104 complex systems. We found that the present method could successfully identify the correct ligand binding mode for 72% (75/104) of the complexes for current test set. The possible reasons for the failure of the method in the other cases were investigated in detail, to enable future improvements.

Uncovering abnormal changes in logP after fluorination using molecular dynamics simulations.

The fluorination-induced changes in the logP (1-octanol/water partition coefficient) of ligands were examined by molecular dynamics simulations. The protocol and force field parameters were first evaluated by calculating the logP values for n-alkanes, and their monofluorinated and monochlorinated analogs. Then, the logP values of several test sets (1-butanol, 3-propyl-1H-indole, and analogs fluorinated at the terminal methyl group) were calculated. The calculated results agree well with experiment, and the root mean square error values are 0.61 and 0.68 log units for the GAFF and GAFF2 force fields, respectively. Finally, the logP estimation was extended to a drug molecule, TAK-438, for which fluorination-induced abnormal logP reduction has been observed experimentally. This abnormal change was qualitatively reproduced by the molecular dynamics simulations. We found that the abnormal logP reduction can be mainly attributed to the effect of fluorination-induced dipole change. Our results suggest that molecular simulation is a useful strategy to predict the fluorination-induced change in logP for drug discovery applications.

Multiscale enhanced sampling of glucokinase: Regulation of the enzymatic reaction via a large scale domain motion.

Enhanced sampling yields a comprehensive structural ensemble or a free energy landscape, which is beyond the capability of a conventional molecular dynamics simulation. Our recently developed multiscale enhanced sampling (MSES) method employs a coarse-grained model coupled with the target physical system for the efficient acceleration of the dynamics. MSES has demonstrated applicability to large protein systems in solution, such as intrinsically disordered proteins and protein-protein and protein-ligand interactions. Here, we applied the MSES simulation to an important drug discovery target, glucokinase (GCK), to elucidate the structural basis of the positive cooperativity of the enzymatic reaction at an atomistic resolution. MSES enabled us to compare two sets of the free energy landscapes of GCK, for the glucose-bound and glucose-unbound forms, and thus demonstrated the drastic change of the free energy surface depending on the glucose concentration. In the glucose-bound form, we found two distinct basins separated by a high energy barrier originating from the domain motion and the folding/unfolding of the α13 helix. By contrast, in the glucose-unbound form, a single flat basin extended to the open and super-open states. These features illustrated the two distinct phases achieving the cooperativity, the fast reaction cycle staying in the closed state at a high glucose concentration and the slow cycle primarily in the open/super-open state at a low concentration. The weighted ensemble simulations revealed the kinetics of the structural changes in GCK with the synergetic use of the MSES results; the rate constant of the transition between the closed state and the open/super-open states, kC/O = 1.1 ms-1, is on the same order as the experimental catalytic rate, kcat = 0.22 ms-1. Finally, we discuss the pharmacological activities of GCK activators (small molecular drugs modulating the GCK activity) in terms of the slight changes in the domain motion, depending on their chemical structures as regulators. The present study demonstrated the capability of the enhanced sampling and the associated kinetic calculations for understanding the atomistic structural dynamics of protein systems in physiological environments.

Exploring the Stability of Ligand Binding Modes to Proteins by Molecular Dynamics Simulations: A Cross-docking Study.

Docking has become an indispensable approach in drug discovery research to predict the binding mode of a ligand. One great challenge in docking is to efficiently refine the correct pose from various putative docking poses through scoring functions. We recently examined the stability of self-docking poses under molecular dynamics (MD) simulations and showed that equilibrium MD simulations have some capability to discriminate between correct and decoy poses. Here, we have extended our previous work to cross-docking studies for practical applications. Three target proteins (thrombin, heat shock protein 90-alpha, and cyclin-dependent kinase 2) of pharmaceutical interest were selected. Three comparable poses (one correct pose and two decoys) for each ligand were then selected from the docking poses. To obtain the docking poses for the three target proteins, we used three different protocols, namely: normal docking, induced fit docking (IFD), and IFD against the homology model. Finally, five parallel MD equilibrium runs were performed on each pose for the statistical analysis. The results showed that the correct poses were generally more stable than the decoy poses under MD. The discrimination capability of MD depends on the strategy. The safest way was to judge a pose as being stable if any one run among five parallel runs was stable under MD. In this case, 95% of the correct poses were retained under MD, and about 25-44% of the decoys could be excluded by the simulations for all cases. On the other hand, if we judge a pose as being stable when any two or three runs were stable, with the risk of incorrectly excluding some correct poses, approximately 31-53% or 39-56% of the two decoys could be excluded by MD, respectively. Our results suggest that simple equilibrium simulations can serve as an effective filter to exclude decoy poses that cannot be distinguished by docking scores from the computationally expensive free-energy calculations.

Exploring the stability of ligand binding modes to proteins by molecular dynamics simulations.

The binding mode prediction is of great importance to structure-based drug design. The discrimination of various binding poses of ligand generated by docking is a great challenge not only to docking score functions but also to the relatively expensive free energy calculation methods. Here we systematically analyzed the stability of various ligand poses under molecular dynamics (MD) simulation. First, a data set of 120 complexes was built based on the typical physicochemical properties of drug-like ligands. Three potential binding poses (one correct pose and two decoys) were selected for each ligand from self-docking in addition to the experimental pose. Then, five independent MD simulations for each pose were performed with different initial velocities for the statistical analysis. Finally, the stabilities of ligand poses under MD were evaluated and compared with the native one from crystal structure. We found that about 94% of the native poses were maintained stable during the simulations, which suggests that MD simulations are accurate enough to judge most experimental binding poses as stable properly. Interestingly, incorrect decoy poses were maintained much less and 38-44% of decoys could be excluded just by performing equilibrium MD simulations, though 56-62% of decoys were stable. The computationally-heavy binding free energy calculation can be performed only for these survived poses.

Discovery of a Novel Series of Pyrazolo[1,5-a]pyrimidine-Based Phosphodiesterase 2A Inhibitors Structurally Different from N-((1S)-1-(3-Fluoro-4-(trifluoromethoxy)phenyl)-2-methoxyethyl)-7-methoxy-2-oxo-2,3-dihydropyrido[2,3-b]pyrazine-4(1H)-carboxamide (TAK-915), for the Treatment of Cognitive Disorders.

It has been hypothesized that selective inhibition of phosphodiesterase (PDE) 2A could potentially be a novel approach to treat cognitive impairment in neuropsychiatric and neurodegenerative disorders through augmentation of cyclic nucleotide signaling pathways in brain regions associated with learning and memory. Following our earlier work, this article describes a drug design strategy for a new series of lead compounds structurally distinct from our clinical candidate 2 (TAK-915), and subsequent medicinal chemistry efforts to optimize potency, selectivity over other PDE families, and other preclinical properties including in vitro phototoxicity and in vivo rat plasma clearance. These efforts resulted in the discovery of N-((1S)-2-hydroxy-2-methyl-1-(4-(trifluoromethoxy)phenyl)propyl)-6-methyl-5-(3-methyl-1H-1,2,4-triazol-1-yl)pyrazolo[1,5-a]pyrimidine-3-carboxamide (20), which robustly increased 3',5'-cyclic guanosine monophosphate (cGMP) levels in the rat brain following an oral dose, and moreover, attenuated MK-801-induced episodic memory deficits in a passive avoidance task in rats. These data provide further support to the potential therapeutic utility of PDE2A inhibitors in enhancing cognitive performance.

Design and synthesis of potent and selective pyridazin-4(1H)-one-based PDE10A inhibitors interacting with Tyr683 in the PDE10A selectivity pocket.

Utilizing structure-based drug design techniques, we designed and synthesized phosphodiesterase 10A (PDE10A) inhibitors based on pyridazin-4(1H)-one. These compounds can interact with Tyr683 in the PDE10A selectivity pocket. Pyridazin-4(1H)-one derivative 1 was linked with a benzimidazole group through an alkyl spacer to interact with the OH of Tyr683 and fill the PDE10A selectivity pocket. After optimizing the linker length, we identified 1-(cyclopropylmethyl)-5-[3-(1-methyl-1H-benzimidazol-2-yl)propoxy]-3-(1-phenyl-1H-pyrazol-5-yl)pyridazin-4(1H)-one (16f) as having highly potent PDE10A inhibitory activity (IC50=0.76nM) and perfect selectivity against other PDEs (>13,000-fold, IC50=>10,000nM). The crystal structure of 16f bound to PDE10A revealed that the benzimidazole moiety was located deep within the PDE10A selectivity pocket and interacted with Tyr683. Additionally, a bidentate interaction existed between the 5-alkoxypyridazin-4(1H)-one moiety and the conserved Gln716 present in all PDEs.

Design and synthesis of a novel 2-oxindole scaffold as a highly potent and brain-penetrant phosphodiesterase 10A inhibitor.

Highly potent and brain-penetrant phosphodiesterase 10A (PDE10A) inhibitors based on the 2-oxindole scaffold were designed and synthesized. (2-Oxo-1,3-oxazolidin-3-yl)phenyl derivative 1 showed the high P-glycoprotein (P-gp) efflux (efflux ratio (ER)=6.2) despite the potent PDE10A inhibitory activity (IC50=0.94 nM). We performed an optimization study to improve both the P-gp efflux ratio and PDE10A inhibitory activity by utilizing structure-based drug design (SBDD) techniques based on the X-ray crystal structure with PDE10A. Finally, 1-(cyclopropylmethyl)-4-fluoro-5-[5-methoxy-4-oxo-3-(1-phenyl-1H-pyrazol-5-yl)pyridazin-1(4H)-yl]-3,3-dimethyl-1,3-dihydro-2H-indol-2-one (19e) was identified with improved P-gp efflux (ER=1.4) and an excellent PDE10A inhibitory activity (IC50=0.080 nM). Compound 19e also exhibited satisfactory brain penetration, and suppressed PCP-induced hyperlocomotion with a minimum effective dose of 0.3mg/kg by oral administration in mice.

Two-dimensional replica-exchange method for predicting protein-ligand binding structures.

We have developed a two-dimensional replica-exchange method for the prediction of protein-ligand binding structures. The first dimension is the umbrella sampling along the reaction coordinate, which is the distance between a protein binding pocket and a ligand. The second dimension is the solute tempering, in which the interaction between a ligand and a protein and water is weakened. The second dimension is introduced to make a ligand follow the umbrella potential more easily and enhance the binding events, which should improve the sampling efficiency. As test cases, we applied our method to two protein-ligand complex systems (MDM2 and HSP 90-alpha). Starting from the configuration in which the protein and the ligand are far away from each other in each system, our method predicted the ligand binding structures in excellent agreement with the experimental data from Protein Data Bank much faster with the improved sampling efficiency than the replica-exchange umbrella sampling method that we have previously developed.

Peptide conformational preferences in osmolyte solutions: transfer free energies of decaalanine.

The nature in which the protecting osmolyte trimethylamine N-oxide (TMAO) and the denaturing osmolyte urea affect protein stability is investigated, simulating a decaalanine peptide model in multiple conformations of the denatured ensemble. Binary solutions of both osmolytes and mixed osmolyte solutions at physiologically relevant concentrations of 2:1 (urea:TMAO) are studied using standard molecular dynamics simulations and solvation free energy calculations. Component analysis reveals the differences in the importance of the van der Waals (vdW) and electrostatic interactions for protecting and denaturing osmolytes. We find that urea denaturation governed by transfer free energy differences is dominated by vdW attractions, whereas TMAO exerts its effect by causing unfavorable electrostatic interactions both in the binary solution and mixed osmolyte solution. Analysis of the results showed no evidence in the ternary solution of disruption of the correlations among the peptide and osmolytes, nor of significant changes in the strength of the water hydrogen bond network.

Ab initio prediction of protein-ligand binding structures by replica-exchange umbrella sampling simulations.

We have developed a prediction method for the binding structures of ligands with proteins. Our method consists of three steps. First, replica-exchange umbrella sampling simulations are performed along the distance between a putative binding site of a protein and a ligand as the reaction coordinate. Second, we obtain the potential of mean force (PMF) of the unbiased system using the weighted histogram analysis method and determine the distance that corresponds to the global minimum of PMF. Third, structures that have this global-minimum distance and energy values around the average potential energy are collected and analyzed using the principal component analysis. We predict the binding structure as the global-minimum free energy state on the free energy landscapes along the two major principal component axes. As test cases, we applied our method to five protein-ligand complex systems. Starting from the configuration in which the protein and the ligand are far away from each other in each system, our method predicted the ligand binding structures in excellent agreement with the experimental data from Protein Data Bank.

Trimethylamine N-oxide influence on the backbone of proteins: an oligoglycine model.

The study of organic osmolytes has been pivotal in demonstrating the role of solvent effects on the protein backbone in the folding process. Although a thermodynamic description of the interactions between the protein backbone and osmolyte has been well defined, the structural analysis of the effect of osmolyte on the protein backbone has been incomplete. Therefore, we have performed simulations of a peptide backbone model, glycine(15), in protecting osmolyte trimethylamine N-oxide (TMAO) solution, in order to determine the effect of the solution structure on the conformation of the peptide backbone. We show that the models chosen show that the ensemble of backbone structures shifts toward a more collapsed state in TMAO solution as compared with pure water solution. The collapse is consistent with preferential exclusion of the osmolyte caused by unfavorable interactions between osmolyte and peptide backbone. The exclusion is caused by strong triplet correlations of osmolyte, water, and peptide backbone. This provides a clear mechanism showing that even a modest concentration of TMAO forces the protein backbone to adopt a more collapsed structure in the absence of side chain effects.

Dependency of ligand free energy landscapes on charge parameters and solvent models.

We performed replica-exchange molecular dynamics (REMD) simulations of six ligands to examine the dependency of their free energy landscapes on charge parameters and solvent models. Six different charge parameter sets for each ligand were first generated by RESP and AM1-BCC methods using three different conformations independently. RESP charges showed some conformational dependency. On the other hand, AM1-BCC charges did not show conformational dependency and well reproduced the overall trend of RESP charges. The free energy landscapes obtained from the REMD simulations of ligands in vacuum, Generalized-Born (GB), and TIP3P solutions were then analyzed. We found that even small charge differences can produce qualitatively different landscapes in vacuum condition, but the differences tend to be much smaller under GB and TIP3P conditions. The simulations in the GB model well reproduced the landscapes in the TIP3P model using only a fraction of the computational cost. The protein-bound ligand conformations were rarely the global minimum states, but similar conformations were found to exist in aqueous solution without proteins in regions close to the global minimum, local minimum or intermediate states.

Analysis of helix-helix interactions of bacteriorhodopsin by replica-exchange simulations.

We performed long-time replica-exchange Monte Carlo simulations of bacteriorhodopsin transmembrane helices, which made it possible that wide conformational space was sampled. Using only the helix-helix interactions and starting from random initial configurations, we obtained the nativelike helix arrangement successfully and predicted a part of the configurations (three helices out of seven) precisely. By the principal component analysis we classified low-energy structures into some clusters of similar structures, and we showed that the above nativelike three-helix configuration is reproduced properly in most clusters and that not only the van der Waals interactions but also the electrostatic interactions contributed to the stabilization of the native structures.

Discovery of 3-Benzyl-1-( trans-4-((5-cyanopyridin-2-yl)amino)cyclohexyl)-1-arylurea Derivatives as Novel and Selective Cyclin-Dependent Kinase 12 (CDK12) Inhibitors.

Cyclin-dependent kinase 12 (CDK12) plays a key role in the coordination of transcription with elongation and mRNA processing. CDK12 mutations found in tumors and CDK12 inhibition sensitize cancer cells to DNA-damaging reagents and DNA-repair inhibitors. This suggests that CDK12 inhibitors are potential therapeutics for cancer that may cause synthetic lethality. Here, we report the discovery of 3-benzyl-1-( trans-4-((5-cyanopyridin-2-yl)amino)cyclohexyl)-1-arylurea derivatives as novel and selective CDK12 inhibitors. Structure-activity relationship studies of a HTS hit, structure-based drug design, and conformation-oriented design using the Cambridge Structural Database afforded the optimized compound 2, which exhibited not only potent CDK12 (and CDK13) inhibitory activity and excellent selectivity but also good physicochemical properties. Furthermore, 2 inhibited the phosphorylation of Ser2 in the C-terminal domain of RNA polymerase II and induced growth inhibition in SK-BR-3 cells. Therefore, 2 represents an excellent chemical probe for functional studies of CDK12 and could be a promising lead compound for drug discovery.

Preferential solvation in urea solutions at different concentrations: properties from simulation studies.

We performed molecular dynamics simulations of urea solutions at different concentrations with two urea models (OPLS and KBFF) to examine the structures responsible for the thermodynamic solution properties. Our simulation results showed that hydrogen-bonding properties such as the average number of hydrogen bonds and their lifetime distributions were nearly constant at all concentrations between infinite dilution and the solubility limit. This implies that the characterization of urea-water solutions in the molarity concentration scale as nearly ideal is a result of facile local hydrogen bonding rather than a global property. Thus, urea concentration does not influence the local propensity for hydrogen bonds, only how they are satisfied. By comparison, the KBFF model of urea donated fewer hydrogen bonds than OPLS. We found that the KBFF urea model in TIP3P water better reproduced the experimental density and diffusion constant data. Preferential solvation analysis showed that there were weak urea-urea and water-water associations in OPLS solution at short distances, but there were no strong associations. We divided urea molecules into large, medium, and small clusters to examine fluctuation properties and found that any particular urea molecule did not stay in the same cluster for a long time. We found neither persistent nor large clusters.

Discovery of Clinical Candidate N-((1S)-1-(3-Fluoro-4-(trifluoromethoxy)phenyl)-2-methoxyethyl)-7-methoxy-2-oxo-2,3-dihydropyrido[2,3-b]pyrazine-4(1H)-carboxamide (TAK-915): A Highly Potent, Selective, and Brain-Penetrating Phosphodiesterase 2A Inhibitor for the Treatment of Cognitive Disorders.

Phosphodiesterase (PDE) 2A inhibitors have emerged as a novel mechanism with potential therapeutic option to ameliorate cognitive dysfunction in schizophrenia or Alzheimer's disease through upregulation of cyclic nucleotides in the brain and thereby achieve potentiation of cyclic nucleotide signaling pathways. This article details the expedited optimization of our recently disclosed pyrazolo[1,5-a]pyrimidine lead compound 4b, leading to the discovery of clinical candidate 36 (TAK-915), which demonstrates an appropriate combination of potency, PDE selectivity, and favorable pharmacokinetic (PK) properties, including brain penetration. Successful identification of 36 was realized through application of structure-based drug design (SBDD) to further improve potency and PDE selectivity, coupled with prospective design focused on physicochemical properties to deliver brain penetration. Oral administration of 36 demonstrated significant elevation of 3',5'-cyclic guanosine monophosphate (cGMP) levels in mouse brains and improved cognitive performance in a novel object recognition task in rats. Consequently, compound 36 was advanced into human clinical trials.

Prediction of Ligand Binding Affinity by the Combination of Replica-Exchange Method and Double-Decoupling Method.

A prediction method for ligand binding affinities to proteins is proposed. We first predict the structures of protein-ligand complex by the replica-exchange umbrella sampling or its extension. We then calculate ligand binding affinities based on these predicted ligand-protein bound structures by the double-decoupling method. As a test of the effectiveness of the proposed method, we applied it to the system of the oncoprotein MDM2 and a ligand. The value of the predicted binding affinity turned out to be in good agreement with that from experiments.

Discovery of 1-[2-fluoro-4-(1H-pyrazol-1-yl)phenyl]-5-methoxy-3-(1-phenyl-1H-pyrazol-5-yl)pyridazin-4(1H)-one (TAK-063), a highly potent, selective, and orally active phosphodiesterase 10A (PDE10A) inhibitor.

A novel series of pyridazinone-based phosphodiesterase 10A (PDE10A) inhibitors were synthesized. Our optimization efforts using structure-based drug design (SBDD) techniques on the basis of the X-ray crystal structure of PDE10A in complex with hit compound 1 (IC50 = 23 nM; 110-fold selectivity over other PDEs) led to the identification of 1-[2-fluoro-4-(1H-pyrazol-1-yl)phenyl]-5-methoxy-3-(1-phenyl-1H-pyrazol-5-yl)pyridazin-4(1H)-one (27h). Compound 27h has potent inhibitory activity (IC50 = 0.30 nM), excellent selectivity (>15000-fold selectivity over other PDEs), and favorable pharmacokinetics, including high brain penetration, in mice. Oral administration of compound 27h to mice elevated striatal 3',5'-cyclic adenosine monophosphate (cAMP) and 3',5'-cyclic guanosine monophosphate (cGMP) levels at 0.3 mg/kg and showed potent suppression of phencyclidine (PCP)-induced hyperlocomotion at a minimum effective dose (MED) of 0.3 mg/kg. Compound 27h (TAK-063) is currently being evaluated in clinical trials for the treatment of schizophrenia.

Prediction of Protein-Ligand Binding Structures by Replica-Exchange Umbrella Sampling Simulations: Application to Kinase Systems.

We have applied our prediction method, which is based on the replica-exchange umbrella sampling for protein-ligand binding structures, to two kinase systems (p38 and JNK3) with two different ligand molecules for each kinase. Starting from configurations in which the protein and the ligand are far away from each other, our method predicted the ligand binding structures in excellent agreement with the experimental data from PDB in all four cases, which suggests the general applicability of our method to kinase systems. In addition, the protein flexibility was shown to be essential to predict the correct binding structure for one of the systems, where dihydroquinolinone was bound to p38 alpha kinase (PDB ID: 1OVE ).