Ligand-Based and Structure-Based Investigation for Alzheimer's Disease from Traditional Chinese Medicine

Alzheimer's disease is a neurodegenerative disease that was conventionally thought to be related to the sedimentation of beta-amyloids, but drugs designed according to this hypothesis have generally failed. That FKBP52 can reduce the accumulation of tau proteins, and that Tacrolimus can reduce the pathological changes of tau proteins are new directions away from the long held amyloid-beta-centric concept. Therefore, the screening of traditional Chinese medicine compounds for those with higher affinity towards FKBP52 than Tacrolimus may be a new direction for treating Alzheimer's disease. This study utilizes ligand-based and structure-based methods as the foundation. By utilizing dock scores and the predicted pIC50 from SVM, MLR, and Bayesian Network, several TCM compounds were selected for further analysis of their protein-ligand interactions. Daphnetoxin has higher affinity and complex structure stability than Tacrolimus; Lythrancine II exhibits the most identical trends in FKBP52 interactions as Tacrolimus, and 20-O-(2′E,4′E-decadienoyl)ingenol may be further modified at its hydrocarbon chain to promote interaction with FKBP52. In addition, we observed the residue Tyr113 of FKBP52 may play a key role in protein-ligand interaction. Our results indicate that Daphnetoxin, 20-O-(2′E,4′E-decadienoyl)ingenol, and Lythrancine II may be starting points for further modification as a new type of non-amyloid-beta-centric drug for Alzheimer's disease.


Introduction
Alzheimer's disease (AD) is the most common form of neurodegenerative disease [1][2][3] with symptoms ranging from intellectual deterioration, cognitive impairment [4,5] to abnormal behavior, personality changes, depression, and sleep disorders [6,7]. Neuronal loss, senile plaques in the cerebral cortex, and neurofibrillary tangles in the temporal and frontal lobes are prominent in the brains of AD patients. Deterioration in the limbic system, including pathological changes in the hippocampus and atrophy of pyramidal cells and the amygdale, has also been reported [8,9]. Formation of -amyloid plaques [10][11][12][13] and neurofibrillary tangles due to abnormal phosphorylation of tau proteins [14] have been linked to AD.
Pioneer research has discovered that the concentration of FKBP52, a FK506-binding protein, is of high density in brain [15] and 40 times higher within the central nervous system than in the immune system [16]. FKBP52 is not only related to immune functions but is also important for CNS protective properties. Interestingly, they can easily bind with highly phosphorylated tau proteins [16], thereby reducing the accumulation of tau proteins [16][17][18]. FKBP52 was initially discovered as a cochaperone of steroid receptor heterocomplexes [19,20] and is a member of the FK506binding protein (FKBP) of immunophilins. FKBP52, which possesses a chaperone function, has a PPIase domain also called FK506-binding domain (FKBD), composed of the first 138 amino acids from the N-terminal, exhibits peptidylprolyl isomerase (PPIase) activity, and plays an important role in regulating tau proteins.
Tacrolimus (FK506) is an immunosuppressant typically used to reduce graft rejection during organ transplants [21]. Research shows that FK506 can also reduce pathological 2 Evidence-Based Complementary and Alternative Medicine changes in tau proteins [22] and promote neuroregeneration [23,24]. In addition, FK506 can disrupt steroid receptor heterocomplexes, leading to the release of FKBP52 from the complex and enabling it to function in the nervous system [24][25][26].
In recent years, most beta-amyloid-centric drug developments have been unsuccessful during Phase III of clinical trials [27,28]. Currently, the personalized medicine and biomedicine defined an important knowledge for diagnosis and treatment analysis [29,30], such as rare diseases [31,32], regional disease [33], signal pathway [34][35][36], and gene association [37,38]. The traditional Chinese medicine is an important medicine culture in Asia defined as a personalized medicine, system biology, and biomedicine in medicine practices widely applied as treatments for cancer [39,40], cardiovascular disease [41,42], diabetes [43], virus disease [44], inflammation [45], and metabolic disorders [46], indicating its large therapeutic potential for various diseases [47]. In this moment, the cloud computation and system biology could help the screening for TCM application [48,49]. We have already designed many therapies to different targets [48,[50][51][52][53][54][55][56][57] and have successfully developed method for designing drug in disordered protein [58,59]. The aim of this research is to screen for potential non-amyloid-beta-centric leads from traditional Chinese medicine targeting FKBP52 and to investigate its mechanisms with the hope of providing important insights for forwarding Alzheimer's disease.

SVM/MLR/Bayesian
Correlation of experimental and predicted pIC 50 using the constructed SVM, MLR, and Bayesian Network models is illustrated in Figure 1. 2 values indicate good prediction accuracy of the constructed models. Table 1   Observed activity (pIC 50 )

Bayesian Network
Training set Test set 95% prediction bands interaction sites for hydrophobic interactions and H-bonds, respectively. Cross comparison reveals that TCM candidates form interactions with all Tacrolimus-interacting residues except Glu85. In addition, each TCM candidate is anchored by interactions with residues not recorded for Tacrolimus.
The higher number of interactions observed for TCM candidates might be associated with higher dock scores and suggests higher binding affinity to FKBP52 than to Tacrolimus.    Table S2. H-bond trajectories depicting time-dependent bond distance variations are illustrated in Figure 6. Since optimum H-bonds are formed between 2.3 and 3.2Å and most bond distances formed between Tacrolimus and FKBP52 residues exceed this distance, we conclude that, during MD, the protein-ligand complex stability was maintained by interactions with Glu110 and Ala 112. Daphnetoxin forms several discontinuous H-bonds with residues of FKBP52. However, the H-bond occupancies at Tyr113 with O16, O18, and H52 of Daphnetoxin are 31.47%, 2.00%, and 9.89%, respectively. It suggests that the Tyr113 might play an important role in protein-ligand complex stability. Lythrancine II and residues of FKBP52 form two stable Hbonds: Ile87 forms a stable H-bond at O36 from 6 to 17.5 ns. Tyr113 and O19 on Lythrancine II form a stable H-bond from 0 to 17.5 ns. The two residues are important for protein-ligand complex stability.

Torsion MD Simulation.
Analysis of torsion helps us understand the stability of ligand and protein binding ( Figure 7, Figure S1; see Supplementary Material available online at http://dx.doi.org/10.1155/2014/364819). Torsion (A1) of Tacrolimus changes from significant unstable rotations (0-7.5 ns) to stable rotations (7.5-17 ns), and change to unstable rotations after 17 ns. Comparing to H-bond simulation results, we can assume that Glu110 forms a more stable Hbond. Torsion (A4) is more stable from 2 to 6 ns. During this period, H-bond simulation shows the formation of an H-bond between FKBP52 and Ser118. Torsion (A5) shows some changes, but the small torsion fluctuations at (A2), (A3), (A6), and (A7) show that they have higher stability. Torsion (B8) of Daphnetoxin shows transient stability. Torsion (B9) is stable with small fluctuations. Torsion (B10) has smaller torsion fluctuations from 0 to 10 ns; torsion (B11) also has very small changes. During this period, they separately form Evidence-Based Complementary and Alternative Medicine       Figure 8 we can observe more significant changes at residue 51. In Tacrolimus, residue 51 is a bend, but in Daphnetoxin, it turns into a turn at 0-1 ns and 12-18 ns. Similarly in 20-O-(2 E,4 Edecadienoyl)ingenol, the residue becomes a turn after 1 ns and after 11 ns. In Lythrancine II, it first appears as a turn from 0 to 2 ns, then turns to bend, and then fluctuates between bends and turns. Comparisons to RMSFs show that RMSF of residue 51 is also larger. In addition, we observed that the SASA values for the four protein complexes (Figure 9(a)) and the ligands (Figure 9(b)) during MD is relatively stable, indicating no significant changes in the protein structure.

Radius of Gyration (Rg) in MD Simulation.
Radius of gyration (Rg) enables one to assess the compactness changes of a ligand-protein complex. Average and maximum values for Tacrolimus, Daphnetoxin, 20-O-(2 E,4 Edecadienoyl)ingenol, and Lythrancine II are listed in Table 3. No significant changes were observed for Rg complex (Figure 9(c)) and Rg ligand (Figure 9(d)), implying a sustained stability and compactness of the complexes. However, fluctuations in Rg complex and Rg ligand were recorded for 20-O-(2 E,4 E-decadienoyl)ingenol, suggesting a loss of compactness for its complex. Comparing to H-bond simulation results in Figure 6, the Rg complex of Tacrolimus presents maxima at 7 ns, and Gly84, Tyr113, and Ser118 trajectories also have the same trade in Figure 6(a). The Rg complex of Daphnetoxin presents maxima at 3.5 ns, and Tyr113 with H52 at 3.5 ns also presents maxima. Moreover, the Rg complex at 6 ns and 11 to 14 ns have higher values, and they correspond to Ile87, Glu85, and Tyr113 trajectories in Figure 6(b). The minima of Rg complex at 6.5 ns in 20-O-(2 E,4 E-decadienoyl)ingenol, similar to that of Tyr113, Lys121 trajectories (Figure 6(c)). In addition, the Rg complex presents maxima at 9 ns, and Tyr113 trajectory matches it. The Rg complex of Lythrancine II presents two maxima at 11 ns and 15 ns, corresponding to Asp68 and Ser118 trajectories and Tyr113 trajectory, respectively (Figure 6(d)). In conclusion, when Rg complex is higher, the compactness of ligand-protein complex is lower causing the interactions between ligand and protein to be more weak. to find the conformation with the highest frequency during the end of MD. According to our analysis results, the representative conformation for each protein-ligand complex was represented by the conformations at 19 Comparing the conformation at 0 ns and the representative conformation for each complex is shown in Figure 11.
Tacrolimus initially form an H-bond with Tyr113, but this Hbond is lost following stabilization. Moreover, the movement of Tacrolimus into the interior of the protein complex (Figure 11(a)) and the higher ligand RMSD (Figure 4) were observed. Therefore, the significant change of Tacrolimus might cause the H-bond to disappear and impact the interaction between Tacrolimus and FKBP52 deeply. The complex RMSD of Daphnetoxin is higher, but the ligand RMSD reach    dynamic equilibrium with small fluctuations (Figure 4), it might cause Daphnetoxin to form two H-bonds with Tyr113, and this was maintained throughout MD ( Figure 11(b)). The conformation of 20-O-(2 E,4 E-decadienoyl)ingenol and protein shows great differences between 0 ns and 19.58 ns at MD simulation. Especially, the position and the side chain of ligand have obvious changes. Comparing to Figure 4, its complex RMSD and ligand RMSD are higher and show larger fluctuation. As a result, 20-O-(2 E,4 E-decadienoyl)ingenol forms H-bonds with Tyr113 and Lys121 at the beginning of MD, but only the bond with Lys121 was maintained following stabilization (Figure 11(c)). Lythrancine II forms H-bonds with Asp68, Arg73, and Tyr113 at the beginning of MD, but these bonds are lost after stabilization, even though its complex RMSD and ligand RMSD are not higher. It might be attributed to the huge conformation change of the ligand. Similar to Tacrolimus, Lythrancine II was observed to be more embedded in the interior of the protein complex ( Figure 11(d)).

Discussion
Discrepancies between docking and MD simulation results were observed. During docking, Tacrolimus, Daphnetoxin, 20-O-(2 E,4 E-decadienoyl)ingenol, and Lythrancine II were predicted to form H-bonds with Tyr113 of FKBP52. Nonetheless, occupancies of the H-bond with Tyr113 varied greatly. Similarly, RMSF also shows that Tyr113 has great differences among the four complexes. In Table S2, we can see that the occupancies of Tyr113 in the Tacrolimus complex are only 0.3% and 1.6%. The occupancies of Tyr113 in the Daphnetoxin complex are 31.47% and 2.00%; therefore the H-bond can still be observed when stabilized. The occupancies of Tyr113 in the 20-O-(2 E,4 E-decadienoyl)ingenol complex are only 1.8% and 0.6%. But the bonds with Lys121 were observed during docking and were maintained during MD simulation with occupancies of 9.09% and 8.99%. The occupancies of Tyr113 with Lythrancine II complex are 9.29% and 3.5%, but H-bonds were not observed in the stabilized complex.
From other MD simulation analysis we can make the following summarizations. Despite the larger protein complex RMSD of Daphnetoxin compared to Tacrolimus, the H-bond with Tyr113 was presented during both docking and MD simulation, indicating a better affinity towards FKBP52 than Tacrolimus. Summarizing the results from RMSF, SASA, and Rg, the protein complex of Daphnetoxin is stable during MD simulation.

Conclusion
We aimed to utilize structure-based and ligand-based methods to screen for high affinity lead compounds for FKBP52 as alternatives for Tacrolimus   corrected for H-atoms using the Prepare Ligands module in Discovery Studio Client v2.5 (DS2.5; Accelrys Inc., San Diego, CA). Ligand structures and activity data from Gopalakrishnan's study [61] were used to construct quantitative structureactivity relationship (QSAR) prediction models for FKBP52. TCM ligands used for virtual screening were downloaded from TCM Database@Taiwan [62]. Tacrolimus (FK506), the clinically used immunosuppressant, was used as the control. applied to adjust ionization states of the downloaded ligands. Downloaded ligands and Tacrolimus were docked to the PPIase domain of FKBP52 under a forcefield of Chemistry at HARvard Molecular Mechanics (CHARMm) using the LigandFit module. Ligands were evaluated based on structural compatibility to the PPIase domain and dock score was selected as the scoring function.

Support Vector Machine (SVM) and Multiple Linear Regression (MLR) Prediction
Models. The 37 ligands [61] were randomly divided into a training set of 30 compounds and a test set of 7 compounds. Molecular properties of each ligand were calculated with Calculate Molecular Properties and the ten representative descriptors most related to bioactivity were determined by Genetic Function Approximation (GFA). Different prediction models were constructed with the representative descriptors and the strength of each model ranked by a square correlation coefficient ( 2 ). Descriptors in the highest 2 model were used to construct a nonlinear SVM model with LibSVM and a linear MLR model with MATLAB. The constructed models were validated with the test sets and applied to predict bioactivity of selected TCM candidates.

Bayesian Network.
Training/test set groups and representative descriptors determined by GFA were also used to construct the Bayesian Network model. According to distribution characteristics, descriptors and pIC 50 were discretized into a maximum of five categories. Linear regression analysis for each pIC 50 category in the training dataset was then applied. Banjo [63] was used to discover relationships among the representative descriptors and pIC 50 values. The algorithm for predicting pIC 50 was written in MATLAB codes integrating Banjo and Bayes Net Toolbox (BNT; https://code.google.com/p/bnt/). Following validation, the algorithm was applied to predict bioactivity of selected TCM candidates.

Candidate Selection Criteria.
All TCM ligands with dock scores lower than Tacrolimus were eliminated. Next, a consensus voting system was used for candidate selection. The three highest scoring candidates in dock score, SVM, or MLR were given a score of "1"; candidates with the three highest Bayesian Network scores were given a score of "2. " Sum of scores was calculated and the three ligands with the highest scores were selected as candidates for further analysis.

Molecular Dynamics (MD) Simulation.
Candidates and Tacrolimus were prepared with SwissParam (http://swissparam.ch/) [64] prior to MD. MD simulation was conducted with GROMACS 4.0.7 under the forcefield of CHARMm27. Distance of protein to box boundaries was 1.2 nm; solvate TIP3P water model was then added to the system. Complex charges were neutralized with sodium and chloride ions using 0.145 M salt model. Simulation was conducted at 310 K under a pressure of 1 bar. Each complex was minimized with 5,000 steps of Steepest Descent, and the final minimized structure was used as the initial structure for MD simulation. Electrostatic interactions were calculated with Particle-Mesh-Ewald (PME) [65]. Equilibration protocol was used to restrain and relax protein-ligand position; first-order kinetics was started from 300 K. Minimized system was used to simulate a five-thousand ps configuration production. MD simulation was conducted for 20 ns with time steps of 2 fs under PME.