Molecular perception of interactions between bis(7)tacrine and cystamine-tacrine dimer with cholinesterases as the promising proposed agents for the treatment of Alzheimer’s disease

The infamous chronic neurodegenerative disease, Alzheimer’s, that starts with short-term memory loss and eventually leads to gradual bodily function decline which has been attributed to the deficiency in brain neurotransmitters, acetylcholine, and butylcholine. As a matter of fact, design of compounds that can inhibit cholinesterases activities (acetylcholinesterase and butylcholinesterase) has been introduced as an efficient method to treat Alzheimer’s. Among proposed compounds, bis(7)tacrine (B7T) is recognized as a noteworthy suppressor for Alzheimer’s disease. Recently a new analog of B7T, cystamine-tacrine dimer is offered as an agent to detain Alzheimer’s complications, even better than the parent compound. In this study, classical molecular dynamic simulations have been employed to take a closer look into the modes of interactions between the mentioned ligands and both cholinesterase enzymes. According to our obtained results, the structural differences in the target enzymes active sites result in different modes of interactions and inhibition potencies of the ligands against both enzymes. The obtained information can help to investigate those favorable fragments in the studied ligands skeletons that have raised the potency of the analog in comparison with the parent compound to design more potent multi target ligands to heal Alzheimer’s disease.


Introduction
In the twenty-first century, Alzheimer's disease has become a global concern in a way that extensive studies have been and are being done throughout the world to improve the disease treatment (Ballard et al., 2011;Ferri et al., 2005). Unfortunately, the percentage of people with Alzheimer's disease is rapidly increasing which makes high spiritual and material costs to be imposed on societies (Ferri et al., 2005;Fratiglioni, Ronchi, & Agüero-torres, 1999;Wancata, Musalek, Alexandrowicz, & Krautgartner, 2003;Wimo, Winblad, & Jönsson, 2010). Despite the considerable researches that have been performed by now, definitive factors that express Alzheimer's disease are not still well known (Dong, Duan, Hu, & Zhao, 2012;Jefferies et al., 2013). As a matter of fact, some hypotheses have been proposed and extensively studied including the cholinergic, amyloid cascade, and tau hypothesis (Tumiatti et al., 2010). The earlier approach to treat Alzheimer's disease was based on the cholinergic hypothesis (Bartus, Dean, Beer, & Lippa, 1982;Contestabile, 2011) which suggested that the reduced levels of acetylcholine and butylcholine (Tabet, 2006) in the synaptic space may lead to Alzheimer's disease. Currently available drugs are produced based on the mentioned hypothesis (Huang et al., 2012;Lawrence & Sahakian, 1998). Recently, some other factors such as oxidative stress and metal ions imbalance, particularly copper, have been suggested as significant factors that induce Alzheimer's pathogenesis (Anand, Gill, & Mahdi, 2014;Galimberti & Scarpini, 2012;Markesbery, 1997;Palmer, 2002).
Since several factors simultaneously involve the occurrence and progression of Alzheimer's disease, using one multipurpose ligand (multitarget-directed ligands (MTDLs)) in order to control several factors at the same time would definitely help a better treatment of the disease (Cavalli et al., 2008;Espinoza-Fonseca, 2006;Morphy & Rankovic, 2005;Youdim & Buccafusco, 2005).
Acetylcholinesterase (AChE) enzyme is not only responsible for acetylcholine break down in body but also plays an important role in the accumulation of amyloid-beta monomers, abnormal phosphorylation of tau protein, expression of APP95 protein, inflammatory processes, and cerebral blood flow modulation (Alvarez, Opazo, Alarcón, Garrido, & Inestrosa, 1997;Silman & Sussman, 2005;Soreq & Seidman, 2001). Butyrylcholinesterase (BChE) is another important member of cholinesterase family with an almost identical amino acid sequence (65%) to that of AChE (Campiani et al., 2005;Lane, Potkin, & Enz, 2006) and also a conformation and active site very similar to those of AChE (Suárez & Field, 2005;Vellom et al., 1993). Coincide with the disease progression, BChE levels in the brain increase as well. Although the rate of acetylcholine hydrolysis by BChE is low, its considerable impacts in reducing acetylcholine levels are non-venial (Mesulam, Guillozet, Shaw, & Quinn, 2002;. The active sites of the enzymes consist of three amino acids: serine, histidine, and glutamate that are called the catalytic triad (Nicolet, Lockridge, Masson, Fontecilla-Camps, & Nachon, 2003;Wiesner, Kříž, Kuča, Jun, & Koča, 2007). The catalytic sites (CASs) of both AChE and BChE are formed of narrow gorges located in the depth of 20 Å where the hydrolysis reaction takes place (Cokugras, 2003;Sussman et al., 1991). In AChE, acylbinding pocket is composed of Phe288 and Phe290 residues (residues numbers are from Torpedo californica AChE (TcAChE)) in the catalytic triad vicinity while in human BChE (huBChE), the equivalent residues are Leu286 and Val288. Other considerable amino acids that participate in AChE choline binding site (anionic site) are Trp84 and Phe330 that play a pivotal role in the hydrolysis process (Xu et al., 2008). In BChE, tryptophan amino acid is conserved (Trp82), but Phe330 is replaced by Ala328. These aromatic residues properly orient and stabilize the substrate in the gorge through cation-π interactions with acetylcholine quaternary ammonium group (Taylor & Radii, 1994). Two glycine and one alanine are constitutive amino acids in oxyanion hole near the choline binding site, but oxyanion hole in BChE consists of highly conserved N-H dipoles (Bajda et al., 2013). AChE peripheral anionic site (PAS) is composed of Tyr70, Asp72, Tyr121, Trp279, and Tyr334 while in that of BChE, the PAS site includes Asp70 and Tyr332 residues (Lane et al., 2006). Presence of less aromatic amino acids in BChE active pocket ordains a more spacious site for this enzyme (Taylor & Radii, 1994). Overall, these differences would dictate specific binding of the substrate and/or different ligands in the target enzymes CAS (Bajda et al., 2013).
Studies have shown that interactions of amyloid-beta monomers with the PAS site lead to a formation of an aggregated amyloid-beta form that is toxic for brain, while its monomeric form is relatively nontoxic. This process is called AChE-induced Aβ aggregation. Therefore, simultaneous inhibition of the active site and PAS with suitable dual inhibitors can prevent acetylcholine degradation and accumulation of beta amyloid monomers (Alvarez et al., 1997;Carvajal & Inestrosa, 2011).
Tacrine was the first approved drug by FDA to treat Alzheimer's disease in 1993, although it was long synthesized by Adrien Albert in 1940 (Summers, 2006). Since the introduction of tacrine, various tacrine hybrid derivatives have been designed and synthesized (Tumiatti et al., 2010). In 1996, bis(7)tacrine was introduced by Pang et al. as an effective MTDL (Pang, Quriam, Jelacic, Hong, & Brimijoin, 1996). Inhibition of AChE by bis (7)tacrine is 1000 times more than tacrine, moreover its special structure enables simultaneous inhibition of the active and PASs. Furthermore, bis(7)tacrine has beneficial interactions with the β-secretase enzyme and NMDA and GABA A receptors (Li et al., 2009;Minarini et al., 2012;Pang, 2001;Rydberg et al., 2006;Shu et al., 2012).
In 2012, Minarini et al. introduced a new compound called cystamine-tacrine dimer. Their experimental studies have shown that cystamine-tacrine dimer has lower toxicity than bis(7)tacrine. This compound is able to inhibit the cholinesterase enzymes, self-and AChE-induced Aβ aggregation similar to bis (7)tacrine. In addition to these characteristics, it can protect the neuroblastoma SH-SY5Y cell line against H 2 O 2 -induced oxidative injury by activating the extracellular signal-regulated kinase 1 and 2 (ERK1/2) and Akt/protein kinase B (PKB) pathways (Minarini et al., 2012).
The role of the molecular dynamic simulations (MD simulations) for a better understanding of biological processes is quite distinctive for any researcher in this field (Boehr, Nussinov, & Wright, 2009;Hansson, Oostenbrink, & van Gunsteren, 2002). MD simulations enable the explanation of the data from protein-ligand interactions that has been obtained experimentally, so that it helps to design more potent compounds. Additionally, the time and cost could be saved in the drug design field (Carlson & Mccammon, 2000;Durrant & McCammon, 2011;Karplus & Kuriyan, 2005).
This study is an attempt to perform a deep analysis on the interactions of both cholinesterases with tacrine as the oldest drug and bis(7)tacrine and cystamine-tacrine dimer as the promising proposed agents for the treatment of Alzheimer's disease with the aid MD simulations. Flexibility of the enzymes and stability of the ligands in the active sites have been thoroughly studied. The different interactions of both target cholinesterases with the target ligands have been demonstrated and the roles of the important residues in the active sites have been investigated. On the other hand, the MD simulations results have been used to interpret structure-affinity relationships. The obtained observations can increase our understanding of the key binding interactions which taking place at the active sites of AChE and BChE. On the other hand, these data can be applied to design better MTDLs.

Preparation of the initial configurations
The original crystallographic structures of TcAChE and hBuChE (human butyrylcholinesterase) were taken from the protein data bank with PDB ID of 2CKM and 4BDS, respectively (Nachon et al., 2013;Rydberg et al., 2006). Because AChE can also perform its biological task in a globular monomer form (Tripathi & Srivastava, 2008), the asymmetric unit of it was used in docking procedure and MD simulations.
All ligands, heteroatoms and other crystallographic agents were removed from the original proteins structures. The structures were then monitored for the missing residues and fixed with the aid of loop modeling server (Ko et al., 2011;Lee, Lee, Park, Coutsias, & Seok, 2010). The protonation state of the acidic and/or basic residues were decided by estimating the corresponding residues pKa values with the application of PROPKA web server and also considering their surrounding environment (Li, Robertson, & Jensen, 2005). Since the protonation states of the titratable residues have a great impact on the enzyme performance (Wiesner, Kříž, Kuča, Jun, & Koča, 2010), in order to have more accurate calculations and close to that of biological conditions, we have considered two protonated states for AChE to find the optimal mode according to a report published by Wiesner et al. (2010).
Lys, Arg, Asp, and Glu were considered charged for the first set of calculations (the state that is commonly applied for the mentioned residues in computational approaches). The protonated states of His amino acids were assigned considering their surrounding amino acids, so that the maximum number of hydrogen bonds between these residues and the neighboring ones were feasible. The second set of calculations was performed while Lys, Arg, Asp, and Glu were left charged except for Glu199 and Glu445 that were kept protonated according to their predicted pKa values. The protonation states of histidine amino acids were considered similar to those of the first set of calculations. Above-mentioned calculations were also performed for BChE. Accordingly, in the first set of calculations, all titratable residues were considered charged, and in the second set, all were again kept charged except for Glu197 and Glu443 that were considered protonated.
All ligands structures were optimized based on Hartree-Fock level of theory with 6-31G ** basis set and the corresponding atom charges were calculated with the aid of Gaussian 03 package (Frisch et al., 2004).

Docking procedure
The optimized structures of the proteins and ligands were used as input for the subsequent docking studies. The target ligands structures ( Figure 1) were prepared employing AutodockTools package (ADT; version 1.5.6) (Morris et al., 2010). Polar hydrogens were added to the proteins structures and partial atomic charges were assigned by Kollman-united charges method. The ligands optimized structures were applied when non-polar hydrogens had been merged to their structures by ADT. The atomic charges were assigned according to Gasteiger-Marsili charge, and flexible torsions in the ligands configurations were considered as well. Docking calculations were done by Autodock Vina (Trott & Olson, 2009) with default parameters except exhaustiveness value which was set to 100. The grid box of 42 × 42 × 42 Å (x, y, and z) was assigned on the enzymes binding pocket. In order to evaluate the docking procedure, the respective crystallographic ligands of the target enzymes (with PDB IDs 2CKM, 1ACJ, and 4BDS) were elicited and re-docked over the active site that testified the accuracy of the docking calculations.

Classical MD simulations
An appropriate protein-ligand complex structure with a productive orientation and also the lowest docking binding energy was selected as an input for the upcoming MD simulations. MD simulations were carried out using GROMACS 4.6.5 (Pronk et al., 2013) with Amber03 force field (Duan et al., 2003). The ligands topology files were prepared based on the published original paper for the force field. The ligands partial atomic charges were calculated at the B3LYP/cc-pVTZ//HF/6-31G ** level with a dielectric constant of four to mimic Amber03 force field condition and based on RESP method (Bayly, Cieplak, Cornell, & Kollman, 1993;Comell, Cieplak, Bayly, & Kollman, 1993). The topological descriptions of the ligands were built using acpype-antechamber (Sousa & Vranken, 2012) and the general amber force field (Wang, Wolf, Caldwell, Kollman, & Case, 2004).
The desired structure was then immersed in a periodic water box of cubic shape with a minimal distance of 14 Å from any edge of the box. The box was solvated using the TIP3P solvation model. The concentration of 150 mM of NaCl (physiological conditions) was added to neutralize the simulation system. Particle mesh Ewald algorithm was used to calculate the electrostatic interactions when all bonds were constrained using LINCS algorithm. A twin range cutoff was used for long-range interactions including 12 Å for Van der Waals and 10 Å for electrostatic interactions. The system was minimized with the steepest descent algorithm followed by a second energy minimization step with the use of conjugate gradient algorithm. The system went through NVT ensemble MD simulation at constant temperature of 300 K for 500 ps and then continued employing NPT ensemble in a periodic boundary condition, at a constant temperature of 300 K and constant pressure of 1 bar with the aid of V-rescale thermostat and Parrinello-Rahman barostat, with a coupling time of τ t = .1 ps and τ p = .5 ps, respectively. The conditions of position restraints for heavy atoms were applied in both stages. MD simulations were carried out, NPT for 20 ns with timesteps of 2 fs.
Graphical representation software Chimera 1.8 software (Pettersen et al., 2004) and Discovery studio v3.5 (Accelrys Software Inc.) (2013) were applied to observe and depict the interactions between the ligands and the active site residues of the enzymes. Grace was used to represent the plots of classical MD simulations.

Results and discussions
The probable modes of interactions between the ligands and the target enzymes were obtained with the aid of docking methodology. The docking procedure was evaluated by docking tacrine and bis (7)tacrine over AChE and also tacrine over BChE. A comparison between the obtained modes of interactions and ligands orientations in the corresponding binding pockets with those of the crystallographic evidence confirmed the accuracy of the docking calculations. For ligands bis(7)tacrine and cystamine-tacrine dimmer, π-π stacking electrostatic interactions are inevitable between the ligands proximal tacrine ring fragment and amino groups of Trp84 and Trp82 residues in AChE and BChE active pockets, respectively ( Figure 2). Six of the fourteen bulky aromatic residues at positions Tyr70, Tyr121, Trp279, Phe288, Phe290, and Phe330 in AChE are occupied by nonaromatic residues of Asn68, Gln119, Ala277, Val286, Val288, and Ala328 in BChE (Figure 3). Presence of more aromatic residues in AChE entrance tunnel enables the enzyme to establish extra hydrophobic and electrostatic interactions especially through peripheral anionic part.
The studied compounds docked structures, orientations, and binding modes in the target enzymes' active pocket as well as their calculated docking binding Figure 2. Interactions between enzymes and ligands after docking calculations (a) superimposition of AChE-bis(7)tacrine (green) and AChE-cystamine-tacrine dimer (orange) complexes, (b) superimposition of BChE-bis(7)tacrine (green) and BChE-cystamine-tacrine dimer (orange) complexes, and (c) BChE-tacrine complex (yellow). energies and experimental IC 50 values were thoroughly studied and compared. Tacrine best pose of interaction with BChE evidences π-π stacking electrostatic interactions with Trp82 residue, but in that of AChE, tacrine is docked in two distinct positions with close docking binding energies. In one of the modes, tacrine establishes π-π stacking electrostatic interactions with Tyr70, Tyr121, and Trp279 and also forms hydrogen bonds with Ile275 (Figure 4(a)), and in the other one, the ligand is stuck in a groove structured by Trp84, Phe330, Tyr334, and Trp432 through π-π stacking electrostatic and hydrogen bond formation with Asp72 and Tyr334 (Figure 4(b)). The docking binding energy of tacrine in the latter mode is slightly less than that of the former position.
Further, a comparison was done in between the experimentally obtained IC 50 values of the ligands and the calculated docking binding energies which are listed in Table 1. According to these results and the IC 50 values for tacrine in complex with AChE and BChE, it is inferred that due to the presence of two binding sites for tacrine in AChE active pocket structure (the choline and PAS sites) rather than the only binding site in that of BChE (the choline site), tacrine inhibits BChE specifically while it behaves selectively against AChE which probably results in higher consumption of the inhibitor to stop the enzyme activity.
Due to rigidity of the docking technique, classical MD simulations was applied in order to monitor the ligands behavior and stability in the active pocket while both proteins and the ligands have chances to get into their best conformations through which they can compose the best modes of interactions. The energy components, temperature, pressure, density, and volume of the simulation systems have been examined and assured to have appropriate stability throughout the 20 ns of MD simulations procedure in order to ensure the equilibration of the systems during the study.
The protonation states of residues Glu199 and Glu445 (with pKa values of 8.66 and 7.31, respectively) in AChE and residues Glu197 and Glu443 (with pKa values of 6.01 and 6.66, respectively) in BChE active pockets were ensured prior to the main simulations   operation. Therefore, the enzymes stability was examined by computing the enzymes behavior in different protonation states for 20 ns of MD simulations. According to the pKa values predicted by PROPKa server, the mentioned residues are expected to be uncharged in AChE pocket and charged in that of BChE. According to the obtained results of the performed MD simulations for AChE, the system reaches the steady state earlier with lower backbone root mean square deviations (RMSDs) when the target residues are considered uncharged in comparison with the other one in which the mentioned residues are considered charged (Supplementary data, Figure S1 (a)). Additionally, it was observed that the protonation state has a large impact on the root mean square fluctuations (RMSFs) of AChE binding site residues, especially Trp84, Ser200, Glu327, Tyr334, and His440 (Supplementary data, Figure S2 (a)). The protonated state of residues Glu199 and Glu445 decreases the CAS residues mobility especially His440 which can affect the overall mechanism of action. According to the report by Zhang et al., the presence of hydrogen bonds between residue His440 and two other residues, Ser200 and Glu327, plays a significant role in AChE catalytic activity (Zhang, Kua, & McCammon, 2002). To evaluate the desired hydrogen bonds stability, the distance between the hydrogen bonds donor and acceptor groups of the hydrogen bond participant residues (O γ of Ser200 and N ε of His440 imidazole ring, and O ε1 of Glu327, and the protonated N γ in His440 imidazole ring) was plotted during the simulations time (Supplementary data, Figure S3 (a) and (b)). A comparison between the obtained results of the two protonation states shows that the maintenance of the desired hydrogen bonds is only possible when Glu199 and Glu445 are protonated that is in agreement with the results reported by Wiesner et al. (2010). As a matter of fact, Glu199 and Glu445 in AChE structure was considered protonated for the upcoming simulations. All mentioned computations and analysis were also repeated for that of BChE enzyme that was contrary to the results obtained for AChE. In fact, according to the backbone RMSD plots of the two uncharged and charged states for BChE, it is concluded that the charged state of Glu197 and Glu443 residues retain the enzyme stability (Supplementary data, Figure S1 (b)). Further, the active site residues RMSFs decrease when Glu197 and Glu443 residues are considered charged (Supplementary data, Figure S2 (b)). The distances between the hydrogen bonds donor and acceptor groups of hydrogen bond participant residues (O γ of Ser198 and N ε of His438 imidazole ring, and O ε1 of Glu325, and the protonated N γ in His438 imidazole ring) were monitored as well which shows that the conservancy of the purpose hydrogen bonds are only possible in the charged state of the examined residues (Supplementary data (Figure S3 (c) and (d)). Therefore, the mentioned Glu residues were considered charged in BChE pocket during the subsequent MD simulations.
MD simulations were conducted for 20 ns for the best poses of the ligands-protein complexes. Backbone RMSDs of the enzymes structures were graphed against the simulations time to ensure the systems convergence into an equilibrated state. As can be seen in Figure 5(a), backbone RMSD plot of AChE in complex with bis(7)tacrine shows the lowest deviations in comparison with the two other ligands. For AChE in complex with cystamine-tacrine dimer, the system reaches the steady state after a gradual increase in the backbone RMSD up to .15 nm and then stays constant during the rest of the simulation. AChE backbone RMSD in complex with tacrine shows the most flexibility. Tacrine only fills the choline binding site of AChE, and as a result, the PAS site remains vacant. A comparison between the backbone RMSD plots of AChE in its free and the ones in the complex with the ligands forms shows that the ligands binding restrict AChE conformational changes to some extent. The backbone RMSD plots of the free BChE and those in complex with the desired ligands are also presented in Figure 5(c). In comparison with the observed behaviors for AChE in complex with the ligands, the difference in the backbone RMSD plots of BChE in complex with the ligands is not significant and even not very different to that of the free BChE enzyme.
To observe the ligand stability throughout the simulations time, the RMSDs of the ligands were also graphed ( Figure 5(b) and (d)). The significant deviations that are seen in tacrine RMSD plots represent rotation of the ligand in both enzymes choline binding sites. This phenomenon usually occurs for small ligands in commodious sites where ligands have enough space to bustle (Bagherzadeh et al., 2015). Tacrine rotates horizontally in the active site of AChE with 180°rotation which leads to the displacement of its aromatic and nonaromatic rings (Figure 6(a)). Considering Figure 7(c), this rotation facilitates the hydrogen bond formation between the NH 2 group of tacrine and oxygen atom of His440 (Supplementary data, Figure S4 (c)). Due to the presence of different residues in the choline binding site of BChE, tacrine rotates in an irregular vertical manner (Figure 6(b)) in a way that only the NH 2 group is displaced to aid a hydrogen bond formation with the oxygen atom of His438 while the position of the rings are almost constant. This different behavior can be attributed to the presence of two aromatic residues, Trp84 and Phe330 in the choline binding site of AChE, so that tacrine rings are trapped in between while the absence of Phe residue that is replaced by a non-aromatic residue, Ala328, in BChE pocket disturbs the formation of firm electrostatic interactions in a mode that is observed for AChE-tacrine complex. Instead, the presence of Met437, Trp430, and Ala328 hydrophobic residues facilitates appropriate implantation of tacrine within the choline binding site of BChE (Supplementary data, Figure S4 (f)) which is another reason for the specific binding of tacrine in BChE binding pocket and the considerable reduction in the IC 50 value (Nachon et al., 2013).
AChE and BChE active sites are well filled by bis (7)tacrine, and the modes of interactions between the ligand and the key residues are identical to those that were obtained with the use of docking procedure. The almost stable RMSD plots of the enzymes backbone and bis(7)tacrine also confirm this issue. In fact bis (7)tacrine is firmed in AChE active site through π-π stacking electrostatic interactions between tetrahydroanthracene fragments of the ligand and residues Tyr70, Trp84, Trp279, Phe330, and Trp432. Further, a water-mediated hydrogen bond is formed between Trp84 (nitrogen atom of indole fragment) and the NH 2 group of tacrine fragment of bis (7)tacrine during MD simulations (Figure 7(a)). Perceiving bis(7)tacrine in complex with BChE, fewer π-π stacking electrostatic interactions are observed for the ligand in the active site (between one of the tetrahydroanthracene fragments and Trp82) which is referred to the anatomy of the pocket that has been discussed previously. Due to larger solvent accessible surface area (SASA) for BChE in complex with bis(7)tacrine than that for AChE (246.562 and 234.514 nm S −2 N −1 , respectively) and lower number of aromatic residues in BChE binding pocket, more water molecules diffuse into BChE pocket that results in making further watermediated hydrogen bonds interactions like those that are observed between Asn68, Ala277, Ala328, and Gly439 and nitrogen atoms of the tetrahydroanthracene fragments, and the NH 2 groups (Figure 7(d)). SASA plots of the studied systems as a function of simulations time (Supplementary data, Figure S5) further confirm the noticeable distinction in the two enzymes SASAs. Also, the plots show that the enzyme-ligand complexes during MD simulations are relatively stable, indicating no considerable changes in the enzymes structure. Figure 5. RMSD plots of the enzymes backbone and the ligands (a) Backbone RMSD plots of AChE in complex with bis(7)tacrine (violet), cystamine-tacrine dimer (green), and tacrine (Orange) and the free (cyan); (b) RMSD plots of the ligands in complex with AChE, bis(7)tacrine (violet), cystamine-tacrine dimer (green), and tacrine (Orange); (c) Backbone RMSD plots of BChE in complex with bis(7)tacrine (magenta), cystamine-tacrine dimer (green), and tacrine (blue) and the free BChE (cyan); and (d) RMSD of the ligands in complex with BChE; bis(7)tacrine (magenta), cystamine-tacrine dimer (green), and tacrine (blue).  Modes of interaction of cystamine-tacrine dimer are different from those that were obtained by the docking calculations, especially for the ligand in complex with AChE. The dynamic nature of the classical MD simulations permits both the ligands and proteins to fit into their best conformations where they can interact best. These deviations from docking poses can be attributed to alternations in C 1 -S 1 -S 2 -C 2 dihedral angle (Figure 8) that has undergone severe changes in the minimization and MD production steps. In general, disulfides show a distinct preference for dihedral angles approaching 90°. As it is shown in Figure 8, C 1 -S 1 -S 2 -C 2 dihedral angle is 172°after docking calculation over AChE while this value changes to 124.43°and 86.47°after minimization and MD production steps, respectively. In addition, the hydrogen bond formation between one of the sulfur atoms of the bridging disulfide bond and residue Gln74 persuades partial deviance of cystamine-tacrine dimer from the choline binding site of AChE (Figure 7(b)). But because of Gln74 tendency to establish hydrogen bonds with Tyr34 and Asp72, this hydrogen bond is temporary exchanged between AChE residues and the ligand during MD simulations (Figure 9(a)).
The C 1 -S 1 -S 2 -C 2 dihedral angle of cystamine-tacrine dimer in complex with BChE is 115.76°after docking calculations and changes to 92.19°after minimizations step and then reaches 92.30°at the end of MD simulations time. These changes help to establish a stable and strong hydrogen bond between cystamine-tacrine dimmer sulfur atom and residue Thr120 (Figure 9(b)).
In order to understand the effects of different ligands on residues fluctuation throughout the simulations time, RMSF of the key residues have been computed and depicted in Figure 10. The key residues (the CAS and PAS sites) show small degree of fluctuations that indicate their appropriate bindings and interactions with the ligands in the active sites. The least fluctuations are observed for the residues in AChE-bis(7)tacrine and BChE-cystamine-tacrine dimer complexes.
Additionally, fluctuations of the omega loops (Ω loop) are shown in Supplementary data, Figure S6. The Ω loop is a disulfide link between two cystein residues located over AChE and BChE active sites, including Cys67 to Cys94 in AChE and Cys65 to Cys92 in BChE (Wiesner et al., 2007;Cokugras, 2003). AChE Ω loop reveals its lowest fluctuations in the presence of bis(7)tacrine (Supplementary data, Figure S6 (a)). The compatible interactions between the bis(7)tacrine and the indole ring of Trp84 and Tyr70 play a significant role in limiting the mobility of the Ω loop. The shifts in Gln74 orientation toward sulfur atoms of cystamine-tacrine dimer (that results in the notable fluctuation seen for the Figure 8. Superimposition of the changes in C 1 -S 1 -S 2 -C 2 dihedral angle of cystamine-tacrine dimer in complex with AChE; the docked configuration (green), configuration after the energy minimization step (magenta), and the configuration after 20 ns of simulations (cyan). residue) as well as the ligand C 1 -S 1 -S 2 -C 2 dihedral angle conversion dictate a conformational rearrangement in the region to establish an efficient interaction. Also, π-π stacking electrostatic interactions between the ligand, Trp279, and Trp432 in the adjacent loop and α-helice, respectively, help better stabilization of the ligand which amplifies its interactions with the Ω loop as well. The Ω loop residues that directly interact with the ligand and their neighbors incur significant changes in the fluctuations. Although presence of the ligands reduces the fluctuations of the most Ω loop residues, the fluctuations of Phe75, Pro76, and Gly77 have increased in complex with cystamine-tacrine dimer. The position of AChE Ω loop in complex with cystamine-tacrine dimer and bis(7)tacrine are compared in Figure 11(a). Gln74 tends to establish a hydrogen bond with sulfur atom of cystamine-tacrine dimer, thus this residue and its neighbor residues have shifted toward cystamine-tacrine dimer. Additionally, Asp72 residue has the highest fluctuations in the absence of the ligand (Figure 10(a)). The key roles of Asp72 are steering of cationic ligands to the entrance of the gorge and helping to maintain the ligands within the gorge (Hosea et al., 1996;Tara et al., 1998). Drastic reduction that is seen in the fluctuations of this residue in the presence of the neutral ligands is related to strength and number of hydrogen bonds between the residue and its surroundings residues (Gln74, Phe75, and Ser81).
RMSFs of the Ω loop residues of BChE are shown in Supplementary data Figure S6 (b). Fluctuations of residues Asn68 to Phe76 in the Ω loop are declined as a result of their direct and/or indirect interactions with the ligands. Like Asp72 in AChE, Asp70 and Tyr332 in BChE show the highest fluctuations in the absence of the ligands except that Tyr332 is not an omega loop participant residue. These residues are located in the PAS site of BChE and play a crucial role in the initial binding of charged substrates and also the enzyme activation control (Masson et al., 1999;Masson, Froment, Bartels, & Lockridge, 1996). The studies of Masson et al. have shown that although these residues are not directly involved in the binding of the neutral substrate but control the functional architecture of BChE active site gorge. As can be seen in Figure 10(b) and Supplementary data, Figure S6 (b), significant reduction in the fluctuations of these residues are observed in the presence of the ligands. Hydrogen bond formation between Tyr332 and the sulfur atom of cystamine-tacrine dimer has caused a sharp reduction in the fluctuations of this residue. Tyr332 also interacts with tacrine fragment through electrostatic interactions (Figure 11(b)).   (7)tacrine and cystamine-tacrine dimer, (a) AChE in complex with bis(7)tacrine (salmon) and cystamine-tacrine dimer (cyan) and (b) BChE in complex with bis(7)tacrine (green) and cystamine-tacrine dimer (orange). Figure 12. Superimposition of the average structures of (a) AChE in complex with bis(7)tacrine, (b) AChE in complex with cystamine-tacrine dimer, (c) AChE in complex with tacrin, over the free AChe, (d) BChE in complex with bis(7)tacrine, (e) BChE in complex with cystamine-tacrine dimer, and (f) BChE in complex with tacrin, over the free BChE. Lone enzymes are presented in yellow and the ones in complex in magneta for AChE and cyan for BChE. Further, in order to depict the extent of the active pocket residues geometry rearrangement that have been induced to the active pocket of the enzymes as a result of the ligands presence, an average structure was obtained from the equilibration state for each complex and was then superimposed over the corresponding free enzyme average structure (Figure 12). Considering the free AChE structure and those in complex with bis(7)tacrine and cystamine-tacrine dimer, it is observed that bis(7)tacrine implantation in the active site changes the orientations of the residues, and consequently, the general conformation of the choline and PAS sites more significantly than that of cystamine-tacrine dimer (Figure 12(a) and (b)). Considering tacrine, it is interesting that although there are no direct interactions between the ligand and the PAS site, significant alternations are seen in Trp279 and Tyr70 orientations (Figure 12(c)). According to the free BChE structure superimposition over those in complex with the ligands it can be seen that while the active pocket residues orientations do not need to go under significant rearrangement to interact with bis(7)tacrine, notable alternations are seen in the conformation of the active site residues in order to better interact with cystamine-tacrine dimer (Figures 12(d) and (e)). Further, as it is expected, tacrine does not induce remarkable variations to the overall structure of the active site (Figures 12(f)). These results better justify the reported IC 50 values for each of the ligands as well.
The reported crystallographic structure of AChE in complex with bis(7)tacrine was superimposed over the structure obtained after 20 ns of MD simulations study. According to Figure 13(a), a consistent agreement is observed between the active site residues and ligands orientations and also mode of interactions. On the other hand, superimposition of the crystallographic structure of tacrine in complex with BChE over the structure obtained after 20 ns of MD simulations confirms the ligands' firm binding in the active pocket. It is worth noting that the ligand oscillates between the two already mentioned orientations, so that the only observed difference in Figure 13(b) is the displacement of the NH 2 group of the ligand while there is no alternation in tacrine rings position. As already discussed, the rotation aids the hydrogen bond formation between the NH 2 group of tacrine and oxygen atom of His440. Further, BChE-tacrine complex structure is probably crystallographed in one of the probable orientations while the other one that is obtained through MD simulations is possible as well.

Conclusions
Inhibition of AChE and BChE enzymes is one of the main objectives to design promising drugs to treat Alzheimer's disease. Although the overall structures of the two enzymes are almost identical but the slight differences in their residues sequences ordain distinct ordinances for the same ligands in the active sites of the enzymes, and as a consequence, different modes of interactions are observed. This alternation is more serious in case of tacrine where the ligand has a choice to select either the PAS or the CAS sites in AChE active pocket to interact with, even preferably one more than the other. Meanwhile, the presence of only one hydrophobic binding pocket in BChE active site that consists of aromatic residues, the choline binding site, results in considerably lower the IC 50 value for BChEtacrine complex. Also, lack of conductive interactions Figure 13. Superimposition of (a) crystallographic structure of AChE-bis(7)tacrine complex (the protein in green and the ligand in orange) over the structure obtained after 20 ns of MD simulations (the protein in cyan and the ligand in yellow) and (b) crystallographic structure of BChE-tacrine complex (the protein in green and the ligand in orange) over the structure obtained after 20 ns of MD simulations (the protein in cyan and the ligand in yellow). between bis(7)tacrine and the PAS site of BChE increases the IC 50 value in comparison with the reported IC 50 value for this ligand in complex with AChE. Integration of bis(7)tacrine and cystamine as an antioxidant has introduced cystamine-tacrine dimer with a very good inhibition performance that simultaneously affects several pathogens of Alzheimer's disease. The characteristic of C-S-S-C dihedral angle induces different orientations of tacrine fragment in cystamine-tacrine dimer body within the choline binding sites of both enzymes along with substantial changes in the active site residues orientations. Sulfur atom of cystamine-tacrine dimer can establish hydrogen bonds with the active site residues of both enzymes. This hydrogen bond is stronger and even more stable in BChE-cystamine-tacrine dimer complex than that in AChE-cystamine-tacrine dimer complex. More stable hydrogen bonds and appropriate implantation of tacrine fragment in the choling binding site of BChE decrease the IC 50 value for cystamine-tacrine dimer in complex with BChE in comparison with the ligand in complex with AChE. Having a clear understanding of the interactions between these enzymes and the already proposed inhibitors with the help of MD simulations is an effective aid to clarify those structural features that induce favorite behaviors and desirable efficiency for a compound to inhibit/activate a specific activity that fastens the process of designing new more effective and less toxic ones in this context.