New Series of Double-Modified Colchicine Derivatives: Synthesis, Cytotoxic Effect and Molecular Docking

Colchicine is a well-known anticancer compound showing antimitotic effect on cells. Its high cytotoxic activity against different cancer cell lines has been demonstrated many times. In this paper we report the syntheses and spectroscopic analyses of novel colchicine derivatives obtained by structural modifications at C7 (carbon-nitrogen single bond) and C10 (methylamino group) positions. All the obtained compounds have been tested in vitro to determine their cytotoxicity toward A549, MCF-7, LoVo, LoVo/DX, and BALB/3T3 cell lines. The majority of obtained derivatives exhibited higher cytotoxicity than colchicine, doxorubicin and cisplatin against the tested cancerous cell lines. Additionally, most of the presented derivatives were able to overcome the resistance of LoVo/DX cells. Additionally, their mode of binding to β-tubulin was evaluated in silico. Molecular docking studies showed that apart from the initial amides 1 and 2, compound 14, which had the best antiproliferative activity (IC50 = 0.1–1.6 nM), stood out also in terms of its predicted binding energy and probably binds best into the active site of βI-tubulin isotype.


Synthesis of 2 and 3
In vitro antiproliferative activity Molecular docking studies LC-MS chromatograms and mass spectra, 1 H NMR and 13 C NMR spectra of compounds 2-17. Figure S1. The LC-MS chromatogram and mass spectra of 2. Figure S2. The 1 H NMR spectrum of 2 in CDCl3. Figure S3. The 13 C NMR spectrum of 2 in CDCl3. Figure S4. The LC-MS chromatogram and mass spectra of 3. Figure S5. The 1 H NMR spectrum of 3 in CDCl3. Figure S6. The 13 C NMR spectrum of 3 in CDCl3. Figure S7. The LC-MS chromatogram and mass spectra of 4. Figure S8. The LC-MS chromatogram and mass spectra of 5. Figure S9. The 1 H NMR spectrum of 5 in CDCl3. Figure S10. The 13 C NMR spectrum of 5 in CDCl3. Figure S11. The LC-MS chromatogram and mass spectra of 6. Figure S12. The 1 H NMR spectrum of 6 in CDCl3. Figure S13. The 13 C NMR spectrum of 6 in CDCl3. Figure S14. The LC-MS chromatogram and mass spectra of 7. Figure S15. The 1 H NMR spectrum of 7 in CDCl3. Figure S16. The 13 C NMR spectrum of 7 in CDCl3. Figure S17. The LC-MS chromatogram and mass spectra of 8. Figure S18. The 1 H NMR spectrum of 8 in CDCl3. Figure S19. The 13 C NMR spectrum of 8 in CDCl3. Figure S20. The LC-MS chromatogram and mass spectra of 9. Figure S21. The 1 H NMR spectrum of 9 in CDCl3. Figure S22. The 13 C NMR spectrum of 9 in CDCl3. Figure S23. The LC-MS chromatogram and mass spectra of 10. Figure S24. The 1 H NMR spectrum of 10 in CDCl3. Figure S25. The LC-MS chromatogram and mass spectra of 11. Figure S26. The 1 H NMR spectrum of 11 in CD2Cl2. Figure S27. The 13 C NMR spectrum of 11 in CD2Cl2. Figure S28. The LC-MS chromatogram and mass spectra of 12. Figure S29. The 1 H NMR spectrum of 12 in CDCl3. Figure S30. The 13 C NMR spectrum of 12 in CDCl3. Figure S31. The LC-MS chromatogram and mass spectra of 13. Figure S32. The 1 H NMR spectrum of 13 in CDCl3. Figure S33. The 13 C NMR spectrum of 13 in CDCl3. Figure S34. The LC-MS chromatogram and mass spectra of 14. Figure S35. The 1 H NMR spectrum of 14 in CDCl3. Figure S36. The 13 C NMR spectrum of 14 in CDCl3. Figure S37. The LC-MS chromatogram and mass spectra of 15. Figure S38. The 1 H NMR spectrum of 15 in CDCl3. Figure S39. The 13 C NMR spectrum of 15 in CDCl3. Figure S40. The LC-MS chromatogram and mass spectra of 16. Figure S41. The 1 H NMR spectrum of 16 in CD3CN. Figure S42. The 13 C NMR spectrum of 16 in CD3CN. Figure S43. The LC-MS chromatogram and mass spectra of 17. Figure S44. The 1 H NMR spectrum of 17 in CDCl3. Figure S45. The 13 C NMR spectrum of 17 in CDCl3.
In 3D illustrations the interacting residues predicted from pairwise per-residue binding free energy decomposition calculations (E < -

General
All solvents, substrates and reagents were obtained from TriMen Chemicals (Poland) or Sigma Aldrich and were used without further purification. Spectral grade solvents were stored over 3 Å molecular sieves for several days. TLC analysis was performed using pre-coated glass plates (0.2 mm thickness, GF-254, pore size 60 Å) from Agela Technologies and spots were visualized by UV-light. Products were purified by flash chromatography using high-purity grade silica gel (pore size 60 Å, 230 -400 mesh particle size, 40-63 μm particle size) from SiliCycle Inc. Solvents were removed using a rotary evaporator.
Electrospray ionization (ESI) mass spectra were obtained on a Waters Alliance 2695 separation module with a PDA 2996 UV detector and Waters Micromass ZQ 2000 mass detector equipped with Restek Ultra Biphenyl 50 x 3 mm, 3 μm column eluted with 0.3 mL/min flow of 3-100% gradient (over 6 min) of acetonitrile in water.

Compound 2
To a solution of 1 (1.0 equiv.) in EtOH a methylamine (solution 33% in EtOH, 10.0 equiv.) was added. The mixture was stirred at reflux for 24 h and then concentrated under reduced pressure to dryness. The residue was purified using column flash chromatography (silica gel; DCM/MeOH, 40/1 v/v) and next lyophilized from dioxane to give the pure product 2 as a yellow solid with a yield of 80%. To a solution of compound 2 (1.0 equiv.) in dioxane, 2M HCl (10.0 equiv.) was added and the mixture was stirred at reflux. Reaction progress was monitored by LC-MS. Then the reaction mixture was neutralized with 4M NaOH to pH~10 and extracted four times with EtOAc. The organic layers were combined, washed with brine, dried over Na2SO4, filtered off and evaporated under reduced pressure. The residue was purified using column flash chromatography (silica gel; DCM/MeOH, 20/1 v/v) and next lyophilized from dioxane to give the pure product 3 as a yellow solid with a yield of 73%. ESI

Cell viability assays
Twenty-four hours before adding the tested compounds, all cell lines were seeded in 384-well plates (Sarstedt, Germany) in appropriate media with 1.0x10 3 cells per well for A549 cell line, 1.5x10 3 cells per well for MCF-7 cell line and 2.0x10 3 cells per well for LoVo, LoVo/DX and BALB/3T3 cell lines. All cell lines were exposed to each tested agent at different concentrations in the range 100-0.001 μg/mL for 72 h. The cells were also exposed to the reference drug cisplatin (Teva Pharmaceuticals, Poland) and doxorubicin (Accord Healthcare Limited, UK). Additionally, all cell lines were exposed to DMSO (solvent used for tested compounds) (POCh, Poland) at concentrations corresponding to those present in tested agents dilutions. After 72 h sulforhodamine B assay (SRB) was performed [2].

SRB
After 72 h of incubation with the tested compounds, the cells were fixed in situ by gently adding of 30 μL per well of cold 50% trichloroacetic acid TCA (POCh, Poland) and were incubated at room temperature for one hour. Then the wells were washed four times with water and air dried. Next, 25 μL of 0.1% solution of sulforhodamine B (Sigma-Aldrich, Germany) in 1% acetic acid (POCh, 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 70μl of 10 mM Tris base (Sigma-Aldrich, Germany). Absorbance of each solution was read off from a Synergy H4 Hybrid Multi-Mode Microplate Reader (BioTek Instruments, USA) at the 540 nm wavelength.
Results are presented as mean IC50 (concentration of the tested compound, that inhibits cell proliferation by 50%) ± standard deviation. IC50 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 [3]. Compounds at each concentration were tested in triplicates in single experiment and each experiment was repeated at least three times independently.

Ligand preparation
The ligand structures were prepared using Ligprep from the Schrödinger suite [4]. Conformations and tautomeric states were assigned to the ligands by following the ligand preparation protocol implemented in Schrödinger suite with default settings. LigPrep generates variants of the same ligand with different tautomeric, stereochemical, and ionization properties.

Tubulin model
The tubulin crystal structures available in the PDB are those for bovine protein. The bovine tubulin structure of tubulin (PDB ID: 1SA0) [5] was used as a template to construct the homology model of human αβ-tubulin isotypes (βI (UniProtKb: P07437), which is the most abundant isotype in most tumors using the Molecular Operating Environment (MOE) software package [6]. The sequence corresponding to the gene TUBA1A (UniProt ID: Q71U36) was chosen as a reference sequence for human tubulin, whereas the 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.

Molecular dynamics simulations
The missing hydrogens for heavy atoms were added using the tLEAP module of AMBER 14 with the AMBER14SB force field [7]. The protonation states of all ionizable residues were determined at pH = 7 using the MOE program [6]. Each protein model was solvated in a 12 Å box of TIP3P water. In order to bring the salt concentration to the physiological value of 0.15 M, 93 Na + ions and 57 Cl − ions were added. Minimization of the structure was carried out in two steps, using the steepest descent and conjugate gradient methods successively. At first, minimization was made in 2 ps on solvent atoms only, by restraining the protein-ligand complex. Next, minimization was run without the restraint in 10 ps. After minimization, the molecular dynamics (MD) simulations were carried out in three steps: heating, density equilibration, and production. At first, each solvated system was heated to 298 K for 50 ps, with weak restraints on all backbone atoms. Next, density equilibration was carried out for 50 ps of constant pressure equilibration at 298 K, with weak restraints. Finally, MD production runs were performed for all systems for 70 ns. The root-mean-square deviation (RMSD) of both the entire tubulin structure and the colchicine-binding site were found to reach a plateau after 40 ns. Clustering analysis of the last 30 ns of the generated MD trajectory was carried out using the Amber's CPPTRAJ program [8] to identify representative conformations of the tubulin dimer. Clustering was made via the hierarchical agglomerative approach 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. On the basis of the clustering analysis, three representative structures of the tubulin dimer were found. The docking was performed on all the three representative structures and the one with the highest docking score was selected, which was the largest cluster (about 70% of the simulation) conformation of the tubulin structure. During the modeling, the 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.

Docking simulations
We used the AutoDock Vina [9] and DOCK (http://dock.compbio.ucsf.edu/) programs to predict the binding pose of the ligands under flexible ligand and rigid receptor conditions. Dockbox package was used to facilitate preparation of docking inputs and post-processing of the docking results [10]. Docking simulations performed with a cubic box (size 30.0 Å) were centered at the middle of binding pockets and the docking was run separately on the tubulin structure. Every generated pose was energy-minimized using Amber14 by keeping the protein fixed and was re-scored using the MOE's GBVI/WSA dG scoring function [6]. No constraints were applied in the docking studies. For each compound/protein-structure pair, the pose with the best score was identified and used as an initial configuration for molecular mechanics Gibbs-Boltzmann surface area MM/GBSA computations.
Binding energy calculations using MM/GBSA method The MM/GBSA technique is used to calculate the free energy associated with the binding of double modified colchicine derivatives [11]. This method combines molecular mechanics with continuum solvation models. We performed MM-GBSA integrated in Amber. The binding free energy is estimated as: where G is the average free energy of the complex, protein, and ligand, are calculated according to the equation: where EMM are determined with the SANDER program and represent the internal energy (bond, angle, and dihedral), van der Waals and electrostatic interactions (see equation (3)). TS is the entropy contribution estimated using normal mode (nmode) analysis.
The solvation free energy can be calculated as the sum of polar and nonpolar contributions. The polar parts are obtained by using the generalized-born (GB) model-resulting in the MM/GBSA method, whereas the nonpolar terms are estimated from a linear relation (equation 4) to the solvent accessible surface area (SASA).
In the present study, a 2 ns-duration MD trajectory was run in TIP3P water using Amber14, for every top pose generated at the end of the docking step. It is worth noting that to assess the performance of MM/GBSA methodology [12], we evaluated the prediction accuracy of this method by various simulation protocols including 1 ns MD production calculations using PDBbind data set. Too long an MD simulation could be prejudicious for the overall success of the MM/GBSA method. According to this study and the common practice to calculate binding energies using MM/GBSA, we have decided to run MD production simulation for 2 ns. The 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, which corresponds to a modified GB model equivalent to model II in reference.

Free energy decomposition analysis
The interaction between inhibitors and each residue were computed using the MM/GBSA decomposition process by the mm_gbsa program in AMBER 12.0. The binding interaction of each inhibitorresidue pair includes three energy terms: van der Waals contribution (ΔEvdw), electrostatic contribution (ΔEele), and solvation contribution (ΔGGB + ΔGSA), in which ΔEvdw and ΔEele are van der Waals and electrostatic interactions between the inhibitor and each protein residue that could be computed by the Sander program in AMBER 12.0. The polar contribution of desolvation (ΔGGB) was calculated using the generalized Born (GB) model. The nonpolar contribution of desolvation (ΔGSA) was computed based on SASA determined with the ICOSA method. All energy components were calculated using 300 snapshots