Homology Modeling of Tubulin Isotypes to Investigate MT-Tau Interactions

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.


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].
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].

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.

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, Drugreceptor interactions, dynamics of protein folding, protein aggregations, protein DOI: http://dx.doi.org/10.5772/intechopen.95792 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.
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, S. no.

Parameter Equation Component
1 RMSD  β-bridge, bend, turn, coils, π-helix and 3 10 -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].

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.

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 (ΔG bind ) of tubulin and TauR2 was calculated by using the following Eq. (1),

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.

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.
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)

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 Figure 6.
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 noncovalent 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.
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 (R g ) and this analysis is discussed in the next section.

Compactness of tubulin-TauR2 complexes
The radius of gyration (R g ) indicates the level of compactness of the protein system which is helpful in getting an insight into the stability of the proteinprotein complex. It also helps to understand folding or unfolding of protein structure during the simulation. The R g values for all the studied tubulin-TauR2 complex ranges from 38.8-40.5 Å (Figure 7A). The complex βIII/α/βIII-TauR2 shows stable R g value for the entire simulation period however other complexes 6CVN-TauR2 6CVN*-TauR2, βI/α/βI-TauR2, βIIb/α/βIIb-TauR2 show variations in their R g values. The absence of C-terminal tail region in the complex 6CVN-TauR2 leads to the less R g values when compared to other tubulin-TauR2 ( Figure 7A). Figure 7B represents the R g values of only TauR2 in different tubulin-TauR2 complexes. The R g 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 R g value of ~18 Å and βI tubulin subunits have largest R g value of 22.5 Å as shown in Figure 7B. On the other hand, TauR2 bound to 6CVN shows uninterrupted decline in R g 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 R g 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 R g 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].

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.

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.

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 3 10 -helix during the simulation. The TauR2 in βI/α/βI-TauR2 complex does not form either α-helix or transient 3 10 -helix as shown in Figure 12C. The   TauR2 in βIIb/α/βIIb-TauR2 complex shows the formation of short-lived α-helix and 3 10 -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.  Table 5.

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   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. 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.

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.