The structural, functional, and dynamic effect of Tau tubulin kinase1 upon a mutation: A neuro‐degenerative hotspot

Alzheimer's disease (AD) is a progressive disorder that causes brain cells to degenerate and die. AD is one of the common causes of dementia that leads to a decline in thinking, behavioral and social skills that disrupts a person's ability to function independently. Tau‐tubulin kinase1 (TTBK1) is a crucial disease regulating AD protein, which is majorly responsible for the phosphorylation and accumulation of tau protein at specific Serine/Threonine residues found in paired helical filaments, suggesting its role in tauopathy. TTBK1 involvement in many diseases and the restricted expression of TTBK1 to the central nervous system (CNS) makes TTBK1 an attractive therapeutic target for tauopathies. The genetic variations in TTBK1 are primarily involved in the TTBK1 pathogenesis. This study highlighted the destabilizing, damaging and deleterious effect of the mutation R142Q on TTBK1 structure through computational predictions and molecular dynamics simulations. The protein deviation, fluctuations, conformational dynamics, solvent accessibility, hydrogen bonding, and the residue‐residue mapping confirmed the mutant effect to cause structural aberrations, suggesting overall destabilization due to the protein mutation. The presence of well‐defined free energy minima was observed in TTBK1‐wild type, as opposed to that in the R142Q mutant, reflecting structural deterioration. The overall findings from the study reveal that the presence of R142Q mutation on TTBK1 is responsible for the structural instability, leading to disruption of its biological functions. The mutation could be used as future diagnostic markers in treating AD.


| INTRODUCTION
The hyperphosphorylation of tau is one of the signature hotspots of Alzheimer's disease (AD). The results of tau hyperphosphorylation or its accumulation results in aggregation, leading to reduced affinity to microtubules that leads to the destruction of tau-associated cellular activities, which abnormally affects axonal growth, vesicle transport, and signal propagation mediated by microtubules. 1,2 More than 20 kinases are involved in phosphorylation of AD sites associated with tau formation and subsequently neurodegeneration. Tau-tubulin kinase1 (TTBK1) belongs to the casein kinase1 superfamily. It is a highly expressed protein observed in the granular cell layer, the hippocampus, the midbrain cerebellum Purkinje cells and the substantia nigra. 3,4 The protein can phosphorylate tau and tubulin at ten different sites that are majorly associated with AD.
TTBK1 phosphorylates tau at Ser198, Ser199, Ser202, and Ser422 2 and can express highly in many tissues, including liver, pancreas, skeletal muscle, kidney, heart, and found mainly in the cerebellum of the brain. 5 TTBK1 is upregulated in the frontal cortex of patients with AD. 6 Apart from AD, the protein TTBK1 is involved in cancer progression, TDP-43 accumulation and transporter stimulation, thereby it is considered to be an important drug target ( Figure 1). 7,8 The transgenic mouse lines overexpressing TTBK1 demonstrated significant age-dependent memory impairments that are also dependent on the dosage of TTBK1 6,9 Interestingly, TTBK1 overexpression induced increase in neurodegeneration has been observed in different animal model systems including mice. 8,10,11 Moreover, genetic variations of TTBK1 (SNPs: rs2651206, rs10807287, and rs7764257) are associated with late-onset Alzheimer's disease (LOAD) observed in two large cohorts of Spanish and Chinese populations, 12,13 indicating the importance of the TTBK1 gene in the development of tauopathy and AD pathogenesis.
At the molecular level, the structure of TTBK1 resembles other protein kinases due to an enriched β strand N-terminal domain and an α-helical C-terminal domain. An extended "hinge" region is present at residues 108-111, connecting both terminals. The P-loop located at residues 40-49 is also a part of the N-terminal domain, while Asp-Phe-Gly, also called as DFG motif, is situated at positions 176-178, and a flexible activation loop (residues 178-202) is located at the fragments of the C-terminal. 14 TTBK1 mutants have been The pictorial representation of the workflow to study the effect of TTBK1 R142Q mutation on the protein stability, using prediction algorithms and MD simulations reported in the literature; however, the precise structural and functional relationship of the mutants remains to be studied. 5,9,14 By analyzing the 1000 genome browser, we were able to identify over 180 missense variations in TTBK1. The prediction of pathogenicity potential of these variations by different in-silico tools has suggested R118Q (R142Q in PDB ID: 4BTK) as one of the most pathogenic and deleterious mutations. The significant structural and dynamic impact of R142Q mutation is studied here. Molecular dynamics (MD) simulations are reliable techniques that can provide the evidence to check the stability of proteins in biological environments through computational validations. [15][16][17] Hence, we relied on MD simulations using GROMACS to check the dynamical behavior of mutant R142Q on TTBK1 protein. In the current study, we executed several parameters, namely, root mean square deviation (RMSD), root mean square fluctuation (RMSF), energy parameters, solvent accessible surface area (SASA), intramolecular hydrogen bond, radius of gyration (Rg), secondary structure element (SSE), principal component analysis (PCA), free energy landscape (FEL), density distribution, dynamic cross-correlation matrices (DCCM) and residue frustration which collectively confirmed the destabilization of the TTBK1 native structure ( Figure 1). 18,19 The present study revealed that the selected TTBK1 R142Q mutant exerts a destabilization effect on the protein, leading to structural and conformational destabilization, thereby negatively regulating the protein activity in association with AD.

| Hardware and software
The MD simulations were carried out on a highperformance computing IBM server with 128GB RAM in each node and four CUDA-enabled NVIDIA (Model: Nvidia Tesla V100) graphics processing units with 16GB RAM each. UniProt was used to retrieve the sequence information, the PDB database 20 was used to retrieve the 3D structures and PyMOL 21 was utilized for visualization. Bioinformatics tools and web servers were used to check the effect of the mutation on protein structure. MD simulations were carried out using GROningen Machine for Chemical Simulations (GROMACS) 5.13.3 software package. 22

| Prediction algorithms
The online tools sorting intolerant from tolerant (SIFT), 23 POLYPHEN, 24 I-mutant, 25 PROVEAN, 26 SNP&GO, 27 SNAP2, 28 and Pmut tools 29 were used to predict the functional impact of selected mutations on TTBK1 protein. 30 Further, the destabilizing effects of R142Q at the substituted position 142 was also examined by predicted changes in folding Gibbs free energy (ΔΔG) by DUET, 31 mCSM, 32 ENCoM, 33 DynaMut, 34 SDM, 35 Mutpred, 36 Fathmm, 37 and MuPro web 38 servers provided the threshold scores effectively. 39 From these prediction algorithm examinations, the destabilizing, the damaging and deleterious role of R142Q on TTBK1 was examined.

| Wild type (WT) and R142Q structures
The crystal structure 4BTK of human TTBK1 consisting of 337 amino acid residues was retrieved from PDB and was used as WT in the study. We mapped the deleterious mutation R142Q on TTBK1 at the corresponding position and was mutated using in silico mutagenesis plug-in embedded in SWISS-PDB Viewer. 40

| MD simulations
MD simulations were carried out by GROMACS 5.18.3 software package. 22 The determination of charge states for the ionizable residues was performed using H ++ calculation to understand the protein's protonation state. 41 The WT and R142Q topology parameter files were created by the CHARMM27 force field included with CMAP correction. 42 The intermolecular (nonbonded) potential, represented as the sum of Lennard-Jones force-based switching with a cutoff distance range of 8-10 Å, pairwise Coulomb interaction and the long-range electrostatic force were determined by particle mesh Ewald (PME) approach. 43 The real-space cutoff was set to 1.2 nm for the PME calculations. 44 The velocity Verlet algorithm was applied for the numerical integrations and the initial atomic velocities were generated with a Maxwellian distribution at the given absolute temperature. Then the system immersed with the default TIP3P water model, and protein was placed at the center of the cubic grid box (1.0 nm 3 ). 45 The neutralization was performed to make the concentration of the system to 0.15 M. The neutralized system was then subjected to energy minimization using the Steepest Descent and Conjugate Gradient algorithms utilizing a convergence criterion of 0.005 kcal mol-1 Å-1. The twostandard equilibration phase was carried out separately constant number of particles, system volume and temperature (NVT) and constant number of particles, system pressure and temperature (NPT) ensemble conditions such as constant volume and constant pressure for each complex similar time scale with LINear Constraint Solver to restrain the bonds involving hydrogen atoms. Nosé-Hoover thermostat and Parinello-Rahman barostat was applied to maintain the system's temperature and pressure, respectively. 46,47 The system was maintained constant at 1 bar and 300 K, with a coupling time of τP = 2 ps and τT = 0.1 ps, respectively. The periodic boundary condition used for integrating the equation of motion by applying the leap-frog algorithm with a 2fs time step. Finally, to make the system ready for production, the fixing of constraints was achieved with the relaxation of the grid box with water and counterions. 48 The run of 100 ns was performed and all the trajectories were recorded at every 4 ps for further validation. The calculations of RMSD, RMSF, SASA, Rg, intramolecular hydrogen bond, FEL and density distributions were analyzed. 49,50 The MD analysis, with the parameters mentioned above, was conducted using GROMACS toolsets. The energy plots and graphs were marked and visualized using the QtGrace visualization tool.

| Principal component analysis
PCA is a widely used statistical technique that reduces complexity in analyzing the data and extracts the rigorous motion information in simulations, which is essentially correlated for biological function. 51,52 PCA works on the variance/covariance matrix constructed based on the trajectories after removing the rotational and translational movements with their projection and the first two principal components using the Essential Dynamics method. The dynamic protein motions of WT and R142Q in the subspace were identified by Cartesian trajectory coordinates, projecting the most important eigenvectors from the complete analysis. Principal components were selected as reaction coordinates and the free energy of the state, Gα, was calculated using: where, k is the Boltzmann constant, T is the temperature of simulation. P(qα) corresponds to the estimate of probability density function extracted from the histogram of the MD data that was constructed using joint probability distribution (PC1 and PC2 as reaction coordinates). Pmax(q) represents the probability of the most probable state.

| Screening and destabilizing effect of R142Q
The functional effects of R142Q on TTBK1 protein were predicted using SIFT, PolyPhen, I-mutant 2.0, PRO-VEAN, PhD SNP online servers (Table 1). SIFT analysis provided the tolerance index score of R142Q between 0 and 1, revealing its deleterious effect on the amino acid substitution at the respective positions on the protein, ranked based on the aligned sequences. PolyPhen predicted R142Q mutation to be deleterious with a score of 0.97. I-mutant predicted R142Q as a destabilizing mutation (ΔG = −1.54). The PROVEAN scores also confirmed that the R142Q mutation confers a damaging effect on the native structure. The PHD-SNP method also projected that the selected R142Q mutation falls under disease-related variants based on the local sequence environment. The screening results indicate the deleterious and pathogenic effects of the R142Q on the 3D structure of TTBK1 by altering its stability and function. The impact on stability due to changes in R142Q of TTBK1 was further evaluated using ΔΔG calculation with DUET, mCSM, EnCOM, DynaMut, Fathmm, MutPred, MuPro, and SDM web servers ( Table 1). The results depicted that the mutation renders destabilizing effects and deleterious impact on the stability of the TTBK1 native structure (Figure 2A).

| Protein deviation and RMSF
The MD trajectories of WT and R142Q structures were analyzed by RMSD calculations. The comparative structural deviation analysis projected that WT-TTBK1 reached equilibrium at~5000 ps and remained unchanged throughout the simulation. The MD simulation of the WT noticeably displayed a steady RMSD of ∼0.19 nm throughout the simulation, whereas R142Q shows a fluctuation of ∼0.25 nm ( Figure 2B). However, for the R142Q mutant structure, equilibrium at~9000 ps and large deviations were observed from 30000 ps to 90000 ps run time, clearly showing nonsignificant changes over the Cα-atom. The further investigation proceeded with the calculations of the probability distribution function (PDF) of RMSD for WT and R142Q ( Figure 2C and Table S1). The conformers stretched a range of 0.1-0.35 nm for WT and R142Q of TTBK1. The PDF analysis portrayed the significant effect of R142Q mutation on the conformation of TTBK1 protein. It was also evident from the analysis that a single amino acid change affected the protein's conformations due to the substitution of a charged residue Arginine at the position 142 by an uncharged polar Glutamine. Summarily, our analysis results also revealed a maximum effect on the flexibility of TTBK1 native structure upon the mutation. The R142Q presented deleterious effects upon the amino acid alteration on TTBK1 protein, found to be damaging and destabilizing for TTBK1, altering its structure and presumably affecting its pathophysiology. The overall flexibility of the WT and R142Q structures were analyzed with Cα-RMSF fluctuations. The negative drift of the graph indicated the increased movement of atoms on TTBK1 by R142Q that are predicted to cause instability of the protein compared to WT ( Figure 2D and Table S2). The analysis revealed that the fluctuations were lower at the mutated residue and the internal motion of the residue was averaged to get the comparative graph between R142Q and WT. Thus, it is observed that R142Q brings about higher fluctuations mainly in the region around ∼90-95, ∼200-230, and ∼260-270, which is the consequence of the positioned mutation at position 142 (∼0.08 nm). The results revealed a disturbance in fluctuation with a 0.08 nm range which further imply that the substitution R142Q imparts rigidity of TTBK1, leading to the most significant degree of flexibility loss to the WT structure.

| SASA and Rg analysis
The SASA and the compactness of structures were monitored for WT and R142Q structures of TTBK1 through the MD simulation. The SASA analysis revealed that R142Q had a higher range of SASA values, whereas WT showed a less SASA range that explains the protein compactness. The PDF of SASA revealed higher drifts for R142Q with a range of 173.96 nm 2 , whereas the same for WT was 173.77 nm 2 . ( Figure 3A) and (Table S1).
The decrease in SASA for R142Q was noticed when compared to WT, indicating that the mutation significantly alters the compactness of TTBK1.
Rg analysis was carried out to investigate conformational changes in WT and R142Q further. Rg reflects compactness of a protein by the mass-weight root mean square distance of a collection of atoms from their common center of mass. The analysis displayed that the Rg of WT TTBK1 presents a smooth trend and a slight decrease from~2 to~1.95 nm and a small increment till 100,000 ps ( Figure 3B) and (Table S1). Contrastingly, the Rg of R142Q showed higher fluctuations. The mutant is confirmed to show a loss of Rg between the atoms with an average Rg of 1.98 nm when compared to that of WT with an average Rg of 1.99 nm ( Figure S1). Thus, the SASA analysis and Rg data reflect the overall conformational stability of WT TTBK1 and R142Q. However, the TTBK1 mutant displayed increased Rg towards the end of the simulation, leading to structural destabilization, loss of shape and protein compactness, in contrast to the WT.

| Secondary structure, H-bond, and PCA analysis
The change in the secondary structure during the MD simulation analysis was examined further to understand the SSE, contributing to the shape of the protein 3D structure. From the SSE analysis, it was evident that mutant R142Q projected deformation from helix to coil, sheets to coil, bends to coil and helix to bend, revealing a loss of structural stability, consequentially leading to differential changes in the TTBK1 structure, when compared to WT ( Figure 4A) and R142Q (Table S3). To a large extent, the stability of a protein is determined by the number of hydrogen bonds formed between its atoms. Hence, we calculated the number of intramolecular hydrogen bonds formed in the WT and R142Q. The overall analysis showed R142Q has 220.25, whereas WT has 219.31 H-bonds ( Figure 4B and Table S1). The WT displayed a higher number of hydrogen bonds, with no alterations in the intramolecular bond numbers, revealing strong contacts between the residues, whereas, in contrast, the mutant displayed a smaller number of hydrogen bonds leading to loss of contacts between the atoms. Further, the simulation trajectories were examined for the essential subspace where the protein dynamics are studied by eigenvectors associated with the eigenvalues. The dynamic behavior of WT and R142Q, on the basis of the comparison of the clusters, revealed that the clusters are well defined in the WT structures by covering the minimum region whereas mutant occupied maximum regions. We further investigated eigen-vector1 and eigenvector2, which were plotted to understand the atomic fluctuation, which indicates the nature of motion ( Figure 4C) and (Figures S2 and S3). The results revealed that R142Q showed maximum fluctuation in comparison to WT throughout the MD simulations. The comparative study of WT and R142Q revealed conformational changes in R142Q mutant due to significant disturbance in internal atomic motion, causing loss of TTBK1 stability.

| FEL, density distribution, and DCCMs analysis
The free energy changes of WT and R142Q were explored with the MD simulations using FEL analysis. The comparison studies indicated stable global free energy minima confined within one basin in WT with a range of 11.6 kJ/mol. However, R142Q has wider basins of sample conformations with numerous meta-stable conformations associated with multiple energy minima with a range of 12.8 kJ/mol ( Figure 5A [i-ii] and Figure S4). The results revealed that the R142Q triggered structural destabilization of TTBK1 and also marked significant differences in the folding behavior leading to unsteady conformation and stability loss.
The atomic density distribution analysis was further investigated to understand the atomic orientation and distribution changes plotted using densmap. The partial density area of WT was observed to be stable with minimum values of 3.46 nm −3 whereas a high density was noticed on mutation R142Q with a range of 3.68 nm −3 ( Figure 5B [i-ii]). The comparison of density distribution between WT and R142Q confirmed a significant impact on the structural transition of TTBK1, leading to stability loss upon the amino acid substitution.
The conformational motions of the WT and R142Q of TTBK1 were monitored by the DCCMs analysis. 53 The analysis was highlighted with two colour representations, where red indicates a positive correlation of residue movements in the same direction and blue indicates anticorrelation of residue movements in the opposite direction ( Figure 5C [i-ii]). The results revealed that the R142Q mutant presented anticorrelation movements of the residues whereas, WT showed a strong correlation. Furthermore, the overall comparative DCCM revealed a significant loss of contacts among the residues, which indicated stabilization of TTBK1 structure upon R142Q mutation.

| Diagonalized covariance matrix and frustration analysis
The WT and R142Q mutant structures were further investigated for diagonalized covariance matrix to understand the positional fluctuations of Cα-atoms on TTBK1.
The atomic behavior of the residual motions was differentiated with red and blue color representation, where red implies small fluctuations between the atoms and blue imply large fluctuations. The amplitude and the intensity of WT magnified with a value of 0.056 nm 2 , whereas higher differences were observed in R142Q ranged with 0.157 nm 2 ( Figure S5). The comparison of residual displacement between WT and R142Q revealed the structural and molecular deformation of TTBK1 on the native structure. A porcupine plot was generated to explain the differences in motions and magnitude between the WT and R142Q ( Figure 6A [i-ii]). The direction of the arrow indicates the direction of the eigenvector and the magnitude of the corresponding value. The larger deviation patterns were observed in R142Q due to the unfolded structure of TTBK1, which is crucial for proper function. The porcupine plot of the first eigenvector also presented the movements and magnitude of Cα-atoms obtained through PCA analysis of WT and R142Q. The residue frustrations of WT and R142Q were compared and resulted with a well global and local energy minima ( Figure S6). The results showed maximum frustrations in mutations, noticed from re-sidues~100-150,~150-200, and~350 for R142Q, leading  Table S3). The residue frustration analysis presented a significant unsteady shift in the patterns at substitution positions of R142Q with inclined loss of contacts, further leading to the deterioration of the TTBK1 structure.

| DISCUSSION
Several kinases phosphorylate tau proteins. Most of these kinases phosphorylate AD sites, out of which TTBK1 and TTBK2 are most common. 8,48 The genetic variations of the TTBK1 have been associated with LOAD, observed in various communities. 12,13 Since protein mutations bring about changes in its folding, dynamics, stability and functions, some of these mutations often manifest in pathological consequences. [54][55][56] In this study, we provide insights toward the effects of point mutation, R142Q on the folding, dynamics and stability of TTBK1, which may be responsible for the development of neurodegenerative diseases like AD, amyotrophic lateral sclerosis or spinocerebellar ataxia Type 11. 1,2 The mutational analysis predicts the mutation R142Q to confer a significant damaging and deleterious effect on the native TTBK1 structure, leading to destabilization of the protein. This destabilization is also corroborated from the free energy calculation (ΔΔG) scores which also implicated a steep loss of protein stability, presumably promoting the functional loss of TTBK1 protein. Performing molecular simulations of TTBK1 WT and R142Q strongly portrayed the molecular effects of the TTBK1, monitored by various parameters like RMSD, RMSF, Rg, SASA, SSE, and FEL analysis. The R142Q mutant caused inclined flexibility and deviations, evident from the RMSF and RMSD analysis and its comparison with the WT. The interchange of charges can potentially cause a loss of a helix upon the point mutation at respective amino acid residue position 142, which also corroborates with outcomes of the stability predictions. Based on the obtained results, it was also evident that high flexible conformations were induced due to R142Q, which may facilitate entry of the drugs/inhibitors/ligands in the binding cleft of the TTBK1. The large oscillations in Rg and low SASA values displayed that the TTBK1 native structure might be undergoing a significant structural transition which is also supported by conformational dynamics, residual motion and atomic density map analysis. During the MD simulations on WT, the intramolecular hydrogen bond formation displayed a strong relation in the bond formation, whereas R142Q demonstrated a smaller number of H-bonds formation throughout the simulation.
Furthermore, the corresponding eigenvalues' results also projected the level of variation and dynamic nature of TTBK1 in the simulation system, where basins were entirely restricted on WT and the wide basins with multiple wells were observed in R142Q, indicating loss of stability. The secondary structure changes in R142Q displayed high structural plasticity, reflecting a high degree of flexibility of TTBK1 protein compared to the WT. The analysis of conformational dynamics, frustration, residual motion and magnitude demonstrates the loss of native protein structure in the presence of mutant when compared to the WT. Overall, the MD results illustrated that the mutation R142Q significantly affects the structural and dynamic behavior of TTBK1 to cause structural disruptions.
Summarily, from the analysis of computational algorithm predictions and dynamics, it may be concluded that R142Q imparts destabilization of the TTBK1 structure. The substitution of amino acid residue Arginine (R), a positively charged, polar, hydrophilic, to amino acid residue Glutamine (Q), a polar, non-charged, probably induces electrostatic changes inducing conformational changes in the TTBK1 structure.

| CONCLUSION
The combination of computational algorithms and MD simulation studies facilitated our understanding of the effect caused by the mutant on TTBK1 contributing to the diseased state of AD. The implementations of MD simulations give a molecular basis of the changes induced by the important TTBK1 mutation. The R142Q mutant displayed a destabilizing, damaging and deleterious effect on the TTBK1 protein structure-function compared to WT. To conclude, the study provided a comprehensive view of destabilizing mechanisms of the R142Q mutant on TTBK1, which gives valuable insights for future investigations and experimental evaluations for designing AD drugs/inhibitors.