Synthesis and Biological Evaluation of Novel Triple-Modified Colchicine Derivatives as Potent Tubulin-Targeting Anticancer Agents

Specific modifications of colchicine followed by synthesis of its analogues have been tested in vitro with the objective of lowering colchicine toxicity. Our previous studies have clearly shown the anticancer potential of double-modified colchicine derivatives in C-7 and C-10 positions. Here, a series of novel triple-modified colchicine derivatives is reported. They have been obtained following a four-step strategy. In vitro cytotoxicity of these compounds has been evaluated against four human tumor cell lines (A549, MCF-7, LoVo, and LoVo/DX). Additionally, the mode of binding of the synthesized compounds was evaluated in silico using molecular docking to a 3D structure of β-tubulin based on crystallographic data from the Protein Data Bank and homology methodology. Binding free energy estimates, binding poses, and MlogP values of the compounds were obtained. All triple-modified colchicine derivatives were shown to be active at nanomolar concentrations against three of the investigated cancer cell lines (A549, MCF-7, LoVo). Four of them also showed higher potency against tumor cells over normal cells as confirmed by their high selectivity index values. A vast majority of the synthesized derivatives exhibited several times higher cytotoxicity than colchicine, doxorubicin, and cisplatin.


Introduction
Colchicine (1, Figure 1) is one of the oldest therapeutic substances known to mankind. Although its medical properties have been known for centuries, the drug was approved for clinical use by the United States Food and Drug Administration only in 2009 [1]. Clinically, colchicine is approved and used for the treatment of familial Mediterranean fever, Behcet's disease, acute gout, chondrocalcinosis, and other types of microcrystalline arthritis [2][3][4][5][6]. Other therapeutic indications include primary biliary cirrhosis, psoriasis, amyloidosis, various forms of dermatitis, relapsing polychondritis, necrotizing vasculitis, Sweet's syndrome, leukocytoclastic vasculitis, and cardiovascular diseases, such as particular pericarditis, atrial fibrillation caused by inflammation, and ischemic episodes [7][8][9][10][11][12]. Colchicine is of Many efforts have been made to develop colchicine analogs with the goal of reducing toxicity and increasing bioavailability. Additionally, diverse structure-activity relationship studies have been performed to elucidate the structural features required for tubulin binding. Substitution of methoxyl group in the C-10 position by a thiomethyl group increases molecular stability and thiocolchicine binds to tubulin more rapidly than colchicine itself [16][17][18][19][20][21]. In 2011 Hiromitsu Takayama's research group [22] published results of their studies on C-4 halogen substituted colchicine derivatives. Some of them exhibited more potent cytotoxicity for tumor cells compared to 1. In our previous studies, we prepared double-modified, in the C-7 and C-10 positions, colchicine derivatives and evaluated their cytotoxicity [23]. Encouraged by these results, we decided to further develop the concept of diversified urethane substituents in the C-7 position and combine it with bromination in the C-4 position. We synthesized a series of triple-modified 4-bromo-N-deacetyl-7-carbamatethiocolchicines maintaining the stereochemistry of the C-7 center, which is critical for antimitotic activity [15,21].

General
All precursors for the synthesis and solvents were obtained from Sigma-Aldrich (Merck KGaA, Saint Louis, MO, USA) and were used as received without further purification. Spectral grade solvent (CDCl3) was stored over 3 Å molecular sieves for several days. Thin layer chromatography (TLC) was carried out on precoated plates (TLC silica gel 60 F254, Aluminum Plates Merck, Merck KGaA, Saint Louis, MO, USA) and spots were detected by illumination with a UV lamp. All the solvents used in flash chromatography were of HPLC grade (CHROMASOLV from Sigma-Aldrich, Merck KGaA, Saint Louis, MO, USA) and were used as received. The elemental analysis of compounds was carried out on Vario ELIII (Elementar, Langenselbold, Germany).

Spectroscopic Measurements
The 1 H, 13 C spectra were recorded on a Varian VNMR-S 400 MHz spectrometer (Varian, Inc., Palo Alto, CA, USA). 1 H-NMR measurements of 2-8 (0.07 mol dm −3 ) in CDCl3 were carried out at the operating frequency 402.64 MHz. The error of the chemical shift value was 0.01 ppm. The 13 C-NMR spectra were recorded at the operating frequency 101. 25 MHz. The error of chemical shift value was 0.1 ppm. All spectra were locked to deuterium resonance of CDCl3. The 1 H-and 13 C-NMR spectra are shown in the Supplementary Materials. Many efforts have been made to develop colchicine analogs with the goal of reducing toxicity and increasing bioavailability. Additionally, diverse structure-activity relationship studies have been performed to elucidate the structural features required for tubulin binding. Substitution of methoxyl group in the C-10 position by a thiomethyl group increases molecular stability and thiocolchicine binds to tubulin more rapidly than colchicine itself [16][17][18][19][20][21]. In 2011 Hiromitsu Takayama's research group [22] published results of their studies on C-4 halogen substituted colchicine derivatives. Some of them exhibited more potent cytotoxicity for tumor cells compared to 1. In our previous studies, we prepared double-modified, in the C-7 and C-10 positions, colchicine derivatives and evaluated their cytotoxicity [23]. Encouraged by these results, we decided to further develop the concept of diversified urethane substituents in the C-7 position and combine it with bromination in the C-4 position. We synthesized a series of triple-modified 4-bromo-N-deacetyl-7-carbamatethiocolchicines maintaining the stereochemistry of the C-7 center, which is critical for antimitotic activity [15,21].

General
All precursors for the synthesis and solvents were obtained from Sigma-Aldrich (Merck KGaA, Saint Louis, MO, USA) and were used as received without further purification. Spectral grade solvent (CDCl 3 ) was stored over 3 Å molecular sieves for several days. Thin layer chromatography (TLC) was carried out on precoated plates (TLC silica gel 60 F254, Aluminum Plates Merck, Merck KGaA, Saint Louis, MO, USA) and spots were detected by illumination with a UV lamp. All the solvents used in flash chromatography were of HPLC grade (CHROMASOLV from Sigma-Aldrich, Merck KGaA, Saint Louis, MO, USA) and were used as received. The elemental analysis of compounds was carried out on Vario ELIII (Elementar, Langenselbold, Germany).
The FT-IR spectra of 2-8 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 × 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, and 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.

Synthesis of 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 the 72 h. Reaction time was determined by TLC. The reaction was quenched with saturated aqueous Na2S2O3. The whole mixture was extracted four times with CH2Cl2, and the combined organic layers were dried over MgSO4, filtered, and evaporated under reduced pressure. The residue was purified by CombiFlash ® (EtOAc/MeOH, increasing concentration gradient) to give 2 (MW = 478.3 g/mol, Figure 2) as amorphous yellow solid with yield 95% (569 mg) [22]. 1 13

Synthesis of 3
To a mixture of 2 (500 mg, 1.01 mmol) in MeOH/water (1/1, v/v, 5 mL), the sodium methanethiolate (solution 21% in H2O, 0.79 mL, 2.1 mmol) was added. The mixture was stirred in at RT for 72 h. Reaction time was determined by TLC. After that time the reaction mixture was quenched by the addition of water (150 mL). The whole mixture was extracted four times with CH2Cl2, and the combined organic layers were dried over MgSO4, 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 ( Figure 3, MW = 494.4 g/mol) as

Synthesis of 3
To a mixture of 2 (500 mg, 1.01 mmol) in MeOH/water (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. Reaction time was determined by TLC. After that time 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 ( Figure 3, MW = 494.4 g/mol) as amorphous yellow solid amorphous yellow solid with yield 75% (388 mg) [24]. 1

Synthesis of 4
Compound 4 ( Figure 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. Reaction time was determined 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 CH2Cl2, and the combined organic layers were dried over MgSO4, filtered, and evaporated under reduced pressure. The residue was purified by CombiFlash ® (EtOAc/MeOH, increasing concentration gradient) to give 4 (MW = 452.4 g/mol) as amorphous brownish solid with yield 80% (366 mg) [25]. 1

Synthesis of 4
Compound 4 ( Figure 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. Reaction time was determined 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 (MW = 452.4 g/mol) as amorphous brownish solid with yield 80% (366 mg) [25]. 1

Synthesis of 4
Compound 4 ( Figure 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. Reaction time was determined 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 CH2Cl2, and the combined organic layers were dried over MgSO4, filtered, and evaporated under reduced pressure. The residue was purified by CombiFlash ® (EtOAc/MeOH, increasing concentration gradient) to give 4 (MW = 452.4 g/mol) as amorphous brownish solid with yield 80% (366 mg) [25]. 1     were added: Et 3 N (1 mL, 7 mmol), and triphosgene (69 mg, 0.23 mmol). The mixture was first stirred at 0 • C temperature for 20 min and then for the next 20 min at RT. After that time respective alcohol (11 mmol) was added and the mixture was stirred at RT for the next 48 h. Reaction time was determined by TLC. The solution was filtered to remove triethylamine hydrochloride. The THF was evaporated and the residue was quenched by the addition of CH 2 Cl 2 (100 mL) and was washed sequentially with a solution of HCl(aq) (0.5 M) and then with water. The organic layer was evaporated to dryness under reduced pressure and purified by CombiFlash ® (hexane/ethyl acetate, increasing concentration gradient) to give respective compounds as amorphous yellow solids (5)(6)(7)(8)(9)(10)(11)(12).  Compounds 5-12 were obtained directly from compound 4. To a solution of compound 4 (100 mg, 0.22 mmol) in tetrahydrofuran (THF, 5 mL) cooled to the 0 °C temperature, the following compounds were added: Et3N (1 mL, 7 mmol), and triphosgene (69 mg, 0.23 mmol). The mixture was first stirred at 0 °C temperature for 20 min and then for the next 20 min at RT. After that time respective alcohol (11 mmol) was added and the mixture was stirred at RT for the next 48 h. Reaction time was determined by TLC. The solution was filtered to remove triethylamine hydrochloride. The THF was evaporated and the residue was quenched by the addition of CH2Cl2 (100 mL) and was washed sequentially with a solution of HCl(aq) (0.5 M) and then with water. The organic layer was evaporated to dryness under reduced pressure and purified by CombiFlash ® (hexane/ethyl acetate, increasing concentration gradient) to give respective compounds as amorphous yellow solids (5)(6)(7)(8)(9)(10)(11)(12).

Compound 6
Amorphous yellowish brown solid, yield 47 mg, 42%, MW = 509.1 g/mol ( Figure 6). 1  Compounds 5-12 were obtained directly from compound 4. To a solution of compound 4 (100 mg, 0.22 mmol) in tetrahydrofuran (THF, 5 mL) cooled to the 0 °C temperature, the following compounds were added: Et3N (1 mL, 7 mmol), and triphosgene (69 mg, 0.23 mmol). The mixture was first stirred at 0 °C temperature for 20 min and then for the next 20 min at RT. After that time respective alcohol (11 mmol) was added and the mixture was stirred at RT for the next 48 h. Reaction time was determined by TLC. The solution was filtered to remove triethylamine hydrochloride. The THF was evaporated and the residue was quenched by the addition of CH2Cl2 (100 mL) and was washed sequentially with a solution of HCl(aq) (0.5 M) and then with water. The organic layer was evaporated to dryness under reduced pressure and purified by CombiFlash ® (hexane/ethyl acetate, increasing concentration gradient) to give respective compounds as amorphous yellow solids (5)(6)(7)(8)(9)(10)(11)(12).

Compound 11
Amorphous brownish solid, yield 54 mg, 44%. MW = 551.1 g/mol ( Figure 11). 1 13 13 3285, 2936, 1718, 1607, 1546, 1463, 1410, 1348, 1324, 1249, 1137, 1081, 1020  Compound 12 Amorphous orangish solid, yield 54 mg, 42%. MW = 584.5 g/mol ( Figure 12). 1    All cell lines were cultured during the entire experiment in a humid atmosphere at 37 °C and 5% CO2. Cells were tested for mycoplasma contamination by mycoplasma detection kit for conventional PCR: Venor GeM Classic (Minerva Biolabs GmbH, Berlin, Germany) and negative results was obtained. The procedure was repeated every year or in the case of less frequently used lines: after thawing.  . All culture media contained antibiotics: 100 U/mL penicillin and 100 µg/mL streptomycin (Polfa-Tarchomin, Warsaw, Poland). All cell lines were cultured during the entire experiment in a 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 was obtained. The procedure was repeated every year or in the case of less frequently used lines: after thawing.

The Antiproliferative Assays In Vitro
Twenty-four h before adding the tested compounds, all cell lines were seeded in 96-well plates (Sarstedt, Nümbrecht, Germany) in appropriate media with 10 4 cells per well. All cell lines were exposed to each tested agent at four different concentrations in the range 100-0.01 µg/mL for 72 h. Cells were also exposed to the reference drug cisplatin (Teva Pharmaceuticals Polska, Warsaw, Poland) and doxorubicin (Accord Healthcare Limited, Middlesex, UK). Additionally, all cell lines were exposed to DMSO (solvent used for tested compounds) (POCh, Gliwice, Poland) at concentrations corresponding to these present in tested agents' dilutions. After 72 h sulforhodamine B assay (SRB) was performed [26].

SRB
After 72 h of incubation with the tested compounds, cells were fixed in situ by gently adding 50 µL per well of cold 50% trichloroacetic acid TCA (POCh, Gliwice, Poland) and were incubated at 4 • C for one hour. Following, wells were washed four times with water and air dried. Next, 50 µL of 0.1% solution of sulforhodamine B (Sigma-Aldrich, Merck KGaA, Saint Louis, MO, USA) in 1% acetic acid (POCh, Gliwice, Poland) were added to each well and plates were incubated at room temperature for 0.5 h. After incubation time, unbound dye was removed by washing plates four times with 1% acetic acid whereas stain bound to cells was solubilized with 10 mM Tris base (Sigma-Aldrich, Merck KGaA, Saint Louis, MO, USA). Absorbance of each solution was read at Synergy H4 Hybrid Multi-Mode Microplate Reader (BioTek Instruments, Inc., Winooski, VT, USA) at the 540 nm wavelength.
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-2014, http://www.cheburator. nevozhay.com, freely available) for each experiment [27]. Compounds at each concentration were tested in triplicates in single experiment and each experiment was repeated at least three times independently. Results are summarized in Table 1. 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).

Molecular Docking Simulations
A combination of different computational methods was used to explore ligand-tubulin interactions. The ligands structures first minimized then fully optimized based on the RHF/cc-pVDZ level of theory in GAMESS-US version 2010-10-01. Since there is no crystal structure available for human βI tubulin (TBB5_HUMAN), we obtained its sequence from UniProt (ID: Q13509). We used the tubulin structure 1SA0.pdb as a template to construct the homology model for βI tubulin using MOE2015. We then docked the colchicine library to the protein using the Autodock4 program [28] under flexible ligand and rigid receptor conditions (Table 2). To verify the computed binding free energies, we employed two accurate but computationally expensive methods, namely MM/PBSA and MM/GBSA. The IC 50 value is defined as the concentration of a compound at which 50% growth inhibition is observed. The SI (selectivity index) 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 (resistance index) indicates how many times a resistant subline is chemo-resistant 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 drug-resistance.

Chemistry
In the present paper, we report the synthesis and spectroscopic analysis of a series of eight novel triple-modified colchicine derivatives (5)(6)(7)(8)(9)(10)(11)(12). As shown in Scheme 1, 4-bromocolchicine (2), 4-bromothiocolchicine (3), and 4-bromo-N-deacetylthiocolchicine (4) were prepared according to previously described methods [22,24,25]. Compound 4 became the starting material for the synthesis of triple-modified colchicine derivatives. Compounds 5-12 were readily available from 4 by treatment with triphosgene as an activating reagent in the presence of triethylamine and the excess of respective alcohol or glycol in dry THF with yields from 30% to 48%. The structures of all products 2-12 were determined using the ESI-MS, FT-IR, 1 H and 13 C-NMR methods. Spectra are shown in Supplementary  Materials (Figures S1-S23) and the results are also presented in the Materials and Methods section. All the spectroscopic and mass spectrometry data presented in Supplementary Materials confirm the structures of the studied compounds. In the 13 C-NMR spectra of the 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 introduction of thiometyl group in C-10 positions shifts of the 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 the 4 signals for the acetyl group, which were observed at 170.0 and 22.9 ppm in compound 3, had disappeared. In carbamates (5)(6)(7)(8)(9)(10)(11)(12), shifts of the signal for the carbamate carbon atom were observed in the range 153.5-156.0 ppm.

In Vitro Determination of Drug-Induced Inhibition of Human Cancer Cell Line Growth
The eight triple-modified colchicine derivatives (5)(6)(7)(8)(9)(10)(11)(12), three other colchicine derivatives (2)(3)(4), and starting material (1) were evaluated for their in vitro antiproliferative effect on normal and cancerous cells. Each compound was tested on four human cancer cell lines, including one cell line displaying various levels of drug resistance, namely human lung adenocarcinoma (A549), human breast adenocarcinoma (MCF-7), human colon adenocarcinoma cell line (LoVo) and doxorubicin-resistant subline (LoVo/DX). The antiproliferative effect was also studied on the normal murine embryonic fibroblast cell line (BALB/3T3) for better evaluation of cytotoxic activity of the compounds studied. The mean values IC50 ± SD of the tested compounds are summarized in Table 1.
To evaluate the agents' activity against the cells with MDR (multidrug resistance) phenotype, one drug resistant cancer cell line, i.e., LoVo/DX, was tested and the resistance index values (RI) were calculated (see Table 1). The RI values indicate how many times more resistant the subline was in In the 13 C-NMR spectra of the 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 introduction of thiometyl group in C-10 positions shifts of the 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 the 4 signals for the acetyl group, which were observed at 170.0 and 22.9 ppm in compound 3, had disappeared. In carbamates (5-12), shifts of the signal for the carbamate carbon atom were observed in the range 153.5-156.0 ppm.

In Vitro Determination of Drug-Induced Inhibition of Human Cancer Cell Line Growth
The eight triple-modified colchicine derivatives (5)(6)(7)(8)(9)(10)(11)(12), three other colchicine derivatives (2)(3)(4), and starting material (1) were evaluated for their in vitro antiproliferative effect on normal and cancerous cells. Each compound was tested on four human cancer cell lines, including one cell line displaying various levels of drug resistance, namely human lung adenocarcinoma (A549), human breast adenocarcinoma (MCF-7), human colon adenocarcinoma cell line (LoVo) and doxorubicin-resistant subline (LoVo/DX). The antiproliferative effect was also studied on the normal murine embryonic fibroblast cell line (BALB/3T3) for better evaluation of cytotoxic activity of the compounds studied. The mean values IC 50 ± SD of the tested compounds are summarized in Table 1. To evaluate the agents' activity against the cells with MDR (multidrug resistance) phenotype, one drug resistant cancer cell line, i.e., LoVo/DX, was tested and the resistance index values (RI) were calculated (see Table 1).
The RI values indicate how many times more resistant the subline was in comparison to its parental cell line.
The data presented in Table 1 show that all of the studied compounds including unmodified colchicine and the two reference anticancer drugs less effectively inhibited the proliferation of the resistant LoVo/DX cell line than the sensitive LoVo cell line. However, all of the colchicine derivatives The values of the selectivity index (SI) were calculated to evaluate the toxicity of compounds against normal cells. The SI values calculated for A549, MCF-7, and LoVo cell lines were especially high (SI ≥ 4.3) for compounds 6-9 (except SI of 8 for the MCF-7 cell line). These values were much higher than the SI values of commonly used drugs, such as doxorubicin and cisplatin. High SI values result from large differences between the cytotoxicity against cancer and normal cells and this means that cancer cells will be killed at a higher rate than normal ones. Moreover, compounds 6-9 showed very high antiproliferative activity against A549, MCF-7, and LoVo cell lines, which is expressed by very low IC 50 values (IC 50 = 0.010-0.030; 0.013-0.027 (without 8), and 0.007-0.018 for A549, MCF-7, and LoVo, respectively).

Molecular Docking: In Silico Determination of Drug-Induced Inhibition of βI Tubulin
βI tubulin, one of the subunits of microtubules in the cytoskeleton of eukaryotic cells, is a well-known and well-studied target for chemotherapeutic drugs selected to inhibit the growth and proliferation of cancer cells. The computational evaluation of binding energies between drug candidates and βI tubulin performed by docking is a fast and inexpensive prediction method to investigate and rank the ability to arrest cancer cell proliferation by new compounds, which inhibit microtubule assembly. Here, the binding energies of new colchicine derivatives were calculated by docking the referred compounds to the colchicine-binding pocket of βI tubulin.
A combination of different computational methods was used to explore ligand-tubulin interactions. The ligand structures were first energy-minimized, then fully optimized based on the RHF/cc-pVDZ level of theory implemented in the software package GAMESS-US, version 2010-10-01. Since there is no crystal structure for human βI tubulin (UniProt ID: P07437) available in the Protein Data Bank (PDB), the bovine tubulin structure 1SA0.pdb was used as a template to construct the homology model for human βI tubulin using the software package MOE2015. Autodock4 program under flexible ligand and rigid receptor conditions was also used to dock the small library of investigated colchicine derivatives to the target protein structure (see Table 2). AutoDock4 software is designed to predict how ligands bind to a receptor of a known 3D structure and it consists of two main modules: (1) autodock, which performs the docking of the ligand to a set of grids describing the target protein; and (2) auto grid, which pre-calculates these grids.
Based on our in silico calculations, 11, 7, 6, and 3 show the strongest binding energies of −9.25, −9.20, −8.99, and −8.60 kcal/mol, respectively. Note that two Met 259 and Lys 352 residues in the colchicine binding site of βI tubulin interact with C-20 (sidechain donor) and oxygen of carbonyl (sidechain acceptor) on ring C of the new colchicine derivatives, respectively. The diagrams depicting these interactions can be found in Table 2.
The evaluation of IC 50 values in a cell-based assay involves interactions of the tested compounds with all tubulin expressed by the cell and βI tubulin isotype is not the only isotype of tubulin expressed, although it is expected to be one of the most abundant ones. Even though both normal and cancer cells in humans contain the same tubulin isotypes, their expression levels and distribution of tubulin isotypes in each of the cell lines are different. In particular, the most abundant isotype in most tumors is βI isotype (TUBB) followed by, βIVb (TUBB4B), βIIa (TUBB2A), βV (TUBB6), and βIII (TUBB3), with 47%, 38%, 8.9%, 3.1%, and 2.2%, respectively, and then βIVa (TUBB4), βIIb (TUBB2B), and βVI (TUBB1) with levels below 0.5% of the total expressed β tubulins. Since The expression levels of tubulin isotypes in each individual cell lines are unique, the binding energies for each of these isotypes would differ and affect the overall response to cytotoxic agents and make the computational prediction fairly complex [29]. To quantify the assumption that isotype expression levels correlate with cytotoxicity of the compounds acting on these isotypes of tubulin, the same docking simulation method was applied between the novel colchicine derivatives and βIIa (UniProt ID: Q13885), βIIb (UniProt ID: Q9BVA1), βIII (UniProt ID: Q13509), βIVa (UniProt ID: P04350), βIVb (UniProt ID: P68371), and βVI (UniProt ID: Q9H4B7) tubulin isotypes.

Linear Regression with Two Independent Variables
Following the docking simulations, in order to determine the level of correlation between experiment and computational simulations, we have calculated the linear regression coefficients between experimental values of IC 50 and computational prediction for different β tubulin isotypes. To have a more realistic correlation coefficient, a logarithmic value of solubility of the compounds, the Moriguchi octanol-water partition coefficients (MLogP), were calculated using a software package called ADMET Predictor 8.0 (ADMET Predictor, Simulations Plus, Lancaster, CA, USA) and taken into account for predicting the biopharmaceutic properties like permeability and the understanding of transport mechanisms of the drugs in vivo [30]. The partition coefficient is also a useful factor in estimating and comparing the distribution of the drugs within the cells, organs, and the body.
A very good correlation of 0.66 and 0.84 involving log IC 50 for LoVo and LoVo/DX cell lines, respectively, was the result of linear regression analysis with two independent variables, namely the compounds' MlogP values and the other being the binding free energies of our compounds and the tubulin βI isotype. These two cell lines have also a good correlation of 0.72 and 0.80 with the binding free energies of the compounds to the βIVb isotype, respectively. MlogP and the experimental IC 50 values of 4-bromocolchicine based series are listed in Table 3, Supplementary Materials Figures S24  and S25. For the A549 cell line, the regression coefficients found were 0.43, 0.43, and 0.51 with binding free energies for βI, βIVb, and βIVa isotypes, respectively. The reported regression coefficient for the A549 cell line is still acceptable, while its value for the MCF-7 and BALB/3T3 cell lines is very low (see Table 3). The reason for the latter discrepancy between computation and experiment may be due to additional biological factors, such as the upregulation of MDR proteins that act as efflux pumps and prevent the drugs from exerting their cytotoxic action. Another possibility could involve off-target interactions. To further investigate the accuracy of the reported linear regressions between our computational results obtained using the docking method and the experimentally generated IC 50 values, more accurate but also more time-consuming computational methods were applied to calculate the binding free energies between the novel colchicine derivatives and β tubulin isotypes. This analysis is described below.

MM/PBSA and MM/GBSA: In Silico Determination of Drug-Induced Inhibition of β Tubulin Isotypes
The molecular operating environment (MOE) was used to build a 3D model for the α-βI tubulin heterodimer [31]. Human tubulin protein sequences were obtained from UniProt [32]. The sequence corresponding to the gene TUBA1A (UniProt ID: Q71U36) was chosen as a reference sequence for human α-tubulin whereas gene TUBB associated to βI isoform (UniProt ID: P07437) was chosen for human β-tubulin. Homology modeling was performed using MOE 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. As a template, the crystallographic structure of α-βIIB tubulin isotype complexed with colchicine (PDB ID: 1SA0) was used [33]. 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 eventually protonated at neutral pH and minimized using a MOE's built-in protocol.
In order to equilibrate our model and get representative conformations, molecular dynamics (MDs) simulations were run using Amber14 [34]. Amber's antechamber utility was applied to generate MD parameters-e.g., partial charges, force constants, etc.-for the four cofactors from the Gasteiger charge method. Amber's tleap program was applied to solvate the system in TIP3P water. Minimization of the structure was carried out in two steps, both using the steepest descent and conjugate gradient methods successively. First, minimization was done during 2 ps on solvent atoms only, by restraining the protein-ligand complex. Next, minimization was run without the restraint during 10 ps. 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. The root-mean-square deviation (RMSD) of both the entire tubulin structure and the colchicine binding site were found to plateau after 40 ns (see Figure 13).
Minimization of the structure was carried out in two steps, both using the steepest descent and conjugate gradient methods successively. First, minimization was done during 2 ps on solvent atoms only, by restraining the protein-ligand complex. Next, minimization was run without the restraint during 10 ps. 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. The root-mean-square deviation (RMSD) of both the entire tubulin structure and the colchicine binding site were found to plateau after 40 ns (see Figure 13).  Clustering analysis of the last 30 ns of the generated MD trajectory was carried out using Amber's CPPTRAJ program [35] to identify representative conformations of the tubulin dimer. Clustering was done via a hierarchical agglomerative using the RMSD of atoms in the colchicine binding site as a metric. An RMSD cutoff of 1.0 Å was set to differentiate the clusters. From clustering analysis, three representative structures of the tubulin dimer were found. The structures were further used as a rigid target for the screening of 4-bromocolchicine derivatives.
Docking of 4-bromocolchicine derivatives was performed using the AutoDock Vina [36] program, which makes use of an iterated local search global optimizer as a searching method. The Vina scoring function combines aspects from knowledge-based and empirical potentials. Tested on the same set used for Autodock4, Vina enabled to identify the correct binding pose in 78% of the cases as compared to 53% for Autodock4 [37]. For our docking simulations, a cubic box with size 30.0 Å centered at the center of mass of the bound colchicine was used. 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 three tubulin representative structures obtained from clustering of the MD trajectory. Every generated pose was energy-minimized using Amber14 by keeping the protein fixed and was re-scored using the Vina software. For each compound/protein-structure pair, the pose with the best score was identified and used as an initial configuration for molecular mechanics Poisson-Boltzmann surface area MM/PBSA computations.
The MM/PBSA technique was used to calculate the free energy associated with the binding of 4-bromocolchicine derivatives [38,39]. This method combines molecular mechanics with continuum solvation models. The binding free energy is estimated as where ∆E MM − T∆S can be regarded as the change in the free energy of the system in vacuum (gas phase); it includes the change in the molecular mechanics energy ∆E MM = E MM bound − E MM unbound and the change in the conformational entropy ∆S due to the binding. Since our goal was to compare the binding free energy of similar compounds derived from colchicine, ∆S was not estimated when calculating ∆G bind as each compound is assumed to provide comparable ∆S values. ∆G solv stands for the difference of solvation free energies due to the binding, which is given as where every term on the right-hand side is given as the sum of polar and nonpolar contributions. The polar parts are obtained by solving the Poisson-Boltzmann (PB) equation or by using the generalized-born (GB) model-resulting in the MM/GBSA method, whereas the nonpolar terms are estimated from a linear relation to the solvent accessible surface area (SASA). The values of ∆E MM and ∆G solv are generally computed as ensemble averages requiring a short MD trajectory of the solvated complexed system as input of the MM/PBSA method. In the present case, a 1 ns-duration MD trajectory was run in TIP3P water using Amber14, for every top pose generated at the end of the docking step. The MM/PBSA and MM/GBSA calculations were performed on a subset of 200 frames collected at regular time intervals from the trajectory. For PB calculations, an ionic strength of 0.0 nM (istrng = 0.0) and a solvent probe radius of 1.6 Å (prbrad = 1.6) were used. For GB calculations, the igb parameter was set to 5 that corresponds to a modified GB model equivalent to model II in reference [40]. For each of the tested compounds, the best PB and GB score out of the three trajectories associated with the three representative structures of the tubulin dimer were collected and reported in Table 4.

Linear Regression with Two Independent Variables
Following the MM/PBSA simulations, we have calculated the linear regression coefficients between experiment and computational simulations of different β tubulin isotypes. The result of linear regression with two independent variables between the binding free energies of βIVb isotypes, MlogP and IC 50 for the LoVo and LoVo/DX cell lines is that very good correlation coefficients were found, namely 0.78 and 0.67, respectively. Also, reasonably good correlation of 0.65 and 0.54 was obtained with the binding free energies of βI isotype for these two cell lines, respectively. These correlation coefficients are also consistent with the correlation coefficients calculated using the docking binding energies and these cell lines except for LoVo/DX and βI which were lower but still acceptable. For the A549 cell line, the regression coefficients found were 0.62, 0.58, and 0.48 for the binding free energies with βIII, βIVb, and βI isotypes, respectively. The reported regression coefficient for the A549 cell line was improved a little bit by using the more expensive method but the MM/PBSA method could not improve the correlation coefficient values for the MCF-7 and BALB/3T3 cell lines and they are still very low (Table 5). The LoVo and LoVo/DX cell line correlation with βI and βIVb isotypes 0.60, 0.66, 0.65, and 0.71 respectively. The result also shows a good correlation of 0.7 between binding energies for βV tubulin and LoVo/DX cell line. For, the regression coefficients of 0.55 and 0.51 with free energies of βI, the isotypes which showed for the A549 cell line were in the same range with docking and MM/GBSA results. There were no improvements for the correlation values for the MCF-7 and BALB/3T3 cell lines with the MM/GBSA simulations (Table 6). The ranges of correlation coefficients between the IC 50 values of all the cell lines and the corresponding binding free energies calculated with three different methods are the same. Reassuringly, the best correlation values from all three methods used were found for βI and βIVb tubulin isotypes, which are reported to be most abundant in these types of cells.

Discussion
We synthesized a series of novel triple-modified 4-bromothiocolchicine derivatives (5-12) as well as 4-bromocolchicine (2), 4-bromothiocolchine (3) and 4-bromo-N-deacetylthiocolcicine (4) and evaluated their biological activity according to the in vitro antiproliferative tests as well as in silico studies. Biological activity was evaluated on four human cancer cell lines and a normal murine embryonic fibroblast cell line. The results of our study have clearly showed that the cytotoxicity of almost all colchicine derivatives (2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12) is higher than the corresponding cytotoxicity of commonly used cytostatic agents-doxorubicin and cisplatin against A549, MCF-7, LoVo (except 4 and 12) and LoVo/DX cancer cell lines. The majority of the derivatives also exhibit a higher cytotoxicity than unmodified colchicine. Particularly noteworthy are the compounds 6-9, which show very high antiproliferative activity against A549, MCF-7 and LoVo cell lines, that is expressed by very low IC 50 values at nanomolar concetrations (IC 50 = 0.010-0.030; 0.013-0.027 (without 8) and 0.007-0.018 µM for A549, MCF-7 and LoVo, respectively). Moreover, compounds 6-9 were demonstrated to be less toxic to normal murine fibroblast cells than the currently used anticancer drugs, such as cisplatin and doxorubicin, which is confirmed by particularly high selectivity index (SI) values calculated for A549, MCF-7, and LoVo cell lines are (SI ≥ 4.3, except SI of eight for MCF-7 cell line). The toxicity is a major challenge in designing a potential colchicine-based drug candidate. High SI values result from large differences between the cytotoxicity against cancer and normal cells and this means that cancer cells exposed to the same concentration of the compound will be killed at a higher rate than normal ones, which might lead to a potential colchicine-based drug candidate.
Based on our in silico calculations, 11, 7, 6, and 3 show the strongest binding energies of −9.25, −9.20, −8.99 and −8.60 kcal/mol, respectively. Two Met 259 and Lys 352 residues in the colchicine binding site of βI tubulin interact with C-20 (sidechain donor) and oxygen of carbonyl (sidechain acceptor) on ring C of the new colchicine derivatives, respectively. In order to determine the level of correlation between experiment and computational simulations, we have calculated the linear regression coefficients between IC 50 values and the binding free energies involving the compounds and tubulin isotypes as well as the MlogP coefficients for these compounds. A very good correlation of 0.66 and 0.84 with log IC 50 for LoVo and LoVo/DX cell lines, respectively, has been found. For the A549 cell line, the regression coefficient found is 0.43, still acceptable, while its value for the MCF-7 and BALB/3T3 cell lines is very low. This may be explained by a number of additional effects taking place in living cells compared to the computational simulations that focus only on the binding mode of the compounds to the target. Specifically, off-target interactions involving efflux pumps with different affinities for the individual compounds may explain the observed partial correlation between IC 50 values and binding free energies. Additionally, differences in the solubility values and membrane permeability may have to be accounted for when ranking the various compounds in biological assays and comparing them to computational predictions based on binding affinity alone. We have additionally supported these findings with calculations using two very accurate methods of calculating the binding energies of ligands to proteins, namely MM/PBSA and MM/GBSA. The results using all three in silico approaches have been found consistent.
In short, these results confirm that a suitable chemical modification of colchicine aimed for an improved binding affinity to human tubulin ßI or ßIVb isotypes is a promising approach to finding highly biologically active and less toxic compounds. Some of the obtained compounds are suitable candidates for further tests (ex vivo, in vivo). The synthesis of new colchicine derivatives with diverse substituents is also a next step in developing structure-activity relationship (SAR) of colchicine-binding site inhibitors.