Structure-based screening and molecular dynamics simulation studies for the identification of potential acetylcholinesterase inhibitors

ABSTRACT Alzheimer’s disease (AD) is an irreversible, progressive neurodegenerative disease characterised by the appearance of amyloid β plaques and neurofibrillary tangles in the brain. The loss of cholinergic neurons is considered to be one of the contributing factors for cognitive and memory deficits in the disease. Donepezil, the most successful and prescribed drug to treat the symptoms of AD, is acetylcholinesterase (AChE) inhibitor that improves the memory and brain’s cognitive functions. In this study, the pharmacophoric features of donepezil were mapped from its cocrystal with AChE and were used to screen the Zinc15 database. The obtained hits were subjected to molecular properties and PAINS filter. The homology modelling and molecular dynamics (MD) were performed to prepare selected AChE protein model (PDB id 4EY7) for advanced studies. The virtual screening and precision docking led to the identification of five compounds. The ADMET property prediction and free energy calculations were carried out to obtain three final compounds. The Alanine scanning and MD study of the compounds viz. ZINC000013719534, ZINC000035551243 and ZINC000035596918 produced stable complexes of the ligands. The identified virtual leads have the potential for better AChE inhibition. Highlights Pharmacophore mapping of donepezil was developed from PDB id 4EY7. Pharmacophore-based screening yielded 601 compounds with RMSD < 0.5 Å. PAINS and Drug likeliness filtration resulted in 199 compounds for further virtual screening and precision docking studies. Free binding energy and ADMET studies provided three active compounds. The molecular dynamics study indicated that the compounds ZINC000013719534, ZINC000035551243 and ZINC000035596918 formed stable interactions with the protein.


Introduction
Alzheimer's disease (AD) is a multifactorial neurodegenerative disorder involving various contributing factors responsible for its progression [1]. An estimate suggests that about 47 million people suffered from dementia in the year of 2015 worldwide, and more than 152 million will suffer from one of the forms of dementia by 2050, if no new cure or treatment is in place [2]. The amyloid plaque and neurofibrillary tangle of hyperphosphorylated tau are the classical hallmarks of the disease. AD affects memory, learning and behaviour in a time-dependent manner leading to the vegetative state. Currently approved drugs for AD include acetylcholinesterase (AChE) inhibitors -Donepezil, Galantamine and Rivastigmine and N-methyl-Daspartate (NMDA) receptor antagonist -Memantine [3]. The other therapeutics which are in clinical development and has gained focus in recent years include Amyloid based therapies, Tau based therapies, secretase inhibitors and modulators (βand γ-secretase), reducing oxidative stress, metal chelation, antiinflammatory therapies and targeting Ca 2+ dysregulation [1,4,5].
AChE is a carboxylesterase, which is a member of the family of protein with α/β-hydrolase fold [6]. It is expressed in various forms, including globular forms (G1, G2 and G4) and asymmetric forms, i.e. attached collagen-like coiled structures (A4, A8, A12) [7]. The primary role of AChE is the hydrolysis of the synaptic acetylcholine (Ach) and regulation of the cholinergic neurotransmission in the body. Structurally, it consists of a 'large central mixed β-sheets' surrounded by '15 α-helices' [8]. The active site is located at the bottom of a narrow gorge of AChE consisting of the esteratic site (Ser203, Glu334 and His447) and anionic site (AS) (Trp86). The Trp86 residue forms cation-π interaction with quaternary nitrogen of the Ach along with Phe330. Another site referred as Peripheral anionic site (PAS) is 20 Å from the catalytic centre and consists of Tyr72, Asp74, Tyr124, Trp286 and Tyr341 residues. The aromatic residue's side chains create 40% surface of the gorge and are located on the loop providing greater conformational flexibility [9].
Peter Davies linked memory loss with a cholinergic deficit in AD patients in 1976. The hypothesis suggested that either decrease in the production of Ach or increased AchE activity may be the primary cause of AD [10]. The use of AChE inhibitors is a pivotal part of the current therapy for symptomatic treatment of AD and is quite successful. The other hypothesis is based on the misfolded protein Amyloid-β (Aβ) that gets deposited in the brain leading to apoptosis and neurotoxicity. Aβ forms various neurotoxic oligomeric species, may be due to the mutation in amyloidogenic pathway components, i.e. presenilin 1 and 2, APP and ApoE4 [1,11]. AChE is also found in non-cholinergic neurons and other tissues of the body, including hematopoietic and neoplastic cells [12]. It is suggested that AChE plays a crucial role in neuritogenesis, synaptogenesis, amyloidosis, dopamine neuron activation, regulation of apoptosis, nerve regeneration and hematopoiesis and lymphocyte activation [12,13].
The in vitro and in vivo studies supported the relationship between APP processing and cholinergic activation via muscarinic and nicotinic receptors [14,15]. Tacrine decreased the secretion of sAPPα in PC12, fibroblast, neuroblastoma and glial cells, and also reduced the level of secreted Aβ, Aβ  and Aβ  in neuroblastoma plasma cell lines, thus highlighting the significance of AChE inhibitors [16]. The studies suggested that Aβ peptide and AChE were localised in the AD brain [17]. The PAS of AChE acts as an adhesive site to amyloidogenic conformer of Aβ, leading to its conformational changes and formation of amyloid fibrils [18]. Thus designing the PAS binding AChE inhibitor(s) may prevent Aβ aggregation as well as enhance the cholinergic transmission. In our previous works, various in silico techniques, i.e. pharmacophore modelling, molecular docking, database mining, quantitative structure-activity relationship and de novo fragment growing strategies, were successfully implemented and led to the identification of potent AChE inhibitors [19][20][21][22].
The pharmacophore modelling, through ligand with cocrystallised protein structure, helps in the quick prediction of hits. The combination of structure-based pharmacophore modelling along with molecular dynamics (MD) study may help in better predictivity of results ( Figure 1). This study comprises of pharmacophoric approach along with other in silico tools to identify improved AChE inhibitors.

Material and methods
The computational tasks were performed on an Intel(R) Xeon (R) CPU E3-1225v5@ 3.30 GHz processor, RAM 32.0 GB system with Nvidia 'Quadro P600' GPU running on a Linux 64 operating system.

Pharmacophore hypothesis
The pharmacophore hypothesis was developed by using the cocrystallised donepezil molecule with AChE (PDB id-4EY7) [23]. The features, i.e. hydrogen bond acceptor, hydrogen bond donor, aromatic group and hydrophobic centre, were used to generate pharmacophore in Ligand Scout 4.1 [24,25].
The Zinc15 library, with 123,073,955 conformers of 13,147,339 compounds (accessed on 17-05-2018), was used to search virtual hits. The search was restricted to one conformer per molecule and one molecule per hit [27].

Molecular property filters
The PAINS (pan assay interference compounds) were identified using Knime (ver. 3.5.3) to remove unwanted compounds that may interfere in enzyme inhibition assays. The workflow developed by Saubern et al. was used in the study [28]. Several factors may have impact on bioavailability and blood-brain permeability of the molecules and were considered during the identification of hits. The molecular property filters including molecular weight <450, logP <5, H-bond donor <3, H-bond acceptor <7, number of rotatable bonds <8, total number of H-bonds <8, pKa (neutral or basic) 7.5-10.5 and polar surface area <70 Å were also used [29]. The filters were implemented using the Knime workflow to obtain suitable hits ( Figure S1 in SI) [30].

Homology modelling and model validation
The PDB (PDB id: 4EY7) was selected as a protein model for human AChE and was obtained from the protein data bank (https://www.rcsb.org/) [23]. The protein structure was carefully examined in Chimera-1.13.1rc to identify the missing sequences and corresponding structural elements [31]. The missing elements were appended using homology modelling, and the model was developed by SWISS-MODEL accessible via the ExPASy web server through user template mode [32]. The selected PDB was applied as a template, while the target sequence (access code: P22303) was obtained from UniProt (https://www.uniprot.org). Further, the quality of the homology model was evaluated using the PROCHECK [33], RAM-PAGE [34], Molprobity [35], GMQE, QMEAN and QMEANDisco [36]. The visualisation of the model was performed with Chimera.

Protein preparation
The PDB2PQR server (http://nbcr-222.ucsd.edu/pdb2pqr_2.1. 1) was employed to assign protonation states of amino acid residues of the developed homology model at pH 7.4 using the AMBER force field. It also optimised the hydrogen bonding network of the protein. The PDB was subjected to a 5 ns MD simulation run in order to prepare protein for molecular docking in AMBER2018. Amber ff14SB force field was used to prepare topology of protein, and the AM1-BCC charge scheme was used to assign the atomic charges and force field parameters for donepezil in Antechamber tool. The final coordinate and topology of the complex were built by using tleap. The complex was hydrated with TIP3PBOX water using a box size of 89 × 91 × 103 under periodic boundary conditions (PBC) with 8 Å non-bonded cut-off and Particle Mesh Ewald (PME) for the long-range electrostatics. The complex was minimised by an exhaustive protocol including convergence criterion that energy minimisation would halt when rms reached lower than dmrs, i.e. 1 × 10 −4 kcal/mol-Å (Table T1 in SI). It was followed by heating of the system for 100 ps to raise the temperature to 310.15 K as the isothermal-isobaric (NPT) ensemble using Langevin dynamics. Further, the density equilibration and final equilibration were carried out for 100 ps each as isothermal-isobaric (NPT) ensemble at a temperature of 310.15 K. The final production run of 5 ns was carried out as isothermal-isobaric (NPT) ensemble at 310.15 K and 1 atm [37,38]. The last frame of the trajectory of the MD run was converted to PDB and further into pdbqt.

Ligand preparation
The identified hits were subjected to energy minimisation in Open Babel 2.4.0. The ligand was minimised at physiological pH using General Amber Force Field (GAFF) with 2500 as maximum steps [39].

Grid generation and validation
Protein-Ligand Interaction Profiler (PLIP) was used to identify the residues involved in non-covalent interactions with donepezil in the selected PDB file. The identified residues were taken as reference points to build a grid box around the active site. The autogrid 4.0 was used to calculate grid maps of interaction energies with various atom types present in the ligands (A, C, HD, NA, N, OA, S, Br, Cl and I). The grid size was set to 74 × 70 × 82 with a grid point spacing of 0.375 Å. Further, the developed grid was validated by redocking donepezil with randomised coordinates. The redocked pose was evaluated by determining RMSD corresponding to the heavy atoms with the co-crystallised structure of donepezil.

Virtual screening and molecular docking
Virtual screening was performed through docking on Auto-Dock 4.2, which was followed by precision docking of selected compounds. The Lamarckian Genetic Algorithm (LGA) was employed for conformational search in both the studies. The LGA parameters used for virtual screening and precision dockings are summarised in Table T3 (SI) [40,41]. The virtual screening results were processed through vstools_v0.16, a python script. Post-docking analysis and visualisation were performed by Autodock tools 1.5.6 and Discovery Studio visualiser.

MM-GBSA/MM-PBSA assay
The MM-GBSA and MM-PBSA of the selected compounds were performed using python script, i.e. MMPBSA.py of Amber18 [42]. The topology and parameter files of a protein-ligand complex, protein and specific ligands were built through tleap in the gas phase. Further, the proteinligand complexes were hydrated as a cubic box using TIP3P water molecules with a cutoff distance of 12 Å between any atom initially present in solute and the edge of the periodic box. The hydrated system was subjected to energy minimisation, heating, equilibration followed by MD run of 1 ns. The parameters for the processes are described in Tables T4  and T5 (SI). The obtained MD simulations were subjected to free energy calculations. The first 50 frames were considered for calculations with a salt concentration of 0.1 mM. The GB model developed by Onufriev et al. was used for MM-GBSA calculations [43]. Further, the total non-polar solvation free energy in MM-PBSA was modelled as linearly proportional to the solvent accessible surface area (SASA) in a single term. The entropy calculation was neglected in both cases [44].

Alanine scanning
The alanine scanning for the compounds selected on the basis of free energy calculations was carried. The residues interacting with the ligands were mutated, and their contributions to the binding energy of ligands were calculated using the same simulation run as used for free energy calculations. The individual residues were mutated to alanine followed by building the topology and coordinate files of the respective mutated PDBs by tleap. The MMGBS.py, a python script, was used to perform virtual alanine scanning for the first 50 frames of each simulation at a salt concentration of 0.1 mM for both MM-GBSA and MM-PBSA calculations [45].

Molecular dynamics
MD simulations were carried out to evaluate chemistry at the molecular level. MD was performed for ACHE, AChE-ZINC000035596918, AChE-ZINC000035551243 and AChE-ZINC000013179534 complexes through Amber18. The protein-ligand complexes were soaked in TIP3P hydrated cubic box with a cutoff distance of 12 Å. The residual charges of the system were neutralised by using counter ions placed in the system through the Monte Carlo method. The PME was used to handle long-range electrostatic interactions. The MD systems were subjected to energy minimisation followed by heating, density equilibration and equilibration under periodic boundary conditions. The detailed experimental protocols are summarised in Tables T1 and T6 (SI). The final production phase of 50 ns was carried out at 310.15 K as an NPT ensemble. The graphical representation of the 3D models was developed using VMD (Visual Molecular Dynamics) [44].

Pharmacophore hypothesis
Donepezil is a reversible AChE inhibitor producing mixed as well as a non-competitive mode of inhibition with an IC 50 of 5.7 nM [46]. It bears indanone and N-benzylpiperidine moieties with a selective action on AChE. Further, it was observed that 5,6-dimethoxylindanone moiety helped in enhancing the potency [46]. The disease-modifying effect of donepezil was observed in terms of neuroprotection against various insults caused by exogenous glutamate and Aβ. It attains peak plasma level within 3-4 h of administration, and most of the drug is excreted unchanged in the urine with a few drug-drug interactions. Thus it is a suitable scaffold for the identification of new AChE inhibitors [47]. Ligand scout identified proteinligand interactions of donepezil with AChE derived pharmacophore features, which was used as a filter to find initial hits. We obtained six pharmacophore features, i.e. two aromatic rings (A) and four hydrogen bond acceptors (HA) (Figure 2).

Database screening using pharmit
The obtained pharmacophore was queried against Zinc15 database, through Pharmit server, to identify the virtual hits. Six hundred and one molecules with diverse structures containing various heterocyclic rings including thiazoles, benzothiazoles, piperazines, pyrimidines, pyridines, quinazolines, pyrazoles, piperidines, etc., were obtained.

Molecular property filters
It is a common practice to pre-filter extensive databases, based on suitable molecular properties, in order to reduce the computational efforts and eliminate unnecessary compounds. PAINS may display difficulties that might produce assay artefacts, false-positive results and sometimes cause off-target effects that interfere with testing methods including aggregation of compounds in the assay medium, auto-fluorescence, cysteine oxidation, membrane disruption, chelation and reactivity towards protein. Thus there is a high risk of failure associated with such compounds for clinical investigation due to the misleading results [48]. PAINS compound filtration was performed by using SMARTS pattern search and 571 compounds showed the absence of PAINS SMARTS patterns.
The blood-brain barrier (BBB) has tight junctions between endothelial cells and drugs have to cross multiple membranes through transmembranous diffusion. It is observed that the CNS active drugs have physicochemical property values in short ranges in comparison to others. They display lower molecular weight, less flexibility and polarity with high values of cLogP and smaller molecular volume [29]. Therefore, the compounds were further subjected to physicochemical property filter to obtain 199 compounds.

Homology modelling and model validation
The careful inspection of human AChE PDB (PDB id-4EY7) indicated the absence of residues 259-264 and 495-497. Therefore, homology model was developed to fill the missing secondary structures using SWISS-MODEL. Ramachandran plot was obtained from two web servers: PROCHECK (https:// servicesn.mbi.ucla.edu/ PROCHECK) and RAMPAGE (http://mordred.bioc.cam.ac.uk/~rapper/rampage.php) ( Figure  3). Quality of the model was appropriate with more than 90% of residues in the most favoured region ( Table 1).
The developed model displayed a MolProbity score of 0.93 indicating the acceptable quality of the model, displayed only a single steric clash between Arg18 and Asp61 [49].
GMQE Score for the model was 0.99 indicating the closeness of the model and template. Further, the model had QMEAN score of 0.88 and favourable torsion angle potential. The Z score of the model had a standard deviation between 0 and 1, was an indicator of its resemblence to experimentally obtained protein structures in the real world ( Figure S2 in SI). Further, QMEANDisCo indicated all the AChE active site and tunnel residues had good local quality scores ( Figure S3 in SI). Hence, most of the validation parameters indicated confidence over the developed model of human AChE.

Protein preparation
The PDB2PQR server was used to prepare protein structure at physiological pH 7.4. It helped in adding any of the missing atoms, missing nomenclature, assigning atomic radii and charges to the residues according to the Amber force field [50]. The PROPKA utility was used to assign the correct protonation state to amino acid residues. The charges on the amino acid residues at particular pH are essential for determining its binding affinity with the ligands and the structural organisation of protein.
Proper inclusion of flexibility during protein structure preparation is essential for the success of virtual screening and docking studies. Thus energy minimisation, heating at 310.15 K, followed by NVT and NPT equilibration and a mini MD run of 5 ns was performed to enable the side chain of residues to be relaxed and acquire suitable conformation at physiological conditions provided. The production end frame was converted to PDB and was used to perform virtual screening and molecular docking. MD run of 5 ns displayed a steady trajectory with a mean RMSD deviation of 1.412 ( Figure S4 in SI).
Further, the prepared protein was converted to pdbqt to utilise it as an input file by AutoDock 4.2. The AutoDock atom types were assigned by merging non-polar hydrogen, since it used United Atom Force-Field (UAFF) of AMBER that had polar hydrogens only. Further, the merging of nonpolar hydrogens with heavy atoms helped in the reduction of computational cost [51].

Ligand preparation
The selected compounds were subjected to energy minimisation in order to rectify the geometry, remove the steric clashes and overlapping atoms of the compounds. The energy minimisation was performed on Open babel software. The total potential energy of the system depends upon the energy contributions due to bonded and non-bonded interactions. The local energy minimum was identified using the steepest descent algorithm.

Grid generation and validation
PLIP server was used for the identification of essential interactions of donepezil with protein residues present in the selected PDB. The interactions with a distance of less than 5 Å were selected for the grid generation. Donepezil displayed hydrophobic interactions with Trp86, Tyr337, Phe338 and Tyr341 along with a hydrogen bonding with Phe295. It also showed π-π interactions with Trp86, Trp236 and Tyr341 (Table T2 in SI).
The autogrid4 was used to generate pre-calculated atomic affinity potentials for each type of atom in order to improve the computational efficiency of the docking procedures. Further, the selected grid box dimensions and grid point spacing criteria were validated. The accuracy of docking protocol can be determined by the ability to efficiently reproduce the pose similar to the co-crystallised ligand conformation with least deviation. The selected grid parameters displayed acceptable average RMSD deviations of 0.967 and 1.465 Å between redocked donepezil pose with homology model and MD derived donepezil poses respectively, over 28 pairs of heavy atoms (Figure 4). The average RMSD deviation in the two cases was quite low and within the acceptable limit of less than 2 Å. Thus it was inferred that the generated grid would produce reproducible and reliable results, which could be considered for virtual screening and precision docking [52].

Virtual screening and molecular docking
The docking-based virtual screening is one of the typical methods used for the discovery of new inhibitors and has resulted in the successful identification of potent compounds that may have resemblance with the co-crystallised structure or as novel leads [53][54][55]. In this study, LGA was employed and ligands were treated as flexible entities [56]. The studies were performed using AutoDock 4.2, which used a semiempirical free energy force field in its scoring function. The scoring function constituted contributions of various energies including van der Waals, electrostatic, hydrogen bond, desolvation and torsional penalty involved in protein-ligand binding [51].
The scoring function is a fast and reliable method of sorting out and reducing the number of compounds from a large database.
During the virtual screening of 198 ligands, we obtained 21 compounds with binding energy ranging from −11.94 to −7.69 kcal/mol. The virtual screening results displayed a near normal distribution in terms of binding energies of compounds with a mean binding energy of −10.41 kcal/mol. The compounds displaying a cluster size of 6 and above were only considered for precision docking ( Figure 5). The precision docking was carried out by using extended LGA parameters with a population size of 100. The cluster size of 50 and above with an RMSD tolerance limit of 2 Å was used to filter the compounds in order to ascertain the binding site of each ligand. Five ligands with a cluster size greater than 50 were obtained and selected for further study. The order of precedence was found to be ZINC000169753041, ZINC000013719534, ZINC000035551243, ZINC000035596918 and ZINC000014996252, in terms of better binding energy in precision docking. The detailed docking results are reported in Table  2. Compound ZINC000013719534 displayed two hydrogen bond interactions with Phe296 and Arg295. The naphtho[2,3-d] [1,3]dioxol-5 (8H)-one ring of the ligand interacted with Trp286 and Tyr341 through π-π interactions, while Tyr124 interacted with terminal tolyl ring. The N-ethyl group displayed π-alkyl interactions with Trp286 and Tyr341. Further, piperidine ring interacted with Phe338 and Tyr341 and naphtho[2,3d] [1,3]dioxol-5(8H)-one ring with Leu289. Gly122 and Gln291 displayed weak van der Waals interactions. The precision docking of ZINC000169753041 with AChE displayed hydrogen bonding interaction with one of the PAS residue, i.e. Tyr72. The methoxy group as well as pyrazole hydrogen also formed hydrogen bond interactions with Thr75 and Arg295. The PAS residues Trp286 and Tyr341 interacted with the two divergent phenyl rings of the molecule and also with pyrazole ring. Tyr124, Ser293 and Tyr341 formed van der Waals interactions with the ligand.

ADMET property
ADMET profile of a compound is an important aspect for its clinical efficacy. PreADMET server provides a user-friendly web interface to predict various parameters [57]. The absorption of the drug is a limiting factor for its bioavailability and response. The predicted Human intestinal absorption (HIA) for all the selected compounds was greater than 95%. ZINC000035596918 displayed maximum HIA of 99.26%. The BBB permeability is another critical criteria to identify CNS active molecules. Ma et al. developed QSAR-based BBB permeability model and observed that the permeability depended on polar surface area, octanol/water partition coefficient and Balaban Index of the molecules [58]. The molecules with BB (C brain /C blood ) range from 0.1 to 2.0 indicate moderate BBB penetration. The identified compounds displayed BB (C brain / C blood ) ranging 0.2-0.6, i.e. moderate BBB penetration. The compounds were found to be an inhibitor of CYP2C19, and CYP3A4 and non-substrate of CYP2D6.
Further, moderate plasma protein binding was observed in the case of all the selected molecules. The rodent carcinogenicity model based on data of the National Toxicology Program and US FDA of mice and rats for two years was used to predict the carcinogenicity. The selected compounds were found to be non-carcinogenic in rat and mouse (Table 3).

MM-GBSA/MM-PBSA assay
Free energy calculations are a group of methods that are used to predict ligand binding energy by considering various atomiclevel interactions responsible for the affinity of ligand toward protein [59].

Protein + inhibitor
DG binding of complex  (Tables  4 and 5). The results of precision docking, however, showed deviation from these results as compound ZINC000169753041 was found to have better binding energy than others, which was not observed in case of both free energy calculations. The reason for the observation might be noninclusion of the water molecules in the docking, hence, making free energy calculations a superior method for the calculation of binding energies than the scoring functions used to evaluate docking poses.

Alanine scanning
Alanine scanning is a crucial method to determine the contribution of specific residues in binding with the ligand. In this study, the virtual alanine scanning for three compounds, i.e. ZINC000013719534, ZINC000035551243 and ZINC000035596918, was performed. The alanine scanning for ZINC000013719534 displayed that Trp286 and Tyr341 played a pivotal role in its binding due to higher losses in the binding energy of −4.78 and −4.49 in case of MM-GBSA and −3.82 and −2.46 in case of MM-PBSA respectively. The docking pose also reflected similar results, as both residues displayed π-π and π-alkyl interactions with the ligand. The contribution from Tyr124 (π-π) and Phe338 (π-alkyl) was the next significant contribution. Further, the contributions due to hydrogen bonding of Arg296 ranked third significant contribution towards free binding energies. The contribution made by Leu289 through π-alkyl interaction was also important. However, the contributions of Phe295 and Gln291 in the form of hydrogen bonding and carbon-hydrogen bonding were found to be insignificant. This may be because these residues were part of flexible loops and had higher positional fluctuation during MD. Virtual alanine study of ZINC000035551243 indicated that the interactions with Ser125 and Phe338 may not have a significant impact on ligand binding in both MM-GBSA and MM-PBSA. Further, Phe297 having π-alkyl interaction displayed a significant loss of binding energy in case of MM-GBSA but not in MM-PBSA. The two crucial interactions with AS and PAS residues were with Trp86 and Tyr341, respectively. Further, docking study also displayed that the binding of the ligand was contributed by van der Waals interactions with Gln71, Tyr72, Asp74, Thr75, Asn87, Pro88, Gly121, Gly122, Tyr124, Gly126, Trp286, Ser293, Val294, Phe295, Arg296, Try337 and Gly342. The hydrogen bond interaction of Tyr124 and ZINC000035596918 had significant contribution in ligand binding with a loss of −3.06 and −1.56 kcal/mol in MM-GBSA and MM-PBSA respectively. Tyr337, Phe338 and Tyr341 displayed significant binding contribution, which was reflected in alanine scanning. Val294 and Phe295 displayed significant contribution in MM-GBSA, and no significant contribution was made by Glu292. The contributions of Gly121 and Gly122 could not be studied due to the limitations of Amber18 (Figure 7).

Molecular dynamics
MD is a state-of-the-art technique to study the motion and trajectory of the molecules in the presence of other molecules along with the various interactions within a system. It helps to learn the conformational changes in the molecules, solvation of the molecules, structural features of the protein and drugprotein interactions. In this study, three protein-ligand complexes of ZINC000013719534, ZINC000035551243 and ZINC000035596918 were studied through MD simulations. A detailed analysis of the trajectories was carried out by using cpptraj [60].

Pre-MD phase analysis
Three selected protein-ligand complexes were prepared for MD through various stages including system building, energy minimisation, temperature equilibration, density equilibration and a small NPT ensemble before MD run. The prepared systems were subjected to energy minimisation. The energy states of protein-ligand complexes of ZINC000013719534, ZINC000035551243 and ZINC000035596918 were found to be −193540, −193390 and −193610 to −289910, −289690 and −289840 kcal/mol respectively, which indicated that the energy minimisation protocol had substantially minimised the complexes. It was further subjected to heating, which indicated that the temperature of 310.15 K was attained at 11, 7 and 6 ps and density of system attained to 1 g/cm 3 at 40, 41 and 43 ps for protein-ligand complexes of ZINC000013719534, ZINC000035551243 and ZINC000035596918 respectively. A short MD simulation run of 200 ps was used to check the stability and accuracy of the MD run. The protein-ligand complexes of ZINC000013719534, ZINC000035551243 and ZINC000035596918 displayed mean temperatures of 310.12, 309.94 and 310.09 K with the average total potential energy of −213775, −213775 and −213773 kcal/mol and mean RMSD of 0.788, 0.706 and 0.777 Ǻ respectively. These simulations attained a stable RMSD at 96, 100 and 132 ps for complexes of ZINC000013719534, ZINC000035551243 and ZINC000035596918 respectively ( Figures S5, S6 and S7 in SI).

RMSD analysis
RMSD is a measure of the average displacement changes in an atom's position for a particular frame with respect to a reference frame. Smaller the deviation, better is the stability of the  than standard donepezil-protein complex, indicating the stability of protein-ligand complexes (Figure 8).
The docking study of protein-ZINC000035551243 complex displayed hydrogen bonding with Trp86, i.e. AS which was also indicated in RMSF plot through low fluctuation in comparison to AChE. The other carbon-based non-conventional hydrogen bond formed by Ser125 residue also indicated similar RMSF pattern in that region. The regions of 283-298 and 331-342 amino acid residues displayed π-π and π-alkyl interactions with lower RMSF. ZINC000035596918 interacted with Gly121, Gly122 and Tyr124 through hydrogen bond, leading to their subsequent stabilisation as indicated by the low RMSF fluctuations. The residues Glu292 and Phe295 formed hydrogen bonds, Val294, Tyr337, Phe338 and Tyr341 formed π-π and π-alkyl interaction with ligand resulting in the stabilisation of protein in these regions. Further, the standard donepezil-AChE complex has displayed higher fluctuation than that of the selected ligands, which further highlighted the fact that the three selected ligands outperformed donepezil in terms of stability.
The ligand RMSF analysis also displayed some interesting observations. The fused benzo[d] [1,3] dioxolyl moiety had maximum RMSF deviation in ZINC000013719534, while in case of ZINC000035551243 and ZINC000035596918 parafluoro phenyl and benzo[d] [1,3]dioxolyl moieties respectively were having maximum RMSF. The indene moiety of donepezil also displayed fluctuation. The ligand RMSF monitoring indicated that the solvent-exposed regions of all ligands were less stable ( Figure 9). for AChE, ZINC000013719534, ZINC000035551243 ZINC000035596918 and donepezil complexes of AChE respectively. It was observed that there was a decrease in SASA values after ligand binding in all the cases, which were highest for donepezil followed by ZINC000035551243. This indicated that the selected ligands were binding in the tunnel around PAS, resulting in decreased SASA in contrast to donepezil ( Figure 10).     the crucial residues were contributing protein-ligand stability through hydrogen bonding ( Figure 11).

Secondary structure changes on ligand binding
The secondary structure analysis measures the changes induced by ligand binding in protein structure with time [61]. The binding of ZINC000013719534 to AChE led to no significant change in the secondary structure of the protein and thus have quite stable RMSD values. The binding of ZINC000035551243 and ZINC000035596918 led to increase in α-helix and a decrease in turns. The bends and 3-10 helix were also decreased in the case of ZINC000035596918 (Figure 12).

Conclusion
Computational drug design has played a crucial role in the drug design and understanding of protein functions. MD methods, though computationally expensive, have significantly reduced the cost of drug discovery to a great extent. AD is multifaceted disorder characterised by loss of cholinergic neurons and release of AChE monomers. The AChE monomers act as nucleating agent and aggravate the formation of Aβ plaques.
The released AChE increases the metabolic turnover of Ach and produces a cholinergic deficit in frontal cortex and hippocampus, leading to memory loss and cognitive impairment. Therefore, targeting AChE provides a multifactorial advantage in AD. The structure-based pharmacophore model was developed from donepezil, and pharmacophore-based screening of ZINC 15 database of 12, 125, 125 compounds was performed. The data-mining led to the identification of 601 compounds, which were further subjected to subsequent drug-likeliness and PAINS filters by removing unnecessary compounds. The crystal structure of AChE complexed with donepezil was obtained from PDB (PDB id. 4EY7). The missing residues were modelled through homology modelling and the model displayed GMQE, and QMEAN Z-scores of 0.99 and 0.88 respectively. The modelled protein still had bad contacts and incorrect side chain(s) positions, so energy minimisation and MD of 5 ns were performed to obtain a relatively stable protein conformation before docking studies. Through virtual screening and precision docking, 21 and 5 ligands were obtained respectively. The in silico pharmacokinetics and toxicity parameters for all the compounds were determined in order to ascertain molecules with potential clinical candidature. Free binding energy calculations were performed on five molecules; three of them displayed better binding free energy than donepezil. The MD study displayed that the stable protein backbone RMSD of AChE by three ligands, i.e. ZINC000013719534, ZINC000035551243 and ZINC000035596918 in comparison to donepezil. ZINC000013719534 displayed stable ligand RMSD and was most stable among the three compounds. The residue position fluctuation results were also in line with RMSD study. ZINC000013719534 and ZINC000035551243 displayed lower protein and ligand RMSF due to greater hydrogen bonding contacts with PAS for longer time duration than donepezil. Thus these compounds are proposed as the virtual leads and are expected to display predicted behaviour in experimental studies.