Interaction of microtubule depolymerizing agent indanocine with different human αβ tubulin isotypes

Tubulin isotypes are known to regulate the stability and dynamics of microtubules, and are also involved in the development of resistance against microtubule-targeted cancer drugs. Indanocine, a potent microtubule depolymerizing agent, is highly active against multidrug-resistant (MDR) cancer cells without affecting normal cells. It is known to disrupt microtubule dynamics in cells and induce apoptotic cell death. Indanocine is reported to bind to tubulin at the colchicine site i.e. at the interface of αβ tubulin heterodimer. However, it’s precise binding mode, involved molecular interactions and the binding affinities with different αβ-tubulin isotypes present in MDR cells are not well understood. Here, the binding affinities of human αβ-tubulin isotypes with indanocine were examined, employing the molecular modeling approach i.e. docking, molecular dynamics simulation and binding energy calculations. Multiple sequence analysis suggests that the amino acid sequences are different in the indanocine binding pockets of βI, βIIa, βIII and βVI isotypes. However, such differences are not observed in the amino acid sequences of βIVa, βIVb, and βV tubulin isotypes at indanocine binding pockets. Docking and molecular dynamics simulation results show that indanocine prefers the interface binding pocket of αβIIa, αβIII, αβIVb, αβV, and αβVI tubulin isotypes; whereas it is expelled from the interface binding pocket of αβIVa and αβI-tubulin isotypes. Further, binding free energy calculations show that αβVI has the highest binding affinity and αβI has the lowest binding affinity for indanocine among all β-tubulin isotypes. The binding free energy decreases in the order of αβVI > αβIVb > αβIIa > αβIII > αβV > αβIVa > αβI. Thus, our study provides a significant understanding of involved molecular interactions of indanocine with tubulin isotypes, which may help to design potent indanocine analogues for specific tubulin isotypes in MDR cells in future.


Introduction
Microtubules are dynamic cytoskeleton filamentous proteins; they play essential roles in cell division, cell movement, and intracellular transport [1].They are polymers of α/β-tubulin heterodimers.These α/β-tubulin are encoded by multiple genes which are expressed tissue-specifically e.g.βI is ubiquitous, βIII is expressed in neuronal and testicular cells, βIVa in neuronal and glial cells, βVI is observed in the erythroid cells and platelets etc. [2,3].In humans, ten β-tubulin and seven α-tubulin isotypes exist [4,5] and these isotypes show considerable difference at the C-terminal end.Tubulin isotypes composition plays an essential role in regulating microtubule dynamics [6][7][8].The essential roles of microtubules during the cell division make them important and attractive targets to design new anticancer agents.Anticancer agents are generally classified into microtubule stabilizing agents (MSA) and microtubule destabilizing agents (MDA).The MSAs prefer to bind at the 'taxol site' (e.g.paclitaxel, epothilone) whereas MDAs prefer to bind at the 'colchicine' and 'vinca' site (e.g.colchicine, indanocine, vinblastine), leading to cell death due to apoptosis in both the cases [5].
A major difficulty with the effectiveness of microtubule-targeting agents arises from the emergence of drug resistance, which is mainly due to the mutation in β-tubulin protein and an increased expression of the P-glycoprotein pump [9].In addition, an over-expression of βtubulin isotypes in cancerous cells also plays a crucial role in drug resistance, as they show lesser binding affinities for numerous anti-mitotic agents [10][11][12][13].Different tubulin isotypes are overexpressed in cancerous cells; particularly overexpression βI, βII, βIII, βIV, and βVtubulin isotypes are associated with multidrug-resistant cancer [14][15][16][17].Furthermore, it has also been observed that βII, βIII and βIV tubulin isotypes show differential binding affinities for a variety of anticancer drugs e.g.taxol, colchicine, DAMA-colchicine, and nocodazole [10,12,13,18].Therefore, these drug-resistant tubulin isotypes have been highlighted as interesting targets for designing new anticancer agents.
Indanocine, a synthetic indanone, is a microtubule depolymerizing agent with potent anti-proliferative activity [19].Indanocine acts against multidrug-resistant (MDR) cancer cells and kills non-dividing and quiescent cells [19], but it does not affect the normal nonproliferating cells.Indanocine affects the microtubule dynamicity at very low concentration and inhibits the migration of metastatic cancer cells [20].It prefers to bind at the interface of αβ tubulin heterodimer i.e. at the colchicine binding site [21].Indanocine is a flexible molecule in which the indanone group and the dimethylphenol group are connected by a single bond (Fig 1B).
Indanocine binds to the αβ tubulin heterodimer in a reversible manner, and it binds to tubulin at a faster rate than colchicine [21].Indanocine is a potent microtubule de-polymerizing agent and it acts against multidrug-resistant (MDR) cancer cells without affecting the normal cells.However, its precise binding mode involved molecular interactions and the binding affinities with different αβ-tubulin isotypes present in MDR cells are not well understood.Here, the binding affinities of human αβ-tubulin isotypes with indanocine were examined, employing molecular modeling approach.
To the best of our knowledge, structures of human tubulin isotypes bound with indanocine have not been determined using either X-Ray Crystallography or NMR techniques.Hence, the crystal structure of αβ-tubulin heterodimer (1SA0.pdb)from Protein Database was used as a template or reference.Bovine 1SA0.pdb has the βII tubulin which is identical in sequence to human βII tubulin [10].Therefore, crystal structure 1SA0.pdb was used as a template to build the 3D model of human tubulin isotypes.The sequence alignments of seven different human β-tubulin isotypes and the template sequence were performed to pinpoint the difference in residues at the indanocine binding pocket using Clustal Omega tool of EMBL-EBI [22].

Homology modeling of human αβ-tubulin isotypes
In this study, the colchicine-bound crystal structure of αβ-tubulin heterodimer (source code:1SA0.pdb)[23] was used as the template structure to build the 3D structures of seven human tubulin isotypes using homology modeling technique.Here, word template implies an initial structure used to build the desired model structure in the absence of required crystal structure, using homology modeling approach.The crystal structure of αβ-tubulin heterodimer '1SA0.pdb' is from bovine source with a resolution of 3.58Å [23].Previously, we have shown that human βIIa, βIII, βIVa are identical i.e. 100% similar to bovine βIIa, βIII, and βIVa [10].Similarly, we have also checked the similarity between other tubulin isotypes, and have found that the human βIVb, βV, and βVI are also identical i.e. 100% similar to bovine βIVb (Q3MHM5), βV(Q2KJD0), βVI(Q2HJ81) whereas human βI and bovine βI (E1BJK2) show 88.89% similarity.
We considered chains A and B from the crystal structure 1SA0.pdb for homology modeling whose missing residues of β-tubulin (37 to 47) and α-tubulin (1, 275-284) were built using MODELLER9v18 [24].This refined structure of αβ-tubulin heterodimers obtained from MODELLER9v18 was used as template and will be referred to as tubulin 1SA0 hereafter.The homology models of seven human αβ-tubulin isotype heterodimers such as αβI, αβIIa, αβIII, αβIVa, αβIVb, αβV, and αβVI were built using the template structure 1SA0.pdbthrough Here, H7 denotes the α-helix number 7, T7 stands for T loop number 7, and B9 implies the β-sheet number 9. The GTP in α-tubulin, GDP in β-tubulin are shown using space-fill models.The white, grey, red, blue and golden yellow colors represent carbon, hydrogen, oxygen, nitrogen and phosphorous atoms, respectively.(B) The structure of indanocine has a dimethoxyaniline group (labeled as X ring) and a dimethylphenol group (labeled as Y ring).The carbon, oxygen, nitrogen, and hydrogen atoms of indanocine are shown in green, red, blue and white color, respectively.https://doi.org/10.1371/journal.pone.0194934.g001 MODELLER9v18 and the best models were selected on the basis of their DOPE score [24].The C-terminal ends of tubulin isotypes were not included in our study as they are not present at the interface of αβ-tubulin.Therefore, we do not expect C-terminal tails to play a direct role in drug binding at the interface.C-terminal ends of the different isotypes have been modeled by Luchko et al. [25] for conformational analysis.A recent study of interaction of different human tubulin isotypes with drug DAMA-Colchicine [10] shows that C-terminal end does not qualitatively affect the binding of DAMA-Colchicine with tubulin isotypes.The stereochemical quality of αβ-tubulin models was evaluated using PROCHECK [26] and Verify-3D [27] to ensure the reliability of the homology models, whose details are given in the Supplementary S1 Text.Subsequent energy minimization was performed on tubulin 1SA0 and seven different αβ-tubulin isotype heterodimers using 5,000 steps of steepest descent method; out of which 3,000 steps were of conjugate gradient method using AMBER12 software [28].For energy minimization, the parameters of guanosine triphosphate (GTP), and guanosine diphosphate (GDP) and Mg 2+ were obtained from the AMBER database [29,30].These energy minimized structures of αβ-tubulin isotypes were then used for the docking of indanocine using AutoDock4.2[31].
In this study, we kept the α-tubulin as constant and varied the beta-tubulin isotypes as accurate combinations of different αβ-tubulin isotypes is not well known experimentally.Indanocine works on various multi-drug resistant cancer cell types where over-expression of different beta-tubulin isotypes (as compared to the α-tubulin isotypes) leads to drug-resistance.Hence, we did not use all possible combinations of different αβ-tubulin isotypes in this study.

Molecular docking of indanocine with αβ-tubulin isotypes
To identify the interactions of tubulin 1SA0 and different human αβ-tubulin isotypes with indanocine, molecular docking was performed using AutoDock4.2[31].For molecular docking, energy-minimized 3D atomic coordinates of indanocine were generated using the PRODRG server [32].Since indanocine was suggested to bind at the interface of αβ tubulin [21], an autogrid was used to outline the putative binding pocket around the interface of αβ tubulin [31].The Gasteiger charges were added to αβ tubulin using AutoDock4.2[31].Here, we used local docking methodology to delineate the binding mode of indanocine with tubulin [33,34].
A grid box of size 60Å×60Å×60Å with a spacing of 0.375Å was prepared at the αβ tubulin interface i.e. the putative indanocine binding site.The Lamarckian Genetic Algorithm (LGA) was used for molecular docking with default parameters [31].Here, a total of 50 independent flexible ligand dockings were conducted, each composed of 100 LGA runs, which yielded a total of 5,000 conformations.They were subsequently clustered using an all-atom RMSD cutoff of 4Å; which were then analyzed considering cluster size and binding free energy calculated by a scoring function of AutoDock4.2[31,35].The lowest binding free energy docked conformation of indanocine was selected for further hydrogen bonding interactions analysis and for molecular dynamics simulations.

Molecular dynamics simulation
Molecular dynamics simulations were performed for indanocine-docked complexes with tubulin 1SA0 and seven different human αβ tubulin isotypes i.e. αβI, αβIIa, αβIII, αβIVa, αβIVb, αβV, αβVI, and αβVII using the SANDER module of AMBER12 [28].The AMBER ff99SB force field was applied for protein, and the parameters for guanosine triphosphate (GTP), guanosine diphosphate (GDP) and Mg 2+ were taken from the AMBER database [29,30].Parameters for indanocine were generated by using the 'Antechamber' module of AMBER12 [10,36].The implicit 'Generalized Born/Surface Area (GB/SA)' model was used to represent the solvent effect by using the parameters described by Tsui [37] to explore the interactions of protein-ligand [10,36].The molecular dynamics simulations steps such as minimization, heating, equilibration and production run were performed using the same parameters as in our earlier studies [10].The trajectories of molecular dynamics simulations were visualized and analyzed using VMD [38] and PyMol [39].VMD was employed to produce the molecular dynamics simulation movies by setting the Trajectory Smoothing Window size for protein to 5, and for GTP, GDP and indanocine to 3.

Binding energy calculations
The αβ tubulin isotype-indanocine binding free energy calculations were estimated using MM-GBSA approaches using AMBER12 [28].The binding free energy was calculated using 10,000 frames from the last 2ns of molecular dynamics trajectories with an interval of 5 for each system using mmpbsa module of AMBER12 similar to our earlier study [10].The entropy calculations are computationally expensive and hence omitted in this study, as done in an earlier study [10].The need of explicit calculation of the entropy can be avoided in this study as we are comparing the relative trend of binding free energies of different isotypes which are related systems (there is a difference of few residues among them) [10,40].

Sequence analysis and homology modeling of αβ-tubulin isotypes
Multiple sequence analysis of the seven above-mentioned different human β-tubulin isotypes against bovine β II tubulin (PDB code: 1SA0, chain B) as reference sequence was performed using Clustal Omega tools of EMBL-EBI [22].The multiple sequence analysis study shows that human β-tubulin isotypes show residue composition variations at different locations (Fig 2).We further analyzed the residue composition variations at the indanocine binding pocket of different β-tubulin isotypes.The indanocine binding pocket of βI has five residue changes i.e.Val236-Ile, Cys239-Ser, Ala315-Cys, Val316-Ile, and Thr351-Val, βIIa has a single amino acid change i.e.Val316-Ile, βIII has three residue changes i.e.Cys239-Ser, Ala315-Thr, and Thr351-Val (Fig 2 ), and βVI also has three residue changes similar to βIII tubulin isotypes such as Cys239-Ser, Ala315-Thr, and Thr351-Val.There is no residue composition variation in βIVa, βIVb, and βV at the indanocine binding pocket (Fig 2).We then built homology models of these seven human αβ-tubulin isotypes and performed molecular docking of indanocine and molecular dynamics simulations of αβ tubulin-indanocine complexes to explore the effect of residue composition on the binding interaction of indanocine.
Molecular docking results suggest that the residue composition variation in and around the indanocine binding pocket results in differences in binding energy, conformation, and hydrogen bonding interactions among the different αβ-tubulin isotypes-indanocine complexes (Fig 3A -3H and Table 1).Thus, our docking studies suggest that Lys-252, Lys-350, Cys-239, Val-236, Ala-248, Leu-246 of β-tubulin and Ala-101, Ser-178, Thr-179 and Val-180 play a key role in the stabilization of indanocine at the interface of all αβ-tubulin heterodimer (S9A-S9H Fig) .Further, we calculated the electrostatic contact potential over the tubulin 1SA0 and seven different β-tubulin isotypes-indanocine complex using PyMol [39] (S10A-S10H Fig) .The electrostatic contact potentials show that the dimethylphenol group is immersed inside the cavity of β-tubulin while the dimethoxyaniline group of indanocine is located out of β-tubulin protein.In addition, the effect of residue composition variation in and around the indanocine binding pocket in different β-tubulin isotypes was further elucidated using molecular dynamics simulations and binding free energy calculations.

Molecular dynamics simulation of αβ tubulin isotypes-indanocine complexes
To understand the refined binding mode of αβ-tubulin isotypes with indanocine, we performed molecular dynamics simulations over the lowest energy αβ-tubulin-indanocine docked complexes (Fig 3A -3H) as our starting structure using AMBER12 [28].The primary analysis was done by looking at the molecular dynamics simulation stability (Fig 4) and analysis of molecular dynamics simulated average structures (Fig 5A -5H).
The Root mean square deviations (RMSD) of Cα backbone atoms of a production molecular dynamics simulations was calculated to examine the stability of the molecular dynamics simulation.The RMSD analysis of all the different αβ-tubulin isotypes-indanocine complexes suggests that all αβ-tubulin-indanocine complexes reached their equilibrium conformation after a time period of 20ns and then retained their stability with fluctuations between 2.5-4.5Å(Fig 4).Molecular dynamics simulation results clearly show that indanocine prefers to bind at the αβ-tubulin interface binding pocket in tubulin 1SA0 and αβIIa, αβIII, αβIVb, αβV, and αβVI tubulin isotypes (S1 Movie, S3 and S4 Movies, S6-S8 Movies) respectively, while in case of tubulin isotype αβI (S2 Movie) and αβIVa (S5 Movie), indanocine is expelled from the αβtubulin interface binding cavity.In αβI and αβIVa tubulin isotypes, the T7 loop of β-tubulin moves backward, while the B9 sheet of β-tubulin and T5 loop of α-tubulin also undergo conformational changes which makes ample space at the interface leading to the expulsion of indanocine from the interface of αβ-tubulin isotypes heterodimer (S2 Movie and S5 Movie).Such conformational changes are seen in cases of αβI and αβIVa tubulin isotypes but are not seen on other tubulin isotypes (S1 Movie, S3 and S4 Movies, S6-S8 Movies).The residues present in the B9 sheet, H7 helix, T7 loop and H8 helix of β-tubulin and T5 loop of α-tubulin have important contributions in the binding of indanocine at the interface in other αβ-tubulin isotypes (S1 Movie, S3 and S4 Movies, S6-S8 Movies).A detailed analysis of residues involved in the hydrogen bonding interactions with indanocine is discussed in the next section.

Binding energy calculations
As reported earlier [10], the binding free energies for different αβ-tubulin isotypes with indanocine were calculated ignoring the entropic contribution to the binding free energy (Table 3).The estimated binding free energies (ΔE bind ) of tubulin 1SA0 and different αβI, αβIIa, αβIII, αβIVa, αβIVb, αβV, and αβVI tubulin isotypes with indanocine are -49.90,-41.39, -44.03, -43.47, -41.50, -44.57, -42.97, and -50.70 kcal/mol, respectively (Table 3).The αβVI has the highest binding free energy for indanocine, whereas αβI has the lowest binding free energy among the other αβ-tubulin isotypes.The binding free energy decreases is in the order of αβVI > αβIIb > αβIVb > αβIIa > αβIII > αβV > αβIVa > αβI.The lower binding free energy of indanocine for αβI-tubulin isotype is due to maximum residue changes in the binding pocket of βI-tubulin such as Val236-Ile, Cys239-Ser, Ala315-Cys, Val316-Ile, and Thr351-Val as compared to other tubulin isotypes (Fig 2).However, indanocine exhibits a good binding affinity for αβVI, αβIIb, αβIVb, αβIIa, αβIII, and αβV (Table 3).The van der Waal (ΔE vdw ) and electrostatic (ΔE ele ) interactions are important for the binding of the protein-ligand complex.Here, van der Waal interactions make the highest contribution towards the binding free energy (Table 3), while the solvation energy (E sol ) is unfavorable for binding of ligand.The αβI-tubulin isotype shows the lowest van der Waals interaction energy in comparison to the other αβ-tubulin isotype-indanocine complexes.The net binding free energy is decided by a competition between E gas and E sol , and is the lowest for αβI-tubulin isotype (Table 3).

Conclusion
In this study, the binding affinity of indanocine with tubulin 1SA0 and seven human tubulin isotypes αβI, αβIIa, αβIII, αβIVa, αβIVb, αβV, and αβVI was investigated using sequence analysis, homology modeling, molecular docking, molecular dynamics simulations and binding free energy calculations.The residue compositions were found to be different at the indanocine binding pocket of human βI, βIIa, βIII and βVI tubulin isotypes, whereas no such differences were found in the βIVa, βIVb and βV tubulin isotypes.
Further, molecular docking results show that indanocine prefers to bind at the interface of all αβ-tubulin isotypes i.e. at the colchicine binding site, as observed in the previous experimental study [21].Indanocine shows different binding mode and binding energy for different αβ-tubulin isotypes; this might be due to the residue composition changes in and around the binding pocket of β-tubulin isotypes.The residues in the H7-Helix (Cys-239, Ile-236), T7-loop (Leu-246, Ala-248), H8-helix (Lys-252, Asn-256) and B9-sheet (Lys-350) of β-tubulin and T5-loop (Ser-178, Thr-179, Val-180) of α-tubulin are involved in the hydrogen bonding interactions with indanocine.
Molecular dynamics simulations were performed on αβ-tubulin isotype-indanocine docked complexes, to further investigate the effect of residue composition differences on the binding of indanocine.Our molecular dynamics simulations results show that indanocine is completely adopted inside the binding pocket of tubulin 1SA0 and αβIIa, αβIII, αβIVb, αβV and αβVI tubulin isotypes, whereas it is expelled from the interface of αβI, and αβIVa-tubulin isotype.Here, the T7-loop of β-tubulin moves backward; meanwhile the B9 sheet of β-tubulin and T5-loop of αtubulin shows conformational change.This leads to making an ample space at the interface of αβI, and αβIVa tubulin isotypes which is helpful to expel indanocine from the interface cavity.Whereas in case of other αβ-tubulin isotypes-indanocine complexes, the T7 loop moves forward which helps to adopt indanocine at the binding pocket.Further, binding free energy calculations show that the tubulin isotypes αβIIa, αβIII, αβIVa, αβIVb, αβV and αβVI have the highest binding free energy and αβI-tubulin isotype has the lowest binding free energy for indanocine.One of the reasons behind the less binding free energy of αβI-tubulin isotype toward indanocine might be due to maximum residue changes at the binding site.Thus, our present computational study provides a detailed understanding of the molecular interactions of human αβ-tubulin isotypes with indanocine and provides insight for designing superior indanocine analogues with isotype specificity.These superior analogues can be valuable in the treating patients with advanced carcinomas which exhibit tubulin isotype specificity or can be helpful in developing personalized medicines for cancer patients.

Fig 1 .
Fig 1.The representative structure of αβ-tubulin dimer and indanocine.(A) α/β-tubulin subunits.α subunit is shown in blue color and β subunit is shown in olive green color.Indanocine binding site i.e. 'colchicine site' is at the interface of α-tubulin and β-tubulin heterodimer (red dotted circle); while the 'taxol site' is present on only β-tubulin over the H7 helix (blue dotted circle).The indanocine binding pocket consists of H7 helix (shown in orange), T7 loop (shown in cyan), H8 helix (shown in magenta) and T5 loop of α-tubulin (shown in cyan).Here, H7 denotes the α-helix number 7, T7 stands for T loop number 7, and B9 implies the β-sheet number 9. The GTP in α-tubulin, GDP in β-tubulin are shown using space-fill models.The white, grey, red, blue and golden yellow colors represent carbon, hydrogen, oxygen, nitrogen and phosphorous atoms, respectively.(B) The structure of indanocine has a dimethoxyaniline group (labeled as X ring) and a dimethylphenol group (labeled as Y ring).The carbon, oxygen, nitrogen, and hydrogen atoms of indanocine are shown in green, red, blue and white color, respectively.

Fig 2 .
Fig 2. Multiple sequence alignment of human β-tubulin isotypes.The isotype βI shows variations of residues Val236-Ile, Cys239-Ser, Ala315-Cys, Val316-Ile and Thr351-Val, βIIa shows a change of Val316-Ile, βIII and βVI shows a change of Cys239-Ser, Ala315-Thr, and Thr351-Val at the indanocine binding pocket.Residue variations in the indanocine binding pocket are shown in red.Symbol ' Ã ' denotes positions of amino acid which have a single, fully conserved amino acid residue; the symbol ': ' denotes conservation between groups of strongly similar properties of amino acid; the symbol '.' denotes the conservation between groups of weakly similar properties of amino acids, and the symbol '-' denotes gaps inserted to maximize sequence alignment [22].https://doi.org/10.1371/journal.pone.0194934.g002