Computational Design of Apolipoprotein E4 Inhibitors for Alzheimer's Disease Therapy from Traditional Chinese Medicine

Apolipoprotein E4 (Apo E4) is the major genetic risk factor in the causation of Alzheimer's disease (AD). In this study we utilize virtual screening of the world's largest traditional Chinese medicine (TCM) database and investigate potential compounds for the inhibition of ApoE4. We present the top three TCM candidates: Solapalmitine, Isodesacetyluvaricin, and Budmunchiamine L5 for further investigation. Dynamics analysis and molecular dynamics (MD) simulation were used to simulate protein-ligand complexes for observing the interactions and protein variations. Budmunchiamine L5 did not have the highest score from virtual screening; however, the dynamics pose is similar to the initial docking pose after MD simulation. Trajectory analysis reveals that Budmunchiamine L5 was stable over all simulation times. The migration distance of Budmunchiamine L5 illustrates that docked ligands are not variable from the initial docked site. Interestingly, Arg158 was observed to form H-bonds with Budmunchiamine L5 in the docking pose and MD snapshot, which indicates that the TCM compounds could stably bind to ApoE4. Our results show that Budmunchiamine L5 has good absorption, blood brain barrier (BBB) penetration, and less toxicity according to absorption, distribution, metabolism, excretion, and toxicity (ADMET) prediction and could, therefore, be safely used for developing novel ApoE4 inhibitors.


Introduction
Alzheimer's disease (AD) is the most common harmful neurological disorder affecting patients over the age of 65 [1]. The major neuropathological hallmarks of AD are neurofibrillary tangles and beta amyloid plaques in the entorhinal cortex and hippocampus [2]. Deposition of -amyloid is associated with neuronal loss and cognitive dysfunction [3], which affects the ability to study and think. There is also long-term memory loss during the life span. -amyloid is generated from cleavage of the amyloid precursor protein (APP) by -secretase and -secretase [4][5][6][7][8]. Several genes linked to familial AD have been identified in many studies, but there are still more risk-factor genes that remain to be identified [9]. The variations of apolipoprotein E (AOPE) genes as potential genetic risk factors for AD have been identified by the linkage studies of Pericak-Vance et al. [10]. The ApoE gene is located on chromosome 19q13.2, and an increase in type 4 allele of this gene is associated with several chronic neurodegenerative diseases [11].

Data Sets.
We used 61,000 TCM compounds for database screening, which were obtained from TCM Database@Taiwan [57], with the drug-likeness of all compounds being predicted by ADMET prediction. The crystal structure of ApoE4 was taken from PDB database (PDB code:1GS9), the missing atoms and loops were cleaned up by Prepare Protein module under Accelrys Discovery Studio 2.5.5.9350 (DS 2.5) [58], and all residues were protonated under pH 7.4 conditions. We also employed disorder predict tool (PONDR-FIT) [59] to predict unfolded regions on ApoE4 sequence for structure validation.

Docking
Analysis. The LibDock program [60] of DS 2.5 was used to define protein site features referred to polar and nonpolar features, with a sphere of 35Å radius used as the binding area. Different rigid ligand conformations were placed into the binding area, and all ligand conformations were minimized using the CHARMm force field.  Protein gyrate (nm)  Figure 5: Plots of (a) protein RMSD, (b) ligand MSD, and (c) radius of gyration from ApoE4 with docked ligand or no ligand (apo) with a simulation time of 5000 ps.

Molecular Dynamics
Simulation. The molecular dynamic simulation was performed with GROMACS 4.5.5 package [63] for protein-ligand complexes simulation and the charmm27 force field was used in the simulation system. For box definition, distance of real space cut-off was set to 1.2 nm. The particle mesh Ewald (PME) method was regarded as coulomb type for treating electrostatics, and the cut-off distance of defining van der Waals (VDW) residues was set at 1.4 nm. In pair potentials versus many-body potentials [64][65][66][67], the potential functions representing the nonbonded energy of VDW and electrostatics using the following: All bonds were constrained with the linear constraint solver (LINCS) algorithm for fixing all bond lengths. The TIP3P model was employed for water simulation. Topology files and parameters of small compounds in protein-ligand complexes were generated for GROMACS simulation by the SwissParam web server; the concentration of NaCl model was set as 0.145 M in the solvent system. The Steepest Descent algorithm performed the 5,000 cycle steps used for energy minimization. This was followed by equilibration performed under position restraints for 1 ns under constant temperature dynamics (NVT type) conditions. Following this step, all production dynamics simulations were performed for 5000 ps under constant pressure and temperature dynamics (NPT type). The temperature in all of the simulation systems was set as 310 K. MD conformations were saved every 20 ps for trajectory analysis.

Docking Analysis.
We used the LibDock score to select potent TCM compounds which have high affinity with ApoE4. The results of the docking score are listed in Table 1 and all candidates are ranked by LibDock Score. We also used ADMET prediction to evaluate the drug-likeness of TCM candidates. All TCM candidates had a good absorption prediction for metabolism. In toxicity evaluation no TCM  candidate displayed CYP2D6 inhibiting and hepatotoxicity, suggesting that these ligands have no toxicity in the liver. Blood brain barrier (BBB) penetration showed that these ligands have good penetration and may be suitable for central nerve system therapy. We regard the top three TCM candidates, Solapalmitine, Isodesacetyluvaricin, and Budmunchiamine L5, as potential compounds and further analyze the binding poses in ApoE4 protein structure. Scaffold of the top three TCM compounds is shown in Figure 1, and the docking poses of each ligand are displayed in Figure 2. For H-bond analysis, Solapalmitine had no H-bonds with any residue of ApoE4. The docking pose of Isodesacetyluvaricin revealed H-bonds with Lys157 and ARG158. For Budmunchiamine L5, the docked ligand formed two H-bonds with Glu77. In the 2D diagram of docking poses by Ligplot (Figure 3), Solapalmitine did not reveal any H-bonds, while the surrounding residues, Asp153, Gln156, Gln27, Tyr74, Gln24, Lys157, Val161, Leu93, Tyr162, Glu96, and Arg158, displayed hydrophobic interactions with the ligand. The 2D diagram of Isodesacetyluvaricin was similar to the 3D docking poses, which form H-bond with Arg158; hydrophobic interaction residues included Arg150, Asp153, Asp154, Tyr162, GLu96, Val161, Lys157, Leu93, and Arg158. For Budmunchiamine L5, the 2D diagram is corrected with 3D docking pose, and only Glu77 can form H-bonds with the docked ligand; this ligand also interacted with Ala73, Arg25, Glu70, Tyr74, Glu77, Gln24, Leu159, Gln156, Ala160, Glu27, Trp26, and Gln163 by hydrophobic forces.
In comparison with the docking study, PONDR-FIT was used to predict the disorder region among all residues on ApoE4 ( Figure 4). According to the drug design of protein structure [68], ordered folding regions do not influence ligand binding [69,70]. Most of binding residues for the top three compounds were below 0.5, suggesting that the structure of ApoE4 will not affect ligands binding to the active site. To further analyze the variation of the docked protein structure, we utilized MD simulation to generate a dynamic structure for binding analysis of the top three TCM compounds.

Stability Analysis.
Complexes of ApoE4 with docked ligands were performed by MD simulation at 5000 ps, and ApoE4 with no ligand (Apo protein) were regarded as the control for comparison. Each plot of the root mean square deviation (RMSD), mean square displacement (MSD), and radius of gyration ( ) is displayed in Figure 5. The Apo protein of ApoE4 changed significantly from 0.20 to 0.27 nm after 4000 ps, indicating that ApoE4 with ligand was more stable during MD simulation. The MSD analysis was used to calculate the migration distance among all simulation times. The MSD values for Budmunchiamine L5 increased from 1500 to 5000 ps, which displayed a similar movement distance to Isodesacetyluvaricin at final simulation time. Solapalmitine had an increased MSD value with a simulation time of 5000 ps, suggesting that Solapalmitine is less stable than the other two TCM candidates. The gyration assay was used to analyze the compactness of protein structure.
values of all protein-ligand complexes and Apo protein had similar fluctuations, indicating all structures tended to become stable after MD simulation. For total energy analysis, no significantly increased values were observed among all simulation times (Figure 6). The total energy of all systems remained in −876000 kJ/mol. These results suggest that all structures of the complexes tend to become constant after the initial simulation time.

Residues Fluctuation and Distance
Analysis. Root mean squared fluctuation (RMSF) was carried out to analyze the fluctuation of residues on ApoE4 protein (Figure 7). It is obvious that residues of Apo proteinfrom 70 to 100 exhibit substantial fluctuation, but the three candidates remain stable.
The ligand binding region is included in this region, but the docked residues are not flexible due to the largest fluctuations being exhibited at terminal residues, and these regions are far from the docked residues. The results suggest that the docked ligand could bind stably to ApoE4. The matrices of distance maps for residue-residue distance calculations over 5000 ps are shown in Figure 8. The results display that all complexes with docked ligands are the same as Apo protein, suggesting that the conformations do change among all MD simulations.

Clustering Analysis for Snapshot Selection.
A cluster algorithm was employed to select the most stable conformation over all simulation times. All MD snapshots with docked ligands were grouped into two or four individual clusters ( Figure 9) with each cluster including similar conformations. The middle structure was chosen from each late group as a standard snapshot for the binding analysis of the top three TCM compounds: Solapalmitine (4040 ps), Isodesacetyluvaricin (4240 ps), and Budmunchiamine L5 (4340 ps) ( Table 2). In the binding analysis from the standard snapshot ( Figure 10), we found that Budmunchiamine L5 exhibited H-bond interaction with Glu77, but the H-bonds disappeared between Isodesacetyluvaricin and Arg158. Solapalmitine continued to have no H-bond interactions between residues. The data shows that Glu77 is a crucial residue for Budmunchiamine L5 binding. We also employed CAVER 3.0 software [71] to analyze the migrated ligand tunnels in ApoE4 (Figure 11), and the ligand pathway analysis was used to predict in previous studies [40]. The prediction of ApoE4 showed lower number of channels than Apo form of ApoE4; the results illustrated that the TCM candidates could form stable binding conformation to interact with ApoE4 with all simulation time.

Conclusion
Solapalmitine is the top candidate by LibDock score but displays significantly increasing MSD values due to unstable binding with ApoE4 over 5000 ps simulation time. Isodesacetyluvaricin has H-bond with Arg158; unfortunately, the H-bond is missing during MD simulation. The LibDock score of Budmunchiamine L5 is not the highest score from the TCM database screening; however, the dynamics simulation shows that the docked ligand complex of ApoE4 is stable. In snapshot analysis, Budmunchiamine L5 still forms an Hbond with Glu77; the binding pose is the same as the initial docking pose, suggesting that Budmunchiamine L5 does not change over all simulation times and stably binds to ApoE4. In terms of ADMET analysis, Budmunchiamine L5 has good absorption, BBB penetration, and less toxicity in the human liver and may therefore be regarded as a safe lead drug for designing a novel ApoE4 inhibitor for AD therapy.