Introduction

Tubulin, a structural protein \(\alpha/\beta\) heterodimer, is the building block of microtubules (MT). Each subunit \(\alpha\) and \(\beta\) is expressed in the cell by multiple genes which lead to different isotypes (up to 8 for the human \(\beta\) subunit). (ref) Most of the time, isotypes are highly conserved during evolution suggesting a functional importance. Tubulins are arranged in a head-to-tail fashion to form 13 protofilaments constituting the MTs. Mts exhibit a dynamic instability which is essential for their biological activity since MTs are cytoskeletal components that are critically involved in several key cellular process such as mitosis, intracellular transport or cell motility.(ref)

As a consequence of this key role, tubulin is a well known target for some of the most successful anti-tumor drugs, such as taxanes, vinca alkaloïds or colchicin (refs). After cellular uptake, microtubule network will be the target to distrupt chromatid separation thus preventing the meiosis. Cytotoxic agents will interfere with MT dynamic by inhibiting the polymerisation or depolymerisation of dimers of MT, leading to inhibition of the cancer cell division and proliferation (ref). Paclitaxel (Taxol), for example, binds to MTs and stabilizes the MTs already formed preventing the depolymerisation (ref). At the opposite, vinca alkaloïds (vinblastine) and colchicin prevent the polymerisation (ref).

Despite the fact that most clinical tubulin binding drugs target the \(\beta\) subunit of the dimer, an important issue is left unanswered regarding the tubulin isotypes which should be specifically targeted in cancer chemotherapy (ref). Indeed, recent studies show that glioblastomas, the most common primary brain tumors, exhibit an abnormal \(\beta\) subunit expression with an proliferation of the class III \(\beta\)-tubulin, which is associated with highly malignant tumor phenotype. This over expression is also related with a paclitaxel resistance (ref). The selectivity of this particular isotype could be a key to a successful drug design process.

Until now, most of the studies on the tubulin do not take into account the isotype, mainly because the major changes occurs in the C-terminal tails (CTT) of the subunits (ref). These CTTs are still missing parts of the crystal structures found in the databases, thus making it complicated to study. Considering recent studies, CTTs exhibit primoridial roles for the MTs functionality, such as formation of MT, interaction with microtubule associated proteins and kinesin, acting as a lasso (ref). It is also a hotspot for posttranslational modifications (polyglutamylation, tyrosination, phosphorylation, etc) (ref).

We present a modelisation process of two complete human tubulins of \(\alpha1\beta1\) and \(\alpha1\beta3\) isotypes. Furthermore, we investigated the C-termini of \(\beta\) subunit structural organisation and interaction with the core parts of the tubulin, using multiple molecular dynamics simulations. Our results show... To be continued.

Material and methods

Human complete tubulin modelization

For this study, we have modeled two different tubulin isotypes, the \(\alpha1\beta1\) and the \(\alpha1\beta3\). Each subunit was modeled independently with the same template, a double dimer in complex with the stathmin (PDB entry : 1SA0), chosen for high identity and similarity (cf table \ref{table_identity}), most recently reviewed and more complete among all crystal structures (despite being a curved conformation). It is also the best structure of tubulin already integrated in the template database of the Modeller program which we used for the modelization process (v9.14). (The tubulin type template already present in the Modeller v9.14 database are 1SA0, 1TUB and FtZ like domain models, the prokaryote counterpart of the tubulin). Modeller process is based on the classic comparative modeling method which consists of four sequential steps : template selection, template-target alignment, model building and finally model evalutaion. For the template selection, Modeller use the 3D-template matching method (pariwise alignement with known structure only, present in an associated database). The next step consist of an alignment between the query sequence and the structure template. In our case, the similarity belong to the best possible template alignement (indentity \(>\) 75%, similarity \(>\) 90%), so it do not need any manual rearrangements regarding gaps. The model building depend on the global alignment quality and the overall completeness but in modeller the method used is a modelization by satisfaction of spatial restraints. The template is used to generate a set of spatial contraints and restraints to serve as a guide for the new model construction with the alignment. In the case of Modeller, the main constraints are dihedral angles and distance, assuming aligned residues in the target are submitted to the same contraints than in the template. An two steps optimisation is then performed by adding stereochemical constaints in the model (bond lengths, bond angles, non-bonded contacts, etc) in the first place followed by a minimization of all the violations of the restraints. For the tails building, Modeller proceed by doing an Ab initio method which can be seen as a protein folding problem for small peptide. It is a conformational search in a given environment guided by an energy function. Loops in the model are treated with the same method. The last step of the modelization process is the model evaluation. Modeller use multiple built-in scoring functions : molPDF, DOPE and GA341 scores (ref). MolPDF and DOPE score are energies, then the lower they are, the better is the model, and GA341 take a value between 0.0 and 1.0 (1.0 is the best value). (ajouter tableau des scores pour les modèles réalisés ?). The human tubulin sequences were obtained on the UniProtKB website and can be found with the accession number Q71U36 for the \(\alpha\) subunit, Q9H4B7 and Q13509 for the \(\beta1\) and \(\beta3\) subunits respectively. After the modelization of each part independently, we reassembled them as a straight tubulin dimer with a spatial alignment on the best resolution crystal for this conformation (PDB entry : 1JFF, the PDB we used in our previous study).

\label{table_identity}

Molecular dynamics simulation

After this step, we submitted all models previously obtained to classic molecular dynamics simulation (five models for each isotypes modeled) in order to relax the structure and to explore the structural variability of the C-termini newly modeled. The simulation were made with Gromacs 5.0.4 with the OPLS-AA forcefield in, periodic boundary condition. The procedure consist of two steps of minimization, the first one is an in vacuo during 1000 steps without any constraints using steepest descent algorithm, followed by the addition of a water box of 2 nm around the tubulin filled with TIP3P molecule model and neutralization of the system with randomly added NaCl ions in the box while maintaining a 150mM concentration. The total system contain around 13 700 atoms for the tubulin dimer, 136 000 water molecules, and 300 ions (approx. 60% Na 40% Cl due to the high negative charge of the C-termini).

The second minimization step is done under the same set of parameters during 5000 steps to prevent any water clashes with the tubulin. The next stage before production phase consists of a two-steps equilibration procedure, an NVT equilibration followed by an NPT equilibration. Each stage lasted 100 ps, with an integration step of 2 fs. Temperature was fixed at 300 K using velocity rescale method, all bonds were constrained with the LINCS algorithm, electrostatic interactions were computed with the Particular Mesh Ewald method. For the pressure coupling during the NPT equilibration, we used Parinello-Rahman method at the value of 1 atm.

The production run were then performed with the same set of parameters and algorithms during 100 ns. Trajectories were saved every 10 ps, and the analysis were made on the last 95 ns, considering the first 5 ps as equilibration period.

Trajectory analysis

After obtaining the simulation data, we processed them with the Gromacs 5.0.4 utilities for basics observations and with VMD 1.9 for the graphical part. Most of the more complexes analysis were made with homemade scripts/programs.

For the clustering part, it was a two step process. Due to the high flexibility of the C-termini (mostly the one of the \(\beta\) subunit), which are the parts of interest in this study, classicals methods can not be applied directly. We made multiples analysis regarding the orientation, radius of gyration, contacts with the rest of the dimer to have a preselection of subgroups. These subgroups were then submitted independently to classical hierarchical clustering method (Gromos algorithm) in order to obtain representatives structures of differents states encountered during the simulations. Cutoff for the clustering was 1.5 Å.

Results

Model building

When the modelisation step is over, Modeller provides multiples metrics to asses the models generated. In our case, we took the 5 best rated models by Modeller for each subunit. DOPE score and MolPDF score can be found in table [score_a1, score_b1, score_b3]. As said in the Methods part, these scores are energies, so we want them to be the lowest possible. From the table, we can see that rankings are not always the same for MolPDF and DOPE score, but the changes are minimal for all the models generated. We can deduce that each models should be correct regarding known parts of the tubulin (i.e. not the CTT), without much variability except side chains, the main differences coming from intrinsically disordered region forming the CTT which is quite small compared to the rest of the molecule.

|c|c|c|c|c|c| & model 1 & model 2 & model 3 & model 4 & model 5
DOPE score & -52451.50 & -52086.40 & -52350.84 & -52158.03 & -52277.13
MolPDF score & 3712.60 & 3448.56 & 3578.11 & 3502.33 & 3566.32

\label{score_a1}

|c|c|c|c|c|c| & model 1 & model 2 & model 3 & model 4 & model 5
DOPE score & -50642.49 & -50662.60 & -51099.90 & -51149.93 & -51023.85
MolPDF score & 4442.23 & 4423.80 & 4270.78 & 4591.98 & 4487.63

\label{score_b1}

|c|c|c|c|c|c| & model 1 & model 2 & model 3 & model 4 & model 5
DOPE score & -50821.29 & -50820.81 & -51103.58 & -50769.45 & -50637.14
MolPDF score & 4254.95 & 4408.62 & 4328.90 & 4536.73 & 4429.58

\label{score_b3}

\(\beta\)-subunit C-terminal tail orientation and folding

Contacts \(\beta\)-subunit C-terminal with the tubulin core

Discussion