A Computational Model for Overcoming Drug ... - ACS Publications

Article pubs.acs.org/molecularpharmaceutics

A Computational Model for Overcoming Drug Resistance Using Selective Dual-Inhibitors for Aurora Kinase A and Its T217D Variant Khaled H. Barakat,†,⊥ J. Torin Huzil,‡ Kirk E. Jordan,∥ Constantinos Evangelinos,∥ Michael Houghton,† and Jack Tuszynski*,†,‡,§ †

Li Ka Shing Applied Virology Institute, Department of Medical Microbiology and Immunology, ‡Department of Oncology, and Department of Physics, University of Alberta, Edmonton, Alberta P6G 2M7, Canada ∥ IBM Thomas J. Watson Research Center 1101 Kitchawan Rd, Yorktown Heights, New York 10598, United States ⊥ Department of Engineering, Mathematics, and Physics, Fayoum University, Fayoum, Egypt §

S Supporting Information *

ABSTRACT: The human Aurora kinase-A (AK-A) is an essential mitotic regulator that is frequently overexpressed in several cancers. The recent development of several novel AK-A inhibitors has been driven by the well-established association of this target with cancer development and progression. However, resistance and cross-reactivity with similar kinases demands an improvement in our understanding of key molecular interactions between the Aurora kinase-A substrate binding pocket and potential inhibitors. Here, we describe the implementation of state-of-the-art virtual screening techniques to discover a novel set of Aurora kinase-A ligands that are predicted to strongly bind not only to the wild type protein, but also to the T217D mutation that exhibits resistance to existing inhibitors. Furthermore, a subset of these computationally screened ligands was shown to be more selective toward the mutant variant over the wild type protein. The description of these selective subsets of ligands provides a unique pharmacological tool for the design of new drug regimens aimed at overcoming both kinase cross-reactivity and drug resistance associated with the Aurora kinase-A T217D mutation. KEYWORDS: Aurora kinase-A, mutation, T217D, virtual screening, clustering, dual-inhibitor

■

INTRODUCTION The human Aurora kinase-A (AK-A) is an important mitotic regulator belonging to a small family of Aurora kinases (AKs). The Aurora kinase family is composed of three members commonly known as A, B, and C, each of which localizes at the centrosome during mitosis. Structurally, each AK consists of a highly conserved central catalytic kinase domain and an activation loop, which are flanked by variable N- and Cterminal domains that contain vital regulatory sequences. Since AKs function as key regulators throughout mitosis, they often become overexpressed in human metastases, including breast, urological, skin, and gastrointestinal (GI) tract cancers.1−3 The classification of AK-A as an oncogene has often been proposed, as it is involved in the regulation of centrosome maturation, entry into mitosis, and the formation of the bipolar mitotic spindle. In early investigations, the overexpression of AK-A mRNA was demonstrated as playing an essential role in the transformation of breast tumor cells, leading to the identification of the 20q13 amplicon as a poor prognostic indicator in patients diagnosed with breast cancer.4,5 While the importance of AK-A in the development of breast cancer is clear, its overexpression has also been established as a transformative factor in several other human tumor cell lines. © 2013 American Chemical Society

The observation that aberrant expression of AK-A is correlated with reducing the activity of tumor suppressor genes such as p536 has led to its selection as a promising target for the development of inhibitors that will function as novel chemotherapeutic agents. The development and testing of several novel anticancer dugs that target AKs have yielded promising results in preclinical and early phase I and II clinical trials.7−9 However, while these AK inhibitors (AKI) do offer a new direction in cancer treatment, detailed structural information to understand the mechanisms behind drug sensitivity and resistance is necessary. Structurally, the plasticity of the AK nucleotidebinding pocket has several functional consequences on inhibitor development. This plasticity in the central kinase domain has been implicated in the nonspecific interaction of several AKIs toward unrelated kinases.10−12 Previous results have demonstrated that resistance to chemotherapy treatment is correlated with the overexpression of AK-A through Received: Revised: Accepted: Published: 4572

July 7, 2013 October 1, 2013 October 4, 2013 October 4, 2013 dx.doi.org/10.1021/mp4003893 | Mol. Pharmaceutics 2013, 10, 4572−4589

Molecular Pharmaceutics

■

disruption of the spindle checkpoint normally activated by paclitaxel treatment.13 Cellular resistance to AKI exposure has resulted in the identification of several somatic mutations within the AK-A catalytic (nucleotide binding) domain. These include three single site missense mutations (V174M, S155R, and T217D) and a single nonsense mutation (S361*), which results in the expression of a nonfunctional C-terminally truncated enzyme.14,15 Expression of these mutations results in either constitutively increased kinase activity (V174M), decreased activity due to misregulation (S155R), or a reduction in nucleotide binding affinity (T217D). Of particular interest for the development of inhibitors toward nucleotide binding is the partially drug-resistant mutant, T217D, which has been shown to result in drug resistance, even in the presence of endogenously expressed AK-A.15 The elucidation of the molecular interactions for clinically used AK-A inhibitors is essential to understand the mechanisms behind their specificity. Crystallographic studies have previously identified several binding modes for AKIs and revealed key hydrogen bonds that are formed with Glu-211 and Ala-213 in the kinase domain, while other regions of the AKI have demonstrated binding with the active site surface through noncovalent interactions.16 Computational efforts to identify AKI molecules with superior binding affinities have previously used techniques including three-dimensional common feature pharmacophores to query existing databases for similar molecules.17 While these studies have been successful at classifying active and inactive compounds, selected “active” compounds tended to demonstrate only low micromolar binding affinities.18 Furthermore, molecular dynamics (MD) investigations of AKI binding have attempted to explain inhibitor selectivity toward AK-A over other members of the AK family, suggesting interactions with T217 may be crucial for drug interactions.19 Here, we describe the use of an improved version of the relaxed complex scheme (RCS) to introduce receptor flexibility during computational screening for AK-A inhibitors.20 The methodology described herein has been previously applied to several targets and successfully revealed novel and potent inhibitors for their activities.21−25 Our approach involved three major steps; first, we screened the National Cancer Institute diversity set26 and the panel of DrugBank’s small pharmacological molecules library27 for novel wild type AK-A inhibitors. This intermediate compound database was screened against an ensemble of targets consisting of the primary structural conformation of the nucleotide binding site as it was obtained from the AK-A crystal structure (PDB code: 1MQ4). In addition to screening of the crystal structure, we also examined the five dominant conformations that were obtained from clustering the MD trajectory of the 1MQ4 structure.28 From this set of compounds, the top 300 hits that were obtained from the wild type screening experiments were then redocked to 85 conformations of the mutated structures and ranked for their ability to bind to the T217D variant. The detailed structural information presented in this study is expected to provide the tools to design new AK-A inhibitors, optimize currently available structures for more potency and specificity, and overcome drug resistance mechanisms to qualify them for further preclinical and clinical development.

Article

EXPERIMENTAL SECTION

Preparation of the Initial Aurora A Kinase Structure. The structure of a human AK-A kinase domain, residues 126− 388, was obtained using the PDB identification code 1MQ4.28 Prior to structural refinement, all waters and nonconjugated heteroatoms were removed using the PyMOL Molecular Graphics System, Version 1.5.0.4 Schrödinger, LLC (The PyMOL Molecular Graphics System, Version 1.5.0.4 Schrödinger, LLC). The 1MQ4 structure is missing two amino acids, K258 and R286 that are present on surface loops, which were modeled using MODELER, version 9.12, implementing the model-loop procedure allowing a single adjacent residue on either side of the gap to move.29 A total of 100 short-loop models were ranked using discrete optimized protein energy (DOPE) statistical potentials to obtain an intermediate, lowest energy AKA kinase domain model. Prior to minimization and equilibration, missing hydrogen atoms and incomplete side chains were repaired using the UCSF Chimera packages Dock Prep module.30 Incomplete side chains were modeled using the Dunbrack rotamer library.31 Preparation of the T217D Mutation. To generate structures for the AK-A mutant T217D variant, all 85 representative cluster conformations were structurally aligned to the minimized crystal structure followed by mutational analysis using the Swiss-PdbViewer software (http://spdbv. vital-it.ch/). Rotamer libraries were used to avoid any steric clashes of the newly generated amino acid with its environment and to maximize its hydrogen bonding with the surrounding amino acids. Molecular Dynamics Simulations. The MD simulations were carried out as previously described.22,23 In brief, we performed an initial 110 ns MD simulation on the wild type Apo-enzyme (residues 126−388) at temperature of 310 K after adjusting the protonation states of the different residues at physiological pH (pH 7) conditions using the PDB2PQR software.32 The simulation was performed using the NAMD33 software and the AMBER99SB force field34 in explicit solvent and ions (22 chloride and 17 sodium) environment with a cube of water molecules of 20 Å. The system was minimized, heated with heavy restraints, which were removed gradually in a 100 ps equilibration phase followed by a production run. We then performed an additional 600 MD simulations (2 ns each) for the top 300 hits docked to the wild type and the mutant structures. The same conditions were used in carrying out the MD simulations on the protein−ligand structures after obtaining the parameters for each ligand using the Antechamber module in the AMBER 10 package.35 The root-meansquare deviation (RMSD) and atomic fluctuation (B-factors) for the apoprotein system was calculated for the full 110 ns MD simulation. For each ligand−protein system, we stored snapshots of the trajectories every 2 ps for subsequent MMPBSA binding energy analysis. Extracting Representative MD Structures. To extract dominant conformations from the MD simulation of the AK-A enzyme, we followed the same clustering procedure described in our earlier studies.22,23 Clustering analysis was performed on the 14 residues that constitute the AK-A substrate binding site, including side chains and hydrogen atoms for each residue. These residues comprised L139, G140, K141, G142, K143, G145, N146, L215, G216, T217, E260, N261, A273, and D274. The centroid of each cluster was chosen as the cluster representative structure. We used the average-linkage algorithm 4573

dx.doi.org/10.1021/mp4003893 | Mol. Pharmaceutics 2013, 10, 4572−4589


Article

Automated Clustering of Docked Poses. Here, we implemented the same automated clustering algorithm we used in our previous work.22,23 The objective was not to rely solely on a fixed parameter while clustering enormous numbers of docked poses for different ligands. An example of such parameters that we no longer use in our clustering calculations is the RMSD clustering with a fixed cutoff as implemented in AUTODOCK.45 Conversely, we use an adaptive method to estimate the optimal number of clusters needed to extract almost all the information contained within the different conformations predicted by AUTODOCK (for more details refer to our earlier works22,23). The metric we used in our calculation is the percentage of variance (SSR/SST),37 which is expected to be constant when all useful information is extracted from the examined data. This concept is termed the elbow criterion,37 and we successfully used it to predict the most probable conformation and cluster count in previous work. For each ligand, we use the SOM clustering algorithm and calculate the SSR/SST value at every trail to increase the number of clusters while computing the slope of the SSR/SST versus the number of clusters (i.e., by estimating the derivatives). A cutoff for the number of clusters is reached once the slope reaches almost zero. In this case, the clustering methodology is independent from the structure of the ligand, and we expect different cluster counts for different ligands. Rescoring of Top Hits Using MM-PBSA. We used the molecular mechanics Poisson−Boltzmann surface area (MM− PBSA) method to rescore the top hits predicted by docking. The method is implemented in the AMBER molecular dynamics package, and the free energy of binding is estimated as the summation of the molecular mechanical energies calculated in the gas phase (EMM), the solvation energy as predicted by the Poisson−Boltzmann implicit solvent method and the solute entropy as predicted by normal-mode analysis:

to group similar conformations in the 110 ns trajectory into clusters. The optimal number of clusters was estimated by observing the values of the Davies−Bouldin index (DBI)36 and the percentage of data explained by the data (SSR/SST)37 for different cluster counts ranging from 5 to 500. At the optimal number of clusters, a plateau in the SSR/SST is expected to match a local minimum in the DBI.37 Principal Component Analysis. As described in our earlier studies,22,23 we employed PCA to reduce the system dimensionality toward a finite set of independent variables covering the essential dynamics of the system.38 Similar to the clustering analysis described above, the PCA calculations were performed only on residues that encompass the binding site. For that, we first removed the overall rotations and translations for the different snapshots of the trajectory by fitting the entire MD trajectory to the initial structure. Then, we calculated the covariance matrix for the Cartesian coordinates of the selected set of atoms using: σij = ⟨(ri − ⟨ri⟩)(rj − ⟨rj⟩)⟩

(1)

Using the same analysis, one can predict the completeness of sampling for the MD trajectory. Similar to our earlier work, we used the normalized overlap parameter proposed by Hess39 to estimate the sampling quality: normalized overlap(C1, C2) = 1 −

tr(( C1 −

C2 )2 )

tr(C1) + tr(C2) (2)

where C1 and C2 are the covariant matrices for the two parts of the trajectory, and the symbol tr represents the trace operator. In this case, the whole trajectory was divided into three equal parts, and covariance analysis was performed on each third as described above; then the overlap between each two-thirds is calculated using the normalized overlap parameter. A value close to 0 indicates poor sampling, while a value close to 1 indicates complete sampling. Screened Ligands. The small molecules we screened in this study incorporated two widely used chemical databases, namely, the National Cancer Institute Diversity Set (NCDIS)26 and DrugBank small molecules.27 The former includes ∼2000 diverse structures, while the latter contains ∼1500 FDA approved drugs. Similar to our earlier studies,22,23 here, we used only 1883 structures from the NCI database that could be parametrized for use in AUTODOCK.40,41 In addition, the AKIs; CC3,42 X6D,43 and ZZL,44 were used as positive controls during the screening experiments. Virtual Screening. The overall virtual screening protocol followed the same procedure described in our earlier studies22,23 and used the AUTODOCK45 software along with all its supplemented utilities to prepare and perform the docking simulations. In short, all docking simulations were focused on the pocket formed by residues L139, G140, K141, G142, K143, G145, N146, L215, G216, T217, E260, N261, A273, and D274. The grid spacing for the AUTOGRID maps was 0.375 Å. The Lamarckian genetic algorithm (LGA) method was used with a maximum of 10 million energy evaluations, 100 trails, and an initial population of 400 random individuals. The rest of the parameters were set to the default values. The full list of small molecules was screened against six different rigid AK-A structures comprising five conformations obtained from the clustering analysis (described above) and a minimized version of the AK-A crystal structure.

G = EMM + Gsolv − TSsolute

(3)

The molecular mechanics and solvation terms can be further decomposed into residue contributions as follows: EMM = Eele + Evdw + Eint

(4)

For each ligand, the binding free energy within the AK-A protein is estimated as the difference between the bound and nonbound cases: AK‐A − ligand AK‐A − ligand ΔGo = ΔGgas + ΔGsolv ligand AK‐A + {ΔGsolv + ΔGsolv }

(5)

We used the same procedure of running the MM-PBSA calculations as well as the same parameters for the dielectric constants of the solute and solvent and ionic concentration as we did in our previous studies.22,23

■

RESULTS AND DISCUSSION Improving the selectivity and binding affinity of inhibitors to AK-A is a task well-suited to the use of high-throughput computational screening. Here, we describe the screening of the NCI diversity set26 and the DrugBank set of smallmolecules27 against six different AK-A models. These models represent the collective dynamics of a single AK-A conformational state as extracted from MD simulations on the PDB structure 1MQ4.28 The top 300 hits resulted from the initial screening against the wild type AK-A nucleotide binding site were then docked to a total of 85 MD derived conformations of 4574



Article

docked hits within the active site. During the 110 ns simulation we did not notice any dramatic conformational changes within the protein structure, indicating the stability of the protein toward a specific conformation. Despite the conformational stability of the protein overall structure, the binding site seems to be more flexible. This observation is clearly shown in Figure 2, where the wild type main-chain B-factors (averaged over

the T217D mutant and ranked by their binding energies. A subset of the molecules we describe below show high binding affinities toward the nucleotide binding site of the wild type AK-A and simultaneously are able to overcome steric contacts that are induced in the T217D mutation. Molecular Dynamics Simulations of AK-A. Many kinases, including AK-A, exist in several distinct conformational states. An inhibitor that binds to a particular state may stabilize the enzyme at this conformation and stop its enzymatic activity. Performing a screening exercise of the magnitude described below against all possible conformational states along with their dynamical substates would require enormous and unattainable computational resources. Therefore, in this study, we chose to target the inactive nucleotide bound conformational state of AK-A as it is available in the crystal structure 1MQ4. Other states of the enzyme will be visited in future studies. To construct a representative set of equilibrated AK-A structures for the screening of the above-mentioned chemical libraries, the protein was subjected to a 110 ns MD simulation. The duration of the simulation was predicted using the normalized overlap method and principal component analysis on the MD trajectory (see methodology), and similar to our previously published work,22,23 the objective of this long run was to introduce the conformational dynamics of the active site during subsequent docking analysis. Since all 300 wild-type top hits were used for the T217D mutant screening, we concluded that it was not necessary to generate an ensemble of mutant structures. Alternatively, we can use the same wild type representative structures after inducing the T217D mutation to all of them. The 300 wild-type top hits were then equilibrated within the mutant active site using a 2 ns MD simulation, followed by scoring their binding energies using the MM-PBSA method. The root-mean-square deviations (RMSDs) for the backbone atoms of the wild type enzyme for the 110 ns of the MD simulation show that the protein is stable throughout the simulation (Figure 1). The average RMSD fluctuation of the protein backbone is approximately 1.75 Å. The protein structure dynamics reached equilibration only after 1 ns, which was our main reason for choosing 2 ns as an adequate simulation time for the equilibration and relaxation of the

Figure 2. Average backbone atomic fluctuations (B-factors) for the different residues.

heavy atoms) of the active site residues (L139, G140, K141, G142, K143, G145, N146, L215, G216, T217, E260, N261, A273, and D274) are generally higher than the rest of the residues, especially residues 139 through 146 that appear to play an important role in kinase activity. Principal Component Analysis (PCA). Similar to our previous studies,22,23 we used PCA to reduce the original MD trajectory space of correlated variables to a significantly condensed set of independent variables that encompasses the essential dynamics of the system.38,39 Here, PCA was preformed on the AK-A active site atoms for residues L139, G140, K141, G142, K143, G145, N146, L215, G216, T217, E260, N261, A273, and D274, which comprise the binding site. The whole MD trajectory was projected onto the planes of the three dominant principal components (see Figure 3). As shown in Figure 3, the active site widely fluctuated around a single conformation throughout the MD simulations, indicating the flexibility of the active site and the rigidity of the overall protein structure. The trajectory spanned a single cluster, indicating a favorable folded conformation for the protein. To measure the completeness of conformational sampling of the MD simulation, we divided the whole trajectory into thirds and calculated their normalized overlaps. Table 1 shows a high overlap between the thirds, indicating that the three parts of the trajectory sufficiently sampled the same conformational space and, approximately, covered the dynamical variability available for the binding site. Ensemble-Based Virtual Screening. Accommodating the flexibility of the target protein within docking has been always a major point of research and innovation. Current docking algorithms can efficiently explore the conformational space of small molecules that retain a small number of rotatable bonds.46 Small peptides, however, and large proteins in general possess massive conformational spaces that are considerably unachievable to follow with the currently available computational resources. To deal with this problem, docking algorithms have incorporated several methods to introduce the target flexibility

Figure 1. RMSD of backbone atoms of the 110 ns MD trajectory for the Apo-wild type protein calculated relative to the initial reference structure as a function of simulation time. The upper panels show magnified views of the initial minimization, heating, and equilibration phases up to the first 2 ns and a sample of the MD simulation from 40 to 50 ns. 4575



Article

Figure 3. PCA for the AK-A nucleotide binding site. The whole 110 ns MD trajectory for the Apo-wild type AK-A was projected onto the three principal component planes. The histograms reflect the degree by which the trajectory explored certain regions on the projection with lighter colors indicating more frequently visited areas. The whole trajectory spanned a single but wide cluster, indicating a favorable conformation with a high degree of flexibility.

Table 1. Normalized Overlaps between the Three Thirds of the MD Simulation for the Aurora A Kinase Nucleotide Binding Site apo-enzyme 0.83 0.86 0.81

1st and 2nd 1st and 3rd 2nd and 3rd

in their calculations. These methods include, but not limited to, allowing a selected subset of side chains that are close to the ligand to rotate freely with some confinements to a rotamer library conformations, docking to an average structure, use of soft potential functions, and dock to an ensemble of target structures.46−48 Here, we followed the latter approach by docking the tested small molecules to a selected subset of protein structures that account for backbone dynamics of the target. To create this ensemble of dominant AK-A structures for docking simulations, we used RMSD conformational clustering on the Apo-protein MD trajectory (see methodology). Figure 4 compares the two clustering metrics that we used (i.e., DBI and SSR/SST) for different numbers of clusters. Although the SSR/SST ratio plateaued after 75 clusters, the DBI experienced local minima at cluster counts of 20, 40, 65, 70, 85, 105, 125, 200, 250, and 260. Therefor, we decided to take 85 as the optimal number of clusters. The first five clusters, at this cluster count, comprised almost 85% of the whole trajectory (data not shown). Overall, six different AK-A conformations were used as targets for six independent virtual screening simulations. Each simulation was carried out against the whole collection of ligand structures. The six used structures comprised the five representative structures of the five dominant clusters as well as a minimized version of AK-A crystal structure. Our objective was to reduce the massive number of MD trajectory snapshots without losing information about the conformational dynamics of the binding site. Although, as expected, the overall AK-A protein adopted limited conformational changes, the binding site itself took on different conformations, signifying the importance of accommodating the receptor flexibility in our screening of compounds. AUTODOCK Scoring. Following the clustering of the docked conformations, we used the most populated cluster for each ligand in our subsequent analysis. Only ligands with clusters that are populated with at least 20% of the docking conformations were considered. For each of the six virtual screening simulations that we carried out in this study, the docking pose, for each ligand, with the lowest binding energy as

Figure 4. Estimating the optimal number of clusters for the Apo-wild type AK-A MD trajectory. The SSR/SST started to plateau around a cluster count of 75, while the DBI has minima at cluster counts of 20, 40, 65, 70, 85, 105, 125, 200, 250, and 260. Therefore, 85 was considered as the optimal number of clusters.

predicted by the AUTODOCK scoring function within the most populated cluster was used as the representative conformation for this ligand. For each ligand, we selected only the ligand−protein combination that resulted in the lowest binding energy from the six different docking cases. The ligands were then ranked by their representatives’ binding energies, and top hits were retained for further analysis using MD simulations and MM-PBSA calculations (see below). Ranking of AK-A Top Hits Using the MM-PBSA Scoring Function. Although docking simulations explored all possible conformations of the ligands within the AK-A active 4576



Article

for the entire set of molecules that are being screened. As the binding energies we calculated using the MM/PBSA methodology confirmed the experimental findings in ranking of the used controls, we applied the same technique to analyze the top 300 hits suggested by docking scoring. Table 3 shows the top 20 hits as ranked based on their predicted MM-PBSA binding energies to the wild type protein. Screening of Wild Type Top Hits against the T217D Mutant. To understand the effect of the T217D mutation on the binding characteristics of the top 300 hits we previously identified using docking simulations (see above), we docked them to the 85 representative cluster conformations after inducing the T217D mutation to their structures. Similar to the pose clustering and ranking procedure described above, for each compound, the binding energy and the population of the largest docking cluster provided the criteria for choosing the most preferred protein conformation out of the 85 used representative structure. Each optimal ligand−protein conformation was then used as a starting point for an MD simulation followed by calculating their binding energy using the MM-PBSA method. Ligands with the lowest binding energies in this case are expected to be specific toward the mutant enzyme. The top 20 mutant-specific hits as calculated by the MM-PBSA method are shown in Table 4. Dual (wildtype/mutant) inhibitors are those hits that can bind with low but similar binding affinities to the two protein variants, and they are listed in Table 5. Binding Modes of the Identified Inhibitors. Our calculations revealed several novel bicyclic scaffolds that bind to the AK-A nucleotide binding site and are suitable for further modifications and optimization (see Tables 3−5). As we expected, control ligands showed significant binding in the deep nucleotide binding cleft of the inactive AK-A. Figure 5 describes the binding trends to the wild type protein of three selected hits compared to those of the bound nucleotide and the control CC3 inhibitor. Interestingly, our first example (Figure 5a) resembles the nucleotide in its binding to the active site where its structure mimics the nucleoside while the rest of the molecule forms a bulky tail that completely occupies the depth of the nucleotide binding groove. The second example (Figure 5b) forms a long chain structure that occupies the whole nucleotide site. Part of the molecule protrudes outside the binding groove similarly to the ZZL control. The third compound (Figure 5c) acts oppositely to the first example, where the tail of the compound occupies the mouth of the nucleotide site and the bulky head binds to the location of the base-binding site. Other examples (Table 4) can bind to the active site with unique binding modes that exploit the variations in their scaffolds’ structures and the different conformations adopted by the active site. Most of the top inhibitors we found here form a few hydrogen bonds with the backbone atoms of the active site in the kinase hinge region. Remarkably, we did not notice any water bridging hydrogen bonds mediating the interactions of our hits with the protein. Effect of T217D. Introducing the T217D mutation to the protein structure dramatically affected the binding modes of most of the compounds obtained in the screening of the wild type AK-A protein, including the three controls. As shown in Figure 6, the T217D residue acts as a gate to the nucleotidebinding domain allowing or preventing the compounds to penetrate deeper within the binding cleft. The negative charge on the side chain of D217 also seems to play a significant role in this function. The binding of the inhibitors to the mutated

site, the rigidity of the target protein presented a great hindrance when attempting to reveal the actual binding mode. To solve this problem, we employed in this study the use of MD simulations to predict the optimal binding conformation between the two interacting entities and relax them toward more realistic interactions compared to docking simulations. As we also simulated each ligand−protein complex in an explicit water and ions environment, we were able to gain further insights on the effects of the solvent on the binding modes and whether they bridge any interactions between the ligand and protein (e.g., water-bridged hydrogen bonds). Furthermore, running MD simulations on each system allowed us to use a more accurate scoring method to rank our identified hits and carefully compare their relative binding energies. Similar to our previous studies, here, we used the MM-PBBA method introduced by Kollman et al.49 to predict the binding energies of the top 300 hits we identified against the wild type AK-A binding site. We also used the same methodology to rank the top 300 hits in their binding to the mutated AK-A structures (see below). Over the last decades, the MM-PBSA has gained growing popularity in estimating binding energies, especially, as a postdocking analysis tool to rescore docking results.50 The apparent IC50 values for CC3, X6D, and ZZL were determined to be 42 nM,42 12 nM,43 and 12.5 nM44 at 25 °C, respectively. To compare these IC50 values given as Ki to the estimated binding energies (ΔG), we used the relation: ΔG = RT ln K i

where R is the gas constant, R = 1.987 cal K−1 mol−1 and T is the temperature in Kelvin. Table 2 lists the predicted binding energies for the three controls that we used to their experimental values. Docking Table 2. Relative Ranking of the CC3, X6D, and ZZL Known AK-A Inhibitors Using the MM-PBSA and AUTODOCK Scoring Functions and Compared to Available Experimental Values ranking (kcal/mol) compound

MM/PBSA

AUTODOCK

experimental

CC3 X6D ZZL

−9.4 −10.4 −9.7

−8.4 −9.1 −12.1

−10.13 −10.87 −10.85

scoring function either overestimated the binding energies of the controls (see ZZL for example) or undermined their affinities (CC3). Although docking energies showed that the three compounds bind with high affinities within the active site, the ranking of the controls was not consistent with that found experimentally, giving ZZL the lead. On the other hand, within an error of 1 kcal/mol, the MM-PBSA calculations for the interactions of the three known AK-A inhibitors are in excellent agreement with experimental data. It is noteworthy to mention that the accurate prediction of the relative ranking of a set of screened compounds is the principle outcome of any virtual screening simulation.24,51,52 Nevertheless, as a first approximation to estimate the absolute binding energies, only nonpolar solvation energies (which corresponds to the loss of solvent accessible surface area upon binding) and vibrational contribution to the entropy as predicted by normal-mode analysis is used for a small subset of the MD trajectory.53,54 This approximation is used as it significantly reduces the computational time required, as these calculations are repeated 4577



Article

Table 3. Group A: A List of Inhibitors That Bind Strongly to Both the Wild Type and T217D Mutant Structures (Dual Inhibitors) Ranked by the MM-PBSA Method

4578



Article

Table 3. continued

the DFG motif (residues D274−W276). The newly adopted orientation within the mutated structure performed synergistically to increase its binding affinity. Therefore, we expect that such compounds (Table 4) are more selective for the mutated structure over the wild type protein and can be used as an alternative medication in combating the induced drug resistance during the course of cancer treatment. Analogously

variant induced many conformational changes to the binding site. These local adjustments were facilitated by the inherent plasticity within the nucleotide binding cleft specifically including the p-loop (residues G140−G145). The most noticeable alternating binding mode is for compound 2 (Figure 6c), where D217 pushed the long chained structure to profoundly infiltrate within the nucleotide-binding cleft toward 4579



Article

Table 4. Group B: A List of Inhibitors That Bind Stronger to the Wild Type Structure Ranked by the MM-PBSA Method

4580



Article

Table 4. continued

4581



Article

Table 4. continued

4582



Article

Table 5. Group C: A List of Inhibitors That Bind Stronger to the T217D Mutant Structure Ranked by the MM-PBSA Method

4583



Article

Table 5. continued

to most of the examined structures, the CC3 inhibitor compound (Figure 6b) shifts outside of the binding site. This displacement considerably reduced its binding affinity to the enzyme confirming what was observed experimentally regarding the T217D mutation for AK-A inhibitors. As a result of our analysis of the structural and energetic consequences of the T217D mutation, we classified our top hits into three main groups. The first (Group A) and, perhaps the

most important, represents the dual inhibitors that can simultaneously bind with high affinity to the wild type and to the mutated structures. This group includes very few compounds listed in Table 5. The central feature of this group is that they are small molecules and have the ability to adapt to the two conformations and fit nicely within the nucleotide-binding cleft. An example of this group is compound 1 (Figure 6b) which is predicted to bind similarly in both cases 4584



Article

Figure 5. Alignment of the bound nucleotide (red), the CC3 inhibitor (green), and our hits (gray). The binding location is shown in ribbon representations with different conformations that accommodate the attached ligand. The three compounds occupy the nucleotide binding location within AK-A.

Figure 6. Effect of inducing the T217D mutation on the binding modes of the identified hits (b−f) compared to the CC3 inhibitor (a). The binding sites are shown in ribbon representations for both the wild type (white) and mutant (yellow) with both the threonine residue (red) and the aspartic acid residue (blue) are shown in surface representations. The ligands as bound to the wild type protein are shown in gray, while those within the mutant structure are shown in green.

within the mutated nucleotide-binding cleft compared to the wild type case. Binding Energy Decomposition and Its Implication on Future Redesign of Dual Inhibitors. To gain further insights into the impact of the individual residues on the binding affinities, we decomposed the binding energies into residue contributions. Figure 7 illustrates the total enthalpic contribution per residue for the 12 dual inhibitors listed in Table 3. Here, we focused on the 31 residues closest to all ligands within the binding site and compared their energies for both the wild type and mutant cases. We also broke down these

with comparable binding affinities, which are still higher than the three controls. The second group (Group B) binds stronger to the wild type compared to the mutant. These compounds could not endure the electrostatic energy induced by the negative charge introduced in the side chain of D217. Accordingly, they are shifted away from the binding cleft and bound to the protein with a significantly lower affinity. An example of such compounds is compound 4 (Figure 6e). The last group (Group C) includes the compounds that can bind to the mutated structure tighter than to the wild type. Examples include compounds 2, 3, and 5. These compounds bind deeper 4585



Article

Figure 7. Binding energy decomposition per residue contribution. (A) The wild type structure. (B) The mutant structure.

energies for each residue into their electrostatic, van der Waals, and solvation parts and even divided them further to either side chain or backbone contributions (Supplementary Figures 1−9). Interestingly, and as we expected for dual inhibition, the general behavior of these residues for each ligand and for the two AK-A variants is almost the same. The individual residue contributions ranged from 1 to −3.5 kcal/mol for the wild type structure, while for the mutant AK-A it ranged from 2 to −4 kcal/mol. The most energetically attractive interactions are for

residues L139, V147, K162, and L263. These interactions are mainly hydrophobic (Supplementary Figure 1) as their electrostatic contributions are almost zero expect for Lys162 for the ZINC1250294 (Supplementary Figure 2) where the large attractive electrostatic attractive energy has been balanced by a considerable repulsive electrostatic solvation energy (Supplementary Figure 3). These dominant hydrophobic interactions are mainly attributed to the side chains of these residues rather than to their backbones (Supplementary Figures 4586



Article

can simultaneously bind with high affinity to the wild type and mutated structures; Group B binds stronger to the wild type compared to the mutant; and Group C that can bind to the mutated structure tighter than to the wild type. The results of this study provide researchers with an initial scaffold that targets the inactive form of AK-A. This is a first step toward the development of potent inhibitors of this important target in human cancers, and future studies will investigate additional conformations of both active and mutant AK-A structures.

4−9). A number of residues seem to have a negligible effect on the binding energy for most of the ligands for both the wild type and the mutant. These residues include F144, L149, A160, L164, P214, P259, W277, and G291 (Figure 7). Their vdW interaction is almost zero for all of the compounds (Supplementary Figures 1). Contrary to the attractive residues described above, the repulsive gas phase electrostatic energies are balanced by equal but favorable electrostatic solvation energies (Supplementary Figures 2 and 3). The most unfavorable interactions are concentrated on residues Glu181, Glu211, Glu260, and Asp274 and are linked to their side chains (Supplementary Figures 4−9). Although these residues have attractive solvation electrostatic energies, their vdW interactions are almost neutral, and their gas phase electrostatic energies are much higher than the solvation energies, making their overall interaction less favorable. This observation is clear for ligands ZINC12503294, NSC74702, NSC285656, NSC224314, and NSC17770. Excluding ZINC12503294, NSC224314, and NSC127494, the T217D mutations had favorable interactions with the rest of the ligands. Nevertheless, for these three residues, the unfavorable interaction due to the mutation was compensated for by interaction with K162 (Figure 7). The binding modes described in Figure 6 and the detailed energy decomposition described herein provide the rational tools to modify or even redesign these ligand structures in a way that makes them more potent to both the wild type and the mutant structures. For example, one can focus on the less attractive or neutral residues and modify the compounds in order to have a more favorable vdW interactions with these residues since these interactions are almost neutral as described above. Another option is to add an electrostatic moiety within the ligand close to one or more of these residues to bias the energetic balance toward more favorable gas phase electrostatic energy, which is currently compensated for by the solvation electrostatic energy. Using the data described here, similar strategies can be followed to generate new ligands that can structurally and energetically fit within the interacting residues and bind with high affinities to the wild type and mutant structures.

■

ASSOCIATED CONTENT

S Supporting Information *

vdW, electrostatic, solvation electrostatic, side chain vdW, side chain electrostatic, side chain solvation electrostatic, backbone vdW, backbone electrostatic, and backbone solvation electrostatic energy decomposition per residue contribution. This material is available free of charge via the Internet at http:// pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*Department of Oncology, Faculty of Medicine and Dentistry, University of Alberta, Edmonton, Alberta, T6G 1Z2, Canada. E-mail: [email protected]. Tel.: (780) 432-8906. Fax: (780) 432-8425. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Funding for this research was provided through grants from Canada Excellence Research Chairs (CERC), Canadian Institutes of Health Research (CIHR), Alberta InnovatesHealth Solutions (AIHS), and the Canadian Association of Gastroenterology (CAG). All of the molecular dynamics simulations and virtual screening experiments were produced using the Blue Gene/Q super computer.

■

ABBREVIATIONS AK, Aurora kinase; AK-A, Aurora kinase-A; AKI, Aurora kinase inhibitor; DBI, Davies-Bouldin index; DOPE, discrete optimized protein energy; MD, molecular dynamics; MM-PBSA, molecular mechanics Poisson−Boltzmann surface area; PCA, principal component analysis; PDB, Protein Data Bank; RCS, relaxed complex scheme; RMSD, root-mean-square deviation; SOM, self-organizing maps; VS, virtual screening

■

CONCLUSIONS Human AK-A is an important mitotic regulator belonging to a small family of AKs. Each AK functions as a key regulator throughout mitosis, and they have been shown to become overexpressed in human cancers, making them excellent candidates for the design of novel chemotherapeutic agents.1−3 While encouraging results have been achieved in preclinical and early phase I and II trials for several AKIs, knowledge of the structural features that affect their sensitivity remains limited. Here, we investigated this research question at the atomistic level to discover novel inhibitors for AK-A targeting its nucleotide binding site. We applied state-of-the-art molecular modeling techniques to understand the activity and mode of binding of currently available AK-A inhibitors. We also analyzed the effects induced by one of the most frequent AK-A mutations, namely, the T217D alternation. This knowledge combined with the massive computational power of the Blue Gene/Q super computer allowed us to utilize a rigorous and modified version of the RCS methodology.23 We redocked the top hits from the initial screening to the mutant structure and ranked them for their ability to bind to the T217D variant with a sufficiently high affinity. We have classified our top hits into three main groups: Group A is composed of dual inhibitors that

■

REFERENCES

(1) Fenner, A. Bladder cancer: Aurora kinase inhibitors light up the therapeutic horizon in bladder cancer. Nat. Rev. Urol. 2013, 10, 184. (2) Wissing, M. D.; van Diest, P. J.; van der Wall, E.; Gelderblom, H. Antimitotic agents for the treatment of patients with metastatic castrate-resistant prostate cancer. Expert Opin. Invest. Drugs 2013, 22, 635−661. (3) Xie, L.; Meyskens, F. L., Jr. The pan-Aurora kinase inhibitor, PHA-739358, induces apoptosis and inhibits migration in melanoma cell lines. Melanoma Res. 2013, 23, 102−113. (4) Sen, S.; Zhou, H.; White, R. A. A putative serine/threonine kinase encoding gene BTAK on chromosome 20q13 is amplified and overexpressed in human breast cancer cell lines. Oncogene 1997, 14, 2195−2200. (5) Shi, H.; et al. Single nucleotide polymorphisms in the 20q13 amplicon genes in relation to breast cancer risk and clinical outcome. Breast Cancer Res. Treat. 2011, 130, 905−916.

4587



Article

(6) Liu, Q.; et al. Aurora-A abrogation of p53 DNA binding and transactivation activity by phosphorylation of serine 215. J. Biol. Chem. 2004, 279, 52175−52182. (7) Kollareddy, M.; et al. Aurora kinase inhibitors: progress towards the clinic. Invest. New Drugs 2012, 30, 2411−2432. (8) Fancelli, D.; et al. Potent and selective Aurora inhibitors identified by the expansion of a novel scaffold for protein kinase inhibition. J. Med. Chem. 2005, 48, 3080−3084. (9) Schoffski, P.; et al. Phase I, open-label, multicentre, doseescalation, pharmacokinetic and pharmacodynamic trial of the oral aurora kinase inhibitor PF-03814735 in advanced solid tumours. Eur. J. Cancer 2011, 47, 2256−2264. (10) Aliagas-Martin, I.; et al. A class of 2,4-bisanilinopyrimidine Aurora A inhibitors with unusually high selectivity against Aurora B. J. Med. Chem. 2009, 52, 3300−3307. (11) Zhao, B.; et al. Modulation of kinase-inhibitor interactions by auxiliary protein binding: crystallography studies on Aurora A interactions with VX-680 and with TPX2. Protein Sci. 2008, 17, 1791−1797. (12) Karaman, M. W.; et al. A quantitative analysis of kinase inhibitor selectivity. Nat. Biotechnol. 2008, 26, 127−132. (13) He, W.; et al. AURKA suppression induces DU145 apoptosis and sensitizes DU145 to docetaxel treatment. Am. J. Transl. Res. 2013, 5, 359−367. (14) Greenman, C.; et al. Patterns of somatic mutation in human cancer genomes. Nature 2007, 446, 153−158. (15) Sloane, D. A.; et al. Drug-resistant aurora A mutants for cellular target validation of the small molecule kinase inhibitors MLN8054 and MLN8237. ACS Chem. Biol. 2010, 5, 563−576. (16) Howard, S.; et al. Fragment-based discovery of the pyrazol-4-yl urea (AT9283), a multitargeted kinase inhibitor with potent aurora kinase activity. J. Med. Chem. 2009, 52, 379−388. (17) Morshed, M. N.; et al. Computational approach to the identification of novel Aurora-A inhibitors. Bioorg. Med. Chem. 2011, 19, 907−916. (18) Sardon, T.; et al. Uncovering new substrates for Aurora A kinase. EMBO Rep. 2010, 11, 977−984. (19) Yang, Y.; et al. Molecular dynamics and free energy studies on Aurora kinase A and its mutant bound with MLN8054: insight into molecular mechanism of subtype selectivity. Mol. BioSystems 2012, 8, 3049−3060. (20) Lin, J. H.; Perryman, A. L.; Schames, J. R.; McCammon, J. A. Computational drug design accommodating receptor flexibility: the relaxed complex scheme. J. Am. Chem. Soc. 2002, 124, 5632−5633. (21) Friesen, D. E.; et al. Discovery of small molecule inhibitors that interact with gamma-tubulin. Chem. Biol. Drug Des. 2012, 79, 639− 652. (22) Barakat, K.; Mane, J.; Friesen, D.; Tuszynski, J. Ensemble-based virtual screening reveals dual-inhibitors for the p53-MDM2/MDMX interactions. J. Mol. Graphics Modell. 2010, 28, 555−568. (23) Barakat, K.; Tuszynski, J. Relaxed complex scheme suggests novel inhibitors for the lyase activity of DNA polymerase beta. J. Mol. Graphics Modell. 2011, 29, 702−716. (24) Jordheim, L. P.; et al. Small molecule inhibitors of ERCC1-XPF protein-protein interaction synergize alkylating agents in cancer cells. Mol. Pharmacol. 2013, 84, 12−24. (25) Barakat, K. H.; et al. Virtual screening and biological evaluation of inhibitors targeting the XPA-ERCC1 interaction. PloS One 2012, 7, e51329. (26) http://dtp.nci.nih.gov/branches/dscb/diversity_explanation. html. (27) Wishart, D. S.; et al. DrugBank: a comprehensive resource for in silico drug discovery and exploration. Nucleic Acids Res. 2006, 34, D668−672. (28) Nowakowski, J.; et al. Structures of the cancer-related Aurora-A, FAK, and EphA2 protein kinases from nanovolume crystallography. Structure 2002, 10, 1659−1667. (29) Sali, A.; Blundell, T. L. Comparative protein modelling by satisfaction of spatial restraints. J. Mol. Biol. 1993, 234, 779−815.

(30) Pettersen, E. F.; et al. UCSF Chimeraa visualization system for exploratory research and analysis. J. Comput. Chem. 2004, 25, 1605−1612. (31) Dunbrack, R. L., Jr. Rotamer libraries in the 21st century. Curr. Opin. Struct. Biol. 2002, 12, 431−440. (32) Dolinsky, T. J.; et al. PDB2PQR: expanding and upgrading automated preparation of biomolecular structures for molecular simulations. Nucleic Acids Res. 2007, 35, W522−525. (33) Kalé, L.; et al. NAMD2: Greater Scalability for Parallel Molecular Dynamics. J. Comput. Phys. 1999, 151, 283−312. (34) Hornak, V.; et al. Comparison of multiple Amber force fields and development of improved protein backbone parameters. Proteins 2006, 65, 712−725. (35) Jakalian, A.; Jack, D. B.; Bayly, C. I. Fast, efficient generation of high-quality atomic charges. AM1-BCC model: II. Parameterization and validation. J. Comput. Chem. 2002, 23, 1623−1641. (36) Davies, D. L.; Bouldin, D. W. A cluster separation measure. IEEE Trans. Pattern Anal. Mach. Intell. 1979, 1, 224. (37) Shao, J.; Tanner, S. W.; Thompson, N.; Cheatham, T. E. Clustering Molecular Dynamics Trajectories: 1. Characterizing the Performance of Different Clustering Algorithms. J. Chem. Theory Comput. 2007, 3, 2312. (38) Garcia, A. E. Large-amplitude nonlinear motions in proteins. Phys. Rev. Lett. 1992, 68, 2696−2699. (39) Hess, B. Convergence of sampling in protein simulations. Phys. Rev. E 2002, 65, 031910. (40) Goodsell, D. S.; Morris, G. M.; Olson, A. J. Automated docking of flexible ligands: applications of AutoDock. J. Mol. Recogn. 1996, 9, 1−5. (41) Morris, G. M.; et al. AutoDock4 and AutoDockTools4: Automated docking with selective receptor flexibility. J. Comput. Chem. 2009, 30, 2785−2791. (42) Tari, L. W.; et al. Structural basis for the inhibition of Aurora A kinase by a novel class of high affinity disubstituted pyrimidine inhibitors. Bioorg. Med. Chem. Lett. 2007, 17, 688−691. (43) Bavetsias, V.; et al. Imidazo[4,5-b]pyridine derivatives as inhibitors of Aurora kinases: lead optimization studies toward the identification of an orally bioavailable preclinical development candidate. J. Med. Chem. 2010, 53, 5213−5228. (44) Dodson, C. A.; et al. Crystal structure of an Aurora-A mutant that mimics Aurora-B bound to MLN8054: insights into selectivity and drug design. Biochem. J. 2010, 427, 19−28. (45) Ali, H. I.; Nagamatsu, T.; Akaho, E. Structure-based drug design and AutoDock study of potential protein tyrosine kinase inhibitors. Bioinformation 2011, 5, 368−374. (46) Kroemer, R. T. Structure-based drug design: docking and scoring. Curr. Protein Pept. Sci. 2007, 8, 312−328. (47) Baxter, C. A.; Murray, C. W.; Clark, D. E.; Westhead, D. R.; Eldridge, M. D. Flexible docking using Tabu search and an empirical estimate of binding affinity. Proteins 1998, 33, 367−382. (48) Totrov, M.; Abagyan, R. Flexible protein-ligand docking by global energy optimization in internal coordinates. Proteins 1997, No. Suppl 1, 215−220. (49) Kollman, P. A.; et al. Calculating structures and free energies of complex molecules: combining molecular mechanics and continuum model. Acc. Chem. Res. 2000, 33, 889. (50) Kuhn, B.; Gerber, P.; Schulz-Gasch, T.; Stahl, M. Validation and use of the MM-PBSA approach for drug discovery. J. Med. Chem. 2005, 48, 4040−4048. (51) Leonis, G.; Steinbrecher, T.; Papadopoulos, M. G. A Contribution to the Drug Resistance Mechanism of Darunavir, Amprenavir, Indinavir, and Saquinavir Complexes with HIV-1 Protease Due to Flap Mutation I50V: A Systematic MM-PBSA and Thermodynamic Integration Study. J. Chem. Inf. Model. 2013, 53, 2141−2153. (52) Fu, T.; Jin, Z.; Xiu, Z.; Li, G. Binding free energy estimation for protein-ligand complex based on MM-PBSA with various partial charge models. Curr. Pharm. Des. 2013, 19, 2293−2307. 4588



Article

(53) Liu, J.; et al. Identification of Novel Potential beta-N-Acetyl-DHexosaminidase Inhibitors by Virtual Screening, Molecular Dynamics Simulation and MM-PBSA Calculations. Int. J. Mol. Sci. 2012, 13, 4545−4563. (54) Grazioso, G.; et al. Design of novel alpha7-subtype-preferring nicotinic acetylcholine receptor agonists: application of docking and MM-PBSA computational approaches, synthetic and pharmacological studies. Bioorg. Med. Chem. Lett. 2009, 19, 6353−6357.

4589


A Computational Model for Overcoming Drug ... - ACS Publications

Recommend Documents