Analysis of Biphenyl-Type Inhibitors Targeting the Eg5 α4/α6 Allosteric Pocket.

Eg5 is a mitotic kinesin protein that plays an important role in the formation and maintenance of the bipolar spindle during the mitotic phase. Due to its potentially reduced side effects in cancer therapy, Eg5 is considered to be an attractive target for developing anticancer inhibitors. Herein, we report a computational modeling study involving biphenyl-type inhibitors known to interact with the α4/α6 allosteric pocket of Eg5. Compared to the well-known α2/L5/α3 allosteric inhibitors, biphenyl-type inhibitors show a unique activity profile. In the Eg5-PVZB1194 (a biphenyl-type inhibitor) crystal structure, loop L11, which is located in the entrance of the α4/α6 allosteric-binding pocket, is missing due to crystal-packing effects. To better understand the role of this flexible loop upon biphenyl-type inhibitor-binding, MD simulations were performed to observe the L11 conformations from different states. It was demonstrated that L11 was more stabilized and showed less fluctuation when PVZB1194 was bound to Eg5. Residue Asn287 from L11 forms hydrogen bonding to the sulfone group of PVZB1194, whereby L11 moves inward to the α4/α6 allosteric pocket and moves away from the pocket in absence of the inhibitor. Pharmacophore, three-dimensional (3D)-QSAR, and ADME studies of biphenyl-type inhibitors of Eg5 were also performed. A best pharmacophore model, DDRRH.6, was generated, having correlation coefficients in the 3D-QSAR study of R2 = 0.81 and Q2 = 0.64. Furthermore, docking studies were carried out to observe the interaction between the remaining biphenyl-type inhibitors with Eg5. In addition, on the basis of fragment docking, a structure-based pharmacophore was generated, which shares good overlap of the DHRR features of the pharmacophore model DDHRR.6. The structure-based pharmacophore also contains extra hydrogen-bond acceptors and hydrophobic groups, features which provide possibilities in developing new or improved series of compounds.

Introduction

Kinesins are a superfamily of motor proteins that are involved in a diverse range of physiological functions, such as mitotic spindle assembly, chromosome segregation, vesicular trafficking, and so forth. They are sometimes described as “nanomotors” as they utilize energy from ATP hydrolysis to generate force and transport intracellular cargos along the microtubules (MTs).[1,2] Eg5, also termed KIF11, is a member of the mitotic kinesin-5 family, which is a subgroup of kinesins that only functions during the mitotic phase of cell division. Eg5 acts on mitosis to help bipolar spindle formation, chromosome alignment, and segregation.[3]

In many proliferative tissues, Eg5 is found to be overexpressed. However, it is rarely detected in the nonproliferative ones. Accordingly, the overexpression of Eg5 was found in various cancers, including breast, lung, ovarian, bladder, and pancreatic cancers as well as in leukemia,[4−6] making it a potential target for cancer therapy.[7] Eg5 inhibition prevents separation of centrosomes and formation of mitotic spindle, which leads to formation of “monoasters” consisting of monopolar spindles, and further triggers spindle checkpoint proteins to cause direct mitotic arrest and inhibition of cell division.[8−10] It has been demonstrated that in human xenograft models inhibition of Eg5 causes complete cell death in many cancer-cell lines and displays in vivo antitumor activity.[11,12] As Eg5 only functions during mitosis, inhibitors of Eg5 do not affect the nonproliferative cells. Furthermore, because Eg5 does not exist in the adult peripheral nervous system, inhibitors of Eg5 thus do not induce neuropathic side effects, which are otherwise commonly found for inhibitors of tubulin also developed with the aim to affect mitosis.[7,13]

Ispinesib (SB-743921) and filanesib (ARRY-520) are amongst the inhibitors that specifically target Eg5, which are currently under clinical development.[14−16] Instead of interacting with the ATP-binding site, these inhibitors bind to an extensively studied allosteric site about 10 Å away, which is composed of helix α2, loop L5, and helix α3 (Figure A).[17−19] Biochemical studies indicate that by inhibiting ADP release from Eg5 this set of compounds functions as ATP-uncompetitive inhibitors.[7,14,20] Although they display potent anticancer activity, drug-resistant mutants, D130V and A133D, located in loop L5 have already been identified in cell cultures. Therefore, a new series of inhibitors that target a novel inhibitor-binding pocket could be developed and utilized either alone or together, with known Eg5 inhibitors targeting the α2/L5/α3 region.

Figure 1

(A) Eg5 with bound AMMPNP (PDB ID: 3HQD), showing pockets α2/L5/α3 and α4/L11/α6. (B) PVZB1194-bound Eg5 (blue) superposed to AMMPNP-bound Eg5 (green). Residues Glu129 and Thr104 (in yellow) from PVZB1194-bound Eg5 clash with the nucleotide binding, whereas residues Glu129 and Thr104 (in pink) from nucleotide-bound Eg5 are not in the ligand-binding site.

Biphenyl compounds have been identified as another type of Eg5 inhibitor by several groups.[21,22] This set of compounds was shown to inhibit ispinesib-resistant tumor cells that harbor D130V and A133D mutations with high activity.[23−25] One of the biphenyl-type inhibitors, PVZB1194, was shown to bind to a different allosteric pocket, composed by helices α4 and α6. This pocket is located 15 Å from the nucleotide-binding pocket (Figure A). Once the biphenyl-type inhibitors bind to the α4/α6 allosteric pocket, they are capable of distorting the ATP-binding site through Tyr104 residues, which displace the residues Thr107 and Glu129 into the nucleotide-binding pockets and thereby suppress the binding of ATP and function as ATP-competitive inhibitors[26] (Figure B). Hence, the ATP-binding site of Eg5–PVZB1194 contains neither ADP nor ATP. This is different from the benzimidazole type inhibitor BI8,[27] which binds to the α4/α6 pocket, when Eg5 is in complex with ATP or ADP and thus, similary to conventional allosteric inhibitors, inhibits the release of ADP. Hence, the allosteric effect asserted by biphenyl-type inhibitors, such as PVZB1194, on the binding of ATP is unique. PVZB1194 is furthermore shown to potently inhibit Eg5 in a MT-dependent manner and induce a monoastral phenotype and cause a mitotic arrest.[26] Because the α4/α6 allosteric pocket was only recently detected and biphenyl-type inhibitors targeting this pocket display a unique mechanism in inhibiting ATP binding, it is of interest to investigate in more detail the interaction between PVZB1194 (and other biphenyls) and Eg5. The aim is to assist in developing a new series of inhibitors targeting this pocket. Such compounds can then be applied alone or in combination with inhibitors targeting the α2/L5/α3 allosteric pocket. To gain further insight, molecular dynamics (MD) simulations were applied to investigate the conformational differences upon PVZB1194 binding versus free protein and in particular, loop L11 at the entrance of the α4/α6 pocket, which is missing in the crystal structure due to its high flexibility. A common pharmacophore hypothesis of the biphenyl-type inhibitors was generated and analyzed, and an atom-based three-dimensional (3D)-QSAR model was evaluated to understand the relationship between their biological activity and molecular structure. In addition, physicochemical property calculations and molecular docking studies were also performed, aimed at assisting in the design of novel allosteric Eg5 inhibitors targeting the α4/α6 pocket.

Results and Discussion

MD Simulations

Four regions were missing in the crystal structure of Eg5 complexed to PVZB1194,[26] corresponding to residues 57–58 (L2), 108–126 (regions of the P-loop, N-terminal of α2, and L5), 225–331 (L9), and 270–287 (L11). Of these, L11 is located in the entrance of the α4/α6 allosteric pocket and might be capable of forming hydrogen bonds and other contacts with a bound ligand. When there is a nonhydrolyzable ATP analogue (AMPPNP) bound to Eg5, L11 forms a completely ordered loop, possibly due to interactions provided by the γ-phosphate of AMPPNP having a stabilizing effect on the residues of L11.[26,31] However, as the biphenyl-type inhibitors are competitive for ATP binding, the conformation of L11 in the Eg5–PVZB1194 complex is disordered and missing from the crystal structure. To explore the interaction between L11 and PVZB1194, we mended L11 by using the protein preparation wizard, as outlined above, together with other missing regions and performed MD simulations to sample its conformation with or without the PVZB1194 inhibitor located in the pocket.

The root-mean-square deviation (RMSD) values (Figure A) overlap between the unbound Eg5 proteins than between the Eg5–PVZB1194 complex. The root-mean-square fluctuation (RMSF) values show that with the inhibitor bound L11 is more stable (encircled region in Figure B), with PVZB1194 forming an extra hydrogen bond to Asn287 in L11. PVZB1194 also forms an extra stable hydrogen bond to Glu345 and a weak one to Ser269 (Figure C,D). Therefore, L11 moves toward the α4/α6 pocket in the complex structure, whereas in the unbound protein structure it swings out into a much more open conformation (Figure E). Thus, we can conclude that L11 is stabilized by the ligand binding of residue Asn287 from the L11 hydrogen bonds to sulfone group. The structure–activity relationship data has shown that the sulfone group is important for biological inhibitory activity,[22] and the above observation from the MD simulations fully explains this inhibition. The free-energy analysis of L11 shows how the inhibitor binding influences its conformation. The free energy was calculated by projecting along the distance of the center of mass (COM) of L11 to α6 and α4, respectively. Only one energy minimum is observed for either the free-protein or the complex structure and centers on the distance values between L11 and α6 (α4) of 1.8–2 nm (1.1–1.4 nm) for the free protein and 1.5–1.6 nm (1.5–1.6 nm) when in the complex. Therefore, in the free protein structure, L11 is, relatively speaking, displaced closer to α4, whereas with the ligand-bound structure, it prefers to stay in between the two helices (Figure A,B). In addition, the free-energy minimum is much more narrow and localized for the complex system.

Figure 2

(A) C-α RMSD of the complex and protein during the 500 ns MD simulations. (B) RMSF of the complex and protein during the simulation. (C) Hydrogen bonding between PVZB1194 and Eg5 residues (Asn287, Ser269, and Glu345) during simulation of the complex. (D) Residues of Eg5 forming hydrogen bonding to the PVZB1194 sulfone group and Tyr104 location in Eg5. (E) L11 conformation in the complex (green) and free-protein structure (magenta).

Figure 3

Free-energy analysis (kJ/mol) for (A) free Eg5 protein and (B) Eg5–PVZB-1194 complex structure, as functions of distance of COM of L11 to α4 (y axis) and α6 (x axis).

Tyr104 is located at the bottom of the allosteric pocket (Figure E) and was proposed to be the key residue leading to the ATP-competitive effect of PVZB1194 inhibition(26). After removing PVZB1194 from the pocket, Tyr104, however, remains in the same position as in the complex structure during the whole MD simulation by forming a stable hydrogen bond to the residue Glu270 of L11. Residues Thr107 and Glu129 are displaced into the nucleotide-binding site and thereby suppress the binding of ATP. In the 500 ns MD simulations of both PVZB1194-bound and unbound Eg5 structures, these two residues remain in the same placement as that causing interference with nucleotide binding, although the loop between Thr107 and Glu129 moves slightly in a manner normal for loop fluctuations. Thus, more enhanced sampling methods may be required in this part to observe the larger conformational changes associated with Tyr104 shifting and Thr104 and Glu129 displacement. However, we emphasize that PVZB1194 binding induces significant changes to the nucleotide pocket as compared with the normal undisturbed one such that the placement of ATP into the pocket then becomes essentially impossible during the process of simulations. The inhibitor binding distorts Thr107 and Glu129 to occupy the nucleotide-binding pocket and hinders ATP from entering.

To further explore this aspect, we also built a model with L11 in its existing AMMPNP-bound Eg5 conformation. L11 is shown, for the first time, to feature a completely ordered loop in the crystal structure of Eg5.AMMPNP compared with any other Eg5 crystal structure with ADP binding or Eg5 alone. Here, we define this as an ordered conformation and refer to the previous as a disordered conformation. All other regions were left unchanged, as detailed in the previous section. In the fully ordered conformation, residues of L11 appear to locate in the α4 helix. Next, we superposed this ordered conformation to the PVZB1194-bound Eg5. It is observed that even though L11 is in an ordered conformation PVZB1194 is able to fit into the binding pocket (Figure A,B). We then carried out MD simulations on the Eg5 structure with an ordered L11 conformation, alone or interacting with PVZB1194, to investigate the impact of the ligand on the L11 conformation. We can observe from the RMSF calculations over the whole simulation (Figure C) that L11 is again much more stabilized with PVZB1194 bound to the pocket. Regarding fluctuations of residues Tyr104, Thr107, and Glu129, it was observed that without PVZB1194 binding Tyr104 fluctuates back to the same position as that in AMMPNP-bound Eg5 (Figure D). This is in contrast to the MD simulation with disordered L11 conformation, in which the position of Tyr104 did not change due to hydrogen bonding to Glu270 of L11. From this observation, we propose that L11 plays a role in shifting the Tyr104 position in the process of biphenyl-type inhibitor binding/unbinding. A possible scenario is that when the biphenyl-type inhibitors bind to the α4/α6 allosteric pocket L11 starts to fluctuate, which causes Tyr104 to reorient. This restructuring propagates the conformational changes to the nucleotide-binding pocket and leads to the unbinding of ATP. On the other hand, when ATP binds to the nucleotide-binding pocket, L11 starts to attain an ordered conformation, which shifts Tyr104 back to its normal position and causes the biphenyl-type inhibitors to unbind. In relation to the previous discussion, we concluded that L11 is important for inhibitor binding, and here we furthermore observe that it is necessary for Tyr104 orientation. It is thus possible that when the biphenyl-type inhibitor binds to the α4/α6 allosteric pocket L11 starts to interact with the inhibitor and stabilizes the inhibitor to further propagate the structural changes in the ATP-binding pocket. In both simulations, residues Thr107 and Glu129 remain in the position by which they interfere with nucleotide binding. More enhanced sampling methods might be required to observe the larger conformational changes described above.

Figure 4

(A) AMMPNP-bound Eg5. (B) PVZB1194 bound to the homology model, Eg5, with ordered L11. (C) RMSF of Eg5–PVZB1194 complex and free protein, both with ordered L11 over the 500 ns simulation. The fluctuations of residues 55–65 in the complex are caused by interference from the N-terminal of the protein, located far away from the region of interest. Residues 105–125 fluctuate more in the complex due to the properties of the loop. (D) Tyr104 moves from the PVZB1194-binding position (yellow) to PVZB1194-nonbinding position (blue). AMMPNP-bound Eg5 is colored purple, unbound Eg5 is green.

Pharmacophore Modeling

A set of 66 biphenyl inhibitors with the pIC50 activity range from −1.04 to 2.30 were selected for pharmacophore modeling (Table ).[21,22] The active molecules were defined as the compounds with activity range >1.7, whereas inactive molecules were defined as those with activity range <0. Of the 66 compounds, 31 were classified as inactive, 5 as active, and 30 as intermediate. The defined active and inactive molecules were used to test the specificity of the pharmacophore hypothesis, and the remaining compounds with median activity range (between 0 and 1.7) were used to validate and authenticate the subsequent pharmacophore hypothesis. On the basis of matching by all active ligands, a total of six hypotheses were generated for the biphenyl inhibitors. The inactive ligands were not involved in pharmacophore hypotheses generation but were used to completely remove hypotheses that did not recognize the difference between the active and inactive compounds. The survival, site, vector, volume, and selectivity terms were used to score all generated hypotheses. The outcomes for the six hypotheses are shown in Table ; we can see that the DDHRR.6 hypothesis shows the highest survival (3.785). The DDHRR.6 hypothesis also shows the highest survival minus inactive score, which is the subtraction of the inactive features from the hypothesis.

Table 1

Compound Structures and Their Corresponding IC50 Values

Table 2

Ligand-Based Pharmacophore Hypotheses and Their Scores

ID	survival	survival-inactive	site	vector	volume	selectivity	# matches	energy	activity	inactive
DDHRR.6	3.785	1.854	0.93	0.956	0.898	2.098	5	0	1.745	1.931
DDHRR.5	3.668	1.715	0.89	0.928	0.847	2.096	5	0	1.824	1.953
DDHRR.4	3.631	1.69	0.81	0.952	0.865	2.083	5	0.48	1.745	1.942
DDHRR.3	3.032	1.181	0.6	0.641	0.792	2.092	5	1.736	1.824	1.851
DDHRR.2	2.785	0.972	0.17	0.949	0.662	2.105	5	0.007	1.824	1.813
DDHRR.1	2.676	0.908	0.14	0.905	0.626	2.091	5	0.486	1.824	1.768

The pharmacophore hypothesis DDHRR.6 was mapped on the sets of aligned active and inactive molecules, displayed in Figure A,B, respectively. The components of the DDHRR.6 hypothesis are two hydrogen-bond donors, two aromatic rings, and one hydrophobic feature. The amide groups represent the hydrogen-bond donors, whereas the two phenyl rings on adjacent sides of the compounds represent the two aromatic features, and the hydrophobic feature is mapped onto the trifluoromethyl group. Hence, these groups are crucial for biological activity of this series of compounds and can be further used to screen and design new candidates. By mapping the most active (no. 38) (Figure C) and inactive compound (no. 58) (Figure D) to the pharmacophore features, it can be noticed that although the most inactive compound has four features matched, it has an ethylsulfone instead of the trifluoromethyl substituent, which leads to a dramatic decrease of inhibitory activity, from 5 to 10 904 nM. The reason is that any functional group larger than trifluoromethyl in this position will not fit the pocket of Eg5, which is discussed in detail in the docking section below. Further differences between active and inactive compounds originate from the inactive compounds either lacking a hydrogen-bond donor (Figure E), such as compound 2 with IC50 value of 1000 nM, or having an extra hydrophobic group (Figure F) attached to the phenyl group, such as compound 21 with IC50 value of 6000 nM, in both cases causing less favorable binding to the receptor.

Figure 5

Pharmacophore mapped over (A) the active compounds, (B) the inactive compounds, (C) the most active compound, (D) the most inactive compound, (E) compound 2, (F) compound 21, and (G) structure-based pharmacophore with PVZB1194 placed in the binding pocket.

In addition to the above, a structure-based pharmacophore was generated on the basis of clustering 667 fragments derived from molecules in the medicinal chemistry literature, which were docked into the α4/α6-binding pocket. This eliminates the need for known active inhibitors to generate the pharmacophore models and helps to identify new regions of the active site to target for drug design. By generating a structure-based pharmacophore, it is also possible to further analyze the ligand-based pharmacophore discussed previously. Seven pharmacophore features, AADHHRR, were obtained, that is, one hydrogen-bond donor (D), two ring aromatics (R), two hydrophobes (H), and two hydrogen-bond acceptors (A). By mapping PVZB1194 in the binding pocket to the generated pharmacophore features, it can be seen that the two phenyl groups locate close to the aromatic ring features, the sulfone group to a hydrogen-bond acceptor feature, the trifluoromethyl group to one of the hydrophobic features, and the amide group to a hydrogen-bond donor feature (Figure G). Compared to the ligand-based pharmacophore DDHRR.6, we can see that there is a good overlap of the DHRR features. The extra hydrogen-bond donor in the ligand-based pharmacophore, which was demonstrated to increase the inhibitory activity and incorporates the sulfone group of PVZB1194, overlaps with the extra hydrogen-bond acceptor feature in the structure-based pharmacophore. Taking together, we can constitute a more effective pharmacophore here with ADDHRR. Two more features (A and H) from the structure-based pharmacophore may be possible to explore for targeting a new region of the allosteric site. We note, however, that the hydrogen-bond acceptor (A) feature at Asn289 is relatively distant from the others.

Three-Dimensional-QSAR and Docking

The 66 compounds were randomly divided into a training (70% of the compounds) and test set (30% of the compounds). Using PLS regression, the hypothesis DDHRR.6 with the highest score was explored further in a QSAR study on the basis of the alignment of pharmacophore features to the training and test set compounds. Using only up to four factors of PLS regression, the aligned data set should be well correlated, whereas increase of PLS factors beyond four in general, will not improve model statistics or predictability. In QSAR theory, if the correlation coefficients based on the training set (R2) and test set (Q2) is >0.5, then a satisfactory model is generated. The PLS statistical parameters of the generated 3D-QSAR model are listed in Table . By choosing a PLS factor of 3, the generated QSAR model has values of correlation coefficients R2 = 0.81 and Q2 = 0.64, which implies this QSAR model is acceptable. In addition, the stability of the model is 0.81 (standard stability value >0.3 with a maximum scale of 1) and F value is 60.4 with small P values (1.85 × 10–15); these additional parameters further imply a high degree of validity of the model. The correlation between experimentally observed and predicted activities of the test and training sets are shown in Figure and Table S1, which illustrates that the 3D-QSAR model has good predictive ability.

Table 3

PLS Statistical Parameters of the Selected 3D-QSAR Model

ID	# factors	SD	R2	F	P	stability	RMSE	Q2	Pearson-R
DDHRR.6	3	0.3855	0.8081	60.4	1.85 × 10–15	0.806	0.535	0.6357	0.8509

Figure 6

Fitness graph between the observed (x) and predicted (y) activity for (A) training set and (B) test set compounds.

Contour maps from the 3D-QSAR model were generated to help in understanding the importance of functional groups at specific positions toward biological activity. By analyzing the contour map of the most active compound (Figure ), we can infer that the favorable regions of the hydrogen-bond donor effect are located toward the donor features (D5, D3). In the case of the hydrophobic effect, the favorable regions are seen surrounding the hydrophobic feature (H7), whereas hydrophobically unfavorable areas are identified at the hydrogen bond donor end plus the associated phenyl group or have hydrophobic groups larger than H7. In the case of the electron-withdrawing effect, the favorable regions are seen around the sulfone group and at the position of the phenyl-attached fluorine, indicating that an additional acceptor group in this position could increase the biological activity.

Figure 7

Contour maps for the 3D-QSAR model. (A) Hydrogen-bond donor effect (blue denotes favorable, red denotes unfavorable). (B) Hydrophobic effect (green denotes favorable, purple denotes unfavorable). (C) Electron-withdrawing effect (pink denotes favorable, cyan denotes unfavorable).

On the basis of the free-energy analysis of the MD simulations, an Eg5 structure with the most relaxed L11 conformation was chosen for the docking studies. All 66 ligands were docked into the α4/α6 pocket and the binding modes and interacting residues were analyzed. To simplify presentation of the docking results, we have selected the same four representative ligands as in the above pharmacophore section. The most active compound (entry 38), which has an IC50 value of 5 nM compared to 140 nM for PVZB1194 (entry 12), forms an extra hydrogen bond with Glu345 due to the extra amide group. However, for the most inactive compound, the whole compound shifts due to the large and bulky ethylsulfone group to fit into the bottom of the pocket, whereby one of the amide groups loses the hydrogen bonding to Glu345. For compound 2, no hydrogen bonding was observed, and compound 21 has a flipped docking position due to its side chain being placed in the meta position, which is unfavorable for binding (Figure ).

Figure 8

Eg5 with different compounds docked to the α4/α6 allosteric pocket. (A) Most active compound (38). (B) Most inactive compound (58). (C) Compound 2. (D) Compound 21.

ADME Features

In clinical trials, one of the common failures of drug candidates is poor pharmacokinetics (PKs). Therefore, selecting drug candidates with ideal ADME properties at early stages can help in passing clinical trial studies. With this aim, the relevant descriptors of ADME properties for the tested compounds were examined in this study, among which log Po/w, log S, PCaco, log BB, log MDCK, and PSA were selected. The PK parameters were computed for all 66 compounds (Table ), and although there was a relatively large spread regarding inhibitory activities (from 5 to 10 904 nM), it is still valuable to compare and understand their PK properties. The log  Po/w, which is the octanol/water partition coefficient, is within the range of 1.10–5.70 for all compounds, which suggests that all molecules are ideal. The log S, which is aqueous solubility (S) of the compounds, ranges between −6.43 and −2.49 mol/L. The ability of a compound to interact with the solvent by dipolar or hydrogen-bonding interaction is represented by polar solvent accessible area (PSA) terms. It was observed that compound 22 possesses a PSA of 0 Å, whereas the rest are in the range of 12.54–119.33 Å, which is within the acceptable maximum limit of 120 Å. log BB predicts the brain/blood partition coefficient, and the compounds possess values in the range of −1.82 to 0.51, which again is within the tolerable limit. Furthermore, PCaco indicates the ability to permeate across the gut–blood barrier, and log MDCK is the permeability of the compounds across the Madin–Darby canine kidney cells (nm/s); it was found that all compounds possess ideal PCaco and log MDCK values, which suggests that they can easily be transported across the gut–blood barrier and distal renal epithelia, respectively.

Table 4

Statistical Values of Estimated Physicochemical and PK Parameters along with Recommended Range for Every Propertya,b,c,d,e,f

entry ID	log Po/w (−2.0 to 6.5)	log S (−6.5 to 0.5)	PCaco (<25 poor, >500 great)	log B (−3.0 to 1.2)	PMDCK (<25 poor, >500 great)	PSA (7–200)
1	3.867	–4.22	5470.935	0.205	3105.326	12.535
2	3.843	–3.999	5468.465	0.215	3103.81	12.542
3	4.042	–4.029	5472.442	0.212	3106.251	12.536
4	4.23	–5.499	2063.626	–0.259	1082.502	25.702
5	3.939	–5.657	2064.133	0.083	4862.632	25.696
6	4.079	–5.889	2063.853	0.153	7042.134	25.7
7	4.329	–4.551	2905.022	0.312	10 000	23.127
8	4.037	–4.408	2593.899	0.266	9017.424	26.184
9	4.432	–4.585	3163.218	0.273	10 000	22.569
10	3.789	–4.574	1861.984	0.121	6298.588	36.911
11	3.869	–3.823	3288.081	0.417	10 000	28.123
12	2.524	–4.199	492.638	–0.524	1523.817	64.054
13	3.257	–4.69	1153.88	–0.177	3822.977	50.639
14	3.432	–4.485	1806.434	0.077	6096.523	37.382
15	3.427	–4.89	1516.17	–0.082	5060.716	49.554
16	2.918	–4.483	339.041	–0.77	1010.486	76.383
17	3.259	–4.518	480.75	–0.656	1495.496	72.669
18	2.395	–4.191	264.81	–0.872	783.033	85.659
19	2.22	–3.953	252.681	–0.958	503.584	87.165
20	3.14	–4.549	1084.251	–0.27	2479.281	50.791
21	3.212	–4.705	1215.982	–0.247	2785.946	49.691
22	5.699	–6.433	9906.038	0.192	10 000	0
23	4.591	–4.928	5790.938	0.509	10 000	14.406
24	3.581	–4.878	1390.756	–0.246	3214.652	49.936
25	3.957	–4.276	2573.515	0.189	6151.143	26.311
26	4.266	–4.426	3018.842	0.255	7310.738	22.469
27	2.403	–4.034	488.01	–0.6	1040.702	64.232
28	2.466	–2.921	239.589	–0.474	1203.899	64.406
29	3.862	–5.109	2158.267	0.269	10 000	43.453
30	2.989	–4.694	307.218	–0.944	618.525	81.062
31	2.811	–4.318	328.742	–0.846	684.135	77.294
32	4.322	–5.433	1323.133	–0.244	3083.543	51.981
33	3.593	–4.987	665.065	–0.557	1472.304	66.858
34	3.277	–5.005	571.03	–0.76	1226.418	63.962
35	2.517	–4.032	362.328	–0.8	765.954	73.35
36	2.994	–4.335	769.759	–0.325	1663.525	64.682
37	2.982	–4.593	679.624	–0.3	1818.741	59.53
38	2.989	–4.317	416.324	–0.644	1198.532	74.972
39	2.279	–5.089	166.869	–1.197	329.214	92.091
40	3.039	–4.414	455.018	–0.685	973.676	70.544
41	3.047	–4.514	442.737	–0.788	943.483	81.226
42	2.207	–3.915	172.939	–1.188	341.753	93.738
43	1.985	–3.806	143.303	–1.27	278.826	98.271
44	2.74	–4.413	349.991	–0.913	732.3	84.288
45	3.271	–4.853	394.628	–0.643	1685.705	75.292
46	2.977	–4.567	325.443	–0.774	1033.237	77.318
47	3.178	–4.797	369.846	–0.648	1770.085	74.554
48	2.669	–3.851	1097.759	–0.36	1375.596	50.631
49	3.297	–4.55	1089.11	–0.326	2560.139	58.832
50	1.566	–3.781	201.506	–1.141	350.439	89.605
51	2.454	–2.936	239.938	–0.483	1197.056	64.397
52	1.765	–2.664	239.351	–0.752	268.811	64.421
53	2.409	–2.824	239.555	–0.842	268.957	64.417
54	2.632	–2.99	239.555	–0.833	268.957	64.417
55	2.391	–2.813	269.537	–0.366	1785.925	62.671
56	2.978	–3.528	242.573	–0.488	1926.466	64.281
57	1.142	–2.804	47.829	–1.328	188.732	101.457
58	1.104	–2.493	72.953	–1.428	102.672	97.017
59	3.393	–5.584	302.87	–1.004	906.624	77.43
60	2.678	–5.139	67.958	–1.822	100.001	119.337
61	1.894	–3.839	112.76	–1.409	199.216	99.443
62	1.472	–2.698	323.178	–1.005	151.069	96.299
63	1.874	–3.394	323.178	–0.822	502.438	96.265
64	2.226	–3.923	323.178	–1.141	151.069	94.172
65	2.652	–4.374	323.462	–0.74	945.005	94.921
66	2.53	–4.544	225.988	–0.942	648.421	95.516

log Po/w: predicted octanol/water partition coefficient.

log S: conformation-independent predicted aqueous solubility, log S. S in mol dm–3 is the concentration of the solute in a saturated solution that is in equilibrium with the crystalline solid.

PCaco: predicted apparent Caco-2 cell permeability in nm/s. Caco-2 cells are a model for the gut–blood barrier. QikProp predictions are for nonactive transport.

log BB: predicted brain/blood partition coefficient. Note: QikProp predictions are for orally delivered drugs so, for example, dopamine and serotonin are CNS negative because they are too polar to cross the blood–brain barrier.

PMDCK: predicted apparent MDCK cell permeability in nm/s. MDCK cells are considered to be a good mimic for the blood–brain barrier. QikProp predictions are for nonactive transport.

PSA: Van der Waals surface area of polar nitrogen and oxygen atoms.

Conclusions

MD simulations have been carried out on the Eg5 protein and Eg5–PVZB1194 complex, with ordered or disordered loop L11, to better understand key features of the α4/α6 allosteric pocket. With L11 in a disordered state, the complex structure shows that there are three residues hydrogen bonded to the inhibitor PVZB1194, one of which (Asn287) comes from L11. L11 becomes much more stable and moves toward the pocket and preferably resides between α4 and α6 with inhibitor present, whereas without inhibitor L11 is more flexible and moves outward from the pocket and prefers to stay closer to α4. Residue Tyr104, which plays an important role in generating the ATP-competitive conformation, did not display any large movements during the simulation, even without the inhibitor binding to the pocket due to hydrogen bonding to residue Glu270. Even with an ordered L11 conformation, as observed in the AMMPNP-bound Eg5 structure, PVZB1194 was able to fit into the pocket and cause reduced L11 fluctuation. A ligand-based pharmacophore and atom-based 3D-QSAR study in combination with molecular docking was carried out to identify the key pharmacophore features for obtaining biological activity of biphenyl-based ligands toward the α4/α6 allosteric pocket. The five most active compounds share common pharmacophore features, including two hydrogen-bond donors, two aromatic rings, and one hydrophobic group. Analysis of the 3D-QSAR studies provided further information on how to increase the biological activity of the current compounds. Docking the 66 biphenylic compounds into the allosteric pocket demonstrates the importance of the sulfone amide group in interaction with the receptor. A structure-based pharmacophore was also generated on the basis of extensive fragment docking. The two pharmacophore models show overlapping DHRR features, and two additional features (A, H) were identified that may be explored for additional binding to the active site. From the bioavailability and ADME point of view, essentially all 66 biphenyl compounds would be capable of passing the clinical trial phases. The study herein provides a set of guidelines that will assist in understanding the biphenyl-type inhibitor interactions with Eg5 and in designing new and more potent leads for targeting this specific allosteric pocket.

Methods

MD Simulations and Homology Model

The Eg5–PVZB1194 complex structure was retrieved from the protein databank (PDB ID: 3WPN).[26] The default protein preparation wizard in Schrödinger was used to prepare and refine the structure in the same manner as in our previous studies.[28,29] First, the structure was checked to guarantee the chemical correctness. Hydrogens were added to the structure, and the side chains that were far from the binding pocket and not forming any salt bridges with other side chains were neutralized. Hydroxyl groups, water molecules, and amide groups of Asn and Gln were reoriented to optimize the hydrogen-bonding network. Because no His residues were included in the binding pocket, the appropriate protonation states and orientations of the imidazole rings were automatically determined by the preparation wizard. ACE and NMA residues were added to cap the termini, and missing side chains and missing loops with less than 20 residues were added. In the final step, the OPLS-2005 force field[30] was used to minimize the structure with the heavy atoms constrained. The prepared structure was then used for subsequent homology modeling and MD simulations.

To generate an Eg5 structure with L11 in an ordered state, we carried out homology modeling based on two selected templates, the above prepared Eg5–PVZB1194 structure and the Eg5 bound to the substrate analogue AMMPNP (PDB ID: 3HQD),[31] using Prime (Schrödinger, LLC, NY).[32] The resulting structure was constructed by superposing the two templates and using defined residues from L11 of the AMMPNP-bound Eg5 structure together with the remainder from the above prepared Eg5–PVZB1194 structure.

The molecule PVZB1194 was first minimized at the AM1 semiempirical level,[33] after which the restrained electrostatic potential (RESP) procedure at the HF/6-31G* level was used to calculate the atomic charges. Using the Antechamber module in the AMBER10 program,[34] the GAFF force field parameters[35] and RESP partial charges were assigned to the molecule. The LEaP module was used to add missing hydrogen atoms.

Four conventional 500 ns MD simulations were carried out on the Eg5–PVZB1194 complex structure and Eg5 protein structure only, in which the Amber 99 force field[36,37] and GROMACS software[38] were used. The whole workflow of MD simulations was the same as that in those carried out in our previous studies.[28,29,39] A buffer distance of 10.0 Å was used to solvate the structures with TIP3P water in periodic boxes. To adjust the electroneutrality condition of the systems, 0.1 mol/L of Na+ and Cl– ions were added. Energy minimization with 200 steps of the steepest descent was carried out to remove close contacts in the obtained systems. Following that, to make sure the water molecules would reach more favorable positions, a 2 ns position-restrained simulation with a constant pressure ensemble (NPT) was performed, with a time step of 1 fs, temperature of 298 K, and coupling pressure of 1 bar. For long-range electrostatics, particle-mesh Ewald[40,41] summation was used. The cutoff of Coulomb and Lennard-Jones interactions was set to be 10 Å. The Berendsen coupling algorithm[42] was used to control the temperature and pressure, with the time constants 0.1 ps for temperature and 1.0 ps for pressure coupling. The LINCS algorithm[43] was used to constrain all bond lengths. Following the position-restrained simulation, 500 ns production simulations with NPT ensemble were performed on each system for evaluation and analysis. In this step, the Nosé–Hoover thermostat,[44] with a time constant 0.1 ps, was used to control the temperature and the Parrinello–Rahman barostat,[45] with a time constant 1.0 ps, was used to control the pressure. The other parameters remained the same as those in the position-restrained simulations.

Docking and Structure-Based Pharmacophore Generation

The resulting relaxed Eg5–PVZB1194 complex structure was attained from the final snapshot of the 500 ns MD simulations. On the basis of the obtained structure, the Receptor Grid Generation panel (Schrödinger, LLC, NY) was used to generate grids to specify the position and size of the active site. Grids were defined by centering them on the bound PVZB1194 using the default box size. Both the 66 recently published biphenyl compounds, which were tested for their inhibition effects for Eg5[21,22] (Table ) and a set of 667 fragments from the Schrödinger fragment library were docked into the α4/α6 allosteric site. To grant full flexibility to the ligands, the XP (extra precision)[46] docking function of Glide[47] was used. To select the best docked poses, a post-docking minimization was carried out on the output complexes and the number of docking poses per ligand was reduced from 10 000 to 5. The local optimization feature in Prime was then used to minimize the resulting poses.

Fragment clustering was used for the complex structure of Eg5 with the 667 docked fragments to generate a structure-based pharmacophore. The methodology has been described in detail in the previous works by us and other groups.[28,48] In brief, the pharmacophore features were generated on the basis of the XP descriptor information. Phase v3.0 (Schrödinger, LLC, New York, NY) was used to automatically generate pharmacophore sites on the basis of six possible chemical features, namely hydrogen-bond acceptor (A), hydrogen-bond donor (D), hydrophobe (H), negative ionizable (N), positive ionizable (P), and aromatic ring (R). Quantification and ranking of each feature was followed in the final step and was on the basis of the energy value, which equals the sum of the Glide XP contributions.

The pharmacophore and 3D-QSAR models were generated on the basis of 66 biphenyl compounds.[21,22] Their structures and corresponding IC50 values in nM are shown in Table and Figure S1. Maestro was used to build the molecular structures, and the “Pharmacophore Modeling” module of Phase was used to generate pharmacophore models.[49,50] Before pharmacophore generation, the IC50 value of each compound expressed in μM was converted to the corresponding negative logarithm (pIC50), which were all in the range −1.04 to 2.30; for each molecule, multiple conformers were generated and optimized by performing energy minimization using the OPLS-2005[30] force. In the next step, the ConfGen module[51] was used to explore the conformational space using 100 conformers per rotatable bond and a maximum of 1000 conformers per structure. After the above steps, the five most active compounds with pIC50 >1.7 were used to generate pharmacophore models. The features such as hydrogen-bond acceptor and donor, hydrophobic, negative, positive, and aromatic rings were employed. In our study, pharmacophores with a maximum of five features that match all active ligands were generated and were further scored using default parameters.

Three-Dimensional-QSAR Modeling

To generate the 3D-QSAR models, we can either choose an atom- or a pharmacophore-based approach. Here we choose the atom-based 3D-QSAR model as it includes the entire molecular space of the compounds, whereas in the pharmacophore-based 3D-QSAR the regions beyond the pharmacophore model are not considered.[49,50] In this study, all 66 compounds in the data set were used to develop a 3D-QSAR model. Those compounds were randomly divided into a training set and a test set, which corresponds to 70 and 30% of the total compounds, respectively, using an automated random selection tool in the Phase module. Three maximum partial least square factors in the regression model and 1 Å length of the sides of cubic volume elements were used to generate the 3D-QSAR models. The detailed method for generating pharmacophore and 3D-QSAR models in the Phase module has been outlined previously.[49,50]

Drug Likeness and ADMET Analysis of Biphenyl-Type Compounds

The QikProp module in Schrödinger[52,53] is widely used to calculate ADME/T properties of various compounds and to predict PK and physicochemical parameters that account for “drugability” of a molecule in the drug selection process, especially in preclinical trials. In this study, we calculated the ADME/T parameters of the biphenyl-type inhibitors to fully understand their molecular properties.