# Study on Molecular Recognition betweenEuphorbia Factor L713283 and β-Tubulin via Molecular Simulation Methods

Authors: Shan Chang; Hong-qiu He; Ren Kong; Zhen-jian Xie; Jian-ping Hu
Journal: Journal of Chemistry
 (2015)

---

## Abstract

Euphorbia factor L713283 is a new lathyrane diterpene isolated fromEuphorbia 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.


---

## Body

## 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 [3–13] 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 onEuphorbia 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 [19–21].Recently, Jiao et al. extracted a new diterpenoid moleculeEuphorbia factor L713283 fromEuphorbia lathyris [22], with IC50 value of 248.30
±
6.62 μ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 usedin 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 fromEuphorbia helioscopia, which are all jatrophane-type diterpenoids [25]. However, the possible targets and molecular mechanism of this novel anticancer agent remain ambiguous.Figure 1
The structure of L713283 ligand.To interpret anticancer mechanism of L713283 in the present paper, we employed integratedin 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.Figure 2
Protocol for the determination of drug target and the investigation of the possible anticancer mechanism for L713283 ligand.


### 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 [3–13] 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., 2,3
,
…
,
128) 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-n with the lowest docking energy in the top n 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 [37–40]. 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 [42–44]. 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 Δ
G
b
i
n
d is defined as(1)
Δ
G
b
i
n
d
=
Δ
E
M
M
+
Δ
G
s
o
l
-
T
Δ
S
,where E
M
M is the molecular mechanical energies in vacuo, G
s
o
l is the solvation free energy, T is the temperature, and S is the solute entropy. The solute entropy S is usually estimated by normal mode analysis method [49].E
M
M can be obtained by average of the molecular mechanics energy from the collected snapshots during MD simulation. It is defined as(2)
E
M
M
=
E
i
n
t
e
r
n
a
l
+
E
e
l
e
c
t
r
o
s
t
a
t
i
c
+
E
V
D
W
,where E
i
n
t
e
r
n
a
l includes the bond, angle, and dihedral energies, E
e
l
e
c
t
r
o
s
t
a
t
i
c is electrostatic interaction energy, and E
V
D
W is van der Waals interaction energy.The contribution from the solvent is calculated as(3)
G
s
o
l
=
G
G
B
+
G
S
A
,where G
G
B is the sum of electrostatic salvation energy (polar contribution) and G
S
A is nonelectrostatic salvation component (nonpolar contribution). The polar contribution is calculated using Generalized Born (GB) approximation model [50–52]. The nonpolar contribution is estimated by solvent accessible surface area (SASA) [53], and it is calculated with the expression(4)
G
S
A
=
γ
×
S
A
S
A
+
b
,where γ and b are the default constants. Then, (1) can be written as(5)
Δ
G
b
i
n
d
=
Δ
E
i
n
t
e
r
n
a
l
+
Δ
E
e
l
e
c
t
r
o
s
t
a
t
i
c
+
Δ
E
V
D
W
+
Δ
G
G
B
+
Δ
G
S
A
-
T
Δ
S
.Each 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.


## 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.Figure 2
Protocol for the determination of drug target and the investigation of the possible anticancer mechanism for L713283 ligand.


## 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 [3–13] 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., 2,3
,
…
,
128) 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-n with the lowest docking energy in the top n 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 [37–40]. 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 [42–44]. 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 Δ
G
b
i
n
d is defined as(1)
Δ
G
b
i
n
d
=
Δ
E
M
M
+
Δ
G
s
o
l
-
T
Δ
S
,where E
M
M is the molecular mechanical energies in vacuo, G
s
o
l is the solvation free energy, T is the temperature, and S is the solute entropy. The solute entropy S is usually estimated by normal mode analysis method [49].E
M
M can be obtained by average of the molecular mechanics energy from the collected snapshots during MD simulation. It is defined as(2)
E
M
M
=
E
i
n
t
e
r
n
a
l
+
E
e
l
e
c
t
r
o
s
t
a
t
i
c
+
E
V
D
W
,where E
i
n
t
e
r
n
a
l includes the bond, angle, and dihedral energies, E
e
l
e
c
t
r
o
s
t
a
t
i
c is electrostatic interaction energy, and E
V
D
W is van der Waals interaction energy.The contribution from the solvent is calculated as(3)
G
s
o
l
=
G
G
B
+
G
S
A
,where G
G
B is the sum of electrostatic salvation energy (polar contribution) and G
S
A is nonelectrostatic salvation component (nonpolar contribution). The polar contribution is calculated using Generalized Born (GB) approximation model [50–52]. The nonpolar contribution is estimated by solvent accessible surface area (SASA) [53], and it is calculated with the expression(4)
G
S
A
=
γ
×
S
A
S
A
+
b
,where γ and b are the default constants. Then, (1) can be written as(5)
Δ
G
b
i
n
d
=
Δ
E
i
n
t
e
r
n
a
l
+
Δ
E
e
l
e
c
t
r
o
s
t
a
t
i
c
+
Δ
E
V
D
W
+
Δ
G
G
B
+
Δ
G
S
A
-
T
Δ
S
.Each 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 [3–13], 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 [55–57]. 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.Figure 3
Molecular similarity of L713283 to ten known anticancer drugs.Notably, some previous biological experimental results also indicated thatβ-tubulin might be the target of L713283. Khaleghian et al. extracted the inganen fromEuphorbia 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 fromEuphorbia semiperfoliata (Euphorbiaceae family) [59]. They investigated the interactions between jatrophanes and purified bovine brain tubulin by anin 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 withEuphorbia factor L713283.


### 3.2. Molecular Docking

Figure4 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
The docking result of L713283 withβ-tubulin. (a) The docked conformations were clustered with the RMSD cutoff value of 0.2 nm. (b) The binding modes of both L713283 docked conformations and the crystal structure of Taxol ligand. Taxol ligand is represented using ball-and-stick model in red. The representative docked conformations of L713283 are shown with stick model in green, blue, orange, and purple, respectively.
(a)
(b)Figure4(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%.Figures5(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 0.11
±
0.01 and 0.27
±
0.03 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.Figure 5
Comparative RMSD and RMSF analysis of complex1, complex2, andβ-tubulin systems. (a) The RMSD values of backbone Cα atoms versus simulation time. (b) The RMSD values of L713283 versus simulation time in complex1 and complex2 systems. (c) The RMSF values of Cα atoms in complex1 and β-tubulin systems. (d) The correlation of RMSF between complex1 and β-tubulin systems.
(a)
(b)
(c)
(d)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. Figure5(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.Figure 6
Fluctuation region analysis ofβ-tubulin. (a) The solid ribbon representation of αβ-tubulin system (PDB code: 1JFF). (b) The high (color in red) and low (text with rectangular border, color in blue) fluctuation of β-tubulin.
(a)
(b)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. Table1 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.Table 1
Energy components contributing to the binding free energy (kJ/mol).
Item
Receptor
Ligand
Complex
Deltah
ELEa
−51089.9 ± 385.7
−339.7 ± 3.1
−51489.9 ± 388.1
−60.4 ± 9.2
VDWb
−7806.0 ± 169.7
60.2 ± 19.9
−7968.7 ± 179.4
−222.9 ± 12.3
INTc
36758.4 ± 230.5
547.7 ± 36.5
37306.0 ± 234.3
0.0 ± 0.0
GBSURd
558.7 ± 4.9
23.6 ± 0.2
551.6 ± 4.9
−30.7 ± 0.6
GBe
−25766.8 ± 296.4
−79.6 ± 2.9
−25701.2 ± 301.9
145.2 ± 8.5
TSTOTf
19675.9 ± 71.3
310.6 ± 0.0
19858.4 ± 68.8
−128.0 ± 82.3
Δ
G
b
i
n
d
g
−40.8
aElectrostatic energy. bThe van der Waals energy. cBond, angle, and dihedral energies. dNonpolar solvent free energy. ePolar solvent free energy. fTemperature × entropy. gBinding free energy (ELE + VDW + INT + GBSUR + GB − TSTOT). hDifference (complex − ligand − receptor).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.Figure 7
Energy contribution of each residue to the binding between L713283 andβ-tubulin. The energy distributions are divided into the nonpolar interaction (a), polar interaction (b), and total effect (c).
(a)
(b)
(c)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].Figure 8
The refined binding mode between L713283 andβ-tubulin. (a) The L713283 and key residues are shown in stick model. The four key residue regions are depicted with Cα-wire model in green, blue, orange, and purple. The polar types of each eight amino acids are shown in the corresponding rectangles. (−), (0), and (+) are used to represent the amino acids in negative, nonpolar, and positive charges, respectively. (b) The electrostatic potential map of the active pocket of β-tubulin. L713283 is shown in green stick model. The regions of positive and negative charges are, respectively, colored in blue and red.
(a)
(b)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

Table2 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.Table 2
The dihedral angle information of both rings and side chains for L713283 molecule.
Regions
Dihedral angles
Average value
Standard deviation
Fluctuation scopea
Initial structure
Large ring
C1-C2-C3-C4
−174.79
13.52
7.6%
−172.51
C2-C3-C4-C5
−131.89
14.57
8.1%
−136.26
C3-C4-C5-C11
1.50
7.60
4.2%
3.20
C4-C5-C11-C12
38.74
8.29
4.6%
36.22
C5-C11-C12-C13
−127.44
7.66
4.3%
−128.23
C11-C12-C13-C16
175.50
6.01
3.3%
174.21
C12-C13-C16-C18
−59.91
7.30
4.1%
−60.01
C13-C16-C18-C19
−60.98
6.26
3.5%
−61.29
C16-C18-C19-C1
127.87
4.98
2.7%
128.63
C18-C19-C1-C2
−63.13
7.40
4.1%
−64.52
C19-C1-C2-C3
−33.80
15.41
8.6%
−38.25
Middle ring
C18-C19-C21-C22
35.67
4.57
2.5%
35.29
C19-C21-C22-C23
−33.61
4.39
2.4%
−33.54
C21-C22-C23-C18
42.04
4.63
2.6%
40.52
C22-C23-C18-C19
123.50
4.31
2.4%
123.48
C23-C18-C19-C21
−105.84
5.21
2.9%
−103.24
Side chain
O20-C19-C21-C22
−152.06
4.86
2.7%
−153.91
O25-C23-C18-C19
179.82
6.11
3.4%
177.10
O17-C16-C18-C19
122.90
8.30
4.6%
129.51
C14-C13-C16-C18
−8.41
6.82
3.8%
−10.21
O7-C6-C4-C3
−88.59
5.86
3.2%
−89.54
aThe fluctuation scope is the standard deviation/180.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.Figure 9
Comparative analysis of motion modes for theβ-tubulin and docked complex1 systems. ((a)-(b)) The structures of αβ-tubulin and the β-tubulin subunit are marked with red oval. ((c)-(d)) The L713283-unbound β-tubulin system. ((e)-(f)) The complex1 system. The motion mode is shown with cone model, in which the length is correlative with motion magnitude and the arrow corresponds to motion direction. Protein backbone is represented with Cα stick model in gray, while small molecule L713283 is shown with ball-and-stick model.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.


## 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 [3–13], 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 [55–57]. 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.Figure 3
Molecular similarity of L713283 to ten known anticancer drugs.Notably, some previous biological experimental results also indicated thatβ-tubulin might be the target of L713283. Khaleghian et al. extracted the inganen fromEuphorbia 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 fromEuphorbia semiperfoliata (Euphorbiaceae family) [59]. They investigated the interactions between jatrophanes and purified bovine brain tubulin by anin 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 withEuphorbia factor L713283.


## 3.2. Molecular Docking

Figure4 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
The docking result of L713283 withβ-tubulin. (a) The docked conformations were clustered with the RMSD cutoff value of 0.2 nm. (b) The binding modes of both L713283 docked conformations and the crystal structure of Taxol ligand. Taxol ligand is represented using ball-and-stick model in red. The representative docked conformations of L713283 are shown with stick model in green, blue, orange, and purple, respectively.
(a)
(b)Figure4(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%.Figures5(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 0.11
±
0.01 and 0.27
±
0.03 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.Figure 5
Comparative RMSD and RMSF analysis of complex1, complex2, andβ-tubulin systems. (a) The RMSD values of backbone Cα atoms versus simulation time. (b) The RMSD values of L713283 versus simulation time in complex1 and complex2 systems. (c) The RMSF values of Cα atoms in complex1 and β-tubulin systems. (d) The correlation of RMSF between complex1 and β-tubulin systems.
(a)
(b)
(c)
(d)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. Figure5(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.Figure 6
Fluctuation region analysis ofβ-tubulin. (a) The solid ribbon representation of αβ-tubulin system (PDB code: 1JFF). (b) The high (color in red) and low (text with rectangular border, color in blue) fluctuation of β-tubulin.
(a)
(b)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. Table1 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.Table 1
Energy components contributing to the binding free energy (kJ/mol).
Item
Receptor
Ligand
Complex
Deltah
ELEa
−51089.9 ± 385.7
−339.7 ± 3.1
−51489.9 ± 388.1
−60.4 ± 9.2
VDWb
−7806.0 ± 169.7
60.2 ± 19.9
−7968.7 ± 179.4
−222.9 ± 12.3
INTc
36758.4 ± 230.5
547.7 ± 36.5
37306.0 ± 234.3
0.0 ± 0.0
GBSURd
558.7 ± 4.9
23.6 ± 0.2
551.6 ± 4.9
−30.7 ± 0.6
GBe
−25766.8 ± 296.4
−79.6 ± 2.9
−25701.2 ± 301.9
145.2 ± 8.5
TSTOTf
19675.9 ± 71.3
310.6 ± 0.0
19858.4 ± 68.8
−128.0 ± 82.3
Δ
G
b
i
n
d
g
−40.8
aElectrostatic energy. bThe van der Waals energy. cBond, angle, and dihedral energies. dNonpolar solvent free energy. ePolar solvent free energy. fTemperature × entropy. gBinding free energy (ELE + VDW + INT + GBSUR + GB − TSTOT). hDifference (complex − ligand − receptor).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.Figure 7
Energy contribution of each residue to the binding between L713283 andβ-tubulin. The energy distributions are divided into the nonpolar interaction (a), polar interaction (b), and total effect (c).
(a)
(b)
(c)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].Figure 8
The refined binding mode between L713283 andβ-tubulin. (a) The L713283 and key residues are shown in stick model. The four key residue regions are depicted with Cα-wire model in green, blue, orange, and purple. The polar types of each eight amino acids are shown in the corresponding rectangles. (−), (0), and (+) are used to represent the amino acids in negative, nonpolar, and positive charges, respectively. (b) The electrostatic potential map of the active pocket of β-tubulin. L713283 is shown in green stick model. The regions of positive and negative charges are, respectively, colored in blue and red.
(a)
(b)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

Table2 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.Table 2
The dihedral angle information of both rings and side chains for L713283 molecule.
Regions
Dihedral angles
Average value
Standard deviation
Fluctuation scopea
Initial structure
Large ring
C1-C2-C3-C4
−174.79
13.52
7.6%
−172.51
C2-C3-C4-C5
−131.89
14.57
8.1%
−136.26
C3-C4-C5-C11
1.50
7.60
4.2%
3.20
C4-C5-C11-C12
38.74
8.29
4.6%
36.22
C5-C11-C12-C13
−127.44
7.66
4.3%
−128.23
C11-C12-C13-C16
175.50
6.01
3.3%
174.21
C12-C13-C16-C18
−59.91
7.30
4.1%
−60.01
C13-C16-C18-C19
−60.98
6.26
3.5%
−61.29
C16-C18-C19-C1
127.87
4.98
2.7%
128.63
C18-C19-C1-C2
−63.13
7.40
4.1%
−64.52
C19-C1-C2-C3
−33.80
15.41
8.6%
−38.25
Middle ring
C18-C19-C21-C22
35.67
4.57
2.5%
35.29
C19-C21-C22-C23
−33.61
4.39
2.4%
−33.54
C21-C22-C23-C18
42.04
4.63
2.6%
40.52
C22-C23-C18-C19
123.50
4.31
2.4%
123.48
C23-C18-C19-C21
−105.84
5.21
2.9%
−103.24
Side chain
O20-C19-C21-C22
−152.06
4.86
2.7%
−153.91
O25-C23-C18-C19
179.82
6.11
3.4%
177.10
O17-C16-C18-C19
122.90
8.30
4.6%
129.51
C14-C13-C16-C18
−8.41
6.82
3.8%
−10.21
O7-C6-C4-C3
−88.59
5.86
3.2%
−89.54
aThe fluctuation scope is the standard deviation/180.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.Figure 9
Comparative analysis of motion modes for theβ-tubulin and docked complex1 systems. ((a)-(b)) The structures of αβ-tubulin and the β-tubulin subunit are marked with red oval. ((c)-(d)) The L713283-unbound β-tubulin system. ((e)-(f)) The complex1 system. The motion mode is shown with cone model, in which the length is correlative with motion magnitude and the arrow corresponds to motion direction. Protein backbone is represented with Cα stick model in gray, while small molecule L713283 is shown with ball-and-stick model.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.



---

*Source: 879238-2015-11-23.xml*