Metal based donepezil analogues designed to inhibit human acetylcholinesterase for Alzheimer’s disease

Among neurodegenerative disorders, Alzheimer’s disease (AD) is one of the most common disorders showing slow progressive cognitive decline. Targeting acetylcholinesterase (AChE) is one of the major strategies for AD therapeutics, as cholinergic pathways in the cerebral cortex and basal forebrain are compromised. Herein, we report the design of some copper and other metal based donepezil derivatives, employing density functional theory (DFT). All designed compounds are optimized at the B3LYP/SDD level of theory. Dipole moments, electronic energie, enthalpies, Gibbs free energies, and HOMO-LUMO gaps of these modified compounds are also investigated in the subsequent analysis. The molecules were then subjected to molecular docking analysis with AChE to study the molecular interactions broadly. Ensemble based docking and molecular dynamics (MD) simulations of the best candidates were also performed. Docking and MD simulation reveal that modified drugs are more potent than unmodified donepezil, where Trp86, Tyr337, Phe330 residues play some important roles in drug-receptor interactions. According to ensemble based docking, D9 shows greater binding affinity compared to the parent in most conformations obtained from protein data bank and MD simulation. In addition, it is observed that the π- π stacking with the residues of Trp86, Tyr337, Tyr341, Tyr124 and Trp286 may be required for strong ligand binding. Moreover, ADME/T analysis suggests that modified derivatives are less toxic and have improved pharmacokinetic properties than those of the parent drug. These results further confirm the ability of metal-directed drugs to bind simultaneously to the active sites of AChE and support them as potential candidates for the future treatment of Alzheimer’s disease.


Introduction
In recent decades, Acetylcholinesterase (AChE) has become a major interest in Alzheimer's disease (AD) research. AD, a neural degenerative disorder, is characterized by accumulation of extracellular and vascular amyloid in the brain [1][2][3]. In brief, inhibitors of AChE enhance the a1111111111 a1111111111 a1111111111 a1111111111 a1111111111

Molecular docking analysis
The three-dimensional crystal structure of recombinant human AChE (PDB ID: 4ey7) was retrieved in pdb format from the protein data bank [25]. The model was then subjected to energy minimization using the steepest descent and conjugate gradient technique to eliminate bad contacts of protein atoms. Computations were carried out in vacuo with the GROMOS 96 43B1 parameters set, with implementation using the Swiss-PDB Viewer. For docking analysis, AutodockVina was employed and AutoDock Tools (ADT) of the MGL software package was used to convert pdb into a pdbqt format to input protein and ligands. The size of grid box in AutoDockVina was kept at 58.81735, 61.2066, and 72.8273 respectively for X, Y, Z. Autodock-Vina was implemented through the shell script provided by AutoDockVina developers. The binding affinity of ligand was observed by kcal/mole as a unit for a negative score [26].

Molecular dynamics simulation
To validate the predictions from docking studies, MD simulation was performed using the NAMD [27] software, version 2.9. In this study, the CHARMM force field [28] was utilized, as it is widely applied to describe the macromolecular system. The Transferable Intermolecular Potential3 Points (TIP3P) water model was used by adding Cland/or Na + ions, where the total solvent molecules, 20109, have a density of 1.012 gm/cm 3 . A periodic boundary condition was employed to perform the simulation, where the box size used was 82.4×85.0×98.8Å 3 . Following the steepest descent energy minimization, equilibration of 100 steps was performed by NPT ensemble. Using Langevin Dynamics for constant temperature, full-system periodic electrostatics were maintained using the Particle Mesh Ewald (PME) [29]. Consistently Nose-Hoover Langevin piston [30,31] was used for constant pressure dynamics and SHAKE was used to keep all bonds involving hydrogen atoms at their equilibrium values. Finally, the full system was subjected to MD production run at 300 K for 25 ns in the NVT ensemble. The MD trajectories were saved every 50 ps for analysis.

Ensemble based molecular docking
To further clarify the results of docking predictions, we used an ensemble based docking method, where two different approaches were employed to obtain different conformations from AChE. In the first approach, different crystallographic conformations of AChE were retrieved from protein data bank, PDB IDs: 1b41, 1f8u, 1vzj, 2x8b, 3lii, 4bdt, 4ey6, 4ey8, 4moe, 4pqe, 5foq, 5fpq, 5hf5, 5hf6, 5hf8, 5hf9, 5hfa. In the second approach, conformers were taken from the 25 ns MD simulation (PDB ID: 4ey7) at every 1 ns of the 25 ns MD simulation. Against these conformers, the compounds donepezil, D8, D9 and D10 were subjected for docking using the same protocol discussed above in the methods section.

Pharmacokinetic parameters study
To check the pharmacokinetic parameters and toxicity of the modified compounds and parent compound, the admetSAR server was utilized. We have utilized the admetSAR online database to evaluate the pharmacokinetics parameters related to drug absorption, metabolism and toxicity of the parent drug and its designed analogues [32]. Using structure similarity search methods, admetSAR predicts the latest and most comprehensive manually curated data for diverse chemicals associated with known ADME/T profiles.
For ADMET analysis, the admetSAR program was used in which 96,000 unique compounds with 45 kinds of ADMET-associated properties, proteins, species, or organisms have been carefully curated from a large number of diverse literatures. Although it is quite difficult to verify all of these compounds and to know whether this program included metal-based drugs or not, we used well known Pt-based cisplatin and carboplatin as well as metal-based drugs approved in the FDA and in clinical trials as test candidates to verify our metal-based donepezil drugs.

Strategies and optimization of designed analogue
The new analogues of donepezil used in this study were designed according to the structural properties of the active site of AChE. As described above, among the two binding sites of AChE, the peripheral anionic site plays a significant role in ligand reorganization and allosteric activators [33,34]. The stabilization of the substrates binding on this site is largely π-cation interaction, while choline ester substrate specificity is mediated partly by Phe295 and Phe297 [35]. From detailed analysis of enzyme-inhibitor complexes, it appeared that the indole ring of Trp286 was involved in direct interaction with several inhibitors, showing a number of interaction modes including stacking, aromatic-aromatic, and π-cation, according to the nature of the ligands [36][37][38]. Furthermore, the active site of AChE forms electrostatic interactions with the substrates, as all of the amino acids were distributed with a large dipole moment. Information from the above studies, therefore, motivated us to design new analogues of donepezil, by increasing their electronegativity and the non-covalent interaction capacity between the aromatic rings.
As shown in Fig 1, 2 ]. There were also several additional modifications in D2-D10. D2-D5 were modified by the addition of F (D2), Cl (D3), Br (D4), and I (D5) atoms in the 2,3-dihydroindene ring portion, respectively. In contrast, D6 was designed by corresponding with D5 while modifications occurred only in the attached benzene ring, i.e., benzene ring with CF 3 group. In D7 and D8, attached benzene ring of the parent structure was replaced by naphthalene and anthracene rings, respectively, with no halogen modification; however, replacement of H with F and Cl atoms at the 2,3-dihydroindene ring portion of D8 results new analogues D9 and D10, respectively.
As the conformational features of a molecule critically influences its physical and chemical properties, all of designed compounds along with parent compound, donepezil, were subjected to full geometry optimization using DFT. Table 1 illustrates the stoichiometry, electronic energy, enthalpy, Gibbs free energy and dipole moment of the compounds and the optimized structures are depicted in Fig 2. According to the Table 1, it is clear that modifications on donepezil significantly influenced the structural properties of the compounds in terms of energy, partial charge distribution, and dipole moment. The highest energy, enthalpy, and Gibbs free energy was observed for D10, while D9 showed the highest dipole moment of 13.596 Debye, representing high polarity in nature. It is important to note that incorporation of the-CF 3 group in D5 significantly reduced the dipole moment, as can be seen in D6; however, D1, D7, D8 showed low dipole moments of 13.596 Debye due to the lack of halogens.

Analysis of frontier molecular orbitals
The frontier molecular orbitals are the most important orbitals in a molecule and they are considered to characterize the chemical reactivity and kinetic stability. These frontier molecular orbitals are known as the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO). Table 2 represents the values of orbital energies, along with the two global chemical descriptors, hardness and softness, which are also calculated for all compounds. The highest softness was observed for D8. D8 also showed the lowest HOMO--LUMO gap and hardness, indicating that the molecule is more reactive than other compounds, according to Pearson et al. [39,40]. In Fig 3, the HOMO plot of compound D9 showed that the electrons were localized on the upper part of the piperidine ring, while the LUMO plot showed that the electrons were localized at modified Cu regions only.

Molecular docking analysis
In order to check the binding modes of modified compounds, molecular docking simulations by Autodock Vina were performed. Molecular docking is one of the most common methods used in structure based drug design to analyze the interaction between a small molecule and a protein at the atomic level. Prior to docking, the crystal pose of donepezil was re-docked into the binding site of AChE with specific docking parameters and scoring functions, to check whether the docking software is reliable for the system. The conformation having the lowest negative score was then compared with the crystal pose. The value of the root mean square deviation (RMSD) of the docked conformation with respect to experimental conformation was 1.9659 Å (Fig 4), signifying the reliability of the docking protocol, as the threshold of reliability is 2.0 Å for a good docking protocol.
Afterward, all designed analogues were docked into the same binding site pocket of AChE, using similar optimized docking conditions. The outcomes of the docking analysis showed that all compounds, along with the parent compound, obtain binding affinities ranging from -10.2 to -14.9 kcal/mol. As shown in Table 3, D4, D5, D6 showed low binding affinities compared to parent compound, donepezil, while D1 exhibited high binding affinity. These results indicated that modification of Cu along with a water molecule increased the binding affinity, while addition of halogen groups like Br, I, and CF 3 made some fluctuations in binding affinities; however, modification with naphthalene and anthracene rings increased the binding affinity. As shown with D7 and D8, obtained docking affinities of -13.9 and -14.8 kcal/mol were determined, respectively. The highest binding affinity was observed for the D9 compound. According to the post docking analysis, it was revealed that all compounds, except D6, showed π-alkyl interactions with Tyr337 and Phe338 residues of the PAS in the active site of the enzyme. D6 is positioned to form stabilizing π-alkyl interactions with Trp286, Tyr337, Tyr341 residues. Furthermore, it was also observed that modifications of donepezil increased the π-π interactions with the residues of the active site, while increasing their polarity resulted in the formation of hydrogen bonding interactions. The highest H-bonds were obtained for the D10 compound, forming with Gln291, Ser293, Phe295, Arg296 residues. In contrast, D7, D8, and D9 formed three H-bonds with Tyr72 and Phe295 residues. D8 and D9 showed similar binding conformations, despite having different bonding distances. Along with Trp286, D8, D9, and D10 displayed the maximum π-π interactions with the Trp86 residue denoting the tight binding with the activesite. Reports suggest that Trp286 is considered as the principal component of the PAS, responsible for the accessibility of small molecules to the active site and also in the allosterism, while aromatic interactions with the Trp286 residue modulates the inhibition constants for some AChE inhibitors [41,42]. As the D9 compound showed the highest binding affinity (Fig 5), it was subjected for subsequent MD, along with donepezil, to investigate the dynamic stability of the AChE-inhibitor complex, and also, to ensure the rationality of the sampling strategy. Furthermore, to understand how D9 showed its binding modes with different metals, different metal atoms such as Fe, Co, Zn and Ni were inserted to the same position of D9 where Cu is present and they are renamed as D9-Fe, D9-Co, D9-Zn and D9-Ni, respectively (S1 Fig). These analogues were optimized by DFT and the subsequent molecular docking was performed by the same protocol discussed above in the methods section. Afterward, obtained results were represented in S1 Table. As shown in S1 Table, D9-Fe, D9-Co, D9-Zn, D9-Ni showed low binding affinities compared to D9. As per the post docking analysis, it is shown that D9-Fe, D9-Co, D9-Zn showed the π-alkyl interactions with Tyr337 and Phe338 residues of the PAS of the active site of the enzyme, like D9, respectively, while the Val294 residue only formed π-alkyl interaction with D9-Ni. In addition, it was revealed that all of the modified D9 compounds, except D9-Ni, showed maximum π-π interactions with Trp286 and Trp86. D9-Ni formed major π-π interaction with Trp286 along with the Tyr341 residue, which was also observed in D9-Co and D9-Zn. Furthermore, D9-Co and D9-Zn formed a hydrogen bond with the Phe295 residue while D9-Fe forms H-bonding with both Phe295 and Tyr72 residue, as like the D9 compound (illustrated in S2 Fig). From the different metal based study of D9, analysis finally revealed that D9-Cu performs better binding than other candidates.

Molecular dynamics simulations
In order to understand the binding mechanism, structural behavior, and flexibility of compound D9, we performed MD simulations for 25 ns. The complex of donepezil-protein was also subjected to MD simulation as a reference compound. The atomic RMSDs of the Cα atoms for both protein and the ligand of each complex were calculated and plotted in a time dependent manner (Fig 6). Fig 6A demonstrats the behavior of the protein during the Here X, X1, X2, X3 indicates that, X = Benzyl-4-piperidyl, X1 = 2,3-dihydro-1H-inden-1-one, X2 = Naphthalen-1-ylmethyl-4-piperidyl, X3 = Anthracen-9-ylmethyl-4-piperidyl https://doi.org/10.1371/journal.pone.0211935.t003 simulation; in which, both complexes were observed to achieve equilibrium after 5 ns and fluctuated around 0.5 Å. However, after 20 ns, D9 complex showed lower RMSD and remained afterward. Similar results were also obtained for ligands of each complex, as shown in Fig 6B. As can be seen in the plot, high fluctuation in RMSDs was observed for donepezil, where the high magnitudes were observed at 16 ns to 18 ns. The RMSD results from both protein and ligand indicated that the complexes were stable, suggesting higher stability of D9 in comparison with donepezil. For better understanding on how D9 and donepezil influence the binding mode with AChE, the structural changes of two complexes were examined by means of root mean square fluctuation (RMSF), radius of gyration, and solvent accessible surface area (SASA) of the protein (Fig 7). Fig 7A represents the total SASA of each protein, in which the D9 compound showed decreased SASA after 15 ns of simulation, demonstrating lower compactness of the protein structure. In contrast, the results from the radius of gyration analysis (Fig 7B) described that D9 comparatively produced a higher radius of gyration value than donepezil, denoting loose packing of the protein structure, which eventfully supported the results from the SASA analysis. RMSF values were also calculated from the trajectories, which reflect the flexibility of each residue in the protein [ [39]. According to Fig 7C, it was observed that D9 induced flexibility to some residues in the protein. Highest fluctuations were observed in several regions, ranging from 116-125, 280-290, 310-320, 361-370, and 505-515. Finally, the information of hydrogen bonding interactions formed within the protein, and also between the protein and ligand at the catalytic domain, was collected from the trajectories and represented in Fig  8A. Here, the D9 complex showed maximum intramolecular hydrogen bonds in the donepezil complex, demonstrating the stability of the complex. The intermolecular hydrogen bond analysis between the protein and ligand displayed that donepezil and D9 formed hydrogen bonds with the residues of the catalytic domain (Fig 8B). At the initial step, D9 did not show much interaction; however, after 11 ns, it showed several H-bond contacts. Consequently, donepezil revealed no hydrogen bonding in the docking pose, although it detected several contacts during the simulations. As a corollary, all analyses from the MD simulations suggested that D9 is more stable than donepezil and caused little conformational changes of the protein by undergoing little movement during the MD simulations (Fig 9).

Ensemble based docking
Usually, proteins are flexible macromolecules in nature. However, this property significantly influences ligand binding, especially in molecular reorganization and interactions [44]. Compared to the other program, AutoDock Vina is the most popular docking program to determine the binding pose of the ligand, yet suffers from backbone flexibility in receptors. Therefore, ensemble based molecular docking by AutoDock Vina has been introduced in this study to overcome this limitation. The results obtained are represented in Table 4 and  Table 5. Table 4 and Fig 10A describe the binding affinity of all ligands with different crystallographic conformations of the AChE enzyme. Interestingly, designed compounds showed better results than the standard drug, donepezil. Among these crystal conformers, designed compounds and donepezil produced best docking scores against the 5foq conformer (Fig  10A), and therefore, detailed molecular interactions of this conformer have been investigated and illustrated in Table 6. As can be seen in Table 6, the D9 compound formed two additional hydrophobic interactions with Tyr341 and Trp286 residues followed by π-alkyl and π-π stacked bonds, and also obtained the highest docking score. In case of D8 compound, additional hydrogen bonds were observed with Ser293 and Trp286 residues, while the polar interactions with Tyr72 and Phe295 were seen to disappear. Similarly, loss of hydrogen bonds was  also observed for Gln291, Phe295, and Arg296 residues. All ligands showed better binding affinities against the conformer obtained from MD simulation at 17 ns, as shown in Table 5 and Fig 10B. Detailed molecular interactions were illustrated in Table 6 which represents the breakdown of π-alkyl interactions of Tyr341, Tyr337, and Tyr124 residues. Instead, they formed π-π stacking with the ligands. Also, D9 and D10 compounds showed an additional salt bridge with the Asp74 residue followed by π-cation interactions. It is noteworthy to state that the flexibility of AChE is the major determinant of the binding affinity of the ligands, as evident from ensemble-based docking. π-π stacking with the residues of Trp86, Tyr337, Tyr341, Tyr124, and Trp286 showed a major contribution for strong drug binding and activity. Our study suggested that protein flexibility can give rise to differences in binding affinity and binding interaction of a drug with its target protein.

ADME/T analysis
In order to analyze whether the modified compounds produce any toxicity or altered pharmacokinetic profile, the admetSAR server was utilized. Different pharmacokinetic and pharmacodynamic parameters such as human intestinal absorption, [45] blood-brain barrier penetration, cytochrome P450 inhibition, [46] human ether-a-go-go-related genes inhibition, acute oral Toxicity, and rat acute toxicity were considered. The results are summarized in Table 7. As shown in Table 7, all compounds revealed a positive value (value above the prescribed threshold suggesting good permeability) with high probabilities, in case of blood-brain barrier and human intestinal absorption. Furthermore, modifications of donepezil resulted in a non-inhibitor of P-glycoprotein. The analysis displayed that D2, D4, D5, D6 and D9 were potential compounds of the human ethera-go-go-related gene. All compounds showed a similar oral toxicity profile, while D9 and D2 indicated the highest LD 50 value in rat acute toxicity, demonstrating non-toxic with respect to parent donepezil. In addition, ADME/T prediction of D9-Fe, D9-Co, D9-Zn, and D9-Ni was performed and compared with the D9-Cu analogue. D9 with different metals revealed the same values as D9-Cu. An exception, D9-Zn, showed a negative value in human oral absorption. These results have been summarized in S2 Table. A published review by Mjos et al. 2014 [47] discussed the importance of metal based drugs in the diagnosis of disease and enlisted a number of metalbased drug names which are already approved by the FDA and which have undergone clinical trials (shown in S3 Table). From S3 Table, we also predicted the pharmacokinetic parameters and toxicity of some drugs (data is shown in S4 Table).

Conclusion
In summary, the present study revealed some novel metal directed AChE inhibitors, developed by modifying a known inhibitor, donepezil. Modification with Cu, along with substitution using aromatic rings and halogens increased the dipole moment and π-π interaction capacity of the designed compounds. Furthermore, these modified compounds were more reactive than donepezil, as they showed lower HOMO-LUMO gaps. Molecular interaction analyses of docking simulations revealed similar binding conformations of all compounds at the active site and suggested D9 as a potent inhibitor, which can equally interact with both the CAS (Trp86) and PAS sites (Trp286) of AChE. The structural analysis with subsequent MD simulations demonstrated that D9 formed a stable conformation by creating hydrophobic and aromatic interactions with the active site residues such as Tyr337, Phe295, Tyr72, and Phe338. In addition, π-π stacking interactions with the residues of Trp86, Tyr337, Tyr341, Tyr124, and Trp286 may play a major role for strong drug binding and activity, according to ensemble based docking. Moreover, ADME/T analyses suggested that modified analogues were less toxic and have improved pharmacokinetic profiles than the parent drug. These results further confirmed the ability of Cu and other metal-directed analogues to bind simultaneously to the active sites of AChE and support them as potential candidates for the future treatment of Alzheimer's disease.
Supporting information S1