Open access peer-reviewed chapter - ONLINE FIRST

Homology Modeling of Tubulin Isotypes to Investigate MT-Tau Interactions

By Vishwambhar Vishnu Bhandare

Submitted: July 9th 2020Reviewed: January 4th 2021Published: February 17th 2021

DOI: 10.5772/intechopen.95792

Downloaded: 8


The Homology modeling techniques uses the template structure(s) to model the full-length structure of unknown sequence. It is being used for determining the structure of biological macromolecules, especially proteins. The wide applications of homology modeling approach have helped us to address various challenging problems in the field of biological sciences and drug discovery despite the limitations in using analytical techniques like X-ray, NMR and CryoEM techniques. Here, this chapter emphasize on application of homology modeling in determining MT-Tau interactions which are important in the Alzheimer disease. In Alzheimer diseases, tau detaches from MTs in misfolded shape and forms insoluble aggregates in neurons due to post-translational modifications. MT-tau interactions are largely unknown due to differential expression of neuronal specific tubulin isotypes and intrinsically disordered nature of tau. MTs play crucial roles in important cellular functions including cell division, transport of vesicles, cell signaling, cell motility etc. MTs are composed of different tubulin isotypes which differs mainly at C-terminal tail. In humans, nine β-tubulin isotypes have been reported which are expressed differently in different tissues. Structures for different tubulin isotypes are still lacking due to their complex differential expression pattern and purification. Hence, homology modeling approach allowed us to generate homology models for different neuronal specific tubulin isotypes and study their interactions with tau repeats. It is believed that this study would gain more structural and functional insights to the linked Alzheimer diseases.


  • homology modeling
  • microtubule
  • tubulin isotypes
  • Alzheimer disease
  • molecular modeling

1. Introduction

Bioinformatics is an interdisciplinary science that uses both computational and informational approaches to retrieve, analyze, organize, visualize, store and develop biological data [1]. It is widely applied in the field of life sciences, especially in functional genomics, sequence analysis, proteomics, drug discovery, etc. Prediction of the structure and functions of the genes and proteins have become a fundamental task in the life science researches. The present book chapter involves molecular modeling study to investigate intermolecular interactions between Microtubule (MT) and Tau. Though these interactions are important in Alzheimer’s disease, the detailed knowledge on MT-Tau interactions are still lacking mainly because of two reasons (i) Lack of full-length structure of tau due to its intrinsically disordered nature and (ii) Differential expression of tubulin isotypes in different type of cells in particular brain and neuronal cells. Earlier experimental efforts have been made to elucidate these interactions using solution structure (PDB ID: 2MZ7.pdb) however correct binding mode and atomistic interactions at structural level are poorly understood. Therefore, this chapter focuses on application of molecular modeling techniques in understanding important MT-Tau interactions in the Alzheimers disease. Bioinformatics approaches like sequence analysis, homology modeling. MD simulations and binding energy calculations are employed systematically to address this challenging problem in the field of Alzheimer’s disease.

Tau is intrinsically disordered protein encoded by ‘mapt’ gene located on chromosome 17 [2]. The primary function of the tau protein is to bind and stabilize the microtubule. It is abundantly expressed in the brain and neuronal tissues hence its misregulation is associated with the Alzheimers and other neurodegenerative disorders [3, 4]. Till date about six isoforms of tau are reported in the human central nervous system. The length of these six isoforms varies between 352 to 441 residues [5].

Primary structure of tau contains the projectile domain at N-terminal (residue 1–244) which is composed of the acidic and proline-rich region, and the C-terminal repeat domain which consists of 4 repeats i.e., R1, R2, R3 and R4 (residues 245–441) (Figure 1). The six isoforms of tau mainly differs by the existence of either R3 or R4 repeats at the C-terminal domain [6]. The one of the tau isoforms is referred as longest isoform mostly observed in humans which comprises 4 repeats i.e., R1, R2, R3 and R4. while the shortest isoform of tau has only 3 repeats (R1, R2 and R3). This shortest isoform of tau is reported in the fetus brain and less common in adults [5, 7]. Figure 1A represents the structure of tau repeat region R2 which is bound to the MT composed of β/α/β tubulin subunits [8]. Hereafter, tau repeat R2 will be mentioned as ‘TauR2’ for the simplicity. The tau repeats R1, R2, R3 and R4 prefers to bind at the outer surface of microtubule (MT) to stabilize it (Figure 1A) and regulates MT polymerization [6]. Figure 1B represents domain organization in the tau primary structure and Figure 1C shows the sequence of TauR2 which is reported in the CryoEM model. It is well established that tau primarily helps in the assembly and stabilization of axonal MTs, which contributes to the proper functioning of neuronal cells [9]. However, recent studies have reported new functional role of the tau in addition to the axonal, i.e. labile domain of the MT to promote its assembly [10]. Tau detaches from the MTs and forms abnormal, fibrillar structures of insoluble aggregates due to post-translational modifications in Alzheimer diseases and other neurodegenerative diseases associated with tau [11, 12].

Figure 1.

CryoEM Structure of tubulin subunits bound to TauR2. (A) tubulin subunits bound to TauR2 in CryoEM structure 6CVN.pdb. TauR2 domain binds at the outer surface of the MT. (B) Domain organization in tau, (C) sequence of TauR2. [Source: Bhandare et al, 2019; doi: 10.1038/s41598-019-47249-7].

Full-length structure of tau protein is not yet determined using X-ray crystallographic techniques due to its intrinsically disordered nature. Also, the efforts to find its solution structure using NMR spectroscopy have failed [13]. Thus, the MT-Tau interactions have been studied so far using various biochemical and biophysical techniques [14, 15, 16, 17]. The CryoEM have showed marginal success in determining the structure of tau repeat R2 bound to MT however it shows discontinues density of tau repeats along with each protofilament upon MT binding [8]. Hence, they synthetically developed R1 and R2 repeat of tau and their interactions with MT were examined. These two tau repeats adopts the extended conformation along the crest of protofilament which stabilizes the MT structure by binding to the interface of tubulin dimers [8].

MTs are made from αβ-tubulin heterodimer subunits [18]. In human, seven α-tubulin and nine β-tubulin isotypes are reported showing their tissue-specific expressions. For instance, βI tubulin isotype is ubiquitously expressed in all cells, βII and βIII tubulin isotype are mainly expressed in brain and neuronal cells, and βVI tubulin isotype is expressed in erythroid cells and platelets [19]. The βI tubulin isotype reported to play crucial role in cell viability, βII tubulin isotype is important for neurite growth and βIII tubulin isotype protects nerve cell against free radicals and reactive oxygen species [20]. It has been well known that all β-tubulin isotypes share a significant residue conservation except the C-terminal tail region of MT [21, 22, 23, 24] which is flexible in nature and structurally disordered. The C- tail region of all these isotypes overhang outwards of the MTs. The C-tail shows interactions with various MAPs including tau and regulate MT dynamics [25, 26].

It is well documented that the composition of β-tubulin isotypes (i) affects MT dynamic instability [27, 28], (ii) their interaction with motor proteins [29], (iii) their binding to the anti-drugs [21, 22, 30] and (iv) different MAPs including tau [31, 32]. These tubulin isotypes show tissue specific expression as their relative proportion varies greatly in different type of cells [20, 33, 34]. It is also well established that binding of tau to the MT promote or demote microtubule polymerization [35]. However, the differential binding affinity of tau to the various β-tubulin isotypes expressed in different types of cells is completely unknown. Therefore, we studied relative binding affinity of Tau repeat region R2 with neuronal specific β-tubulin isotypes namely βI, βIIb, and βIII using molecular modeling [36].

2. Methodology

2.1 Sequence analysis and homology modeling of tubulin isotypes

In biological research, sequence analysis plays a major role as it has wide range of applications such as, (i) whole genome sequencing and annotation, (ii) identification of functional elements in the sequence, (iii) gene prediction, (iv) comparative genomics, (v) protein classification, (vi) protein and RNA structure prediction, (vii) evolutionary studies, etc. Protein sequences reveal the evolutionary history and hence, the events occurred during evolutions can be traced from the protein sequences.

The structure 6CVN.pdb is used as template structure for homology modeling of neuronal specific human tubulin isotypes namely βI, βIIb and βIII tubulin (uniprot IDs Q9H4B7, Q9BVA1 and Q136509). The multiple sequence alignment of these sequence was performed using ‘clustal omega’ tool [37]. The multiple sequence alignment reveals that residue variations is mainly at the C-terminal tail when compared to the other regions of the protein. The high-resolution cryo-EM structure of β/α/β-tubulin bound with TauR2 was recently deposited in the RCSB structural databases [8] was used as a template to build βI/α/βI-TauR2, βIIb/α/βIIb-TauR2 and βIII/α/βIII-TauR2 complexes. The structure for C-terminal tail was absent in the template structure (6CVN.pdb), therefore the tail region was modeled using the Modeler 9v20. Hereafter, template structure with modeled C-terminal tail region would be referred as 6CVN* in the further discussion.

Template based homology models for neuronal specific tubulin isotypes βI, βIIb, βIII was build using Modeler 9v20 [38]. The least discrete optimized potential energy (DOPE) score model was selected for further use. The stereo-chemical properties of these modeled subunits were evaluated and validated using the GMQE score [39], verify3D [40], ERRAT score [41] and Ramachandran plot through PROCHECK [42]. The selected subunit models were further used to 6CVN*-TauR2, build βI/α/βI-TauR2, βIIb/α/βIIb-TauR2 and βIII/α/βIII-TauR2. These complexes namely 6CVN-TauR2, 6CVN*-TauR2, βI/α/βI-TauR2, βIIb/α/βIIb-TauR2 and βIII/α/βIII-TauR2 were further subjected for energy minimization to get their least energy state. Here we used Steepest Descent and Conjugate Gradient methods in Gromacs 2018.1 software for minimization [43]. The process of energy minimization is a numerical procedure aimed to find a minimum on the potential energy surface (PES) of the newly modeled conformation which mostly exists at a higher energy level. These minimized models were used as a starting input structures for molecular dynamics simulations to understand the binding mode and binding affinity of TauR2 towards neuronal specific tubulin isotypes βI, βII and βIII.

2.2 Molecular dynamics simulations of tubulin-TauR2 complexes

Molecular dynamics (MD) simulation plays a key role in exploring the structure and function of biological macromolecules [44]. In MD simulations, the dynamic behavior of the molecule is studied as a function of time. Molecular dynamics is being routinely used to address various biological problems such as biomolecular interactions (Protein–protein, protein-DNA/RNA), molecular pathways, Drug-receptor interactions, dynamics of protein folding, protein aggregations, protein structure prediction, etc. Tremendous development in high performance computing and simplicity of the basic MD algorithm has shortened the time required to perform molecular dynamics simulation and hence, studying larger systems became an easier task [45].

All atom molecular dynamics (MD) simulation was performed in explicit solvent on the modeled tubulin-TauR2 complexes (i.e. 6CVN-TauR2, 6CVN*-tau, βI/α/βI-TauR2, βIIb/α/βIIb-TauR2 and βIII/α/βIII-TauR2) using GROMACS 2018.1 [43, 46]. The ‘Amber99SB-ILDN’ force field [47] was chosen for the MD simulation because it is well customized to handle the parameters for GTP,GDP and MG atoms which are the functional players of all the modeled tubulin-TauR2 complexes. The force field parameters for the GDP and GTP molecules were retrieved from the amber parameter database [48, 49]. The ‘xleap’ module of AmberTools was used to generate topology files and initial starting coordinates for all the complexes [50]. All the modeled tau-tubulin complexes were placed at the centre of a cubic shaped solvation box having dimension of 15 Å from the extent of the molecule and TIP3P water model was used for solvation. All the systems were neutralized by adding appropriate number of required counterions. The topology files generated using xleap module of AmberTools were converted to Gromacs compatible topologies with ParmEd tool [51]. Energy minimization was carried out in two steps, In the first step, steepest descent algorithm was used for 50,000 followed by the conjugate gradient [46]. The energy minimized models were equilibrated using canonical ensemble (NVT) followed by isothermal-isobaric ensemble (NPT) for 1 ns. In the NVT equilibration systems were heated to 300 K using V-rescale, a modified Berendsen thermostat [46]. These heated systems were further equilibrated using the Parrinello-Rahman barostat to maintain constant pressure of 1 bar. The production MD simulations were performed for 100 ns without restraining any atoms over all the tubulin-TauR2 complexes using parameters discussed in earlier study [52]. The PME method was used to treat long range electrostatic interactions [53, 54] and covalent bonds involving H-atoms were constrained by using ‘LINCS’ algorithm [55]. The 2 fs time step was set for integrating the newtonian equation during the MD simulation. Similar protocol was adopted to perform MD simulation on three additional systems (i) 6CVN* (without tau), (ii) free tau and (iii) 6CVN*-polyA (as negative control) having 27 amino acids residues. All the MD simulation trajectories were further analyzed by using the GROMACS 2018.1 inbuilt tools [43, 46]. The general parameters explaining the conformational stability such as Root Mean Square Deviation (RMSD), Root Mean Square Fluctuations (RMSF) and Radius of Gyration (Rg) were measured, and the equations used for calculation of these parameters are tabulated in the Table 1.

S. no.ParameterEquationComponent
1RMSDRMSD=i=1Ndi2N𝑑𝑖 is the distance, ‘N’ is number of atoms
2RMSFRMSF=j=1txitjx¯i2t𝑥𝑖 is atom position at time 𝑡, reference position is x¯i
3RgRg=i=1Nxixcom2N𝑥𝑐𝑜𝑚 centre of mass, 𝑥𝑖 is the distance at time 𝑡 from their centre of mass

Table 1.

Equations used to calculate RMSD, RMSF and Rg.

The primary sequence of a protein is a linear chain of amino acids linked by peptide bonds. There is a direct link between the protein sequence, structure and function. The secondary structure of a protein is comprised of coils, α-helices, β-sheets, β-bridge, bend, turn, coils, π-helix and 310-helices. DSSP is an algorithm developed by Kabsch and Sander to extract the secondary structural features based on atomic coordinates [56]. The overall stability of the structure is highly determined by the stable dynamics of these secondary structures and any significant changes in secondary structure attributes to the structural flexibility/fold as well as functional diversity of the protein. Hence, conformational changes in the secondary structure during MD simulation were analyzed using the DSSP programme [56]. The simulation movies over the entire trajectories were generated using the VMD software [57] and publication quality images were generated using the Biovia Discovery studio visualizer [58] and Chimera software [59].

2.3 Calculations of contact surface area (CSA) for tubulin-TauR2 complexes

Solvent Accessible Surface Area (SASA) is used to represent the degree of hydration of a biomolecule. SASA also be especially useful to quantify the stability of the biomolecular complexes in the aqueous medium. The C-terminal tail of the tubulin subunits is highly dynamic in nature and has no definite secondary structure, hence it affects the overall hydrophobic SASA. Therefore, interface of the MT (in this case tubulin trimer made up of β/α/β subunits) where TauR2 binds at the exterior surface has been selected for the calculating the precise CSA. The in-built gromacs tool “gmx sasa” [60] was used to calculate the SASA. In addition, SASA is also calculated for the tubulin subunits and the TauR2.

2.4 Binding affinity of tauR2 towards different neuronal specific tubulin isotypes

The biomolecular recognition pattern mainly depends on the binding ability of the interacting biomolecules. The binding affinity as well as the energy between the two interacting molecules can be calculated using various theoretical approaches like (i) Pathway methods such as Thermodynamic integration (TI) as well as Free energy perturbation (FEP) and (ii) End point methods such as Molecular Mechanics Poission-Boltzman Surface Area (MM/PBSA) and Molecular Mechanics Generalized Born Surface Area (MM/GBSA) [61]. In the present study, MM/PBSA approach was used to calculate relative binding energies of the simulated molecules. This MMPBSA approach is very popular, computationally less expensive, and has better accuracy even for the larger systems [62].

Here, the binding affinity between different neuronal specific tubulin isotypes and TauR2 was estimated by performing relative binding energy calculation similar to earlier studies [63, 64, 65]. The stable trajectory observed in between 70 ns to 100 ns was chosen to perform the binding energy calculations for all the tubulin-TauR2 complexes. The ‘g_mmpbsa’ tool v1.6 was used to perform binding energy calculation using MM/PBSA approach [66]. The parameters for the binding energy calculations were chosen from the earlier similar studies [52, 65, 67, 68, 69]. In the MMPBSA methods binding energy (ΔGbind) of tubulin and TauR2 was calculated by using the following Eq. (1),


Where, the ΔGtubulinTauR2, ΔGtubulinand ΔGTauR2represents the average free energies of the complex (tubulin-TauR2), receptor (tubulin) and ligand (TauR2), respectively. The calculation of the entropic contribution in binding energy is computationally expensive for a larger biomolecular complexes and hence it is omitted as similar to previous studies [21, 22, 70, 71, 72].

3. Results and discussion

In this chapter we employ sequence analysis, homology modeling, MD simulations, and binding energy calculation to (i) gain structural insights to the detailed binding mode, (ii) study atomic level tubulin isotypes-tauR2 interactions and (iii) study relative binding affinity between neuronal specific tubulin isotypes and TauR2.

3.1 Sequence analysis and homology modeling of neuronal specific tubulin isotypes

The residue composition of different β-tubulin isotypes mostly varies at the carboxy-terminal tail region as revealed by the multiple sequence alignment (Figure 2). The βI and βIII tubulin isotypes have longer C-terminal tail regions when compared with the βIIb tubulin isotype. The β-tubulin sequence in the template structure i.e., 6CVN (chain A) and human βIIb tubulin isotypes show 98.65% sequence identity. These sequence variations in the tubulin isotypes are reported to regulate number of protofilaments in the MT and their stability [73]. These β-tubulin isotypes sequences were used to generate three-dimensional homology models using 6CVN as the template. The structures of βI, βIIb, βIII tubulin isotypes were modeled using Modeler 9v20 [38]. The best homology model generated is selected using DOPE score. The DOPE score value for A and C chain of (i) βI subunits are −54299.89 and − 54291.24 (ii) βIIb subunits are −53487.13 and − 53054.42 and (iii) βIII subunits are −53725.86 and − 53054.42. The quality of these models is accessed using GMQE score [74], Verify3D [40], Errat score [41], Z-score [75] and Ramachandran plot [76, 77]. The parameters describing the overall quality of the modeled neuronal specific β-tubulin subunits are shorn in Table 2. The GMQE score provides an estimate of the accuracy of the modeled tertiary structure of neuronal specific β-tubulins. Here, GMQE score for all the modeled β-subunits is 0.98, which represents the accuracy of the generated model. Further, verify3D and ERRAT score also validates the quality of the generated models (Table 1). Ramachandran plots for all the modeled β-tubulin isotypes represents more than 98% of the residues occupy a favored region. The occupancy of amino acid residues in the Ramachandran plot is given in Table 3. These modeled structures of neuronal specific β-tubulin isotypes were used further to build the tubulin and TauR2 complexes such as βI/α/βI-TauR2, βIIb/α/βIIb-TauR2 and βIII/α/βIII-TauR2 using 6CVN.pdb as a template structure. These modeled complexes were used as starting structures to perform MD simulations.

Figure 2.

Multiple sequence analysis of different β-tubulin isotypes. The βI, βIIb, βIII tubulin isotypes and template 6CVN show maximum residue variations mainly at C-terminal tail region. The TauR2 binding regions H12 helix and C-terminal tail region of β-tubulin subunits are shown in hot pink and brown, respectively.

S. No.ChainsGMQEVerify3DErratz-score
1β1 (A) -subunit0.9893.13%81.3212−8.93
2β1 (C) -subunit0.9892.02%83.2569−8.78
3β2b (A) -subunit0.9898.43%87.471−8.65
4β2b (C) -subunit0.9898.43%83.2947−8.37
5β3 (A) -subunit0.9898.44%86.3636−8.54
6β3 (C) -subunit0.9892.67%83.33−8.39

Table 2.

Validation of three-dimensional models generated for βI, βIIb and βIII isotypes chain A and chain C using Swiss model GMQE score, Verify-3D, Errat score.

Regionβ1 tubulinβ2 tubulinβ3 tubulin
Chain AChain CChain AChain CChain AChain C
% of most favored regions98.4 (442)98.9 (444)98.9 (443)98.9 (438)98.9 (443)98.9 (443)
% of additional allowed regions1.3 (6)0.9 (4)0.7 (3)1.1 (5)1.1 (5)0.7 (3)
% of outlier0.2 (1)0.2 (1)0.4 (2)0 (0)0 (0)0.4 (2)

Table 3.

Ramachandran plot showing the percentage of residues in the different regions for tubulin isotypes obtained from the Ramachandran plot using PROCHECK.

3.2 Structural stability of the tubulin-TauR2 complexes

The all atom MD simulations were performed on tubulin-TauR2 complexes namely 6CVN-TauR2, 6CVN*-TauR2, βI/α/βI-TauR2, βIIb/α/βIIb-TauR2, βIII/α/βIII- TauR2 using Gromacs 2018.1 [78]. The stability of these tubulin-TauR2 complexes is accessed by plotting the potential energy during the simulation period, which highlight that all the complexes are well minimized, and simulation trajectories are well converged during the simulation period of 100 ns (Figure 3).

Figure 3.

Energy Minimization Plot Potential energy over the simulation time plotted for 6CVN-TauR2 (black), 6CVN*-TauR2 (orange), βI/α/βI-TauR2 (green), βIIb/α/βIIb-TauR2 (cyan), βIII/α/βIII-TauR2 (violet) are shown.

The parameters describing the stability of tau-tubulin complex such as RMSD (root mean square deviation), RMSF (root mean square fluctuation), and Rg (radius of gyration) was studied. The RMSD values for simulated tubulin-TauR2 complexes, tauR2 and backbone atoms of tubulin trimer without considering disordered C-tail were plotted over the trajectory. This analysis reveals the stability of all the studied complexes throughout the simulation time i. e. 100 ns. Figure 4A and B shows the RMSD plot for studied tubulin-TauR2 complexes and TauR2, respectively. The RMSD for the complex βIII/α/βIII-TauR2 is observed to be relatively more stable than other tubulin-TauR2 complexes. Similarly, structure of TauR2 bound to βIII/α/βIII tubulin trimer expresses stable dynamics during the simulation. The complex 6CVN-TauR2 is stabilized at higher RMSD values, the primary reason for this elevated RMSD value is absence of C-tail region which highlights the importance of C-terminal tail in the stabilizing tubulin-TauR2 complex. Average backbone RMSD value is converged at ~3.5 Å hence represents the equilibration of all above simulated systems (Figure 5). The molecular dynamics simulation movies reveals the stable dynamics of all the simulated systems 6CVN-TauR2, 6CVN*-TauR2, βI/α/βI-TauR2, βIIb/α/βIIb-TauR2 and βIII/α/βIII-TauR2 (,,,, respectively. Further, specificity of TauR2 towards tubulin subunits was accessed by replacing the TauR2 with negative control ‘polyA’ peptide of same length. Interestingly, This system having negative control poly A bound to 6CVN* shows the weak binding during the simulation. These weaker interactions of polyA peptide with tubulin subunits ( represents that tauR2 has specificity towards tubulin subunits. s.

Figure 4.

Stability of the tubulin-TauR2 complex and TauR2. (A) The Root mean square deviation values (RMSD) for tubulin-tauR2 complexes. RMSD values for 6CVN, 6CVN*, βI/α/βI, βIIb/α/βIIb and βIII/α/βIII have been plotted in black, orange, green, cyan and violet, respectively. (B) The Root mean square deviation values for TauR2 shown using same color Scheme as in (A).

Figure 5.

Backbone Root mean square deviation for different tubulin subunits. Color scheme is same as Figure 3.

3.3 Residue fluctuations of tubulin subunits and TauR2

The flexibility of tubulin trimers systems and TauR2 has been studied using RMSF analysis. For this analysis Cα-atom from the backbone was selected to get fluctuations in the overall protein. Figure 6 represents the RMSF for tubulin subunits and TauR2. The residues from the H12 helix of β-tubulin and the C-terminal tail region (residue 400–451) show significant decrease in the RMSF values, as their free dynamics is arrested by tau binding (Figure 6A and B). RMSF values for the tubulin β-subunits in the systems 6CVN*, βIIb/α/βIIb and βIII/α/βIII are lesser than those of 6CVN and βI/α/βI tubulin subunits (Figure 6B). This observation also highlights the binding of TauR2 at the interdimer interface where residual fluctuations are less. However, the part of C-tail region which has no direct contact with TauR2 is highly flexible (Figure 6B). The H12 helix and C-terminal tail region of the tubulin subunits significantly contribute to the non-covalent interactions resulting towards stronger binding of TauR2. Therefore, these intermolecular interactions were analyzed in detail and are discussed in the section ‘Intermolecular interactions between tubulin and tau’. Further, atomic Cα-fluctuations of TauR2 (Figure 6C) was also studied for better understanding its conformational behavior during the MD simulations. It is surprising to observe highest fluctuations at the N- and C-terminal region in TauR2 bound to 6CVN, where the C-terminal tail region is absent (Figure 6C). Interestingly, residual fluctuations expressed by TauR2 bound to βIII/α/βIII-tau complex are much lesser as compared to 6CVN*-TauR2, βI/α/βI-TauR2 and βIIb/α/βIIb-TauR2 complexes (Figure 6C). This also proves that the C-terminal tail region of tubulin subunits plays an important role in the binding of TauR2.

Figure 6.

Root mean square fluctuations (RMSF) of different β/α/β tubulin subunits and TauR2 (A) RMSF of different β/α/β tubulin subunits (B) Magnified view of their C-terminal H12 helix and tail regions (C) RMSF of TauR2 bound with different β/α/β tubulin subunits observed during the simulations5. Color scheme is same as Figure 3.

Overall, RMSF analysis suggests the significance of H12-helix and C-terminal tail region in stabilization of the microtubule by binding of tau repeats (TauR2) and it also reveals the greater affinity of TauR2 towards βIII tubulin isotypes which are overexpressed in neuronal cells and brain. Further compactness of all the tubulin-TauR2 complexes was explored by calculating the radius of gyration (Rg) and this analysis is discussed in the next section.

3.4 Compactness of tubulin-TauR2 complexes

The radius of gyration (Rg) indicates the level of compactness of the protein system which is helpful in getting an insight into the stability of the protein–protein complex. It also helps to understand folding or unfolding of protein structure during the simulation. The Rg values for all the studied tubulin-TauR2 complex ranges from 38.8–40.5 Å (Figure 7A). The complex βIII/α/βIII-TauR2 shows stable Rg value for the entire simulation period however other complexes 6CVN-TauR2 6CVN*-TauR2, βI/α/βI-TauR2, βIIb/α/βIIb-TauR2 show variations in their Rg values. The absence of C-terminal tail region in the complex 6CVN-TauR2 leads to the less Rg values when compared to other tubulin-TauR2 (Figure 7A). Figure 7B represents the Rg values of only TauR2 in different tubulin-TauR2 complexes. The Rg values for TauR2 shows fluctuations between 17.5 to 20 Å in case of 6CVN*, βIIb/α/βIIb, and βIII/α/βIII complexes except for βI/α/βI complex (Figure 7B). The βIII tubulin subunits show Rg value of ~18 Å and βI tubulin subunits have largest Rg value of 22.5 Å as shown in Figure 7B. On the other hand, TauR2 bound to 6CVN shows uninterrupted decline in Rg values from 21.5 Å to 16.5 Å. This analysis also highlights the importance of C-terminal tail region in the stable binding of tau (Figure 7B). It is important to note that βIII tubulin subunits (Figure 8) have Rg values like that of βIII/α/βIII-TauR2 complex (Figure 7A). This highlights that the tubulin subunits composed of βIII tubulin isotype are structurally stable after binding to the TauR2. Thus, calculation of Rg values for tubulin-TauR2 complexes, tubulin subunits and TauR2 reveals (i) structural stability of the βIII/α/βIII-tau complex over other complexes and (ii) importance of the C-terminal tail region in the binding of TauR2. Contact surface area (CSA) and solvent accessible surface area (SASA) was calculated using ‘gmx sasa’ tool of gromacs to understand the exposure of the interface residues of tubulin subunits bound to the TauR2 [46].

Figure 7.

Radius of Gyration (Rg) of different tubulin-TauR2 complexes and TauR2. (A) Rg of 6CVN-TauR2 (black), 6CVN*-TauR2 (orange), βI/α/βI-TauR2 (green), βIIb/α/βIIb-TauR2 (cyan), βIII/α/βIII-TauR2 (violet) (B) Rg for TauR2 in different tubulin-TauR2 complexes. Color scheme same as Figure 3.

Figure 8.

Radius of Gyration for different tubulin isotypes. Color scheme is same as Figure 3.

3.5 Contact surface area (CSA) and solvent accessible surface area for tubulin-TauR2 complexes

The CSA and SASA describes the accessibility of a binding interface and protein surface to the solvent, respectively. It is well documented that, TauR2 binds to the MT exterior surface via C-terminal tail region [8, 79, 80, 81, 82]. Therefore, initially contact surface area (CSA) of the interface between the TauR2 and tubulin trimer, was calculated, without considering flexible C-terminal tail region. The CSA of βIII/α/βIII is very less when compared to other tubulin isotypes (Figure 9A) this represents the tight binding of TauR2 to the βIII/α/βIII tubulin subunits. The higher CSA for βI/α/βI-TauR2 complex indicates weaker binding of the TauR2 to the βI/α/βI tubulin subunits. Furthermore, least SASA in complex βIII/α/βIII-TauR2 represents tight binding of TauR2 to the βIII/α/βIII (Figure 9B). On the other hand, higher hydrophobic SASA of the complex βI/α/βI-TauR2 indicate the exposure of hydrophobic residues which are responsible for loss of native contacts between tubulin and TauR2. The SASA for 6CVN*, βI/α/βI, βIIb/α/βIIb, βIII/α/βIII shows higher SASA values between 4900 and 5400 Å when compared to 6CVN-TauR2 (~4500 Å) due to the presence of C-terminal tail region (Figure 10). To get detailed understanding of the atomic-level interaction between tubulin isotypes and TauR2, further hydrogen bonding interactions were estimated during simulation and in the MD simulated end-structures obtained from trajectory.

Figure 9.

Contact surface area (CSA) and solvent accessible surface area (SASA) of different β/α/β-tubulin subunits andTauR2. (A) CSA for different 6CVN-TauR2 (black), 6CVN*-TauR2 (orange), βI/α/βI-TauR2 (green), βIIb/α/βIIb-TauR2 (cyan), βIII/α/βIII-TauR2 (violet) complexes. (B) hydrophobic SASA for tubulin isotype bound TauR2. Color scheme same as Figure 3.

Figure 10.

Solvent accessible surface area for different tubulin subunits. SASA plotted for 6CVN (black), 6CVN* (orange), βI/α/βI (green), βIIb/α/βIIb (cyan), βIII/α/βIII (violet) are shown.

3.6 Intermolecular interactions between tubulin and TauR2 in tubulin-TauR2 complexes

The total number of hydrogen bonds formed between tubulin isotypes and TauR2 during the MD simulations are calculated using in-built ‘gmx hbond’ command [46]. The cut-off value for the measurement of H-bond was set to 3.4 Å. Consistent H-bond formation was observed throughout the MD simulation in all tubulin-TauR2 complexes. The average number of H-bonds roughly varies between 10 to 20 as shown in Figure 11. The details of atom participating in the hydrogen bonding interactions present between tubulin isotypes and TauR2 in the MD simulation end-structures are listed in Table 4. All the hydrophobic interactions participating in the formation of stable tubulin-TauR2 complexes are listed in Table 5. The βIII/α/βIII-TauR2 complex shows the maximum number of electrostatic interactions when compared to other tubulin-TauR2 complexes (Table 6). Further, to understand the role of TauR2 in stabilizing tubulin subunits, secondary structure analysis on TauR2 was done using DSSP.

Figure 11.

The number of hydrogen bonds formed in between tubulin subunits and TauR2 during MD simulation. Color scheme is same as Figure 3.

SystemAtoms involved in H-bondingDistance (Å)Angle (°)
6CVN-TauR2D: SER16: HG - B: GLU434:OE21.55968170.912
C: LYS392:HZ2 - D: ASP22:OD11.79981155.003
D: SER20:H - B: GLU434:OE11.80025149.621
D: CYS18:H - B: ASP431:OD11.8071147.422
D: GLY19:H - B: ASP431:OD11.81516170.051
D: ASN23:HD21 - C: PHE389:O1.83614168.422
D: ILE5:H - B: GLU415:OE11.87843148.138
B: ARG402:HH22 - D: LYS7:O1.95256131.998
D: LYS21:HZ3 - B:ASP438:O1.96835128.744
C: ARG391: HE - D: SER20:O1.97398162.451
D: ASN6:H - B: GLU415:OE21.97803165.667
D: LYS17:HZ2 - B: ASP424:OD12.01656167.691
6CVN*-TauR2A: SER16: HG - E: ASP431:OD21.56368161.89
A: ASN6:HD21 - G: GLN433:OE11.68122165.335
A: LYS1:HZ1 - G: ASP417:OD11.73491152.719
A: LYS7:HN - E: ALA400:O1.75148168.412
A: SER12: HG - A: ASP10:OD21.77457163.204
E: LYS401:HZ1 - A: ASN6:OD11.77792156.289
A: SER20:HN - E: GLU434:O1.8254168.05
A: ILE4:HN - G: GLN424:OE11.85744159.004
F: ARG391:HH21 - A: SER20: OG1.88141150.648
A: LYS17:HN - E: ASP431:OD21.88237156.579
A: LYS1:HT2 - G: ASP417:OD21.93321163.221
A: LYS25:HZ2 - A: VAL27: OXT1.95566155.991
F: ARG391: HE - A: SER20:O1.95954137.171
A: SER12:HN - A: ASP10:OD21.97395166.966
A: VAL2:HN - G: GLU421:OE21.97952167.491
A: SER20: HG - E: TYR262: OH2.01579158.876
A: CYS18:HN - E: ASP431:OD22.04985143.832
A:ASP22:HN - A: ASP22:OD12.06458123.428
βI/α/βI-TauR2D: SER20: HG - B: GLU434:OE21.71015154.981
D: LYS17:H - B: ASP431:OD21.73879154.124
D: SER16: HG - B: ASP431:OD21.75545159.778
D: VAL2:H - A: GLU421:OE11.79605175.632
D: SER20:H - B: GLU434:OE21.79906165.806
D: LYS21:HZ1 - B: GLU434:O1.82861143.493
B: LYS430:HZ2 - D: VAL14:O1.84508147.656
D: ASN23:HD22 - C: PHE389:O1.87483147.025
D: LYS1:H3 - A: ASP417:OD11.96013159.665
D: LYS7:H - B: ALA400:O2.03926154.297
D: ILE4:H - A: GLN424:OE12.05966141.181
βIIb/α/βIIb-TauR2A: SER16: HG - E: ASP431:OD11.71023174.429
A: LYS17:HZ2 - E: ASP424:OD21.72125162.879
A: LYS1:HT2 - F: GLU421:OE21.73767160.093
A: LYS25:HZ2 - E: GLU445:OE11.7524156.386
A: LYS7:HN - E: ALA400:O1.7726158.577
A: ASN6:HD22 - F: ASP431:OD11.88399166.501
A: LYS17:HN - E: ASP431:OD11.90489145.148
E: ARG402:HH12 - A: LYS7:O1.90591147.259
A: VAL27:HN - E: GLU445:OE11.92907160.86
A: CYS18:HN - E: ASP431:OD12.03765145.506
A: LYS25:HZ3 - E: GLU446:O2.04646165.654
F: GLN424:HE21 - A: LYS1:O2.06218173.384
βIII/α/βIII-TauR2A: SER16: HG - E: ASP431:OD11.57372172.137
A: LYS21:HZ1 - E: GLU443:OE21.74362173.969
F: GLN424:HE21 - A: VAL2:O1.79549177.625
A: GLN15:HE21 - E: GLU443:OE21.81621154.684
A: LYS8:HZ3 - F: GLU433:OE11.84962164.797
A: LYS8:HZ1 - E: ASP396:OD11.85984168.252
G: ARG391:HH11 - A: ASN23:OD11.89423162.494
E: GLY442:HN - A: ILE24:O1.93324171.538
A: LYS17:HN - E: ASP431:OD11.99492142.974
A: CYS18: HG - E: ASP431:OD12.0734155.708

Table 4.

Hydrogen bonding interaction between tubulin subunits and TauR2 after molecular dynamics simulations.

SystemHydrophobic InteractionsDistance (Å)
6CVN-TauR2B: ALA427 - D: LYS174.23972
B: ARG264 - D: CYS184.32305
B: ALA426 - D: VAL144.34456
B: ARG402 - D: ILE45.12692
B: VAL409 - D: ILE45.16258
B: ARG422 - D: LEU115.2126
B: ALA427 - D: VAL145.2724
6CVN*-TauR2E: ALA426 - A: VAL143.78619
E: ALA426 - A: LEU114.1983
E: ALA400 - A: LYS84.29681
G: PHE260 - A: ILE44.3943
A: VAL2 - G:PRO2614.85577
A: LEU9 - A: LEU115.10635
A: LYS21 - E: VAL4375.11425
E: ALA427 - A: VAL145.28218
A: VAL14 - A: LEU115.3507
E: ARG422 - A: LEU95.37331
A: VAL2 - A: ILE45.42169
βI/α/βI-TauR2B: ALA426 - D: LEU114.086
B: ALA426 - D: VAL144.25902
A: PHE425 - D: ILE44.29071
B: ALA427 - D: VAL144.57533
A: PHE260 - D: VAL24.62894
B: TYR262 - D: CYS184.75533
B:PRO263 - D: LYS174.89022
B: ARG402 - D: LYS74.99921
D: CYS18 - B: VAL4355.09522
B: ARG422 - D: LEU95.15149
B: LYS430 - D: VAL145.38523
A: ALA428 - D: ILE45.40848
B: VAL440 - D: LYS215.41789
C: ILE405 - D: ILE245.46457
βIIb/α/βIIb-TauR2E: ALA427 - A: LYS173.91294
E: ALA426 - A: VAL143.93471
E: ALA426 - A: LEU114.19609
A: CYS18 - E: ARG2644.36259
E: ALA400 - A: LYS84.4046
A: CYS18 - E:PRO2635.04073
A: CYS18 - E: ILE2655.21553
G: LYS392 - A: ILE245.27213
E: TYR399 - A: LEU95.32769
βIII/α/βIII-TauR2E: ALA426 - A: LEU113.83867
E: ALA426 - A: VAL144.14469
E: ALA427 - A: LYS174.22472
E: ALA427 - A: VAL144.83218
E: ARG422 - A: LEU114.84097
A: LYS25 - E: VAL4375.06163
A: CYS18 - E: ARG2645.14951
E: ARG422 - A: LEU95.48905

Table 5.

Hydrophobic interactions between different β/α/β-tubulin isotypes and TauR2 after molecular dynamics simulations.

SystemsElectrostatic interactionsDistance (Å)
6CVN-TauR2D: LYS8:NZ - A: ALA430:O4.04432
D: LYS21:NZ - B: SER439:O4.66847
D: LYS25:NZ - B: GLU434:OE24.90999
D: LYS21:NZ - B: GLU434:OE15.38615
6CVN*-TauR2A: LYS1: N - G: GLU421:OE24.85022
A: LYS7:NZ - E: GLU415:OE14.9846
A: LYS25:NZ - E: GLU441:OE15.12197
A: LYS17:NZ - E: ASP424:OD25.27557
βI/α/βI-TauR2D: LYS1: N - A: GLU421:OE12.86182
D: LYS25:NZ - C: GLU412:OE14.31321
D: LYS7:NZ - B: GLU415:OE14.45715
D: LYS21:NZ - B: GLU434:OE24.75044
βIIb/α/βIIb-TauR2A: LYS21:NZ - E: GLU434:OE24.3529
A: LYS8:NZ - E: ASP396:OD25.021
βIII/α/βIII-TauR2A: LYS25:NZ - E: GLU434:OE22.68494
A: LYS25:NZ - E: GLU450:OE22.85019
A: LYS1: N - F: ASP417:OD22.91844
A: LYS1: N - F: GLU421:OE24.28381
A: LYS7:NZ - E: GLU415:OE14.31182
A: LYS21:NZ - E: GLU434:OE14.48844
A: LYS17:NZ - E: ASP424:OD24.70401
G: LYS392:NZ - A: ASP22:OD24.95323
A: LYS21:NZ - E: GLU450:OE25.37409

Table 6.

Electrostatic interactions between different β/α/β-tubulin isotypes and TauR2 after molecular dynamics simulations.

3.7 Conformational changes in TauR2 upon tubulin binding

Tau belongs to the class of intrinsically disordered proteins for which no definitive secondary structure exists. Hence their structure determination is difficult by using existing biophysical techniques like X-ray crystallography and NMR. Previous experimental observations propose that tau repeat undergoes a conformational changes from the disordered to ordered state when it binds to the MT [2, 83, 84, 85, 86]. Hence, the secondary structural changes during the MD simulations in the TauR2 were studied using DSSP [56]. Figure 12 represents conformational changes in the secondary structure of TauR2 upon binding to the tubulin. TauR2 in 6CVN- TauR2 (Figure 12A) and 6CVN*- TauR2 complexes (Figure 12B) show formation of short and transient 310-helix during the simulation. The TauR2 in βI/α/βI- TauR2 complex does not form either α-helix or transient 310-helix as shown in Figure 12C. The TauR2 in βIIb/α/βIIb-TauR2 complex shows the formation of short-lived α-helix and 310-helix (Figure 12D). The terminal residues of TauR2 from Ser293 to Val300 in βIII/α/βIII-TauR2 complex () undergoes turn to α-helix conformational transition (Figure 12E). Therefore, it is proposed that this conformational transition of TauR2 from disordered to ordered state promotes the stable binding of TauR2 with the βIII/α/βIII tubulin subunits.

Figure 12.

The secondary structure changes during MD simulation using DSSP for TauR2. Secondary structure changes observed in (A) 6CVN-TauR2 (B) 6CVN*-TauR2 (C) βI/α/βI-TauR2, (D) βIIb/α/βIIb-TauR2 and (E) βIII/α/βIII-TauR2.

3.8 Relative binding affinity of TauR2 towards neuronal specific tubulin isotypes

The relative binding affinity of TauR2 towards neuronal specific tubulin isotypes (β/α/β) was analyzed by performing MMPBSA calculations for complexes 6CVN-tau, 6CVN*-tau, βI/α/βI-tau, βIIb/α/βIIb-tau and βIII/α/βIII-tau etc. The energy components that govern the binding energy are recorded in Table 7. This analysis reveals that, βIII/α/βIII-tau complex shows most favorable interactions while 6CVN-tau complex is least favorable as supported by the binding energy values listed in Table 7. Thus, it is interesting to note the significance of C-terminal tail of the tubulin subunits in the stable binding of the tau repeat R2 to stabilize this complex. The order of calculated binding energy in between TauR2 with neuronal specific tubulin-TauR2 complexes is βIII/α/βIII > βIIb/α/βIIb >6CVN* > βI/α/βI > 6CVN. The electrostatic interactions in these complexes contribute significantly to the binding energy particularly in βIII/α/βIII-TauR2 and βIIb/α/βIIb-TauR2 complexes when compared to the 6CVN and βI/α/βI tubulin subunits (Table 7). The complex βI/α/βI-TauR2 exhibits higher binding energy leading to its weaker affinity towards TauR2. In addition to βIII/α/βIII-TauR2 the complex βIIb/α/βIIb-TauR2 also exhibits relatively higher affinity towards TauR2 compared to rest other complexes. Further, contribution of the individual residues in the binding energy has been investigated by calculating the decomposition energy for each residue. This analysis reveals that, residues from the H12 helix and C-terminal tail of tubulin subunits shows maximum contribution in the binding energy (Figure 13). The per-residue interactions energy (residue decomposition energy) calculated for various pairs of interacting residues highlights the importance of the interacting residues. The residues from the H12-helices and C-terminal tail region of complex βIII/α/βIII shows maximum contribution (most negative energy) in the non-bonded contacts leading to the stable and tight binding of the tauR2 to the βIII/α/βIII tubulin subunits.

Vdw−97.25 ± 0.55−131.17 ± 0.62−133.38 ± 0.69−124.36 ± 0.60−125.24 ± 0.59
Elec−1232.33 ± 3.43−1534.36 ± 3.12−1423.07 ± 2.77−1659.30 ± 4.68−1768.65 ± 2.77
Polar413.98 ± 4.73349.43 ± 4.18439.27 ± 2.95433.79 ± 4.13505.14 ± 2.92
SASA−12.35 ± 0.06−15.12 ± 0.07−17.05 ± 0.05−15.66 ± 0.06−16.04 ± 0.05
Binding Energy−927.87 ± 3.15−1331.15 ± 3.19−1134.13 ± 1.13−1365.2.26 ± 2.26−1404.7 ± 1.84

Table 7.

The relative binding energy of the tubulin-TauR2 complexes calculated using MMPBSA. All energies are given in kcal/Mol.

Figure 13.

The H12 and C-terminal tail regions show highest energy contribution for the binding of TauR2 in 6CVN*, βI/α/βI, βIIb/α/βIIb and βIII/α/βIII tubulin subunits except in case of 6CVN which does not have C-terminal tail region.

Hence, relative binding energy calculations further support all other MD simulation results highlighting the tight binding of TauR2 to the βIII/α/βIII tubulin isotype which is predominantly expressed in the neuronal cells and brain.

4. Conclusion

MTs are distributed across all types of cells and play an important role in the cellular functions. Structurally MTs are made up of α/β heterodimeric subunits. Large diversity of α and β-tubulin isotypes exists which are differently expressed in different types of cells, this makes MTs unique from one another in relative proportion of isotypes. The much elevated expression levels of βII and βIII tubulin isotypes about 58% and 25% respectively have been reported to in neuronal cells and brain [35]. The present study extensively uses molecular modeling approaches including homology modeling, MD simulation, binding energy to investigate the binding mode and interaction of neuronal specific tubulin isotypes with TauR2.

Extensive analysis on MD simulation trajectory shows a stable complex formation in between different tubulin isotype and TauR2. The stability of these complexes is mainly mediated by the interactions of H12 helix and C-terminal tail of the α/β tubulin isotypes with TauR2. TauR2 shows differential binding affinity towards various neuronal specific β-tubulin isotypes (βI, βII and βIII) the order of binding affinity is ‘βIII> βIIb>βI’. Thus, it is found that TauR2 expresses greater binding affinity with βIII and βIIb tubulin isotypes which are abundantly expressed in neuronal cells and brain. The molecular modeling strategy adopted in this chapter could be potentially used to understand differential binding affinity of other tau repeats such as R1, R3, R4 towards β tubulin isotypes present in other cell lines. The structures for other repeats could be generated using homology modeling and their interactions with neuronal specific tubulin isotypes could also be studied using similar molecular modeling approach. I believe that the knowledge on precise molecular origin of differential binding affinity of tau with β tubulin isotypes present different cell types will pave the way for developing effective treatments against tau related disorders such as Alzheimers, Amyotrophic lateral sclerosis and other tauopathy linked neurodegenerative disorders.

Thus, homology modeling allows us to investigate important biomolecular interactions whose protein structures are not known and/or difficult to get by using biophysical experiments. This, homology modeling, computational tool has made it easier to address various challenging problems in understanding the basic phenomenon’s/pathways in modern biology. Also, it has proven wide successful applications in determining the role of proteins in various genetic diseases, hormonal disorder cancers, neurological disorders and other diseases etc. by exploring protein structure and functions using molecular modeling techniques.


VVB is thankful to IIT Bombay for Institute postdoctoral fellowship. Author is also thankful to Prof. Ambarish Kunwar, Department of Biosciences and Bioengineering, IIT Bombay, Mumbai for fruitful discussion and providing necessary computational resources to perform this research work. Author also sincerely thanks Creative Commons license for giving permission to use data from my manuscript my published research article in Scientific Reports ( A copy of creative commons license can be found at the link

Conflict of interest

The author declares no conflict of interest.

Download for free

chapter PDF

© 2021 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Vishwambhar Vishnu Bhandare (February 17th 2021). Homology Modeling of Tubulin Isotypes to Investigate MT-Tau Interactions [Online First], IntechOpen, DOI: 10.5772/intechopen.95792. Available from:

chapter statistics

8total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us