Synthesis, biological evaluation and molecular docking studies of new amides of 4-bromothiocolchicine as anticancer agents

Colchicine is the major alkaloid isolated from the plant Colchicum autumnale , which shows strong therapeutic effects towards different types of cancer. However, due to the toxicity of colchicine towards normal cells its application is limited. To address this issue we synthesized a series of seven triple-modified 4-bromothiocolchicine analogues with amide moieties. These novel derivatives were active in the nanomolar range against several different cancer cell lines and primary acute lymphoblastic leukemia cells, specifically compounds: 5 – 9 against primary ALL-5 (IC 50 = 5.3 – 14 nM), 5, 7– 9 against A549 (IC 50 = 10 nM), 5, 7 – 9 against MCF-7 (IC 50 = 11 nM), 5 – 9 against LoVo (IC 50 = 7 – 12 nM), and 5, 7 – 9 against LoVo/DX (IC 50 = 48 – 87


Introduction
Natural products and their derivatives play an important role in many modern therapies and as a foundation for the development of new drugs. For example, for cancer drugs developed during the period from the 1940s to the end of 2014, of 175 small molecules approved, 49% were natural products or direct derivatives of natural products [1][2][3][4][5] . Thus, chemical modification of biologically active compounds of natural origin is a proven and efficient approach to novel drug development. Among them, colchicine, an alkaloid and secondary metabolite extracted from plants of the genus Colchicum e.g. Colchicum autumnale, is one of the oldest therapeutic substances known to mankind 6 . Colchicine is still used for the treatment of familial Mediterranean fever, Behcet's disease, acute gout, chondrocalcinosis, and other types of microcrystalline arthritis because of its anti-inflammatory properties [7][8][9][10][11][12][13][14][15][16][17] . Its primary target is β-tubulin where it forms complexes with tubulin dimers and copolymerizes into microtubule structures, suppressing microtubule dynamics and disrupting key cellular processes like mitotic spindle formation [18][19][20] . Compounds which bind to the colchicine site on tubulin and inhibit microtubules (colchicine site inhibitors, CSIs) have been a rich source of new drugs particularly those with tumor vascular disrupting ability [21][22][23] . The parent molecule colchicine has shown very high anticancer activity in vitro, but its clinical potential for human tumors is limited because of its relatively high toxicity, which results from its accumulation in the gastrointestinal tract, as well as neurotoxicity 5,23,24 . Many efforts have been made to optimize the structure of colchicine to generate less toxic and more bioavailable analogues 7, . We have previously reported studies on 4-bromothiocolchicine derivatives with diversified carbamate substituents in the C-7 position 30 . In order to further increase anti-cancer activity and decrease toxicity towards normal cells, in this study we investigated the impact of replacing the carbamate moiety with an amide. Previous independent work studied C-4 halogen substituted colchicine derivatives and 4-chlorocolchicine derivatives bearing an amide moiety in the C-7 position 29 47 . 4-Chlorocolcichine and some of the double-modified derivatives with trifluoroacetyl or propionyl substituents exhibited strong antitumor activities over broad effective dosage ranges in vivo, but their metabolic stabilities were poorer than that of colchicine. In the present study, we synthesized 4-bromothiocolchicine derivatives bearing a thiomethyl group in the C-10 position anticipated to increase molecular stability. Thiocolchicines do not undergo photochemical isomerization, while colchicine is 48 , as well as harsher conditions can be use for deacetylation at the C(7) of thiocolchicines without any undesirables changes in their structures. 49 We report here the synthesis and biological activity in normal and cancer cells of a series of novel triple-modified derivatives, with six different amide substituents and one urea moiety in the C-7 position 50 .
In the 13 C-NMR spectra of 2, a resonance for the C-4 carbon atom of the A aromatic ring was observed at 113.5 ppm, while in 1 it was observed at 107.3 ppm. After the substitution of a thiomethyl group in the C-10 position, shifts of the 13C NMR signal for the C-20 carbon atom in compound 3 were observed at 15.2 ppm, while in unmodified 1 as well as 2 shifts of the signal for the C-20 carbon atom were observed in the range 56.1-56.5 ppm. In the 13 C-NMR spectra of 4, signals for the acetyl group, which were observed at 170.0 and 22.9 ppm in compound 3, had disappeared. In amides (5-10), shifts of the signal for the amide carbon atom were observed in the range 167.0-176.9 ppm while in urea (11) the shift was observed at 155.9 ppm.
*Inhibition of proliferation did not exceed 50% at the highest concentration tested of 10 µM.
SI was calculated for each compound using the formula: SI = IC 50 for normal cell line BALB/3T3 / IC 50 for respective cancerous cell line. A beneficial SI > 1.0 indicates a drug with efficacy against tumor cells greater than the toxicity against normal cells.
The RI indicates how many times a resistant subline is chemoresistant relative to its parental cell line. The RI was calculated for each compound using the formula: RI = IC 50 for LoVoDX / IC 50 for LoVo cell line. When RI is 0-2, the cells are sensitive to the compound tested, RI in the range 2-10 means that the cell shows moderate sensitivity to the drug tested, RI above 10 indicates strong drugresistance. This is of significance since primary ALL cells are not transformed, immortalized or adapted to culture and are thus anticipated to exhibit properties more closely representative of the initial cancer cells than conventional cell lines 54 . The viability curves for ALL-5 cells are presented in Supplementary  Fig. S21.
The selectivity index (SI) is a ratio that measures the window between cytotoxicity and anticancer activity by dividing the given IC 50 values for the normal BALB/3T3 cell line and respective IC 50 values for cancerous cell lines. Higher SI values are desirable since they reflect efficacy with less toxicity. The precursor colchicine derivatives (compounds 2-4) had in general high SI values for each of the non-MDR cell types (Table 1). Data for A549, MCF-7 and LoVo were reported previously (REF) and data for ALL-5 are new and showed a similar trend. The triple-modified colchicine derivatives showed relatively lower SI values on primary ALL-5, A549, MCF-7 and LoVo cell lines (SI = 0.8 -2.0). The exceptions were compound 10 with extremely low SI values (SI = 0.1 -0.2), and compounds 5, 6, 7 and 9 with moderate SI values on primary ALL-5 cells (SI = 2.2, 4.6, 2.8 and 2.2 respectively) and LoVo cell line (SI = 5.3 (6), 2.4 (7) and 2.9 (9)). Importantly however, conventional anticancer agents like doxorubicin and cisplatin typically showed lower SI values than many of the colchicine analogs for a given cell type (Table 1). Overall, the results of Table 1 indicate that triple modified colchicine derivatives have greater potency than the parent colchicine for several different types of cancer cell.
To evaluate activity against the cells with an MDR (multidrug resistance) phenotype, one drug resistant cancer cell line, i.e. LoVo/DX, was tested and the resistance index (RI) values were calculated (Table 1), as described in Materials and Methods. Most of the derivatives were not able to overcome the drug resistance of the LoVo/DX cell line, with RI values ranging from 5.6 to 22.2. One exception was compound 10, with an RI value of 1.8, but this desirable property was offset by a poor SI value (Table 1).

Effects on cell death and mitotic arrest in primary ALL-5 and MCF-7 cells
Since colchicine (1) and its analogs showed favorable activity towards primary ALL-5 cells further mechanistic investigation was conducted. First, we studied whether the compounds induced apoptotic cell death, by assessing DNA content and fragmentation via flow cytometry. We employed propidium iodide staining in order to determine DNA content, and cells with sub-G1 (<2N) DNA were considered to have undergone cell death. ALL-5 cells were subjected to treatment with 1, six analogs characterized by IC 50 values lower than for parent molecule 1, and DX, at concentrations of 5 x IC 50 values (Table  1), for 24, 48 or 72 h. Cells treated with vehicle (0.1% DMSO) for equivalent incubation times served as control. The full set of representative flow cytometric data is included in Supplementary  Fig. S22A. A graphical representation of sub-G1 DNA content derived from the mean of three such experiments is shown in Fig.  1A. For all of the compounds, statistically significant DNA fragmentation was observed after 48 h of treatment and further increased after 72 h. However, for analogs 5, 7, 8 and 9, cell death was more rapidly induced, with significant DNA fragmentation observed as early as 24 h. Since parent molecule 1 is a well-known microtubule targeting agent (MTA) 20 , it was of interest to determine whether the analogs retained this property. Therefore, cells with 4N DNA content, reflecting mitotic arrest, were quantitated after treatment ( Fig. 1B; Supplementary Fig.  S22A). In a population of untreated, asynchronous ALL-5 cells, about 10% have 4N DNA (Fig. 1B). When treated with starting compound 1 and its analogs this amount significantly increased after 24 h (exception 2) and further after 48 h (exception 3 and 7), to subsequently drop down after 72 h due to the growing pool of dead cells (Fig. 1A). Overall, levels of mitotic arrest were low, with maximum extents of about 20%. This is consistent with previous work where it has been shown that in response to microtubule destabilizing dugs, primary ALL-5 cells tend to be susceptible not only in M phase but also in G1 phase. Thus, the above observations are in good agreement with our previously published data, where other two microtubules destabilizers, namely vincristine and eribulin, failed to cause mitotic arrest in primary ALL-5 cells but instead induced death directly in G1 phase when treated with 100 nM or higher concentrations of these drugs 55 . In contrast to ALL-5 cells, MCF-7 breast cancer cells showed only low levels of DNA fragmentation after treatment, but much higher levels of mitotic arrest, with a more rapid onset (Fig. 1C, 1D). This is consistent with MCF-7 cells: having a shorter doubling time versus ALL-5 cells (24 h versus 60 h, respectively) 54 ; exhibiting susceptibility in M phase not interphase; and lacking caspase-3 56 , and thus dying in a manner which does not involve DNA fragmentation. Similar results have been previously reported by us for a series of thiocolchicine urethanes 33 .

A.
B.

The effect of colchicine and 5 on PARP cleavage in primary ALL-5 cells
Colchicine (1) and compound 5 as the most active triplemodified colchicine derivative based on IC 50 value (Table 1) were selected for further studies examining poly (ADP-ribose) polymerase (PARP) cleavage by immunoblotting as an additional marker of apoptotic cell death. The immunoblots are shown in Fig ). Treatment with DX was performed as a positive control, and loss of 116 kDa PARP and generation of the 85 kDa product was observed, as previously reported in both ALL-5 57,58 and HeLa cells 59 . Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was employed as a loading control ( Fig. 2A, lower panel).

Molecular docking studies
Prediction of the binding modes of the 4-bromothiocolchicine derivatives with different tubulin isotypes, namely, -I, -IIa, -III, -IVb and -VI was carried out by first running molecular dynamics (MD) simulations for each isotype during 70 ns in order to generate representative structure(s) of the proteins. Identification of such structures was done by performing RMSDbased clustering analysis on all the structures generated over the last 30 ns of each MD trajectory (see section 4.9). As a second step, docking was performed on every representative structure using Autodock Vina and DOCK 6.5. All the poses generated by both programs were then rescored using AutoDock Vina's scoring function. Table 2 shows the top binding pose of colchicine and the derivatives as predicted by Vina scoring function for -I. Active residues involving non-hydrophobic interactions with the ligand are also specified for each compound. Subsequent to docking, the best Vina scores were collected for each isotype-compound pair. Two-variable linear regressions using the above-mentioned scores as the first variable and the Moriguchi octanol-water partition coefficient (MLogP) as the second variable was carried out in an attempt to fit experimental pIC50 values. Table 3 depicts, for each isotype, the R 2 values computed from such linear regressions using the pIC50 values reported in Table 1. The best correlations were obtained for primary acute lymphoblastic leukemia cells where R 2 values were found to range between 0.359 and 0.525. The highest R 2 (0.525) was produced by using Vina scores computed on the most common tubulin isotype I (see also Fig. 3). The second best yet poor correlation was obtained in the case of MCF-7 cells, where R 2 values were found around 0.2 and the best R 2 value was again obtained in the case of the I isotype. In the cell lines where good correlation was not found, especially LoVo/DX (which is a drug resistant type), additional factors may be playing a role in the compounds' efficacy such as P-glycoprotein based drug resistance, off-target interactions of the compounds or membrane permeability issues which need to be addressed in future studies. Table 2 Binding location and orientation of colchicine and its derivatives on 1 tubulin as predicted by the Autodock Vina software andactive residues interacting with each ligand via hydrogen bonding, acidic, basic or pi-stacking interactions. Interaction diagrams show nearby residues within a 4.0 Å radius . Compounds with no specified active residues (N/A) were found to interact only via hydrophobic interactions.

Conclusions
Seven novel triple-modified colchicine derivatives (5-11) were synthesized with moderate yields. All of the compounds were evaluated for their anti-proliferative activity against several standard cancer cell lines as well as primary acute lymphoblastic leukemia cells. We conclude that the introduction of different amide moieties in the C-7 position had an impact on the biological activity of colchicine derivatives. Several triple modified derivatives were more potent than colchicine when tested against primary cancer cells and established cancer cell lines (Table 1), in many cases by a factor of 10 or greater. These results strongly support the further development of triplemodified bromo-analogues of colchicine as anti-cancer agents. Further studies revealed that both colchicine (1) and its most potent derivatives induced characteristics of apoptotic cell death and induced contrasting effects on the progression of cell cycle in MCF-7 vs. ALL-5 cells. While 1 and studied analogues caused mitotic arrest in MCF-7 cells, such effect was not observed in ALL-5, showing premises of interphase cell death. Based on our molecular docking studies we found reasonably close correlation between the calculated binding affinity for the most comment isotype of tubulin βI with the experimental IC 50 values for most of the cases, in particular for acute lymphoblastic leukemia cells. Other factors may be playing a role in the cell lines, where good correlation was not found. Those include p-glycoprotein based drug resistance, off-target interactions of the compounds or membrane permeability issues. Several of the newly synthesized compounds show excellent potential as a novel cancer chemotherapeutics. They show favorable activity comparing not only to their parent molecule, colchicine but also to standard chemotherapeutic agents such as doxorubicin and cisplatin. However, more studies are warranted in order to increase a clinical potential of this promising class of compounds.

General procedures
All precursors and solvents for the synthesis were obtained from Sigma Aldrich (Merck KGaA, Saint Louis, MO, USA) and were used without further purification. CDCl 3 spectral grade solvent was stored over 3 Å molecular sieves for several days. TLC was performed on precoated plates (TLC silica gel 60 F254, Aluminium Plates Merck, Merck KGaA, Saint Louis, MO, USA) visualized by illumination with an UV lamp. HPLC grade solvents (without further purification) were used for flash chromatography (CHROMASOLV from Sigma Aldrich, Merck KGaA, Saint Louis, MO, USA). The elemental analysis of compounds was performed on Vario ELIII (Elementar, Langenselbold, Germany).

Spectroscopic measurements
The 1 H, 13  The FT-IR spectra of 2-11 in the mid infrared region were recorded in KBr. The spectra were taken with an IFS 113v FT-IR spectrophotometer (Bruker, Karlsruhe, Germany) equipped with a DTGS detector; resolution 2 cm -1 , NSS = 64. The Happ-Genzel apodization function was used.
The ESI (Electrospray Ionisation) mass spectra were recorded also on a Waters/Micromass (Waters Corporation, Manchester, UK) ZQ mass spectrometer equipped with a Harvard Apparatus syringe pump. The samples were prepared in dry acetonitrile (5 x 10-5 mol dm -3 ). The sample was infused into the ESI source using a Harvard pump at a flow rate of 20 ml min -1 . The ESI source potentials were: capillary 3 kV, lens 0.5 kV, extractor 4 V. The standard ESI mass spectra were recorded at the cone voltages: 10 and 30 V. The source temperature was 120 °C and the desolvation temperature was 300 °C. Nitrogen was used as the nebulizing and desolvation gas at flow-rates of 100 dm 3 h -1 . Mass spectra were acquired in the positive ion detection mode with unit mass resolution at a step of 1 m/z unit. The mass range for ESI experiments was from m/z = 100 to m/z = 1000, as well as from m/z = 200 to m/z = 1500.

. 3 . 1 . S y n t h e s i s o f -b r o m o c o l c h i c i n e ( 2 )
A mixture of N-bromosuccinimide (NBS, 279 mg, 1.57 mmol) and 1 (500 mg, 1.25 mmol) in acetonitrile was stirred at RT under nitrogen atmosphere for 72 h. The progress of the reaction was monitored by TLC. The reaction was quenched with saturated aqueous Na 2 S 2 O 3 . The whole mixture was extracted four times with CH 2 Cl 2 , and the combined organic layers were dried over MgSO 4 , filtered, and evaporated under reduced pressure. The residue was purified by CombiFlash ® (EtOAc/MeOH, increasing concentration gradient) to give 2 as an amorphous yellow solid with yield 95% (569 mg). 1

. 3 . 2 . S y n t h e s i s o f -b r o m o t h i o c o l c h i c i n e ( 3 )
To a mixture of 2 (500 mg, 1.01 mmol) in MeOH/H 2 O (1/1, v/v, 5 mL), the sodium methanethiolate (solution 21% in H 2 O, 0.79 mL, 2.1 mmol) was added. The mixture was stirred in at RT for 72 h. The progress of the reaction was controlled by TLC. After 72 h the reaction mixture was quenched by the addition of water (150 mL). The whole mixture was extracted four times with CH 2 Cl 2 , and the combined organic layers were dried over MgSO 4 , filtered, and evaporated under reduced pressure. The residue was purified by CombiFlash ® (hexane/EtOAc (1/1), then EtOAc/MeOH, increasing concentration gradient) to give 3 as amorphous yellow solid with yield 75% (388 mg). 1

. 3 . 3 . S y n t h e s i s o f -b r o m o d e a c e t y l o t h i o c o l c i c i n e ( )
Compound 4 was prepared from 3 by hydrolysis with 2 N HCl. To a solution of compound 3 (500 mg, 1.01 mmol) in MeOH (3 mL), the 2 N HCl solution (5 mL) was added. The mixture was stirred at 90 • C for 72 h. The progress of the reaction was controlled by TLC. After that time the reaction mixture was quenched by the addition of water (100 mL). The whole mixture was extracted four times with CH 2 Cl 2 , and the combined organic layers were dried over MgSO 4 , filtered, and evaporated under reduced pressure. The residue was purified by CombiFlash® (EtOAc/MeOH, increasing concentration gradient) to give 4 as amorphous brownish solid with yield 80% (366 mg). 1

. 3 . . G e n e r a l p r o c e d u r e f o r t h e s y n t h e s i s o f c o l c h i c i n e d e r i v a t i v e s ( 5 -1 1 )
Compounds 5-11 were obtained directly from 4. To a solution of compound 4 (100 mg, 0.22 mmol) in tetrahydrofuran (THF, 5 mL) cooled to 0 • C, the following compounds were added: Et 3 N (2 mL, 14 mmol) and DMAP (catalytic amount). The mixture was first stirred at 0 • C for a few minutes and then the solution of respective acyl chloride (5 -10) or diethylcarbamoyl chloride (11) in THF (0.66 mmol in 2.5 ml) was added dropwise. The mixture was stirred at RT for the next 24 h. The solution was filtered to remove triethylamine hydrochloride. The THF was evaporated and the residue was purified by CombiFlash® (hexane/ethyl acetate, increasing concentration gradient) to give respective compounds as amorphous yellow solids with yield from 18 % to 83 % (5-11). The 1 H and 13 C NMR spectra of the compound 5-11 are shown in the Supplementary Materials.

Cell lines and culturing conditions
Primary ALL-5 cells were derived from the bone marrow of a 37-year old patient as previously described. 54,60 Although these cells can be cultured up to 6 months with no obvious change in their properties 54 , in the present study they were exclusively used at low passage for all experiments, and are thus referred to as primary cells. Primary ALL-5 cells were routinely maintained at 37 °C in a humidified 5% CO 2 incubator in IMDM Modified (HyClone) media supplemented with 10 μg mL -1 cholesterol, 6 mg mL -1 human serum albumin, 2 mM L-glutamine, 2% v/v amphotericin-B/penicillin/streptomycin, 1 μg mL -1 insulin, 200 μg mL -1 apo-transferrin, and 50 μM β-mercaptoethanol, and were subcultured to maintain a density of 1·3 x 10 6 cells mL -1 . Human MCF-7 mammary gland adenocarcinoma cells originally isolated from a 69-year-old Caucasian woman with several characteristics of differentiated mammary epithelium were cultured in Eagle's Minimum Essential Medium (EMEM)  antibiotics: 100 U/ml penicillin and 100 μg/ml streptomycin (Polfa-Tarchomin, Warsaw, Poland). All cell lines were cultured during entire experiment in humid atmosphere at 37 °C and 5% CO 2 . Cells were tested for mycoplasma contamination by mycoplasma detection kit for conventional PCR: Venor GeM Classic (Minerva Biolabs GmbH, Berlin, Germany) and negative results were obtained. The procedure is repeated every year or in the case of less frequently used lines after thawing. Results are presented as mean IC 50 (concentration of the tested compound, that inhibits cell proliferation by 50%) ± standard deviation. IC 50 values were calculated in Cheburator 0.4, Dmitry Nevozhay software (version 1.2.0 software by Dmitry Nevozhay, 2004 http://www.cheburator.nevozhay.com, freely available) for each experiment 53 . Compounds at each concentration were tested in triplicate in individual experiment and each experiment was repeated at least three times independently.

. 5 . 2 . M T T a s s a y
3-(4,5-dimethylthiazol-2-yl)-2,5-dimethylthiazol-2-yl)-2,5diphenyltetrazolium bromide (MTT)-based assay was used to evaluate the effect of drugs on the viability of primary ALL-5 cells 62,63 . Cells (10 5 /well) in 100 µL of complete IMDM Modified medium were seeded in 96-well plates (TPP, Switzerland) and treated with drugs at concentrations up to 10 µM for 120 h with control cells receiving vehicle (0.1% DMSO) alone. After treatment, 10 μL of MTT solution (5 mg/mL) was added to each well, and the plate was incubated at 37°C for 24 h in a humidified 5% CO 2 incubator. Then 100 μL of 10% SDS in 0.01 M HCl was added to each well and the plate was incubated at 37 °C for a further 24 h. Absorbance was recorded at 540 nm using a BioTek Plate Reader. Inhibition of formation of colored MTT formazan was taken as an index of cytotoxicity activity. IC 50 values were determined by non-linear regression analysis using GraphPad Prism 6 for Windows (GraphPad Software). Selectivity index (SI) was calculated by dividing the IC 50 value for BALB/3T3 cells by the IC 50 value for individual cancer cell lines, and resistance index (RI) was calculated by dividing the IC 50 for LoVo/DX cells by the IC 50 for LoVo cells. The Resistance Index (RI) was defined as the ratio of IC 50 for a given compound calculated for resistant cell line to that measured for its parental drug sensitive cell line (Table 1).

DNA content analysis
ALL-5 (1.5 x 10 6 ) and MCF-7 (0.2 x cells 10 6 ) were seeded in 100 mm Petri dishes (Corning, NY) and incubated in the presence of vehicle (0.1% DMSO) or compounds, at concentrations specified in the text, for 24, 48 or 72 h at 37°C in a humidified 5% CO 2 incubator. Cells were then washed with 1 mL phosphate-buffer saline (PBS), fixed with 1-3 mL of 70% ice-cold ethanol and stored at 4 °C prior to flow cytometric analysis. Cells were centrifuged, treated with 500 μL propidium iodide/RNase Staining buffer (BD Biosciences, San Jose, CA, USA) and stored in the dark for 1 h at RT. The stained cells were subjected to a FacsAria IIIu Flow Cytometer (BD Biosciences, San Jose, CA, USA) and data were analyzed using FlowJo software.

Statistical analysis
Unpaired t test with Welch's correction was performed for the significance and p values of <0.05 were considered significant.

Homology modeling and molecular dynamics
A combination of different computational methods was used to investigate in silico the binding modes and binding energies of the colchicine derivatives for different isotypes of the human tubulin dimer. The 3D structures of the human isotypes were built from homology modeling by using the crystallographic structure of the bovine -IIb tubulin isotype complexed with colchicine as a template (PDB ID: 1SA0) 64 . The dimer models generated were made of the IA isotype for -tubulin (UniProt ID: Q71U36) and different isotypes of the -tubulin including I (UniProt ID: P07437), IIa (UniProt ID: Q13885), III (UniProt ID: Q13509), IVb (UniProt ID: P68371), VI (UniProt ID: Q9H4B7) 65 . The latter isotypes of tubulin are the most commonly expressed ones in cancer cells. Homology modeling was performed using the Molecular Operating Environment (MOE) 66 by setting the number of generated models to 10 and by selecting the final model based on MOE's generalized Born/volume integral (GB/VI) scoring function. During the modeling, cofactors including GTP, GDP, colchicine, and the magnesium ion located at the interface between and monomers were kept as part of the environment and included in the refinement step. The final model was protonated at neutral pH and minimized using a built-in protocol in the MOE package.
Subsequent to homology modeling, molecular dynamics (MD) simulations were run on the generated model of each isotype using Amber14 67 . MD parameters, e.g., partial charges and force constants were set via the ff14SB force field for the protein while parameters for GTP, GDP and colchicine were obtained from the GAFF force field. The complex system was further solvated in TIP3P water using the Amber's tLEaP program. Energy minimization of the structure was carried out in two steps, both using the steepest descent and conjugate gradient methods successively. First, 2000 cycles of energy minimization were performed on solvent atoms only, by restraining the proteinligand complex. Next, minimization was run without the restraint for 5000 cycles. The structure was then equilibrated in an NVT ensemble during 20 ps and in an NPT ensemble during 40 ps setting the temperature to 298 K and the pressure to 1 bar. Finally, MD production was run for 70 ns (see Figure S.23).
Clustering analysis of the last 30 ns of the generated MD trajectories was carried out using Amber's CPPTRAJ program 68 to find the representative conformations of each tubulin isotype. Clustering was done via a hierarchical agglomerative algorithm using the RMSD of atoms in the colchicine binding site as a metric. An RMSD cut-off of 1.0 Å was set to differentiate the clusters. For each isotype, the representative structure of each cluster was used as a rigid target for the screening of the colchicine derivatives.

Docking simulations
Docking of the 4-bromothiocolchicine derivatives was performed using two different docking software packages, namely, AutoDock Vina 69 and DOCK6.5 70 . Since those programs are based on different methods for ligand placement and scoring, using both programs simultaneously normally increases the chance of generating the correct ligand pose 71 . Vina includes an iterated local search global optimizer as a searching method and a combination of knowledge-based and empirical potentials as a scoring function [6]. On the other hand, DOCK6.5 is based on the anchor-and-grow algorithm to generate ligand poses and makes use of a force-field-based potential to score them [7]. For our docking simulations, a cubic box with size 30.0 Å centered at the center of mass of the bound colchicine was considered. All cofactors, namely, GTP, GDP, colchicine, and the magnesium ion were removed during docking while the target was kept rigid. For every compound, docking was run separately on each of the tubulin representative structures obtained from clustering. The number of generated ligand poses was set to 10 for both AutoDock Vina and DOCK6.5, meaning that a maximum of 20 ligand poses was produced for every compound/protein structure pair. The ligand poses were eventually rescored using AutoDock Vina's scoring function. For every derivative, the pose with the best Vina score over all representative structures of each tubulin isotype was kept for further analysis, especially to investigate the correlation with experimental pIC 50 values. Besides Vina scores, the Moriguchi octanol-water partition coefficient (MLogP) of every compound was calculated using the ADMET Predictor 8.0 package (ADMET Predictor, Simulations Plus, Lancaster, CA, USA). Both Vina scores and MlogP values were used as inputs to build a two-variable linear regression model for every tubulin isotype.