Structural basis for drug resistance conferred by β-tubulin mutations: a molecular modeling study on native and mutated tubulin complexes with epothilone B

Using molecular modeling, we have investigated the structure and dynamic properties of epothilone B–tubulin complexes with wild-type and mutated tubulin, aimed at identifying the molecular factors involved in the emergence of drug resistance induced by four protein mutations at Phe270Val, Thr274Ile, Arg282Gln, and Gln292Glu. Our results revealed that tubulin mutations render significant changes in the protein conformation in regions involved either in the binding of the ligand or in interdimer contacts that are relevant to the assembly of stable microtubules. In addition, point mutations induce drastic changes in the binding pose of the ligand and in the interaction networks responsible for the epothilone–tubulin association. Large ligand displacements inside the binding pocket and an overall decrease in the strength of drug-receptor polar contacts suggest a looser binding of the ligand in tubulin mutants. These results explain the loss of activity for epothilone B against cancer cells that contain tubulin mutants and provide valuable information to enhance the understanding of the atomic source of epothilones’ activity, which can be helpful to conduct further research on the rational design of more potent therapeutic tubulin-binding agents.


Introduction
Microtubules are cytoskeletal structures composed of αβ-tubulin heterodimers that have a crucial role in cell division, which makes them attractive targets for anticancer drug design. Tubulin-binding agents comprise many chemically diverse compounds that bind to the β-subunit of the tubulin dimer and inhibit cell division either by acting as microtubule-stabilizers (such as taxanes) or microtubule-destabilizer drugs (such as vinca alkaloids) (Jordan, Hadfield, Lawrence, & McGown, 1998). Within microtubule-stabilizers binding agents, epothilones ( Figure 1) have attracted major attention due to their natural origin, aqueous solubility, and superior antitumor properties compared to taxanes (Bollag et al., 1995). In the past years, epothilones have been the subject of extensive investigation and constitute lead compounds in the search for new analogs with improved chemotherapeutic properties. Although these 16-membered macrolides have dissimilar structures compared to Taxol, they bind to a common binding pocket in the intermediate domain of β-tubulin, in a region of association between the M-loop and the helix H7 of the protein. It has been proposed that the overall effect of epothilone binding to tubulin is to lock the relative positions of helix H7 and the M loop in such a way that the formation of straight protofilaments is favored and the assembly equilibrium of tubulin is shifted toward stable microtubules (Diaz, Andreu, & Jiménez-Barbero, 2009).
A major problem concerning epothilones efficacy arises due to the emergence of drug resistance, which can be originated by several mechanisms, such as changes in cellular drug uptake, metabolic drug deactivation, changes in cellular components that interact with the target, or mutations in the drug target (Berrieman, Lind, & Cawkwell, 2004). Experiments have found that point mutants at Phe270Val, Thr274Ile, Arg282Gln, and Gln292Glu affect the ability of epothilones to induce tubulin polymerization as well as inhibit cell growth compared with parental cells (Table 1) (Giannakakou & Snyder, 2008;Giannakakou et al., 1997Giannakakou et al., , 2000He, Yang, & Horwitz, 2001). These mutations are located in the vicinity of the binding pocket and are expected to affect binding affinity of the ligand toward complexation with β-tubulin ( Figure 2). Thr274 is situated in the M-loop of β-tubulin and has been reported to be critical for epothilone association through the formation of intermolecular hydrogen bonds involving the C1-carbonyl and the C3-and C7-hydroxyl groups of the macrocycle (Canales et al., 2014). Arg282 is also located in the M-loop and constitutes a relevant residue for the binding of Taxol; however, this residue does not participate in relevant interactions with the ligand in epothilone-tubulin complexes. Gln292 is not in direct contact with the ligand either in taxol-or epothilone-tubulin complexes; however, it has been proposed that this residue plays a critical role in inter-dimer interactions, which are relevant for the formation and stabilization of microtubules (Verrills et al., 2003). Finally, Phe270 has been recognized as a critical residue for the interaction of the epoxide moiety of epothilones with the β-tubulin receptor (Prota et al., 2013). However, structural details about the mechanism of drug resistance induced by tubulin mutations are still lacking. Identifying the molecular mechanisms that mediate resistance to tubulin-binding agents is vital to improve the efficacy of these agents and conduct further investigation in the search for new more potent tubulin binding agents with application in anticancer therapy.
In this scenario, molecular modeling techniques provide valuable tools to increase the understanding about the epothilone-tubulin association and the role of the reported tubulin mutations in conferring resistance to epothilones. Concerning this aspect, in 2012, Shi et al. reported the use of molecular dynamics and proteinligand docking simulations to investigate the binding mode of epothilone A to wild β-tubulin and mutants at Thr274Ile and Arg282Gln, revealing fewer active residues and smaller interaction energies for the association of epothilone A with tubulin mutants compared to the native receptor (Shi et al., 2012). Also in Natarajan & Senapati (2012), Natarajan et al. performed molecular dynamics simulations to examine the effect of Thr274Ile, Arg282Gln, and Gln292Glu point mutations on the structure of uncomplexed β-tubulin, showing that mutations affect the structure and dynamics of the tubulin dimer. In addition, protein-ligand docking and MM/PBSA analysis were employed to examine the binding of taxol and epothilone A to wild and mutated β-tubulin, revealing that local changes due to mutations perturb the binding region in a significant extent, particularly in the conformation of the M-loop, and reduce the binding affinity of epothilone A and taxol toward complexation with β-tubulin (Natarajan & Senapati, 2012). However, further information is required in order to broaden the molecular understanding of drug resistance, considering the effect of tubulin mutations on the binding of epothilone B, which is the natural ligand with the best antitumor activity and the most promising lead compound in epothilone-based drug development efforts. Recent evidence has suggested that although epothilones A and B possess almost identical structures and occupy a common region in β-tubulin, they do not overlap within the binding site Source: Data retrieved from the reports of Giannakakou et al. (1997Giannakakou et al. ( , 2000 and He et al. (2001).   (Jimenez, Alderete, & Navarrete, 2015). Thus, it is valuable to obtain new structural information concerning mutation-induced resistance to epothilone B and to increase knowledge about the molecular origin of drug resistance within the epothilones series. Our current work focuses on the analysis of the structure and dynamic properties of epothilone B-tubulin complexes obtained from molecular dynamics simulations, considering the native protein dimer and four αβ-tubulin mutants at Phe270Val, Thr274Ile, Arg282Gln, and Gln292Glu. We believe that this work can be helpful to elucidate the molecular mechanistic aspects underlying drug resistance and also enhance our understanding of how epothilones interact with β-tubulin.

Computational aspects
MD simulations for epothilone B-tubulin complexes with native and mutated tubulin were performed with the NAMD software (Phillips et al., 2005) and the all-atom CHARMM27 force field (MacKerell et al., 1998).
The initial structure of epothilone B was obtained from full geometry optimization calculations without any symmetry constraint at the semi-empirical PM3 level, using the Gaussian 03 software (Frisch et al., 2004). The resulting local minima structure was employed as input geometry for all MD simulations. The required topology files and force field parameters for epothilone B were retrieved from the SwissParam web interface at http:// www.swissparam.ch, which provides parameters derived from the Merck Molecular Force Field (MMFF9) for small organic compounds (Zoete, Cuendet, Grosdidier, & Michielin, 2011). Under this approach, epothilone B is treated as single residue units, avoiding the division of the ligands into arbitrarily defined fragments. The validity of this force field parameters was assessed by comparing the final structure of epothilone B in complex with native αβ-tubulin dimer with available experimental information (Canales et al., 2014).
The crystalline structure of epothilone-bound tubulin in zinc-stabilized sheet solved at 2.9Å resolution (PDB code 1TVK) (Nettles et al., 2004) was employed to obtain the initial structure of the native αβ-tubulin dimer. This structure was subjected to mutagenesis in silico using the mutator 1.3 plugin of the VMD software considering the four point mutations at Phe270Val, Thr274Ile, Arg282Gln, and Gln292Glu. Protonation states of ionizable residues in the native and mutated tubulin structures were set to those corresponding to pH 7.0, by using PROPKA 3.0 web interface at http:// propka.ki.ku.dk/ (Olsson, Søndergaard, Rostkowski, & Jensen, 2011). PROPKA 3.0 is a fast empirical method for protein pKa prediction and rationalization, based on the effects of residue desolvation and intra-protein interactions. PROPKA 3.0 uses a PDB file as input to the server and affords high-accuracy pKa values in a fully automatic fashion (Bas, Rogers, & Jensen, 2008). According to PROPKA results, three cationic residues (Lys19, Arg276, and Arg282) and two anionic residues (Glu22 and Asp224) are located in close proximity of the ligand in the native tubulin structure. According to PROPKA 3.0, the side chain of the Glu292 residue in the Gln292Glu mutant is expected to be protonated at pH 7.0 and no changes in the total charge of the protein are predicted. In the case of Arg282Gln mutation, the positively charged Arg282 residue is replaced by the neutral Gln282 residue, which involves a −1 change in the total charge of the protein as a consequence of the point mutation.
The starting structures of epothilone B-tubulin complexes with native and mutated αβ-tubulin were obtained by placing the structure of epothilone B in the binding site occupied by epothilone A in the 1TVK model, such as the structural superposition between these ligands is maximized. This procedure accounts for the fact that the location of the binding site has been properly established, whereas the details of the association are still a matter of controversy. The final geometry of epothilone B in the tubulin complex with the native protein was compared with available experimental information to assess the validity of the simulation (Canales et al., 2014).
TIP3P water model was employed to solvate the systems in a cubic water box of 10 Å of length. The systems were neutralized with Cl − counterions to maintain the overall charge neutrality. Electrostatic interactions were evaluated using the particle-mesh Ewald method. Calculations were performed under periodic boundary conditions in an isothermal-isobaric ensemble (NPT), using NAMD Nosé-Hoover implementations. The cut-off distance for force evaluation was chosen to be 12 Å. MD simulation protocol involved: (a) 10,000 minimization steps, (b) progressive heating from 300 to 400 K, in 10 K increments of 10,000 steps, (c) 200 ps of equilibration at 400 K, (d) progressive cooling from 400 to 300 K, in 10 K increments of 10,000 steps, (e) 10 ns of equilibration at 300 K, and (f) 20 ns of NPT production dynamics at 300 K and 1 bar. Energy minimization was used to remove steric clashes, while the heating-cooling process was employed to enable the system to explore the energy surface, preventing molecules getting stuck in a confined region of the conformational space. At high temperatures, the systems are able to occupy high-energy regions of the conformational space and to pass over high-energy barriers. As temperature falls, the lower energy states become more probable in accordance with the Boltzmann distribution. Equilibration of the system was monitored by checking the protein root mean square deviation (RMSD) along the simulation run. After 10 ns of MD simulation, all systems became equilibrated.
Further 20 ns of production dynamics were used to analyze the structure and conformation of epothilone B-tubulin complexes with native and mutated αβ-tubulin.
The trajectories obtained from production dynamics were analyzed and viewed by using GROMACS and VMD utilities. The g_rmsd utility was used to obtain the RMSD. The g_rmsf application was used to calculate the Root Mean Square Fluctuation (RMSF) per residue along 20 ns trajectories. The NAMD Energy Plugin version 1.4 was employed to obtain average values for the electrostatic (E elec ) and van der Waals (E vdW ) interaction energies (kcal/mol) between the epothilone B and the residues of the native-and mutated-tubulin binding sites along the production dynamics. Tubulin residues were considered to be part of the binding domain if they locate at distances shorter than 4.0 Å from the ligand in more than 50% of the structures sampled from 20 ns of production dynamics. The volume of the binding site cavity of wild-type and mutated structures of epothilone B-tubulin complex was examined through the CASTp tool, using representative structures for each complex. CASTp (http://cast.engr.uic.edu) is an online tool that locates and measures binding pockets and voids on 3D protein structures.
3. Results and discussion 3.1. Structure and dynamics of the tubulin dimer Molecular dynamics simulations have been employed to study the structure and dynamic properties of the wildtype and mutated αβ-tubulin dimer in complex with epothilone B, aimed at providing structural information concerning the source of drug resistance conferred by tubulin mutations. Figure 3 shows the RMSDs for native αand β-tubulin and their mutants at Phe270Val, Thr274Ile, Arg282Gln, and Gln292Glu, obtained from 20 ns of production dynamics of epothilone B-tubulin complexes. RMSD values were calculated for each frame along the trajectory with respect to the average structure of wild-type β-tubulin in complex with epothilone B. According to our results, the protein backbone experiences noticeable displacements compared to wild-type tubulin as a consequence of point mutations, with the exception of the mutant Phe270Val. In this case, the mutated protein displays a very similar RMSD profile compared to wild-type tubulin, in agreement with the less pronounced drug resistance effect induced by the Phe270Val mutation (Table 1). It is relevant to note that the two subunits of the protein dimer exhibit similar displacements, revealing that the structural changes arisen from point mutations spread over an extensive protein region.
Local structural changes due to point mutations have been analyzed from the calculation of RMSFs differences per residue (ΔRMSF), considering the average structure of native αβ-tubulin as reference structure (Figure 4). ΔRMSF reports the degree of motion of each C α around its average position in the β-tubulin mutants compared to the motion observed in the structure of the native protein, thus implying that large ΔRMSF values correspond to protein regions that become more flexible due to tubulin mutations, while regions that are constrained show low ΔRMSF values. According to our results, the most noticeable fluctuations occur in regions that are relevant either for interdimer contacts or for the binding of epothilone B.
Particularly, the M loop experiences a great deal of variation and fluctuates greatly in mutants Thr274Ile, Arg282Gln, and Gln292Glu, which suggest that these mutated residues are critical for retaining the M-loop conformation in a suitable form to enable the effective complexation of the ligand. On the other hand, the mutation Phe270Val is not involved in relevant flexibility  changes in the M-loop, but in regions that participate in lateral interactions along protofilaments. Thus, the aforementioned results suggest that the RMSD and RMSF changes induced by tubulin point mutations can be related to drug resistance by altering the conformation and dynamics of protein regions involved in the effective binding of the ligand and in the assembly of stable microtubules.
A detailed structural comparison between the binding region in ligand-bound β-subunits of wild-type tubulin and its mutants is shown in Figure 5. The figure highlights the conformational changes induced by tubulin mutations in the M loop region, which adopts a more inward orientation in the wild-type and Phe270Val mutant tubulin complexes. A strong intra-protein electrostatic interaction was observed during simulation in the wild-type, and Phe270Val and Thr274Ile mutant tubulin complexes, involving the side chains of Arg276 and Asp224 residues, which remained at average distances of 3.9-4.3 Å in 84-100% of the simulated trajectory ( Figure 5). This interaction was found in a less significant extent in the mutant Gln292Glu (62%) and was not observed in the ligand-bound structure of the mutant Arg282Gln (0%). To check if there could be any induced fit in the epothilone's binding pocket, we have carried out 20 ns of molecular dynamics simulations on the structure of uncomplexed wild-type tubulin, which revealed the absence of the electrostatic interaction between Arg276 and Asp224 in all equilibrated structures, as shown in the Supplementary Information section. This suggests that the M loop must undergo a conformational reorganization upon binding of epothilone B, folding up toward the binding site. Considering that the M loop is not only relevant to the epothilone binding, but also participates in lateral interdimer contacts that confer stability to microtubules, these results suggest that the intra-protein interaction between Arg276 and Asp224 residues is critical to enable the effective Figure 5. Conformational changes in the M loop region induced by tubulin mutations. Notes: The intra-protein interaction between Asp224 and Arg276 residues is highlighted, including the corresponding average distances and occupancy percentages for this interaction.  binding of the ligand and the stabilization of the M-loop in the required conformation to favor the assembly of tubulin into stable microtubules.

Epothilone binding pocket
The binding modes of epothilone B in wild-type and mutated tubulins have been investigated from the analysis of 20 ns of molecular dynamics trajectories. Figure 6 depicts representative protein-superimposed structures showing the time evolution of the ligand's location inside the binding pocket of wild-type and mutated tubulin. According to our results, the ligand adopts different orientations within the binding cavity, thus suggesting a strong influence of point mutations on the binding pose of epothilone B. Large ligand displacements were observed in the case of Phe270Val, Thr274Ile, and Gln292Glu mutants, suggesting a looser binding with the tubulin receptor compared to wild-type tubulin. In addition, the binding pocket of tubulin mutants experience a noticeable volume increase compared to the native protein, which can be related to drug resistance by weakening the drug-receptor interactions responsible for the effective complexation of the ligand (Figure 7). Figure 8. Pictorial summary of the binding site residues identified for epothilone B-tubulin complexes with native and mutated tubulin. Note: Black boxes represent interactions that are found in more than 75% of the molecular dynamics trajectory, whereas gray boxes represent interactions that are observed in between 50 and 75% of the trajectory. A pictorial summary of the binding site residues identified for epothilone-tubulin complexes with native and mutated tubulin is shown in Figure 8. Protein residues were considered as part of the binding site if they remained at distances shorter than 4 Å from the ligand in more than 50% of the molecular dynamics trajectory. Black boxes represent interactions that are found in more than 75% of the simulation, whereas gray boxes represent interactions that are observed in between 50 and 75% of the trajectory. Overall, the simulation results indicate that point mutations render lesser polar contacts with the ligand compared to wild-type tubulin, particularly in the region comprised between Lys19 and Glu27 residues, whereas most apolar interactions are retained in epothilone B complexes with tubulin mutants. Thus, the loss of favorable polar interactions upon mutations can be related to drug resistance by weakening the epothilone B-tubulin association. These results are corroborated by the calculation of average interaction energies between epothilone B and the β-tubulin subunit, as a preliminary source of information concerning the energetic aspects of this association (Table 2). Interaction energies reveal that point mutations induce a significant decrease in the electrostatic term (E elec ) of the drugreceptor association, whereas the van der Waals interaction energy is less affected by tubulin mutations. It is relevant to note that the tubulin mutant at Gln292Glu, which induces the largest drug resistance effect within the series of mutations considered in the present study, experiences the largest binding energy difference compared to the wild-type tubulinparticularly in the electrostatic termrevealing that weaker drug-receptor interactions are related with larger fold-resistance factors.
Representative structures of epothilone B-tubulin complexes with wild-type and mutated tubulin are provided in Figure 9, showing the details of the drugreceptor interactions. In the case of the wild-type tubulin, the epoxide ring of epothilone B is surrounded by the apolar residues Leu228 and Phe270, and the methyl group at C12 is involved in hydrophobic contacts with the side chains of the Leu228 residue. In this model, Thr274 is hydrogen bound to the C7-OH group, in agreement with previous pharmacophore models (Giannakakou et al., 2000), and the C1 carbonyl is in polar contact with the Ser364 residue. The thiazole side chain is deeply buried in a tubulin region surrounded by polar residues Glu22, Asp 26, and Glu27, whereas the His227 residue forms a lateral hydrogen bond to the thiazole side chain of epothilone B, as postulated by electron crystallography studies (Nettles & Downing, 2009). In the case of the mutant Phe270Val, apolar contacts are found involving the methyl groups at C4, C6, and C12 and residues Leu228-Leu273, Ala231, and Leu361, respectively. In this case, the thiazole side chain is found to be exposed to the solvent, whereas the C1carbonyl group is hydrogen bound to Gly277. In the case of the Thr274Ile mutant, the epoxide moiety is located in the region comprised Ala231 and Phe270 residues, whereas the thiazole side chain of the ligand is surrounded by the polar residues Arg276, Asp 224, and Lys216. The C1-carbonyl group interacts with the side chain of the His227 residue, the C3-OH group forms a hydrogen bond to the backbone carbonyl group of Arg276, and the C7-OH group of the ligand interacts with the Arg359 residue. In the case of the mutant Arg282Gln, the epoxide oxygen atom of the ligand is hydrogen bound to the side chain of His227 residue, whereas the C12 methyl group interacts with the Leu228 and Phe270 residues. The C8-C12 region of the ligand interacts with the apolar residues Leu215 and Leu217, whereas the C4 methyl groups are in close contact with the side chain of Leu361. This model preserves the interaction between Thr274 residue and the C5-C7 region of the ligand, as observed in the wild-type tubulin complex, whereas the thiazole side chain is involved in polar contacts with the Arg359 and Asp26 residues. Finally, in the case of Gln292Glu, the thiazole side chain and epoxyde moiety of the ligand are considerably more exposed to the solvent compared to the structure of the wild-type tubulin complex, revealing a less effective binding of the ligand that can be related to the large drug resistance effect induced by the tubulin mutation at Gln292Glu. In this model, the C1-carbonyl oxygen of the ligand binds to the guanidinium group of the Arg276 residue, the C3-OH group interacts with the side chain of Asp224, and the C5-C7 region of the ligand interacts with Arg276 and Ser 275 residues. In addition, the C4 methyl groups are involved in apolar contacts with Leu217, Leu228, and Leu273 residues. Thus, the aforementioned results reveal that point mutations in the tubulin structure induce marked effects on the interaction networks responsible for epothilone's association, which can be related to drug resistance by weakening the affinity of the ligand toward complexation.

Conclusions
Molecular dynamic simulations on epothilone B-tubulin complexes have been employed to retrieve structural information concerning the emergence of drug resistance due to tubulin mutations at Phe270Val, Thr274Ile, Arg282Gln, and Gln292Glu. According to our results, the structure of the ligand-bound protein experiences noticeable conformational and flexibility changes upon mutations, particularly in regions involved in the binding of the ligand or in lateral contacts in protofilaments, which might induce drug resistance by altering the ability of tubulin dimers to bind to the ligand and assembly into stable microtubules. In addition, point mutations bring drastic changes in the binding pose of the ligand and interaction networks responsible for the effective drug-receptor association. A weaker drug binding is predicted for epothilone B complexes with tubulin mutants, particularly due to the loss of relevant polar contacts with the tubulin receptor. The present study provides atomic-level explanations for the drug resistance induced by tubulin mutations, which is expected to be helpful to conduct further work in the search for new more potent tubulin-binding agents that share the epothilones' mechanism of action.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
This work was supported by the FONDECYT [grant number 1120250].

Supplemental data
The supplemental data for this paper is available online at http://dx.doi.10.1080/07391102.2015.1063455.