Abstract

Euphorbia factor L713283 is a new lathyrane diterpene isolated from Euphorbia lathyris and shows strong anticancer activity. By using molecular similarity analysis, β-tubulin was identified as one of the possible targets of L713283. We further investigated the binding modes of L713283 with β-tubulin using molecular docking and molecular dynamics (MD) simulation methods. The results indicated that the binding site between β-tubulin and L713283 was composed of the four regions, that is, residues Phe20~Glu27, Leu225~Thr232, Phe270~Gly277, and Ile356~Met363. MM/GBSA method was used to calculate the binding free energy and determine the key residues for the association of L713283 with β-tubulin. It was found that nonpolar interactions made the major contributions for the binding. In addition, we compared the binding pocket and motion modes of L713283-free and L713283-bound β-tubulin systems. It is proposed that L713283 may bind to β-tubulin and favor the formation of αβ-tubulin dimmer. This work provides possible explanation for molecular mechanism of the anticancer agent L713283, and the strategy used here could benefit the investigation of possible target profile for those bioactive agents with unknown mechanisms.

1. Introduction

Cancer is the second leading cause of death worldwide and has a major impact on global health [1, 2]. Although significant scientific advances in recent years have been made to improve the treatments against cancer, the rapid emergence of drug resistance is the major obstacle for current therapy. Prevalent targets [313] have been identified, such as PTK, PKC, FTase, MEK, TopoI, TopoII, MMP, HSP90, tubulin, and TLMA. Their corresponding representative drugs include STI-571, UCN-01, Zarnestra, PD184352, Camptothecin, Salvicine, RS-130830, Geldanamycin, Taxol, and PIPER, respectively. In the past few years, it is a research hotspot to extract new anticancer drugs originated from plant, which helps to attenuate drug resistance in cancer treatment [14, 15]. Euphorbia lathyris is widely distributed in China and has been used as a traditional folk medicine for the treatment of cancer and warts [16, 17]. Many chemical studies on Euphorbia species were focused on the components of skin-irritant compounds, mainly consisting of ingenane, tigliane, or daphnane skeletons [18]. Among these compounds, a series of diterpenoids exhibiting extremely strong skin irritation have been isolated, which possess tumor-prohibiting effects and show biologic activities such as anticancer and anti-HIV [1921].

Recently, Jiao et al. extracted a new diterpenoid molecule Euphorbia factor L713283 from Euphorbia lathyris [22], with IC50 value of μg/mL for tumor cell inhibition (see Figure 1 and Figure S1 in Supplementary Material available online at http://dx.doi.org/10.1155/2015/879238. For convenience, the atomic sequence number is provided in Figure S1 and used in the following analysis). As known, drug resistance is still the major cause for the mortality of cancer. Thus, there is an urgent need to develop new effective therapies to improve survival rates. The inhibitory activity of L713283 is slightly weak as compared with other anticancer drugs, but it possesses a new chemical scaffold of isolathyrol skeleton characterized by a 5 : 11 : 3 fused ring system [22]. Some diterpenoids similar to L713283 were isolated by Duarte and Ferreira’s group, and their possible biosynthetic pathways were proposed [23]. Duarte et al. further used in vivo and 2D NMR methods to explore the multidrug resistance induced by macrocyclic lathyrane-type diterpenoids [24]. Lu et al. obtained the NMR data of cytotoxic diterpenoids from Euphorbia helioscopia, which are all jatrophane-type diterpenoids [25]. However, the possible targets and molecular mechanism of this novel anticancer agent remain ambiguous.

To interpret anticancer mechanism of L713283 in the present paper, we employed integrated in silico methods including 3D similarity search, structural based docking, and molecular dynamics simulation to explore the possible targets of the compound. These studies show that β-tubulin protein may be the potential target of L713283, providing new insights into the molecular architecture of β-tubulin with L713283.

2. Computational Methods

2.1. Computational Protocol

The crystal structure of L713283 was downloaded from Cambridge Crystallographic Data Centre (CCDC) [26], and then molecular similarity analysis was performed to identify the possible acting target. Furthermore, the binding mode of L713283 with its target protein was obtained by molecular docking experiment. Finally, the possible anticancer mechanism of L713283 was proposed through binding free energy, key residues, and motion mode analysis based on the comparative molecular dynamics (MD) simulation trajectories. Figure 2 shows the protocol on the explorations for the target protein of L713283 and its possible anticancer mechanism.

2.2. Molecular Similarity Analysis

Molecular similarity principle, which is widely used in drug design, states that molecules with similar structure tend to have similar properties [27, 28]. Here, ten widely used representative anticancer molecules [313] acting on different drug targets were chosen for the similarity analysis with L713283. The structures of all the small molecules were created with Chem3D package [29], followed by the energy minimization with Antechamber suite [30]. ShaEP package [31] was employed for the molecular similarity calculation. In detail, L713283 was superimposed on each anticancer molecule using a matching algorithm [31]. Two characteristic scores (i.e., 3D-shape and electrostatic potential) were calculated for comparison. The score range is from 0 to 1, in which 0 and 1 correspond to no similarity and the same molecules, respectively. Finally, the possible anticancer target of L713283 will be predicted through molecular similarity analysis.

2.3. Molecular Docking

The structure of β-tubulin (i.e., the possible target protein of L713283) was built from the X-ray crystal complex structure between αβ-tubulin and its inhibitor Taxol (PDB code: 1JFF) [32]. We only kept chain B of 1JFF from Arg2 to Asp427 (i.e., β-tubulin subunit and total of 426 residues), while the ligand Taxol, ADP, crystal waters, and chain A were deleted. Then, the structure was minimized for 5000 steps with the steepest descent method followed by 5000 steps with the conjugate gradient method [33]. The ultimate minimized system was used as the structure of docking receptor. In addition, the structure of L713283 was considered as the docking ligand after minimization using semiempirical quantum chemistry methods with Antechamber suite [30].

The molecular docking experiment of L713283 with β-tubulin was performed with AutoDock 4.0 package [34]. The program conducts automated docking of a fully flexible ligand to a rigid receptor. The rotations of all the single bonds of the ligand are taken into account during the docking calculation and a total of 17 flexible torsions are defined in this work. The Lamarckian Genetic Algorithm (LGA) [35] is used to search the globally appropriate conformations and the possible binding sites of the ligand, with the semiempirical potential function used as the energy scoring function. Before docking L713283 to β-tubulin, we first redocked Taxol to β-tubulin and compared the docked conformations with the crystal structure. The conformation with lowest potential energy was superimposed to the crystal structure of 1JFF. The superimposition results show that both the binding pocket and molecular pose of the redocked complex model are highly similar to those from 1JFF crystal structure. As shown in Figure S2, the corresponding RMSD of ligand is less than 0.2 nm. This result shows that the docking protocol is suitable for docking the β-tubulin system.

Then, the L713283 was docked to β-tubulin. The coordinate of the grid box center was set to [5.0, −22.0, 14.0], which is approximately the mass center of Taxol. The box grid number was defined as 120 × 120 × 120 with grid spacing 0.0375 nm in each dimension, and the final box size was 4.5 × 4.5 × 4.5 nm3. Other unmentioned parameters of molecular docking were set as default values. Each docking calculation produced 128 complex structures. Clustering analysis was performed on 128 ligand conformations to obtain the representative docking poses. A RMSD cutoff value was set as the threshold of cluster. The RMSD values between the first conformation and the other conformations (i.e., ) were calculated. The conformation is put into the first cluster if the RMSD is less than the threshold value; otherwise, it becomes the conformation in the second cluster. The rest of conformations are then compared to conformations in the first and second clusters. These steps are repeated until all the conformations are fallen into a certain cluster. Finally, the conformation with the lowest potential energy is determined as the representative docked complex structure in the cluster [35, 36]. The analyses of both energy and cluster were performed to gain the final docking results, denoted by complex- with the lowest docking energy in the top cluster. The poses located near the Taxol’s binding site with relative low docking energy and large number in clusters were selected for further investigation in molecular dynamics simulation.

2.4. Molecular Dynamics Simulation

All-atom MD simulation is an important tool in exploring the thermodynamics and kinetics features of biomolecules [3740]. In this work, three independent all-atom MD simulations were carried out for the β-tubulin and complex systems using the AMBER 10 suite of programs [41]. The Amber ff99 all-atom force field was adopted [4244]. In each simulation, the solute was solvated in a truncated periodic box with a 1.2 nm solute-wall distance using the TIP3P water model. Finally, the three simulation systems (complex1, complex2, and β-tubulin) contained 10256, 10256, and 9955 water molecules with a total of 37449, 37449, and 36464 atoms, respectively. Before MD simulations, two stages of minimization were performed for the initial structures. First, the solute atoms were restrained with the force constant 2.09 × 105 kJ·mol−1·nm−2. Water molecules were minimized for 5000 steps with the steepest descent method followed by 5000 steps with the conjugate gradient method. Then, the restraints for the solute were removed. The whole system carried out 5000 steps of steepest descent minimization followed by 5000 steps of the conjugate gradient minimization. The convergence criterion was 4.182 kJ·mol−1·nm−2.

MD simulations were also performed by two stages. First, the solutes were restrained (the restraining force constant was 4.182 × 103 kJ·mol−1·nm−2) and slowly heated up from 0 to 300 K over 1.0 ns. Then, nonrestraint MD simulations at 300 K were performed for 20.0 ns. The SHAKE algorithm [45] was applied to constrain the covalent bonds involving hydrogen atoms and a 1.2 nm cutoff was used for all nonbonded interactions. The MD time step was set as 2 fs, and one snapshot was sampled every 2500 steps. Thus, total 4200 conformations were obtained during each 21.0 ns MD simulation. In these simulations, the trajectories were monitored with the VMD software [46].

2.5. MM/GBSA Methodology and Energy Decomposition

The MM/GBSA method [47, 48] was used to calculate the binding free energy between β-tubulin and L713283. 50 snapshots were extracted from the trajectory of 11.0 to 21.0 ns every 0.2 ns in order to obtain the mean binding free energy. The binding free energy is defined aswhere is the molecular mechanical energies in vacuo, is the solvation free energy, is the temperature, and is the solute entropy. The solute entropy is usually estimated by normal mode analysis method [49].

can be obtained by average of the molecular mechanics energy from the collected snapshots during MD simulation. It is defined aswhere includes the bond, angle, and dihedral energies, is electrostatic interaction energy, and is van der Waals interaction energy.

The contribution from the solvent is calculated aswhere is the sum of electrostatic salvation energy (polar contribution) and is nonelectrostatic salvation component (nonpolar contribution). The polar contribution is calculated using Generalized Born (GB) approximation model [5052]. The nonpolar contribution is estimated by solvent accessible surface area (SASA) [53], and it is calculated with the expressionwhere and are the default constants. Then, (1) can be written asEach energy item is calculated separately by using the conformations extracted from MD trajectory and the binding free energy of forming complex is finally obtained by subtracting the energy terms for complex by those for protein and ligand. The essential idea of energy decomposition is to decompose the energy contribution of each residue to the association of the receptor and the ligand. By using energy decomposition, we could identify the key residues and their contributions for the binding.

3. Results and Discussion

3.1. Molecular Similarity

Based on the literature research, we selected the representative anticancer drugs with different targets. First, the mechanisms of actions of these drugs are clarified and all the targets are validated. Second, the structural information of the targets is available from the Protein Data Bank. These selected anticancer drugs cover most of widely researched cancers [14, 54]. Starting from the knowledge of ten small molecules identified as anticancer drugs [313], we analyzed the molecular similarity between L713283 and these known drugs to understand the possible anticancer mechanisms of L713283. In particular, two important properties, 3D-shape and electrostatic potential (ESP) of L713283, were compared to those of ten drugs. Figure 3 shows the shape similarity and ESP similarity of L713283 to the anticancer drugs. The high value (0.749) of shape similarity between L713283 and Taxol indicates that L713283 shares similar 3D-shape with this inhibitor, and both molecules possess a similar framework of multirings (i.e., 5 : 11 : 3 and 6 : 8 : 6 : 4). L713283 also has a high shape similarity with Geldanamycin with a value of 0.620. The main similar factor is still the multiring framework. It is also shown that L713283 has the relatively high ESP similarity (over 0.7) with most of inhibitors. The high ESP similarity may be due to the existence of plentiful electron-withdrawing functional group close to molecular ring framework, which is similar to that of L713283. For example, the value of ESP similarity (0.769) between L713283 and PIPER is the highest because of the positively charged piperidine group. By combining the two properties, the results show that Taxol has the highest molecular similarity to L713283 with an average similarity of 0.736. Taxol (Paclitaxel) is an important anticancer botanical drug, which blocks the life cycle of tumor cell by inducing the bundling of microtubules in the interphase cells through targeting β-tubulin [5557]. Its anticancer biological activity is exerted by binding to β-tubulin and prohibiting the disassembly of αβ-tubulin dimer. The high molecular similarity of L713283 with Taxol suggests that β-tubulin could also be the target for L713283.

Notably, some previous biological experimental results also indicated that β-tubulin might be the target of L713283. Khaleghian et al. extracted the inganen from Euphorbia tirucalli (Euphorbiaceae family) and found the conformational change of tubulin at the presence of inganen [58]. Both fluorescent spectroscopy and ultraviolet spectroscopy results showed that the inganen molecule may decrease polymerization of tubulin. Miglietta et al. isolated jatrophanes from Euphorbia semiperfoliata (Euphorbiaceae family) [59]. They investigated the interactions between jatrophanes and purified bovine brain tubulin by an in vitro polymerization assay and electron microscopy. As new types of tubulin binders, both inganen and jatrophanes are diterpenoid molecules obtained from Euphorbiaceae family and have highly similar structures with Euphorbia factor L713283.

3.2. Molecular Docking

Figure 4 shows the results of docking experiments of L713283 with β-tubulin. The docked conformations were divided into 11 clusters with the root mean standard deviation (RMSD) threshold of 0.2 nm. The maximum cluster (i.e., cluster 1) has 54 structures accounting for 42.2% of total conformations, whose representative conformational ID is 120. The corresponding minimum binding free energy for this cluster is −7.88 kJ/mol. The second largest cluster (i.e., cluster 2) owns 36 structures accounting for 28.1% of all the conformations, and 45 is the ID of its representative conformation. The minimum binding free energy of this cluster is −7.44 kJ/mol. Both the third (i.e., cluster 3) and fourth (i.e., cluster 4) largest clusters have 12 structures, and the representative conformational IDs are 36 and 11, respectively. The rest seven clusters have very few structures, and then they are shown in Figure S3.

Figure 4(b) shows the representative binding modes of the top four largest clusters and the crystal structure of Taxol ligand. Both the binding sites for clusters 3 and 4 are the flat pocket composed of the two random coils of Gln280~Thr285 and Arg359~Met363. As for clusters 1 and 2, the binding pocket is relatively deep comprising the helix of Asp224~Ala231 and two random coils (i.e., Phe270~Gly277 and Arg359~ Met363). Obviously, the binding sites of clusters 1 and 2 are the same with those regions of β-tubulin recognized by Taxol. We further analyzed the contact residues of β-tubulin for clusters 1 and 2, which mainly included Leu215, Gly223, His227, Pro272, Leu273, Thr274, Arg276, Arg282, Gly360, and Leu361. Here, the residues of the protein in the neighborhood of 0.4 nm of L713283 inhibitors are defined as the contact residues, which contributes to the specific recognition of inhibitors.

The representative conformations in clusters 1 and 2 were selected as the final binding poses and, respectively, named as complex1 (conformational ID = 120) and complex2 (conformational ID = 45). The criteria for binding poses selection are described as follows: (1) the binding site of both complex1 and complex2 structures aligned well with that in the crystal structure of Taxol : β-tubulin [32]; (2) the conformations of clusters 1 and 2 account for over 70% of total conformations and have relative low binding free energy. Thus, the structures of complex1 and complex2 were chosen to perform the following comparative MD simulations.

3.3. Convergence Behavior of Three MD Simulations

By monitoring the evolution of potential energies of the β-tubulin, complex1, and complex2 systems versus time, it is found that the potential energies of the three systems reach equilibrium after 2.0 ns. The mean values of the potential energies for the β-tubulin, complex1, and complex2 systems are −4.41 × 105, −4.53 × 105, and −4.52 × 105 kJ/mol, respectively. Their corresponding standard deviation values are 4161, 4019, and 4090 kJ/mol, which account for the total values less than 1%.

Figures 5(a) and 5(b) show the RMSD evolution of Cα atoms and L713283 for the three MD simulation systems. The RMSD values of the β-tubulin, complex1, and complex2 systems have been calculated by fitting the Cα atoms to the initial minimized structures and converge around 0.198, 0.182, and 0.212 nm after 11.0 ns, respectively. Obviously, the complex1 system is the most stable. As shown in Figure 5(b), the RMSD values for L713283 in complex1 and complex2 systems remain at and  nm, respectively. Roughly, the mobility of L713283 in complex1 is lower than that in complex2. Furthermore, the RMSD values of L713283 in complex2 rise continuously and do not show convergence. The electron density map analysis for the Taxol, β-tubulin crystal structure [60, 61], shows that the bounded Taxol conformation is very similar to that determined independently by energy-based refinement. Considering the systemic stability in MD simulation and the electron crystallographic density result, it is speculated that the structure of complex1 is most likely the final binding mode of L713283 with β-tubulin.

Root mean square fluctuation (RMSF) can be used to describe the conformational change versus the mean structure of the system during a MD simulation. The bigger the value of RMSF is, the higher the flexibility of the system is. Figure 5(c) shows RMSF values in the complex1 and the β-tubulin alone systems. The RMSF distribution of the β-tubulin is similar to that of its complex with L713283 (i.e., complex1). The protein part exhibits lower flexibility after binding with L713283, which is consistent with the result gained from RMSD analysis (see Figure 5(a)). Figure 5(d) gives the correlation result between the RMSD values of β-tubulin and that of complex1. The correlation coefficient is 0.76, which is very high as for a system of 426 residues. The high correlation of RMSF distribution indicates the dependability of the trajectories gained from the two MD simulations.

In order to visualize the flexibility of each residue, high (RMSF > 0.25 nm) and low (RMSF < 0.2 nm) fluctuation regions of β-tubulin are, respectively, colored in red and blue according to the RMSF values in Figure 5(c) of complex1. The rest of the region (0.2 nm < RMSF < 0.25 nm) is defined as medium fluctuation regions and depicted in gray (see Figure 6(b)). As shown in Figure 6(b), the highly flexible regions are mainly composed of seven external random coils, including Ala55~Asp57, Ser95~Ala97, Pro173~Glu181, Lys216~Tyr222, Gly277~Ala283, Val342~Asp348, and Phe384~Ser413. The low flexible regions are composed of two internal α-helices (i.e., Gly17~Ile30 and Ser230~Thr238) and three β-sheets (i.e., Arg2~Ile7, Ile47~Ala54, and Arg58~Ile64). Referring to Figure 5(c), the β-tubulin alone system also exhibited similar flexible distribution.

The anticancer drug Taxol stabilizes αβ-tubulin and reduces its dynamicity, promoting mitotic arrest and cell death [62]. From the above comparative RMSD and RMSF analysis of complex1 and β-tubulin systems, the association with L713283 may stabilize β-tubulin and result in the obvious decline of these two parameters. From this point of view, the influence of the association with L713283 on β-tubulin is similar to that of Taxol ligand.

3.4. Binding Free Energy Calculation between L713283 and β-Tubulin

By combining the free energy calculation and decomposition analysis, important information was obtained, including key residues, driving force, and the charged environment of binding pocket during the association of receptor-ligand. Table 1 lists the energy components contributing to the binding free energy between L713283 and β-tubulin. The binding free energy calculated by MM/GBSA method is −40.8 kJ·mol−1. By analyzing the nonpolar interaction (VDW + GBSUR) and the polar interaction values (ELE + GB), it is found that the associations between L713283 and β-tubulin are mainly driven by the favorable nonpolar interactions. In addition, both the polar interactions (ELE + GB) and the entropy effect (TSTOT) are unfavorable for the association.

In order to investigate the detailed interactions between L713283 and β-tubulin, we used the energy decomposition strategy in MM/GBSA package to analyze the contribution of each key residue to the binding of L713283 with β-tubulin (see Figure 7). The key residues mainly locate in the following four regions: (1) residues Phe20~Glu27, such as Val23 (−8.6 kJ/mol); (2) residues Leu225~Thr232, such as His227, Leu228, and Ala231 (−10.4, −3.4, and −5.9 kJ/mol); (3) residues Phe270~Gly277, such as Phe270, Pro272, Leu273, and Arg276 (−6.6, −6.8, −6.7, and −10.0 kJ/mol); and (4) residues Ile356~Met363, such as Pro358 and Leu361 (−1.3 and −7.2 kJ/mol). Nevertheless, the electrostatic interactions of Asp26 (0.8 kJ/mol, region 1), Glu27 (0.9 kJ/mol, region 1), Asp224 (1.0 kJ/mol, region 2), and Ser275 (2.8 kJ/mol, region 3) are unfavorable for the association of L713283 with β-tubulin. On the whole, these unfavorable factors were balanced out and reversed by the favorable interactions contributed by the key residues nearby.

Subsequently, we analyzed the hydrogen bonds between L713283 and β-tubulin. Here, the formation of hydrogen bond is judged using geometric criterion with the distance of donor-acceptor less than 0.35 nm and the angle of donor-hydrogen-acceptor over 135°. It is found that there is only one stable hydrogen bond between L713283-O34 and Arg276-NH2-1HH2 (see Figure 8(a) and Figure S4). Based on the average structure of complex1 from 11.0 to 21.0 ns, Figure 8 shows the refined binding mode of L713283 with β-tubulin. Among the 32 residues located near the binding pocket, there are 25 nonpolar residues. Figure 8(b) shows the electrostatic potential distribution of the active pocket of β-tubulin. From Figures 8(a) and 8(b), it is obvious that the binding site is a nonpolar pocket. In order to clearly observe the structure and position of L713283 in the pocket, the whole molecule in Figure 8(b) was carried out in a simple rotation operation compared with Figure 8(a). It is shown that the 3-ring possessing more nonpolar character is located slightly away from the molecule, which is mentioned in previous studies [10, 63].

In summary, the calculated binding free energy value between L713283 and β-tubulin via MM/GBSA method is −40.8 kJ·mol−1. Both systemic van der Waals interactions and nonpolar solvent effect are favorable for the formation of L713283 : β-tubulin complex. In addition, the results of energy decomposition and key residues show that the binding site is a nonpolar pocket composed of four nonpolar residue segments.

3.5. Motion Modes of L713283

Table 2 lists the dihedral angle distribution of large ring, middle ring, and side chains for L713283 molecule. All the dihedrals of L713283 keep relatively stable. In particular, the fluctuation scope for the dihedral angles of middle ring is within 3%. In the side chain, C19 is the most stable because of the restrictive hydrogen bond. Unexpectedly, the aromatic ring linked by C23 and the acyl groups linked by C13 and C16 also exhibit 5% or less fluctuation scope. For the large ring, the fluctuation scope of the three dihedral angles (i.e., C19-C1-C2-C3, C1-C2-C3-C4, and C2-C3-C4-C5) is more than 7.5%, while that of the other eight angles remains at 5% or less. The last column of Table 2 lists the initial dihedral parameters of L713283 optimized by semiempirical quantum chemical methods. It is obvious that the optimized values are basically the same as those from MD simulation.

In summary, the dihedral angles of L713283 remain relatively stable during MD simulation, and molecular motion exhibits the integrity characteristic. Obviously, there is only slight change in the conformation of L713283 after its association with β-tubulin. It is consistent to the hydrophobic property and electrostatic environment of the binding pocket. The results from dihedral angle analysis agree well with the lower RMSD values (see Figure 5(b)). Finally, it is found that the C19-C1-C2-C3-C5 chain in the large ring exhibits significant torsion, which mainly results from the initial torsion of C4-C5-C6 ring by monitoring the trajectory of L713283 in complex1 system.

3.6. Comparative Motion Mode Analysis of β-Tubulin

The motion modes of biomacromolecules play key roles in exerting specific biological functions including binding with their ligand [64, 65]. To investigate the possible inhibitory mechanism of L713283, we comparatively analyzed the motion modes of the β-tubulin alone and L713283-bound complex1 systems. The motion directions are based on the average structures from trajectories of 11.0 to 16.0 ns and 16.0 to 21.0 ns, respectively.

To better display the position of β-tubulin in the αβ-tubulin complex, Figures 9(a) and 9(b) show the crystal structure of αβ-tubulin from different viewpoints, respectively. Here, the criterion of setting three-dimensional graphics is the most conducive to exhibiting the mutual contact surface between α-tubulin and β-tubulin. Figures 9(c) and 9(d) show the motion modes of the β-tubulin alone system. The key motion mode of L713283-unbound β-tubulin is the movements parallel to the contact surface by referring to Figures 9(a) and 9(b). Figures 9(e) and 9(f) show the motion modes of the L713283-bound complex1 system. The range of motions is relatively smaller for the β-tubulin system after binding with L713283, which is consistent with the RMSF result from Figure 5(c). In addition, the entrance of hydrophobic pocket of β-tubulin exhibits a slight expansion after the association with L713283, which can increase the contact area with α-tubulin to some extent. Furthermore, the motion direction of β-tubulin is mainly pointing to α-tubulin, which is conducive to the stability of αβ-tubulin aggregate.

Previous studies [62, 66, 67] have shown that the association with Taxol significantly affects the lateral and vertical movements between the contact surfaces of αβ-tubulin. Specifically, both the increment of vertical movements and the decrement of lateral movement can greatly enhance the stability of αβ-tubulin. Considering the above discussion, the anticancer mechanism of L713283 may be similar to that of Taxol. L713283 stabilizes αβ-tubulin, inhibits its depolymerization, and prevents the formation of either α-tubulin or β-tubulin monomer. As is all known, either α-tubulin or β-tubulin monomer is essential to the life cycle of cancer cell [62, 66, 67]. To sum up, the possible anticancer mechanism of L713283 is that it may decrease the existence of α-tubulin and β-tubulin monomer by stabilizing αβ-tubulin. Therefore, the mitosis of cancer cells is reduced, which leads to the decrement in the number of cancer cells.

4. Conclusions

The molecular similarity comparison of L713283 with a series of anticancer drugs shows that L713283 possesses multiring structure and is most similar to Taxol, so it is speculated that the anticancer target of L713283 may be β-tubulin. Referring to experimental information, molecular docking method was used to obtain the binding mode between L713283 and β-tubulin. Subsequent comparative MD simulation and MM/GBSA methods were performed to investigate the binding free energy, conformational stability, and the motion mode differences between L713283-bound and L713283-unbound β-tubulin systems. It is found that β-tubulin exhibits less flexibility after the association with L713283, and this phenomenon is more apparent in the binding region. Meanwhile, the refined L713283 structure from the docked complex system adopts a conformation similar to that determined independently by energy-based semiempirical quantum chemical minimization. Combined with energy decomposition, it is found that the nonpolar interactions are favorable for the association between β-tubulin and L713283. The results of motion mode analysis show that the association with L713283 weakens the lateral movement and enhances the vertical movement between the subunits of αβ-tubulin, which leads to the increased stability of αβ-tubulin complex. It is speculated that the anticancer mechanism of L713283 is similar to that of Taxol through stabilizing the structure of αβ-tubulin and reducing its dynamicity. This work not only provides methodological guidance for exploring the anticancer mechanism for novel inhibitors, but also benefits the application studies of anticancer drug design.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Authors’ Contribution

Shan Chang and Hong-qiu He equally contributed to this study.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (31200990, 81202438, and 11247018), the Key Project of Sichuan Provincial Education Bureau (12ZA066), Leshan Science and Technology Project (14GZD022), and the supercomputer resources at Jiangsu University of Technology.

Supplementary Materials

Figure S1. The structure of L713283 with atomic number.

Figure S2. Superimposition between 1JFF crystal structure and the re-docked model of Taxol molecule with β-tubulin.

Figure S3. The binding modes of β-tubulin with the representative conformation of L713283 in the 5th~11th cluster, which are shown in red, green, blue, yellow, violet, cyan, orange, respectively.

Figure S4. The geometry parameter for the hydrogen bond between L713283 and β-tubulin. (A) Distance of donor atom with acceptor atom. (B) Angle of donor atom, hydrogen atom and acceptor atom.

  1. Supplementary Material