Relaxation and Intermolecular

Subscriber access provided by UNIV OF DURHAM

Pharmaceutical Modeling

Combining Stochastic Deformation/Relaxation and Intermolecular Contacts Analysis for Extracting Pharmacophores from Ligand-Receptor Complexes. Ma'mon M. Hatmal, and Mutasem O Taha J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.7b00708 • Publication Date (Web): 12 Mar 2018 Downloaded from http://pubs.acs.org on March 13, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling 1

Combining Stochastic Deformation/Relaxation and Intermolecular Contacts Analysis for Extracting Pharmacophores from Ligand-Receptor Complexes.

Ma'mon M. Hatmala and Mutasem O. Taha*b

a

Department of Medical Laboratory Sciences, Faculty of Allied Health Sciences, The

Hashemite University, P.O. Box 330127, Zarqa 13133, Jordan. b

Drug Discovery Unit, Department of Pharmaceutical Sciences, Faculty of Pharmacy,

University of Jordan, Amman 11942, Jordan.

*

To whom correspondence should be addressed: Drug Discovery Unit, Department of

Pharmaceutical Sciences, Faculty of Pharmacy, University of Jordan, Amman, Jordan. Tel.: 00962 6 5355000 ext. 23305 Fax: 00962 6 5339649 E-mail address: [email protected] (Mutasem O. Taha, Ph.D.).

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 2 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ABSTRACT We previously combined molecular dynamics (classical or simulated annealing) with ligand-receptor contacts analysis as means to extract valid pharmacophore model(s) from single ligand-receptor complexes. However, molecular dynamics methods are computationally expensive and time consuming. Here we describe a novel method for extracting valid pharmacophore model(s) from a single crystallographic structure within reasonable time scale. The new method is based on ligand-receptor contacts analysis following energy relaxation of predetermined set of randomly deformed complexes generated from the targeted crystallographic structure. Ligand-receptor contacts maintained across many deformed/relaxed structures are assumed to be critical and used to guide pharmacophore development. This methodology was implemented to develop valid pharmacophore models for PI3K-γ, RENIN and JAK1. The resulting pharmacophore models were validated by receiver operating characteristic (ROC) analysis against inhibitors extracted from CHEMBL database. Additionally, we implemented pharmacophores extracted from PI3K-γ to search for new inhibitors from the national cancer institute list of compounds. The process culminated in new PI3K-γ/mTOR inhibitory leads of low micromolar IC50s.


Page 2 of 39

Page 3 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1. INTRODUCTION A pharmacophore is abstract description of steric and electronic features necessary to be present in a particular ligand to allow optimal recognition and binding to specific receptor and to trigger subsequent biological function.1-9 Pharmacophore models are derived either by structure-based or ligand-based approaches.1-17 In structure-based methodologies, pharmacophore models are extracted from crystallographic complexes involving high affinity ligands co-crystallized within targeted proteins.21 However, despite that ligandreceptor binding involves numerous attractive interactions, only few can be defined as pharmacophoric, i.e., critical for complex stability.1-32 It is rather challenging to identify critical binding interactions from single crystallographic complex because of lack of explicit correlation connecting each interaction with bioactivity loss or gain. In fact, structure-based pharmacophores extracted from single crystallographic structures can be too limiting or promiscuous22-25 by including redundant binding features or missing essential binding interactions, respectively. In our quest to achieve valid structure-based pharmacophore models we attempted to discern critical (pharmacophoric) ligand-receptor interactions from redundant ones through two general approaches: (i) QSAR-based approach whereby binding interactions that can explain bioactivity variations among numerous docked ligands are selected as being essential pharmacophoric interactions (e.g., dbCICA).26-31 (ii) Molecular dynamics-based approaches whereby the particular crystallographic complex is exposed to conventional (classical) or drastic (annealing) molecular dynamics simulations to identify resilient ligand-receptor binding interactions that frequently emerge during molecular perturbations (md-LRCA or sa-LRCA).22,23 Such interactions are assumed to represent basic pharmacophoric requirements. Nevertheless, these approaches were either limited



by the need for long lists of bioactive ligands (i.e., db-CICA) or were computationally expensive and time-consuming (md-LRCA ad sa-LRCA).22,23 From our experience in implementing molecular dynamics for gauging the stability and resilience of different ligand-receptor binding interactions, and therefore their pharmacophoric significance, we noticed that simulated annealing (i.e., in sa-LRCA22) yielded superior pharmacophore models compared to classical dynamics (i.e., in md-LRCA23). We believe this is probably because excessive simulated heating (in annealing dynamics) represents dramatic challenge for ligand-protein binding interactions and allow better exploration of the energy surface of the complex. Accordingly, we decided to reproduce a similar situation to that of simulated annealing, i.e., to artificially generate numerous deformed ligand-protein complexes that resemble heated states in simulated annealing setup and allow them to energy-minimize similar to the cooling cycles in annealing dynamics. Interatomic contacts that resist separation in response to deformation/relaxation (D/R) among majority of deformed complexes should represent critical (pharmacophoric) interactions.

The "artificial" deformation of the particular ligand-receptor complex

should be computationally facile compared to the equivalent time-consuming excessive heating step in simulated annealing. Moreover, subsequent energy-minimization, performed via traditional gradient energy minimization, is also simpler than cooling dynamics implemented in simulated annealing. Nevertheless, tandem structural deformationminimization should provide good testing for the robustness of different ligand-receptor interaction and therefore allow useful deductions about the pharmacophoric (essential) interactions within a particular ligand-receptor complex. We tested the validity of this new method, i.e., deformation/relaxation ligand-receptor contacts analysis (dr-LRCA), against three protein complexes: PI3K-γ, renin and JAK1. PI3K-γ is a member of PI3Ks family involved signal transduction of wide range of sig-


Page 4 of 39

Page 5 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


naling pathways involved in cell growth, survival, differentiation and motility.32-34 Frequent activating mutations in PI3K-γ have recently been implicated in several cancers rendering PI3K-γ interesting drug target.34-37 On the other hand, renin is an aspartic protease secreted by the kidneys and involved in regulating blood pressure. Renin inhibitors are useful for the treatment of hypertension.38 Finally, JAK1 is a protein tyrosin kinase involved in certain leukaemias, lymphomas and polycythaemia.39,40 A JAK1 inhibitor has received FDA approval in 2011 for the treatment of polycythaemia.39-42

METHODOLOGY 1.1. Protein Preparation The three-dimensional (3D) structures of PI3K-γ, renin and JAK1 complexed with potent ligands were downloaded from protein data bank (PDB codes: 3lj3, 1bil and 4ivb, respectively). The proteins were hydrated with TIP3 molecules, neutralized with NaCl ions and energy-minimized over 30000 conjugate gradient steps. Then the temperature of the system was raised linearly in the heating phase from 0 to 310 K over 620,000 steps over a time step of 0.0005 ps. Subsequently, the system was equilibrated over 1 ns. Simulations were run using CGenFF force field for ligands (https://cgenff.paramchem.org) and CHARMM22 forcefield for proteins. SHAKE constrains were applied to all hydrogen atoms.22,23

1.2. Complex Deformation/Relaxation Cycles Deformation/relaxation cycles were preformed through the following steps: 1) Each protein complex obtained after equilibration was used to generate 500 different deformed complexes as follows: Each torsion (dihedral) angle within rotable amino acid side chain (except glycine, alanine, and proline) was randomly assigned a value ranging



from 0˚ to 180˚ with 5˚ increments. If two atoms occupied the same space (inter atomic distance ≤ 0.05 Ǻ) in certain deformed structure, the deformation process is repeated with another set of random dihedral angles. 2) To ensure that each deformed complex represents a unique conformation, the complexes were aligned and their RMSD values were calculated (determined for amino acid side chains only, i.e., not protein backbone). Deformed complexes exceeding RMSD threshold of 2.5 Ǻ (against all other conformers) were retained, while those having RMSD ≤ 2.5Ǻ were considered very similar and were therefore discarded except one. Discarded deformed complexes were replaced with new randomly deformed complexes and rechecked for RMSD similarity. This step is done to avoid erroneous frequent ligandreceptor contacts resulting from coincidental emergence of multiple similar deformed protein complexes that converge into close minimized states. Emergence of multiple close conformers should have detrimental effects on the signal-to-noise ratio required to discern essential interactions from redundant ones. 3) Each deformed complex (of 500 complex conformers) is exposed to two distinct relaxation cycles based on conjugate gradient minimization (as implemented in NAMD software, http://www.ks.uiuc.edu/Research/namd/). In the first step, the ligand is minimized over 100 steps while the protein host is kept fixed. This step is done to allow the ligand to adapt into the new conformation of the protein. However, in the second step the whole complex (including the ligand) is minimized for further 300 steps. 4) Energy-minimized complexes are then aligned and their main chain (backbone) RMSD values against the original (un-deformed) crystallographic complex were calculated. Minimized protein conformers of RMSD within 2.5 Ǻ from the un-deformed complex were retained, while those having RMSD value > 2.5Ǻ were considered drastically different from the original structures and were discarded. Discarded complexes are replaced


Page 6 of 39

Page 7 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


with new deformed/relaxed complexes and rechecked for RMSD similarity with the original complex. This step is done to avoid severely deformed conformational substates that are probably un-natural and inaccessible to protein complexes under physiological conditions (probably violate Ramachadran criteria). The RMSD limits in this step are derived from our experience in simulated annealing molecular dynamics performed of variety of protein complexes.23 2.2. Contacts Analysis and Subsequent Pharmacophore Generation and Validation This step was performed as we described earlier in coupling molecular dynamics or simulated annealing with intermolecular contacts analysis (i.e., sa-LRCA and md-LRCA).22,23 Briefly: 500 Deformed/relaxed frames are collected for the particular ligand-receptor complex. The ligand conformation in each frame is evaluated to identify its closest atomic neighbors within the binding pocket. Intermolecular atomic neighbors ≤ 2.5 Å are assigned an intermolecular contact value of "one", otherwise they are assigned a contact value of "zero". Then a two-dimensional binary matrix is generated with row labels corresponding to deformed/relaxed ligand-protein conformer frames and column labels corresponding to binding site atoms. Frequently encountered binding site atoms (i.e., by ligand atoms) among all 500 deformed/minimized protein conformer frames are identified and used to build corresponding binding pharmacophores. Section SM-1 under Supporting Information describes in details how intermolecular contacts analysis is performed and used to generate corresponding pharmacophore models. The resulting pharmacophores were validated using receiver-operating characteristic (ROC) curve analysis, which evaluates the ability of tested pharmacophore to discriminate active compounds from a list composed of known 50 active ligands, i.e., against the targeted enzyme (i.e., PI3K-γ, Renin or JAK1, IC50 ≤ 1.0 µM), and a larger list of decoys (300-750 compounds, IC50 range from 2 to 10,000 µM).43-47 The testing sets were assembled from CHEMBL



database.48 Details of the ROC analysis are shown under Supporting Information Section SM-2.8-17,22,23,26-31,43-47

2.3 In silico screening for new PI3K-γ/mTOR inhibitors 5-Featured PI3K-γ pharmacophores originating from either simulated annealing-contacts analysis (sa-LRCA), molecular dynamics-contacts analysis (md-LRCA), or deformation/relaxation-intermolecular contacts analysis (dr-LRCA) methods were employed as 3D search queries to screen the National Cancer Institute (NCI) list of compounds for new inhibitors. Captured hits were filtered by the SMARTS filter implemented in DiscoverStudio 2.5 to remove compounds of reactive and undesirable functional groups (e.g., alkyl halides, epoxides, aziridine, sulfonate esters, long aliphatic chain, azides, etc.,). The remaining hits were fitted against corresponding capturing pharmacophores using the pharmacophore-ligand mapping protocol in DiscoverStudio 2.5.3-17,22,23,26-32 The highest ranking available hits from the NCI (40 compounds) originate from each method (sa-LRCA, md-LRCA, or dr-LRCA) based on their respective fit values were evaluated against PI3K-γ and mTOR in vitro.

2.4 In vitro testing of captured hits against PKC-θ NCI samples were kindly provided by the National Cancer Institute. Effect of the indicated compounds on the activity of PI3K-γ and mTOR were assessed by SelectScreen Kinase Profiling Service (Invitrogen-Life Technologies, USA). Assays were performed using 10 µM of the tested compounds. Compounds showing inhibitory percentages exceeding 50% at 10 µM were further evaluated at least 4 concentrations ranging from 0.01 to 100 µM to determine their IC50 values. IC50 values were obtained using nonlinear


Page 8 of 39

Page 9 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


regression of the log(concentration) versus inhibition percentages using GraphPad Prism 5.0 (http://www.graphpad.com).22,23

3. RESULTS AND DISCUSSION 3.1. Basic Concept of dr-LRCA. Identification of critical binding interaction(s) within certain ligand-protein complex is rather daunting because many of the observed interactions are virtually redundant.3-17,22-32 For example, the crystallographic structure of inhibitor-renin complex studied herein (See Table 2 below) generates 13-feaured pharmacophore based on interactions perceived from close electronically-complementary groups within the complex. Needless to say that such large pharmacophore is too limiting to act as successful search query towards new bioactive hits. Therefore, 5-featured derivatives from this parent model should act as reasonable alternatives, i.e., as 3D search queries. However, the fact that featurereduction is not supervised by bioactivity trends renders all pharmacophoric features of equivalent weight and any reasonable feature-reduction process should include all possible 5-featured pharmacophore combinations. Therefore, the 13-featured pharmacophore of renin inhibitor should result in 1287 possible unique 5-featured models (corresponding to 

 ). 13!  5!×(13 − 5)!   

This number of pharmacophores is rather extensive to assess and to draw

conclusions about the strengths and significances of individual attractive interactions conceived from the crystallographic structure. Moreover, pharmacophore model(s) extracted from single ligand-protein crystallographic structure might miss essential interactions seen in protein complexes with other ligands, which necessitates similar combinatorial evaluations of crystallographic complexes tying the same protein with other ligands, a luxury not always available.



We previously implemented molecular dynamics approaches (classical or annealing) to perturb ligand-receptor complexes as means to identify critical binding features. Binding interactions that resisted breakage under such challenging circumstances were considered essential and were translated via ligand-protein intermolecular contacts analysis into valid pharmacophoric models.22,23 Nevertheless, molecular dynamics approaches are computationally expensive, which undermines their practicability. Stochastic rotamer libraries of amino acid side chains have been used in docking studies49-52 and in calculating ligand-receptor binding energies.53 The implementation of stochastic selection of amino acid rotamers has the advantage of quick conformational sampling and ability to surmount energy barriers, thus allowing better exploration of complicated energy surface of ligand-protein complexes.54 In fact, stochastic methods can be thought of as facile alternatives to computationally-intensive molecular dynamics methods (classical or simulated annealing) in exploring protein energy minima.55 This encouraged us to couple stochastic sampling of binding site amino acid rotamers with our previously published ligand-protein intermolecular contacts analysis22,23,26-31 to identify critical pharmacophoric binding features. This is the first time to couple stochastic deformation/relaxation methods with intermolecular contacts analysis as means to extract valid pharmacophore models. We termed the new combination deformation/relaxation intermolecular contacts analysis (dr-LRCA). The concept is rather simple: The ability of certain ligand-protein contact to survive across many randomly deformed/relaxed protein structures should point to the resilience of the corresponding surrogate binding interaction(s). Accordingly, top surviving set of contacts should encode for basic pharmacophoric requirements of ligand-protein binding. GRID-based energy probes can be also be used to identify critical attractive interactions


Page 10 of 39

Page 11 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


within protein complexes,24,25 still, they can be computationally costly for numerous protein-ligand frames. However, the stochastic process may artificially over-sample energy minima that are inaccessible under ambient energetic perturbations. Such structural artifacts can lead to artificial frequent intermolecular contacts leading to deducing erroneous pharmacophoric features. Accordingly, we implemented several precautionary steps to minimize the influence of artificial energy minima on intermolecular contacts frequencies. Firstly, the protein-ligand complex was properly hydrated, neutralized with NaCl ions, energy minimized, heated to 310 K and equilibrated using molecular dynamic simulation. It is anticipated that the resulting 'starting' protein conformation closely represents the natural state, which should minimize the possibility of severe conformational deformation and entrapment within artificial potential energy minima. Furthermore, we only allowed amino acid side chains to be involved in stochastic sampling and avoided backbone modifications to avoid drastic unnatural conformational changes. However, we allowed the protein backbone to be involved in subsequent energy minimization to accommodate for minor backbone conformational changes that accompany normal protein perturbations.56-60 Secondly, we compared the amino acid side chains of deformed structures (i.e., before minimization) against each other. Similar deformed structures (i.e., of RMSD ≤ 2.5 Ǻ) were eliminated (i.e., only one is retained). This step should also reduce the possibility of repeated artificial energy minima resulting from similar starting distorted structures, which should enhance the signal-to-noise ratio among frequent contactors. Thirdly, assuming that artificial local minima will probably violate the natural conformational trends of proteins (e.g., Ramachandran criteria), it was decided to exclude any deformed/relaxed protein structure(s) that have its (their) backbone(s) deviating signifi-



cantly from the initial protein backbone (RMSD > 2.5 Ǻ). This step guarantees that resulting deformed/minimized siblings remain reasonably relevant and within the conformational vicinity of the starting protein-ligand complex. As a final precaution to minimize the emergence of unnatural conformational minima, it was decided to generate relatively large number of deformed/relaxed protein-ligand conformational states (500 conformers). This should dilute the influence of artificial discrete minima (and their surrogate intermolecular ligand-protein contacts) while enhance the influence of continuous natural minima that tend to converge. Figures 1-3 compare some randomly selected conformational states resulting from implementing our deformation/relaxation algorithm on ligand-PI3K-γ complex with randomly selected conformational states achieved via simulated annealing and classical molecular dynamics, respectively. On the other hand, figure 4 compares the total energies of the same protein complex calculated for 100 sequential frames generated by the production phase of classical molecular dynamics simulation (over 1.25 ns) and simulated annealing molecular dynamics (100 cycles) compared to 100 structures generated using our deformation/relaxation algorithm. Evaluation of figures 1-3 suggests that the energy substates achieved by the three methods are conformationally related as evident from the close similarity between binding amino acids' conformations generated by all three methods. For example, the hydrogen bonding interactions connecting Val882 and Asp841 with the bound ligand remained intact in most conformations generated by the three methods (Figures 1-3) while other interactions tend to frequently break in the three approaches (e.g., ligand hydrogen bonding to Lys833). Nevertheless, figure 4 shows that our deformation/relaxation protocol accessed more stable conformational energies compared to simulated annealing or classical molecular


Page 12 of 39

Page 13 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


dynamics. Moreover, the deformation/relaxation protocol seems to explore more diverse energy minima compared to both molecular dynamics methods as evident from the conformational energy fluctuations seen among conformers of the deformation/relaxation trait. This trend is probably because the conformers were randomly generated in the deformation/relaxation approach, while they are successively generated in the molecular dynamics methods (i.e., each conformational substate depends on its precursor). Accordingly, one can conclude from figures 1-4 that the three methods explore relevant conformational spaces; however, the lower energy profile achievable by our algorithm indicates its accessibility to more stable protein conformers allowing pharmacophorically relevant intermolecular interactions (i.e., critical binding interactions) to emerge more frequently among generated complex ligand-protein conformers. Incidentally, our stochastic deformation/relaxation conformation-generator is related to Monte Carlo molecular dynamics55,61 albeit computationally simpler and tends to access more drastic "unnatural" conformations more suited to gauge the robustness of pharmacophoric features.

3.2. Case Studies We tested the new approach (i.e., dr-LRCA) against three crystallographic ligand-protein complexes, namely, PI3K-γ (PDB code: 3lj3) and JAK1 (PDB code: 4ivb) and renin (PDB code: 1bil). Moreover, we compared dr-LRCA-based pharmacophores extracted from PI3K-γ with corresponding models obtained via sa-LRCA and md-LRCA.22,23 dr-LRCA modeling runs were performed through 500 virtual structures generated by tandem random deformation-relaxation on each protein-ligand complex (PI3K-γ, JAK1 and renin), while sa-LRCA and md-LRCA were performed on PI3K-γ complex over 50 SAMD cycles (sa-LRCA) and 7 ns (md-LRCA). The structural frames of each de-



formed/relaxed complex (dr-LRCA), at the end of each SAMD cycle (sa-LRCA), or at every 0.0125 ns during MD production (md-LRCA) were collected for contacts analysis whereby the ligand's pose/conformer in each frame is evaluated to identify its closest atomic neighbors within the binding pocket. Intermolecular atomic neighbors at distances ≤ 2.5 Å were assigned a contact value of "one", otherwise they were given a contact value of "zero".22,23,26-30 This threshold was implemented because reversible attractive interactions are most pronounced at this distance. 22,23,26-30 Each binding site atom is given contact frequency (CF) value defined as the percentage of frames in which the ligand contacts this atom (i.e., at distance ≤ 2.5 Å) compared to the number of evaluated frames. Top ranking binding site atoms were considered frequent contactors and were assumed to encode for essential proxy ligand-protein binding interactions. These were used to build corresponding pharmacophore models (as in Supporting Information Section SM-1). Table 1 shows the highest ranking binding site atoms and their corresponding CF values as determined by the three methods (sa-LRCA, md-LRCA and dr-LRCA) in the PI3K-γ case, while table 2 shows the highest ranking binding atoms for renin and JAK1-kinase produced by dr-LRCA method. The two tables also show the corresponding pharmacophores generated based on the highest-ranking binding site contacts. Frames of highest weighted-contacts sums (WCS) were selected to project pharmacophoric features (see point 6 under section SM-1 in Supporting Information).22,23


Page 14 of 39

Page 15 of 39

Journal of Chemical Information and Modeling

15

1 (C) (B) 2 (A) 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 (D) (F) (E) 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 Figure 1: PI3Kγ (PDB ID: 3lj3) binding site residues closely contacting the co-crystallized ligand and how they perturb in some randomly selected protein 42 frames after random deformation-relaxation. Green lines represent hydrogen bonding. (A) Binding site residues in the original crystallographic structure of 43 PI3K-γ. (B), (C), (D), (E), and (F) are snapshots taken for deformed/relaxed protein/ligand complex (number 100, 200, 300, 400, and 500, respectively). 44 ACS Paragon Plus Environment 45 46 47


16

Page 16 of 39

1 (C) (A) (B) 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 (F) (E) 21 (D) 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 Figure 2: PI3Kγ (PDB ID: 3lj3) binding site residues closely contacting the co-crystallized ligand and how they perturb in protein frames at the end of some 43 simulated annealing cycles. Green lines represent hydrogen bonding. (A) Binding site residues in the original crystallographic structure of PI3K-γ. (B), (C), 44 (D), (E), and (F) are snapshots taken after 1, 20, 40, 60, and 100 simulated annealing cycles, respectively. ACS Paragon Plus Environment 45 46 47

Page 17 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47


17

(A)

(B)

(D)

(E)

(C)

(F)

Figure 3: PI3Kγ (PDB ID: 3lj3) binding site residues closely contacting the co-crystallized ligand and how they perturb in some selected protein frames during molecular dynamic simulation. Green lines represent hydrogen bonding. (A) Binding site residues in the original crystallographic structure of PI3K-γ. (B), (C), (D), (E), and (F) are snapshots taken after 2 , 4, 6, 8 and 10 ns of simulation, respectively. ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47

Page 18 of 39

18

Figure 4: Total energy calculated for conformations generated for PI3K-γ/ligand complex (PDB code: 3lj3) via classical molecular dynamics (100 sequential frames over 1.25 ns), simulated annealing molecular dynamics (energy readouts at the end of 100 cycles) and our deformation/relaxation protocol (Randomly selected 100 frames) (MD fluctuations are not clearly seen due to large scale of the y-axis). ACS Paragon Plus Environment

Page 19 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Clearly, the three molecular perturbation methods (i.e., D/R, MD and SAMD) unveiled the redundancy of several interactions apparent in the static complex. For example, Glu880, Tyr867, Met804 and Trp812 disappeared from the vicinity of the bound ligand upon exposure to all three simulation methods suggesting that any putative ligand binding interactions with these amino acids are essentially redundant. On the other hands, new high-frequency binding contacts appeared through the three simulation methods, in fact some of them are common among more than one method. For example, Val882-HN (CFs 97%, 96%, and 96% for D/R, SAMD and MD, respectively), Asp964-HN (CFs 74%, 76%, and 77% for D/R, SAMD and MD, respectively) and Asp841 (CFs 48%, 51%, and 71% for D/R, SAMD and MD, respectively). While Phe961-HE1 is found among the highest contacts of sa-LRCA and md-LRCA (CFs 45% and 48%, respectively) only and was downgraded to lesser frequency in the dr-LRCA case (CF 10%). Nevertheless, some high-ranking surviving contacts are unique for each of the three methods, e.g., LYS833, TRP812 and the combination of THR887/ASP964 uniquely emerged in sa-LRCA, md-LRCA and dr-LRCA, respectively (as in table 1). Top frequently ligand-encountered binding site atoms encoding from 4 to 6 binding interactions (of CF > 30%) were used to construct corresponding pharmacophore models for PI3K-γ, renin and JAK1, as described earlier in md-LRCA and sa-LRCA methods.22,23 (tables 1 and 2). This range of binding features is generally acceptable standard for pharmacophore models. In fact, binding models of lesser features can be promiscuous and those of greater number of features are too limited to serve as search queries towards new bioactive hits.1-23 We will discuss herein how we constructed pharmacophore model for PI3K-γ based on the new deformation/relaxation method (table 1): Emergence of the peptidic NH of Val882 as high contacting atom (CF = 97%) combined with its close proximity to the pyridine nitrogen of the ligand's aza-indole fragment suggested their



mutual hydrogen bonding and prompting us to add hydrogen bond acceptor (HBA) feature originating from the pyridine nitrogen towards the peptidic NH of Val882. Similarly, emergence of the amidic NH of Asp964 (CF = 74%) and carboxylic oxygen of Asp841 (48%) within close vicinity to one of the ligand's resorcinol hydroxyls suggested hydrogen bonding interactions involving these hydroxyls and Asp964 and Asp841 (table 1). We represented these putative interactions with hydrogen-bond donor (HBD) and hydrogen-bond acceptor features (HBA) in the corresponding pharmacophore model as in table 1. Likewise, the close proximity of the alpha-carbon of Asp964 (CF = 47%) to the ligand's resorcinol ring suggested mutual hydrophobic/van der Waals' interactions and prompted us to place hydrophobic feature onto the ligand's aromatic ring, as in table 1. Finally, The high-frequency emergence of the methyl of Thr887 (CF = 49%) and its close proximity to the methyl functionality of ligand's methyl-piperazine suggested mutual hydrophobic attraction that we represented in the corresponding pharmacophore model as hydrophobic feature. We constructed similar pharmacophore models for PI3K-γ based on SAMD (50 cycles) and MD (7 ns), and 8-feaured pharmacophore model based on the crystallographic static structure of the co-crystallized ligand (pharmacophoric features assumed based on distances ≤ 2.5 Å). Likewise, appropriate pharmacophore models were constructed for renin and JAK1 using dr-LRCA and were compared with their corresponding static congeners, as in table 2. Incidentally, the dr-LRCA process filtered the binding features from 13 in the renin static crystallographic complex to 6 only. On the other hand, although the number of pharmacophoreic features remained the same in the JAK1 case after dr-LRCA assessment (i.e., 4 features), one of the features switched from being HBA in the static crystallographic pharmacophore to become HBD upon treatment with dr-LRCA as in table 2. This change is because the deformation/relaxation process shifted the amidic


Page 20 of 39

Page 21 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


carbonyl of Glu966 from the vicinity of the ligand's hydroxyl group and replaced it with the amidic NH of Ser963 (lowest row in table 2). In the static structure the ligand's hydroxyl acts as hydrogen bond donor with the amidic carbonyl of Glu966, while after deformation/relaxation cycles, the same hydroxyl acted as hydrogen bond acceptor and received the hydrogen atom from the amidic NH.

3.3 Validation of Resulting Pharmacophores 3.3.1 Receiver Operating Characteristic (ROC) analysis To validate the overall dr-LRCA procedure, we tested the ability of the resulting pharmacophores to classify corresponding lists of inhibitors (extracted from ChEMBL database)48 into active and inactive via receiver operating characteristic (ROC) analysis (see section SM-2 under Supporting Information).3-17,22,23,26-31,43-47 We also evaluated the ROC performances of PI3K-γ pharmacophores (generated by D/R, SAMD and MD) against a list of mTOR (mammalian target of rapamycin) inhibitors (see section SM-2 under Supporting Information). mTOR and PI3K-γ are close homologues and their joint aberrant activation is commonly observed in several diseases such as cancer, diabetes and Alzheimer's disease.13,62-64 Table 3 shows ROC results of pharmacophores compared to their crystallography-based (static) analogues. Clearly, static pharmacophores corresponding to PI3K-γ and renin failed to capture any hits from the validation lists. This is unsurprising since the two static pharmacophores were too restrictive as they included too many features (≥ 8 features derived based on proximity between co-crystallized ligands and binding pockets atoms, ≤ 2.5 Ǻ). Nevertheless, the static binding pharmacophore derived from JAK1-ligand complex showed reasonable ROC performance probably because it exhibits only 4 binding features.



On the other hand, all three conformational perturbation methods (D/R, SAMD, and MD) refined close contacts and allowed reduction of the number of binding features in corresponding PI3K-γ pharmacophores, as in table 1, yielding superior binding models of better ROC results (Table 3). Interestingly, however, PI3K-γ pharmacophores resulting from the three methods (i.e., dr-LRCA, md-LRCA and sa-LRCA) exhibited measurable differences in their ROC performances which points to the fact that they cover different pharmacophoric spaces. The dr-LRCA approach produced similar ROC enhancement in renin case, as in table 2. However, in JAK1 case the improvement is less obvious albeit noticeable (ROC sensitivity increased from 47 to 92%). The fact that our deformation/relaxation algorithm maintained the basic features of the static pharmacophore in JAK1 case (which was of fair ROC performance), while it filtered out redundant features from static PI3K-γ and renin pharmacophores, underlines the robustness and resilience of dr-LRCA approach: Only redundant features were filtered without compromising essential contacts.. ROC results suggest that dr-LRCA-derived pharmacophores are comparable, and even superior, to those derived from other molecular perturbation/contacts analysis methods (sa-LRCA and md-LRCA), while on the hand, the new deformation/relation method is less time-consuming and computationally simpler compared to the molecular dynamics methods. For example, several days (10 to 20) were required to generate pharmacophore models using sa-LRCA (50 cycles) or md-LRCA (7 ns) from a single crystallographic complex using standard core-i7 computer with 8 processors, while the same procedure requires only 10-30 hours to build valid dr-LRCA-based pharmacophore (from 500 stochastic structures) for the same complex and using similar computational power, which should allow development of numerous pharmacophore models from several protein complexes within reasonable time frame.


Page 22 of 39

Page 23 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


It remains to be mentioned that pharmacophore models generated via dr-LRCA for PI3Kγ and renin were comparable with those generated by supervised ligand-based methods vis-à-vis types of binding features and their 3D arrangement.9,32

3.2.2. Experimental validation of anti-PI3K-γ hits To further validate our dr-LRCA procedure, we decided to mine the national cancer institutes (NCI) list of compounds for new potential anti PI3K-γ/mTOR inhibitors using dr-LRCA-, sa-LRCA- and md-LRCA-based pharmacophores. The highest-ranking 40 available hits (based on their fit values)22,23,29,43,62,65 captured by each pharmacophore model were acquired from the NCI and evaluated against PI3K-γ and mTOR in vitro. Hits that illustrated > 40% inhibition against both targets at 10 µM were further tested at different concentrations to determine their IC50 values. Table 4 show their inhibitory bioactivities against PI3K-γ and mTOR. Figure 5 shows the corresponding dose-response curves of the tested hits. Figures 6, 7 and 8 show the chemical structures of evaluated hits and how they map their corresponding capturing pharmacophores. Clearly from figure 5 and table 3, the dose-response curves of tested hits exhibit reasonable steepness (Hill slope values ≲ 1.5) and excellent correlation coefficients, which strongly suggest their authenticity (i.e., non-promiscuous inhibitors).65-67 The pankinase inhibitor staurosporine was used as standard control inhibitor. Interestingly, one of the active hits was jointly captured by dr-LRCA and md-LRCA pharmacophores against mTOR (e.g., of NCI code 683637), which further proof the robustness of our new methods, and its capability to capture similar hits of other methods.



Page 24 of 39

4. CONCLUSION We introduced a novel method for extracting valid pharmacophore(s) from single crystallographic structure(s) within reasonable time scale. The new method is based on ligandreceptor contacts analysis following energy relaxation of predetermined set of randomly deformed

complexes.

Ligand-receptor

contacts

maintained

across

many

de-

formed/relaxed structures are assumed to be critical and used to guide pharmacophore development. This methodology was used to develop valid pharmacophore models for PI3K-γ, RENIN and JAK1. The resulting pharmacophore models were validated by receiver operating characteristic (ROC) analysis. PI3K-γ pharmacophores were also used to search for new inhibitors from the national cancer institute list of compounds. The process identified new PI3K-γ/mTOR inhibitory leads of low micromolar IC50 values.


Page 25 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1. Top ranking ligand-contacting binding amino acids and corresponding ligand-encounter rates and pharmacophore models based on the three different methods (SAMD, MD, and deformation/relaxation) for the PI3K-γ case. The X,Y and Z coordinates and sphere-sizes are shown in Supporting Information Table S-1. Run

SAMD (50 Cycles)

MD (7 ns)

High ranking Percent contacting Contact atomsc Frequency Val882-HN

96

Asp964-HN

76

Asp841-OD2

51

Phe961-HE1

45

Lys833-HZ3

43

Val882-HN

96

Asp964-HN

77

Asp841-OD1

71

Trp812-HZ3

57

Annotated High ranking Contacting Atoms Within PI3K-γ Binding Pocketd

Phramacophorese

+

+

+ Phe961-HE1

48

Val882-HN

97

Asp964-HN

74

Thr887-HG21

49

Asp841-OD1

48

Asp964-HA

47

+ D/R (500 frames)

Static posea,b

Val882-HN, Val882-CG2, Glu880-O, Tyr867-OH, Phe-965-HN, Asp964-CA, Asp964-CG, Asp964-OD2, Lys833-HZ2, Met804-SD, Trp812-CZ3

a

Based on the crystallographic complex before minimization (PDB code: 3lj3). Contacts ≤ 2.5 Å were considered. The corresponding static pharmacophore included 8 features. cAtom names and number are as encoded in the protein databank, except for hydrogen atoms which were named and numbered according to Discovery Studio (version 2.5.5). dThe frame of highest WCS pose for the particular perturbation method (MD, SAMD or D/R). High frequency contacting atoms were annotated as spheres. eHydrogen bond acceptors are shown as green vectored spheres, Hydrogen bond donors as pink vectored spheres, hydrophobic features as blue spheres

b


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling Page 26 of 39 26 Table 2. Top ranking ligand-contacting binding site amino acids and corresponding ligand-encounter rates and pharmacophore models identified using the deformation/relaxation method for renin and JAK1. The X,Y and Z coordinates and sphere-sizes are shown in Supporting Information Table S-2. Run

Renin

JAK1kinase

a

Percent High ranking Contact c contacting atoms Frequency

Ala218-HB1

98

Gly217-HA1

98

Gly34-HA1

98

Thr295-HG21

95

Ala218-HA

89

Thr77-HG21

82

Asp215-OD2

73

Asp215-OD1

71

Ser963-HN

73

Glu957-O

63

Leu959-HN

58

Asp1021-HB1

44

Ser963-HB1

41

Asn1008-HA

37

Renin Statica,b

Val120-HG21, Asp32-OD1, Gly34-HA2, Tyr-HD1, Tyr-HB2, Gly217-O, Gly217-HA2, Ser76HN, Ser76-HB1, Phe117-HE2, Phe117-HZ, Phe117-HE1, Ala115-HA, Leu114-HB1, Thr77HG1, Thr77-OG1, Thr77-HG22, Thr77-HG21, Ala218-HA, Ala218-HB1, Ile291-HD12, Ala300-HB2, Ser222-OG, Ser222HG, Tyr220-O, Met289-HE2, Ser219-HN, Ser219-OG, His287HE1, Val120-HG21.

JAK1kinase statica,b

Leu881-HB2, Gly882-HA2, Glu966-OE1, Val889-HG11, Val889-HG21, Val889-HB, Asp1012-HB1, Leu1010HD21, Leu959-HN, Glu957-O

Annotated High ranking Contacting Atoms within Binding Pocketsd

Phramacophorese

Based on the crystallographic complex before minimization (PDB code: 1bil and 4ivb, respectively). Contacts ≤ 2.5 Å were considered. The corresponding static pharmacophore included 13 features for Renin and 4 features for JAK1-kinase. cAtom names and number are as encoded in the protein databank, except for hydrogen atoms which were named and numbered according to Discovery Studio (version ACS ParagonSM-1) Plus Environment 2.5.5). dThe frame of highest WCS (as described in supplementary pose for the D/R method. High frequency contacting atoms were annotated as spheres. eHydrogen bond acceptors are shown as green vectored spheres, Hydrogen bond donors as pink vectored spheres, hydrophobic features as blue spheres

b

Page 27 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 3: ROC analysis results for pharmacophores generated from PI3K-γ (mTOR), renin, JAK1 crystallographic structures using the different perturbation methods (see section Supplementary section SM-2 for details about testing sets selection and ROC analysis). Protein

PI3K-γ

Renin

JAK1

Method

Sensitivitya

Specificityb

ROC-AUCc

SAMD (50 cycle)

0.32 (0.34)d

0.87 (0.83)d

0.89 (0.86)d

MD (7 ns)

0.20 (0.24)d

0.90 (0.91)d

0.91 (0.91)d

Deformation/Relaxation (500 frames)

0.06 (0.22)d

0.96 (0.94)d

0.97 (0.95)d

Static

0.00

0.00

--- e


0.10

0.89

0.90

Static

0.00

0.00

--- e


0.92

0.41

0.74

Static

0.47

0.64

0.74

a

True positive rate. True negative rate. c Area under the ROC curve. d Numbers between brackets represent ROC results against mTOR list of inhibitors ligands. e Pharamacophores failed to capture hits b



Page 28 of 39

Table 4: Active hits captured by the 5-featured pharmacophores generated after 50 SAMD cycles, 7ns MD, and DR of mTOR and PI3K-γ Experimental Method

a

Enzyme

mTOR SAMD PI3K-γ MD

mTOR PI3K- γ

mTOR DR

PI3K- γ

a

NCI code

292213 292216 292217 292217 658264 683637 723785 75328 75327 683637 734052 75327 75328

Percent inhibition of at 10 µM 88 47 60 66 46 60 79 92 57 56 54 89 74

IC50 (µM)b

Hill Slopec

2.95 10.51 16.39 2.44 15.89 7.67 2.45 2.00 7.85 9.80 8.85 2.02 2.81

1.60 1.60 0.50 1.00 0.60 1.10 1.00 1.50 1.30 0.87 1.50 1.10 0.70

a

Correlation coefficientd 0.993 0.996 0.839 0.997 0.986 0.992 0.997 0.994 0.999 0.977 0.995 0.885 0.929

Chemical structures shown in Figures 6-8. IC50 values experimentally determined as duplicates for most-active hits determined over concentrations ranging from 0.01 to 100 µM (as in figure 5). c Calculated based on the dose-response curves by GraphPad Prism 5.0. d Correlation coefficient of the corresponding dose-response curves. b


Journal of Chemical Information and Modeling 29 100

Percent inhibition

(B)

80 60 40 20

100 80 60 40 20 0

0 1

2

3

4

1

5

2

60 40 20

Percent inhibition

Percent inhibition

(D)

80

2

3

4

20

1

5

2

(F)

80 60 40 20

100

Percent inhibition

Percent inhibition

100

80 60 40 20 0

1

2

3

4

1

5

2

100 80 60 40 20 0 2

3

5

4

4

5

4

5

4

5

80 60 40 20 0

5

1

2

Log(Conc.) [nM]

3

Log(Conc.) [nM]

(J)

80 60 40 20

100

Percent inhibition

100

Percent inhibition

4

100

Percent inhibition

Percent inhibition

(H)

1

3

Log(Conc.) [nM]

Log(Conc.) [nM]

80 60 40 20 0

0 1

2

3

4

1

5

2

(L)

80 60 40 20 0 1

2

3

4

5

Log(Conc.) [nM]

100

Percent inhibition

100

3

Log(Conc.) [nM]

Log(Conc.) [nM]

Percent inhibition

3

Log(Conc.) [nM]

0

(K)

5

40

Log(Conc.) [nM]

(I)

4

60

0

1

(G)

5

80

0

(E)

4

100

100

(C)

3

Log(Conc.) [nM]

Log(Conc.) [nM]

80 60 40 20 0 1

2

3

Log(Conc.) [nM]

100

(M)

Percent inhibition

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(A)

Percent inhibition

Page 29 of 39

80 60 40 20 0 1

2

3

4

5

Log(Conc.) [nM]

Figure 5: Dose-response curves of different hits (NCI codes) against PI3K-γ or mTOR. (A) 292213 (against mTOR), (B) 292216 (against mTOR), (C) 292217 (against mTOR), (D) 292217 (against PI3K-γ). (E) 658264 (against PI3K-γ), (F) 723785 (against PI3K-γ), (G) 683637 (against mTOR), (H) 75328 (against mTOR), (I) 75327 (against mTOR), (J) 683637 (against mTOR), (K) 734052 (against ACS Paragon Environment mTOR), (L) 75327 (against PI3K-γ), and (M) 75328Plus (against PI3K-γ).


(A)

Page 30 of 39

H2N

(B)

N N OH O

S

O

NH2

(C)

(D) NH2 N N

OH

O

S

O

NH2

(E)

(F) Cl

OH

NH N

OH

N

(G)

(H) Cl

N NH O

O

S O

O

S NH2

Cl

Figure 6: Active hits captured by dr-LRCA pharmacophore model developed for PI3K-γ. (A), (C), (E) and (G) show how NCI hits (codes: 75328, 75327, 683637, and 734052, respectively) fit the corresponding capturing pharmacophore. (B), (D), (F) and (H) show the corresponding chemical structures of captured active hits. ACS Paragon Plus Environment

Page 31 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(A)

(B) OH

O

O

O HO

OH

O

(C)

OH

(D) O

O

O

HO HO

OH

O

(E)

O

(F)

OH

O

CF3

O HN

Cl

S

N NH

O

(G)

(H) OH HO

O O

O

HO OH O

Figure 7: Active hits captured by sa-LRCA pharmacophore model developed for PI3K-γ. (A), (C), (E) and (G) show how NCI hits (codes: 292216, 292217, 658264, and 292213, respectively) map the capturing pharmacophore. (B), (D), (F) and (H) show the corresponding chemical structures of captured active hits.



(A)

Page 32 of 39

(B) Cl OH

H N

OH N

N

(C)

(D)

NH2 N

N

N

N S

N

S

Figure 8: Active hits captured by md-LRCA pharmacophore model developed for PI3K-γ. (A) and (C) show how NCI hits (codes: 683637 and 723785, respectively) map the capturing pharmacophore. (B) and (D) show the corresponding chemical structures of captured active hits.


Page 33 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ACKNOWLEDGMENTS The authors would like to thank The Hashemite University, King Abdullah II Fund for Development (KAFD) and IMAN1 Center – Jordan’s National Supercomputing Center – for providing us with the fund, computational tools and high-performance computing system time to perform our experiments, the authors are deeply indebted to Eng. Zaid Abudayyeh from IMAN1 Center for his help in performing this project. The authors also thank the Deanship of Scientific Research the University of Jordan for generous support.



REFERENCES (1) Kumar, B.V.; Kotla, R.; Buddiga, R.; Roy, J.; Singh, S.S.; Gundla, R.; Ravikumar, M.; Sarma, J.A. Ligand-Based And Structure-Based Approaches In Identifying Ideal Pharmacophore Against C-Jun N-Terminal Kinase-3. J. Mol. Model. 2010, 17, 151163. (2) Du, X.; Li, Y.; Xia, Y.; Ai, S.; Liang, J.; Sang, P.; Ji, X.; Liu, S. Insights Into Protein– Ligand Interactions: Mechanisms, Models, And Methods. Int. J. Mol. Sci. 2016, 17, 144-177. (3) Mortier, J.; Rakers, C.; Bermudez, M.; Murgueitio, M.; Riniker, S.; Wolber, G. The Impact Of Molecular Dynamics On Drug Design: Applications For The Characterization Of Ligand–Macromolecule Complexes. Drug. Discov. Today 2015, 20, 686-702. (4) Sanders, M.; McGuire, R.; Roumen, L.; de Esch, I.; de Vlieg, J.; Klomp, J.; de Graaf, C. From The Protein's Perspective: The Benefits And Challenges Of Protein Structure-Based Pharmacophore Modeling. Med. Chem. Commun. 2012, 3, 28-38. (5) Jacoby, E. Computational Chemogenomics. Wiley Interdiscip Rev. Comput. Mol. Sci. 2011, 1, 57-67. (6) Ermondi, G.; Caron, G. Recognition Forces In Ligand–Protein Complexes: Blending Information From Different Sources. Biochem Pharmacol. 2006, 72, 1633-1645. (7) Kurogi, Y.; Guner, O. Pharmacophore Modeling And Three-Dimensional Database Searching For Drug Design Using Catalyst. Curr. Med. Chem. 2001, 8, 1035-1055. (8) Abuhamdah, S.; Habash, M.; Taha, M.O. Elaborate Ligand-Based Modeling Coupled With QSAR Analysis And In Silico Screening Reveal New Potent Acetylcholinesterase Inhibitors J. Comput. Aided. Mol. Des. 2013, 27, 1075-1092. (9) Al-Nadaf, A.H.; Taha, M. Discovery Of New Renin Inhibitory Leads Via Sequential Pharmacophore Modeling, QSAR Analysis, In Silico Screening And In Vitro Evaluation. J. Mol. Graph. Model. 2011, 29, 843-864. (10) Al-Sha’er, M.A.; VanPatten, S.; Al-Abed, Y.; Taha, M.O. Elaborate Ligand-Based Modeling Reveal New Migration Inhibitory Factor Inhibitors. J. Mol. Graph. Model. 2013, 42, 104-114. (11) Al-Sha’er, M.A.; Khanfar, M.A; Taha, M.O. Discovery Of Novel Urokinase Plasminogen Activator (Upa) Inhibitors Using Ligand-Based Modeling And Virtual Screening Followed By In Vitro Analysis. J. Mol. Model. 2014, 20, 2080-2095. (12) Habash, M.A.; Abdelazeem, A.H.; Taha, M.O. Elaborate Ligand-Based Modeling Reveals New Human Neutrophil Elastase Inhibitors. Med. Chem. Res. 2014, 23, 3876-3896. (13) Khanfar, M.A; AbuKhader, M.M; Alqtaishat, S.; Taha, M.O. Pharmacophore Modeling, Homology Modeling, And In Silico Screening Reveal Mammalian Target Of Rapamycin Inhibitory Activities For Sotalol, Glyburide, Metipranolol, Sulfamethizole, Glipizide, And Pioglitazone. J. Mol. Graph. Model. 2013, 42, 39-49. (14) Shahin, R.; Taha, M.O. Elaborate Ligand-Based Modeling And Subsequent Synthetic Exploration Unveil New Nanomolar Ca2+/Calmodulin-Dependent Protein Kinase II Inhibitory Leads. Bioorg. Med. Chem. 2012, 20, 377-400.


Page 34 of 39

Page 35 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(15) Shahin, R.; AlQtaishat, S.; Taha, M.O. Elaborate Ligand-Based Modeling Reveal New Submicromolar Rho Kinase Inhibitors. J. Comput. Aided. Mol. Des. 2011, 26, 249- 266. (16) Taha, M.O; Qandil, A.M; Al-Haraznah, T.; Abu-Khalaf, R.; Zalloum, H.; Al-Bakri, A.G. Discovery Of New Antifungal Leads Via Pharmacophore Modeling And QSAR Analysis Of Fungal N-Myristoyl Transferase Inhibitors Followed By In Silico Screening. Chem. Biol. Drug Des. 2011, 78, 391-407. (17) Taha, M.O.; Habash, M.; Hatmal, M.M; Abdelazeem, A.H; Qandil, A. Ligand-Based Modeling Followed By In Vitro Bioassay Yielded New Potent Glucokinase Activators. J. Mol. Graph. Model. 2015, 56, 91-102. (18) Morris, G.M.; Olson, A.J.; Goodsell, D.S. Protein-Ligand Docking Methods. Princ. Med. Chem. 2000, 8, 31-48. (19) Kontoyianni, M.; McClellan, L.M.; Sokol, G.S. Evaluation Of Docking Performance: Comparative Data On Docking Algorithms. J. Med. Chem. 2004, 47, 558565. (20) Beier, C.; Zacharias, M. Tackling The Challenges Posed By Target Flexibility In Drug Design. Expert. Opin. Drug. Discov. 2010, 5, 347-359. (21) Langer, T.; Hoffmann, R.D. Pharmacophore Modelling: Applications In Drug Discovery. Expert. Opin. Drug. Discov. 2006, 1, 261-267. (22) Hatmal, M. M.; Taha MO. Simulated Annealing Molecular Dynamics And LigandReceptor Contacts Analysis For Pharmacophore Modeling. Future Med. Chem. 2017, 9, 1141-1159. (23) Hatmal, M. M.; Jaber, S.; Taha M. O. Combining Molecular Dynamics Simulation And Ligand-Receptor Contacts Analysis As A New Approach For Pharmacophore Modeling: Beta-Secretase 1 And Check Point Kinase 1 As Case Studies. J. Comput. Aided. Mol. Des. 2016, 30, 1149-1163. (24) Ortuso, F.; Langer, T.; Alcaro, S. GBPM: GRID-Based Pharmacophore Model: Concept And Application Studies To Protein-Protein Recognition. Bioinformatics 2006, 22, 1449-1455. (25) Alcaro, S.; Artese, A.; Ceccherini-Silberstein, F.; Chiarella, V.; Dimonte, S.; Ortuso, F.; Perno, C. Computational Analysis Of Human Immunodeficiency Virus (HIV) Type-1 Reverse Transcriptase Crystallographic Models Based On Significant Conserved Residues Found In Highly Active Antiretroviral Therapy (HAART)-Treated Patients (Supplementary Material). Curr. Med. Chem. 2010, 17, 290-308. (26) Al-Sha'er, M. A.; Taha, M.O. Application Of Docking-Based Comparative Intermolecular Contacts Analysis For Validating Hsp90α Docking Studies And Subsequent In Silico Screening For Inhibitors. J. Mol. Model. 2012, 18, 4843-4863. (27) Taha, M.O.; Habash, M.; Al-Hadidi, Z.; Al-Bakri, A.; Younis, K.; Sisan, S. DockingBased Comparative Intermolecular Contacts Analysis As New 3D QSAR Concept For Validating Docking Studies And In Silico Screening: NMT and GP Inhibitors as Case Studies. J. Chem. Inf. Model. 2011, 51, 647–669. (28) Taha, M.O.; Habash, M.; Khanfar, M. The Use Of Docking-Based Comparative Intermolecular Contacts Analysis To Identify Optimal Docking Conditions Within



Glucokinase And To Discover Of New GK Activator. J. Comput. Aided. Mol. Des. 2014, 28, 509-547. (29) Abuhammad, A.; Taha, M. Innovative Computer-Aided Methods For The Discovery Of New Kinase Ligands. Future Med. Chem. 2016, 8, 509-526. (30) Jaradat, N. J., Khanfar, M. A.; Habash, M.; Taha, M. O. Combining Docking-Based Comparative Intermolecular Contacts Analysis And K-Nearest Neighbor Correlation For The Discovery Of New Check Point Kinase 1 Inhibitors. J. Comput. Aided. Mol. Des. 2015, 29, 561-581. (31) Habash, M.; Abuhamdah, S.; Younis, K.; Taha, M. O. Docking-Based Comparative Intermolecular Contacts Analysis And In Silico Screening Reveal New Potent Acetylcholinesterase Inhibitors. Med. Chem. Res. 2017, 26, 2768–2784. (32) Taha M. O.; Al-Sha'er M. A.; Khanfar, M. A.; Al-Nadaf, A. H. Discovery Of Nanomolar Phosphoinositide 3 Kinase Gamma (PI3K-γ) Inhibitors Using LigandBasedmodeling And Virtual Screening Followed By In Vitro Analysis. Eur J. Med. Chem. 2014, 84, 454-465. (33) Gyori, D.; Chessa, T.; Hawkins, P. T.; Stephens, L. R. Class (I) Phosphoinositide 3Kinases In The Tumor Microenvironment. Cancers (Basel). 2017, 9, 24. (34) Sadeghi, N.; Gerber, D. Targeting The PI3K Pathway For Cancer Therapy. Future Med. Chem. 2012, 4, 1153-1169. (35) Fruman, D.; Rommel, C. PI3K And Cancer: Lessons, Challenges And Opportunities. Nat. Rev. Drug Discov. 2014, 13, 140-156. (36) Lim, H.; Crowe, P.; Yang, J. Current Clinical Regulation Of PI3K/PTEN/Akt/mTOR Signaling In Treatment Of Human Cancer. J. Cancer Res. Clin. Oncol. 2014, 141, 671-689. (37) Milosevic, Z.; Pesic, M.; Stankovic, T.; Dinic, J.; Milovanovic, Z.; Stojsic, J.; Dzodic, R.; Tanic, N.; Bankovic, J. Targeting RAS-MAPK-ERK And PI3K-AKTmTOR Signal Transduction Pathways To Chemosensitize Anaplastic Thyroid Carcinoma. Transl. Res. 2014, 164, 411-423. (38) Taylor, A.; Pool, J. Clinical Role Of Direct Renin Inhibition In Hypertension. Am. J. Ther. 2012, 19, 204-210. (39) Chen, E.; Staudt, L.; Green, A. Janus Kinase Deregulation In Leukemia And Lymphoma. Immunity 2012, 36, 529-541. (40) Vannucchi, A.; Kiladjian, J.; Griesshammer, M.; Masszi, T.; Durrant, S.; Passamonti, F.; Harrison, C.; Pane, F.; Zachee, P.; Mesa, R.; He, S.; Jones, M.M.; Garrett, W.; Li, J.; Pirron, U.; Habr, D.; Verstovsek, S. Ruxolitinib Versus Standard Therapy For The Treatment Of Polycythemia Vera. N. Engl. J. Med. 2015, 372, 426-435. (41) Kiladjian, J.; Winton, E.; Talpaz, M.; Verstovsek, S. Ruxolitinib For The Treatment Of Patients With Polycythemia Vera. Expert. Rev. Hematol. 2015, 8, 391-401. (42) McKeage, K. Ruxolitinib: A Review In Polycythaemia Vera. Drugs 2015, 75, 17731781. (43) Abuhammad, A.; Taha, M. QSAR Studies In The Discovery Of Novel Type-II Diabetic Therapies. Expert. Opin. Drug Discov. 2015, 11, 197-214.


Page 36 of 39

Page 37 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(44) Al-Aqtash, R.; Zihlif, M.; Hammad, H.; Nassar, Z.; Meliti, J.; Taha, M. LigandBased Computational Modelling Of Platelet-Derived Growth Factor Beta Receptor Leading To New Angiogenesis Inhibitory Leads. Comput. Biol. Chem. 2017, 71, 170179. (45) Triballeau, N.; Acher, F.; Brabet, I.; Pin, J. P.; Bertrand, H. O. Virtual Screening Workflow Development Guided By The “Receiver Operating Characteristic” Curve Approach. Application To High-Throughput Docking On Metabotropic Glutamate Receptor Subtype 4. J. Med. Chem. 2005, 48, 2534-2547. (46) Kirchmair, J.; Distinto, S.; Schuster, D.; Spitzer, G.; Langer, T.; Wolber, G. Enhancing Drug Discovery Through In Silico Screening: Strategies To Increase True Positives Retrieval Rates. CMC 2008, 15, 2040-2053. (47) Kirchmair J, Markt P, Distinto S, Wolber G, Langer T. Evaluation Of The Performance Of 3D Virtual Screening Protocols: RMSD Comparisons, Enrichment Assessments, And Decoy Selection—What Can We Learn From Earlier Mistakes?. J. Comput. Aided Mol. Des. 2008, 22, 213-228. (48) Gaulton, A.; Bellis, L. J.; Bento, A. P.; Chambers, J.; Davies, M.; Hersey, A.; Light, Y.; McGlinchey, S.; Michalovich, D.; Al-Lazikani, B. Chembl: A Large-Scale Bioactivity Database For Drug Discovery. Nucleic Acids Res. 2011, 40, 1100-1107. (49) Jones, G.; Willett, P.; Glen, R.; Leach, A.; Taylor, R. Development and Validation Of A Genetic Algorithm For Flexible Docking 1 1Edited By F. E. Cohen. J. Mol. Biol. Mol. 1997, 267, 727-748. (50) Leach, A. Ligand Docking To Proteins With Discrete Side-Chain Flexibility. J. Mol. Biol. 1994, 235, 345-356. (51) Verdonk, M.; Cole, J.; Hartshorn, M.; Murray, C.; Taylor, R. Improved ProteinLigand Docking Using GOLD. Proteins 2003, 52, 609-623. (52) Bianco, G.; Forli, S.; Goodsell, D.; Olson, A. Covalent Docking Using Autodock: Two-Point Attractor and Flexible Side Chain Methods. Protein Sci. 2015, 25, 295301. (53) Tirado-Rives, J.; Jorgensen, W. Contribution Of Conformer Focusing To The Uncertainty In Predicting Free Energies For Protein−Ligand Binding. J. Med. Chem. 2006, 49, 5880-5884. (54) Meng, X.; Zhang, H.; Mezei, M.; Cui, M. Molecular Docking: A Powerful Approach for Structure-Based Drug Discovery. Curr. Comput. Aided Drug Des. 2011, 7, 146157. (55) Adcock, S.; McCammon, J. Molecular Dynamics: Survey of Methods For Simulating The Activity of Proteins. Chem. Rev. 2006, 106, 1589-1615. (56) Kuznetsov, I. Ordered Conformational Change In The Protein Backbone: Prediction Of Conformationally Variable Positions From Sequence And Low-Resolution Structural Data. Proteins 2008, 72, 74-87. (57) Kuroda, D.; Gray, J. Pushing The Backbone In Protein-Protein Docking. Structure 2016, 24, 1821-1829. (58) Chaudhury, S.; Gray, J. Conformer Selection And Induced Fit In Flexible Backbone Protein–Protein Docking Using Computational And NMR Ensembles. J. Mol. Biol. 2008, 381, 1068-1087. ACS Paragon Plus Environment


(59) Mashiach, E.; Nussinov, R.; Wolfson, H. Fiberdock: Flexible Induced-Fit Backbone Refinement In Molecular Docking. Proteins 2010, 78, 1503-1519. (60) Bonvin, A. Flexible Protein–Protein Docking. Curr Opin Struct Biol. 2006, 16, 194200. (61) Paquet, E.; Viktor, H. Molecular Dynamics, Monte Carlo Simulations, And Langevin Dynamics: A Computational Review. Biomed. Res. Int. 2015, 2015, 1-18. (62) Zoncu, R.; Efeyan, A.; Sabatini, D. mTOR: From Growth Signal Integration To Cancer, Diabetes And Ageing. Nat. Rev. Mol. Cell Biol. 2010, 12, 21-35. (63) Rosner, M.; Hanneder, M.; Siegel, N.; Valli, A.; Fuchs, C.; Hengstschlager, M. The mTOR Pathway And Its Role In Human Genetic Diseases. Mutat. Res. Rev. Mutat. Res. 2008, 659, 284-292. (64) Xu, K.; Liu, P.; Wei, W. mTOR Signaling In Tumorigenesis. Biochim. Biophys. Acta, Rev. Cancer 2014, 1846, 638-654. (65) Taha, M. O.; Dahabiyeh, L. A.; Bustanji, Y.; Zalloum, H.; Saleh, S. Combining Ligand-Based Pharmacophore Modeling, Quantitative Structure−Activity Relationship Analysis And In Silico Screening For The Discovery Of New Potent Hormone Sensitive Lipase Inhibitors. J. Med. Chem. 2008, 51, 6478-6494. (66) Walters, W.; Namchuk, M. A Guide To Drug Discovery: Designing Screens: How To Make Your Hits A Hit. Nat. Rev. Drug. Discov. 2003, 2, 259-266. (67) Shoichet, B. K. Interpreting Steep Dose-Response Curves In Early Inhibitor Discovery. J. Med. Chem. 2006, 49, 7274-7277.


Page 38 of 39

Page 39 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table of Contents Graphic 338x190mm (72 x 72 DPI)


Relaxation and Intermolecular

Recommend Documents