A Steered Molecular Dynamics Simulation - American Chemical Society

Biology, Weizmann Institute of Science, 76100 RehoVot, Israel, and School of Pharmacy,. East China UniVersity of Science and Technology, Shanghai 2002...
2 downloads 0 Views 805KB Size
Download PDF
23730

J. Phys. Chem. B 2005, 109, 23730-23738

Dynamic Mechanism of E2020 Binding to Acetylcholinesterase: A Steered Molecular Dynamics Simulation Chunying Niu,† Yechun Xu,† Yong Xu,† Xiaomin Luo,† Wenhu Duan,† Israel Silman,‡ Joel L. Sussman,‡ Weiliang Zhu,† Kaixian Chen,† Jianhua Shen,*,† and Hualiang Jiang*,†,§ Center for Drug DiscoVery and Design, State Key Laboratory of Drug Research, Shanghai Institute of Materia Medica, Shanghai Institutes for Biological Sciences, Chinese Academy of Sciences, and Graduate School, Chinese Academy of Sciences, 555 Zuchongzhi Road, Shanghai 201203, China, Department of Structural Biology, Weizmann Institute of Science, 76100 RehoVot, Israel, and School of Pharmacy, East China UniVersity of Science and Technology, Shanghai 200237, China ReceiVed: September 17, 2005; In Final Form: September 22, 2005

The unbinding process of E2020 ((R,S)-1-benzyl-4-[(5,6-dimethoxy-1-indanon)-2-yl]-methylpiperidine) leaving from the long active site gorge of Torpedo californica acetylcholinesterase (TcAChE) was studied by using steered molecular dynamics (SMD) simulations on a nanosecond scale with different velocities, and unbinding force profiles were obtained. Different from the unbinding of other AChE inhibitors, such as Huperzine A that undergoes the greatest barrier located at the bottleneck of the gorge, the major resistance preventing E2020 from leaving the gorge is from the peripheral anionic site where E2020 interacts intensively with several aromatic residues (e.g., Tyr70, Tyr121, and Trp279) through its benzene ring and forms a strong direct hydrogen bond and a water bridge with Ser286 via its O24. These interactions cause the largest rupture force, ∼550 pN. It was found that the rotatable bonds of the piperidine ring to the benzene ring and dimethoxyindanone facilitate E2020 to pass the bottleneck through continuous conformation change by rotating those bonds to avoid serious conflict with Tyr121 and Phe330. The aromatic residues lining the gorge wall are the major components contributing to hydrophobic interactions between E2020 and TcAChE. Remarkably, these aromatic residues, acting in three groups as “sender” and “receiver”, compose a “conveyer belt” for E2020 entering and leaving the TcAChE gorge.

Introduction According to the cholinergic hypothesis, Alzheimer’s disease (AD) is associated with the impairment in cholinergic transmission that leads to memory loss.1 The inhibition of acetylcholinesterase (AChE) activity may be one of the most realistic approaches to the symptomatic treatment of AD.2 AChE is responsible for degradation of the neurotransmitter acetylcholine (ACh) in the synaptic cleft of neuromuscular junctions and of neuronal contacts in the central nervous system.3 The threedimensional structure of Torpedo californica AChE (TcAChE)4 provided valuable insight for studying the structure-function relationships of this enzyme. Structurally, there is a narrow active site gorge about 20 Å deep which consists of two separated ligand binding sites, the acylation (or active) site and the peripheral anionic binding site. The acylation site located at the bottom of the gorge contains residues involved in a catalytic triad (His440‚‚‚Glu327‚‚‚Ser200) and the important aromatic residue Trp84. The peripheral anionic binding site is located close to the mouth of the active site gorge. In this segment, an aromatic residue, Trp279, is an important element of this anionic site. At the midway of the deep gorge, a bottleneck formed by the aromatic side chains of Phe330 and Tyr121 has a size permitting a water molecule to permeate. Because the cross section of substrate and inhibitors is much larger than the size of the bottleneck, large-amplitude fluctua* To whom correspondence should be addressed. Phone: 86-2150806600-1210. Fax: 86-21-50807088. E-mail: [email protected]. † Chinese Academy of Sciences. ‡ Weizmann Institute of Science. § East China University of Science and Technology.

tions are necessary for substrate or inhibitors to enter or leave.5-7 These structural assignments were confirmed by biochemical studies, involving site-directed mutagenesis, which attested to the important role of aromatic residues in AChE proposed on the basis of the crystallographic data.8,9 Because of the key role that AChE plays in the nervous system, AChE inhibitors are also used in the treatment of various disorders such as myasthenia gravis and glaucoma. Their use has been proposed not only as a possible therapeutic approach to allay the symptoms of AD but also for the recovery of neuromuscular block in surgery.4,10 Several clinical studies with AChE inhibitors such as physostigmine (PHY) and tacrine (THA)11 have revealed modest improvements in cognitive function in AD patients. However, some factors still continue to restrict the acceptance of these drugs, such as variable oral activity, short duration of action, and side effects.12 E2020 ((R,S)-1-benzyl-4-[(5,6-dimethoxy-1-indanon)-2-yl]-methylpiperidine) is a piperidine derivative, cholinesterase inhibitor (ChEI), known by its trivial name donepezil hydrochloride and marketed as Aricept,13 for the treatment of AD.12 Pharmacological studies have shown that E2020 is a reversible and noncompetitive inhibitor of cholinesterase (ChE) and displays high selectivity for AChE in comparison to butyrylcholinesterase (BChE).14,15 It overcomes the deficits of short duration of PHY and liver toxicity of THA by producing distinct and long-lasting inhibition of brain AChE and increasing brain content of ACh in vivo.12,16 The X-ray crystal structure of the E2020-TcAChE complex revealed a high degree of structural complementarity between the inhibitor and enzyme. E2020 orientates along the active site gorge, extending from the active site, at the bottom

10.1021/jp0552877 CCC: $30.25 © 2005 American Chemical Society Published on Web 11/17/2005

E2020 Binding to Acetylcholinesterase

J. Phys. Chem. B, Vol. 109, No. 49, 2005 23731

Figure 1. Ribbon schematic representations of the solvated E2020-TcAChE model for the SMD simulation (water molecules not shown). During the SMD simulations, E2020 is pulled away from the binding gorge of TcAChE using a harmonic potential symbolized by an artificial spring connected to the center of mass of the indene ring of E2020. This pulling potential moved with constant velocity Vpull (arrow) along the Z axis. A structure diagram of E2020 is shown at the bottom right-hand corner. This picture was rendered in the POV-Ray program.18

near Trp84, to the peripheral binding site, at the top near Trp279; three major functional moieties of E2020, the benzyl group, the piperidine nitrogen, and the dimethoxyindanone moiety, make principal interactions with the active site gorge of TcAChE (Figure 1).17 Different from other AChE inhibitors, such as Huperzine A (HupA), which binds to the acylation site, E2020 is a typical bivalent AChE inhibitor that almost orients along the whole gorge of AChE, capable of binding with both the acylation (active) site and the peripheral “anionic” site simultaneously. In our previous study, taking HupA as an example, we investigated the binding and unbinding processes of acylation site inhibitors to AChE using the steered molecular dynamics (SMD) simulation method.19 The detailed analyses for the structures and interactions in the SMD trajectories detected the important residues associated with HupA binding and unbinding, the flip of the peptide bond between Gly117 and Gly118, and the dynamic effect of the buried water molecules. However, for bivalent AChE inhibitors, several questions are still open. How do bivalent inhibitors enter or leave the active site gorge of AChE? What make bivalent inhibitors cross the narrow bottleneck of AChE? What roles do aromatic residues play in the activity and function of the bivalent inhibitors? To answer these questions, we simulated the unbinding process of E2020 from AChE using the SMD method. SMD is a simulation method that has been widely used to explore the binding and unbinding properties of biomolecules and their responses to external mechanical manipulations at the atomic level.20 In SMD simulations, time-dependent external forces are applied to the ligand to facilitate its unbinding from a protein, as shown in Figure 1. Through the accelerated dissociation process of the ligand, SMD can reveal information about the protein’s flexibility and its response to the dissociation of the ligand. Analyses of interactions between the dissociating ligand and the protein and the relationship between the applied forces and the ligand position could yield important information about the structurefunction relationships of the protein-ligand complex, the binding and unbinding pathways, and possible mechanisms of ligand recognition and inhibition.19,21,22 In the following, we present the results of SMD simulations for E2020 leaving TcAChE in water solution. In this simulation study, 26 pulling

velocities were tested to pull E2020 out from the binding gorge of TcAChE, and corresponding rupture forces and mechanical works under different pulling velocities were also calculated. The subsequent detailed analyses based on a SMD simulation under a velocity of 0.005 Å/ps detected the important residues associated with E2020 unbinding, revealed the roles that aromatic residues play in E2020 association and dissociation, and provided a valid interpretation for pronounced long-lasting inhibition of this inhibitor. Models and Methods The simulation model used in the present study was built up on the basis of the X-ray crystal structure of the E2020TcAChE complex at 2.5 Å resolution taken from the Protein Data Bank (PDB),23 entry code 1EVE.17 The missing atoms of residues Asp2 and His3 were added by using the molecular modeling software Sybyl 6.7.24 The ionization states of some residues in TcAChE were determined.25 Residues His440, Glu443, and Asp392 were neutralized, as pointed out by McCammon et al.,5 and all other ionizable residues were set as in their standard protonation states. All of the simulations were performed with the parallel version of the MD program GROMACS.26 A modification of the GROMOS-87 force field27-29 was applied for the protein. The molecular topology file for E2020 was generated by the program PRODRG30 (http://davapc1.bioch.dundee.ac.uk/programs/prodrg/prodrg.html). The partial atomic charges of E2020 were determined by using the CHelpG method31 implemented in the Gaussian 98 program32 at the level of HF/6-31G*. Before the MD simulations, the E2020-TcAChE complex was solvated with simple point charge (SPC)33 water molecules in a cuboid box, in which the shortest distances between the TcAChE surface and the box walls are larger than 10 Å. To ensure overall neutrality of the system, 43 Na+ and 40 Cl- were added to replace 83 water molecules at physiological concentration in the box. The simulation system is totally composed of 51 879 atoms. The LINCE algorithm34 was used to constrain all covalent bonds to their equilibrium lengths, allowing for a time step of 2 fs. The system was minimized using the steepest-descents algorithm to the tolerance of 100 kJ/(mol‚nm). Afterward, the solvent

23732 J. Phys. Chem. B, Vol. 109, No. 49, 2005

Niu et al.

Figure 2. Time dependence of the RMSD from the crystal structure of the E2020-TcAChE complex for the CR (red) and all atoms (black) of TcAChE in the 5 ns equilibrium MD simulation.

molecules were heated to 300 K using a 25 ps MD simulation with E2020 and protein fixed. Following the equilibration, the protein was relaxed and heated to 300 K using another 25 ps MD simulation, and finally, E2020 was relaxed to heat up to 300 K by using a 5 ps MD simulation. Finally, the dynamic simulation of the entire system was carried out for 5 ns under the normal pressure (1 bar) and room temperature (300 K). The SMD simulations were performed while the equilibrated systems were obtained. E2020 was pulled out of the gorge of TcAChE under different velocities by employing different artificial harmonic potentials, and the force constant (k) was set as 280 pN/Å, as shown by the symbolic “spring” in Figure 1. Considering the structural characteristic of E2020, the harmonic potential was assigned to the center of mass of the indene ring of E2020, allowing E2020 to walk along the previously defined axis that extends from atom Ile444-Cδ to the center of mass of the atomic group of Glu73-CR, Asn280Cβ, Asp285-Cγ, and Leu333-O.19,35 Results and Discussion Equilibrium Simulation. The structure deviation of AChE from the initial X-ray crystal structure is assessed on the basis of root-mean-square deviation (RMSD). It is an important criterion for the convergence of the protein system. The RMSD of the CR and all atoms of TcAChE from its crystal structure versus simulation time are shown in Figure 2. The values of RMSD are below 0.3 nm, and the structure of TcAChE appears to have been stabilized after 1.5 ns of equilibration. At the same time, the energy components and the temperature were inspected and found to have reasonable stability throughout the 5 ns simulation (data not shown). Therefore, we may conclude that the system has been equilibrated after about 1.5 ns, which can be used as the starting point for the further SMD simulations. Comparison of Different Pulling Velocities. In general, the time scale accessible to computer simulation is often much shorter than the natural time scale of the process in studying large systems such as biomolecules.36 Therefore, to study the binding between the protein and ligand, this process should be accelerated and the pulling velocity (Vpull) should be an important parameter in the SMD simulation.19 The magnitude of the pulling velocity will influence the force profile that could be related to essential rupture events such as the breaking of hydrogen bonds or hydrophobic interaction between E2020 and the active gorge of TcAChE. Higher pulling velocity in SMD

Figure 3. (A) Computed unbinding forces as a function of pulling velocity (Vpull). The inset shows the same data on a logarithmic scale. (B) Computed mechanical works as a function of pulling velocity (Vpull). The inset shows the same data on a logarithmic scale. The red points are the values with the pulling velocity 0.005 Å/ps.

simulation may lead to remarkable nonequilibrium effects, which may introduce obvious errors into the simulation results.37 The low-velocity SMD simulation that was carried out on a millisecond time scale can overcome these disadvantages and reproduce actual atomic force microscopy experiments; however, the corresponding computational cost will be very expensive. To obtain information about the dependency of the rupture forces on the pulling velocity and to find an appropriate simulation velocity, 26 different velocities from 0.001 to 0.5 Å/ps were performed in our SMD simulations. E2020, under the largest pulling velocity of 0.5 Å/ps, only needed 50 ps to leave the TcAChE gorge completely. However, it took about 20 ns to pull E2020 out of the gorge under a pulling velocity of 0.001 Å/ps. Figure 3A shows the results of the largest rupture forces as a function of pulling velocity. The computed rupture forces as a function of pulling velocity change from 400 to 1100 pN, showing obvious fluctuation, but in general, the maximum rupture force trends to decrease as the pulling velocity becomes slower. To check the influence of friction, the mechanical work (W(z,V)) done by pulling E2020 out of the gorge from 0 to z (20 Å, the length of the active gorge of AChE) under the action of a force, F(z′,V), was calculated by eq 1

W(z,V) )

∫0zF(z′,V) dz′

(1)

E2020 Binding to Acetylcholinesterase

Figure 4. Rupture force (A), numbers for direct hydrogen bonds (B), atom pairs of hydrophobic interactions (C), and water bridges (D) between E2020 and TcAChE versus time in the process of E2020 leaving the TcAChE binding gorge. The data in parts A and C are smoothed over 10 ps intervals.

where F(z′,V) ) k(Vt - z′) is the external force due to the harmonic constraint at position z′ with pulling velocity V. Figure 3B shows the results of the mechanical works under different pulling velocities. In general, the work increases as the pulling velocity increases, because, as expected, more work is dissipated to conquer larger friction produced by faster pulling velocity of the moving molecule. Despite the fact that the pulling velocity 0.001 Å/ps we used is very slow in the present study (it costs four processors on a SGI Origin3800 2 months for the simulation), values of the rupture force and mechanical work of this velocity do not decrease much in comparison with other velocities. Moreover, the aim of this study is to address the dynamical processes of the E2020-TcAChE unbinding rather than reproduce the accurate binding force or binding free energy. Therefore, in the following, we use the SMD result under a pulling velocity of 0.005 Å/ps, which produces the lowest rupture force and the lowest mechanical works with less computer resources, to discuss the unbinding process of E2020 from TcAChE. Rupture Force Profile. Figure 4A shows the rupture force profile in the leaving trajectory of E2020. Different from the unbinding process of HupA from TcAChE,19 there are four major peaks in the rupture force profile of E2020; among them, the highest peak is located at about 2150 ps with the largest rupture force of ∼550 pN. To elucidate the fluctuation of the rupture forces along with the E2020 moving trajectory, the interactions between E2020 and TcAChE in the process of E2020 leaving the gorge were analyzed. The direct hydrogen bonds (DHBs), water bridges (WBs), and hydrophobic interactions (HIs) between E2020 and TcAChE were considered. The interaction features were analyzed by using the tools of GROMACS26 and the LIGPLOT program.38 The results are shown in Figure 4B-D. To illustrate the most important residues of TcAChE concerning E2020 unbinding, six typical snapshot structures of the E2020-TcAChE complex were isolated from the SMD trajectory, as shown in Figure 5. While moving along the binding gorge, interaction between E2020 and TcAChE changes from one kind to another; for example, old hydrogen bonds (HBs) break and new HBs form continuously. Figure 4B illustrates the variation of the direct HBs (DHBs) formed between E2020 and TcAChE. The number of DHBs is always less than four. Tyr121, Asn280, Phe284, Ser286, Phe288, and Arg289 are the main residues forming

J. Phys. Chem. B, Vol. 109, No. 49, 2005 23733 DHBs during E2020 unbinding. From 0 to 300 ps, E2020 forms DHBs mainly with Tyr121, Phe288, Ser286, and Arg289; from 300 to 1000 ps, only Arg289 forms DHBs with E2020; when the dimethoxyindanone moiety of E2020 begins to move out from the peripheral anionic site at about 1000 ps, E2020 forms new DHBs with Ser286, Asn 280, and Phe284, and this process holds up to about 1100 ps; then, a DHB vacuum appears from 2160 to 2650 ps; at about 2650 ps, a few DHBs are formed with the nitrogen of Phe284; no DHBs are observed between E2020 and TcAChE after 2740 ps. All oxygen and nitrogen atoms of E2020 take part in the formation of hydrogen bonds, but mostly, E2020 forms hydrogen bonds with residues of TcAChE through its carbonyl oxygen atom (O24 in Figure 1) of the dimethoxyindanone group (Figure 5). Figure 4C presents the change for hydrophobic interaction (HI) pairs between E2020 and TcAChE during E2020 moving through the gorge. In general, the HI pairs are decreasing along with E2020 leaving away the active binding site. From 0 to 1000 ps, because E2020 occupies most of the gorge (Figure 5A and B), the hydrophobic interaction is relatively stable, as indicated by the HI pairs. When the dimethoxyindanone moiety of E2020 moves out of the gorge, the HI pairs decrease dramatically. After 3000 ps, the HI pairs gradually disappear. In the whole process of E2020 leaving, the aromatic residues contribute about 70% to hydrophobic interactions. Table 1 lists the main residues in the gorge contributing to the hydrophobic interaction, which indicates that most of the major residues contributing to the HI are aromatic residues. This may be assigned to the fact that the walls of the gorge are lined predominantly by aromatic residue side chains.4,39 The pairs of HIs demonstrated that Trp84 and Trp279 are the two most important residues to the hydrophobic interactions between E2020 and TcAChE; their functions will be discussed below. Both X-ray crystal structures of inhibitor-AChE complexes4,40-42 and theoretical analyses43,44 have demonstrated the importance of water molecules in the ligand-AChE binding; that is, water molecules bridge ligands with AChE. Therefore, more attention was paid to the water bridges (WBs) between E2020 and TcAChE during the SMD simulation. The number of water bridges between E2020 and the residues of TcAChE through a water molecule is displayed in Figure 4D. Different from the WBs formed in the HupA-AChE complex,19 E2020 forms WBs mainly with the peripheral residues of the gorge using two oxygen atoms and one nitrogen atom. From 1000 to 3000 ps, WBs are formed with and break from Ser286, Asn280, Phe288, and Phe284 (Figure 5C and E). As mentioned above, these residues are also involved in the formation of DHBs. Therefore, we conclude that this alternative formation between DHBs and WBs discourages E2020 dissociation from the gorge of TcAChE and, reversely, these DHBs and WBs are possibly the driving force for E2020 binding to TcAChE. The force peaks and troughs in the force profile obtained from the SMD simulation clearly reflect the complexity of the unbinding process of E2020 out of TcAChE, as shown in Figure 4A. There are several major peak regions that correlate with the break of the DHBs, hydrophobic interactions, and WBs. At ∼270 ps, E2020 forms strongly hydrogen bonds with Phe288 and Ser286; at the same time, the hindrance between the benzene ring of E2020 and Phe330 prevents E2020 from going through the bottleneck (a detailed analysis of bottleneck change will be discussed below) (Figure 5A). All of these produce the first force peak. However, the rupture force (∼200 pN) corresponding to the first force peak is less than those corresponding to other high peaks. It is different from the unbinding process of

23734 J. Phys. Chem. B, Vol. 109, No. 49, 2005

Niu et al.

Figure 5. Snapshots isolated from the unbinding SMD trajectory of E2020 at (A) 270 ps, (B) 970 ps, (C) 1615 ps, (D) 2000 ps, (E) 2150 ps, and (F) 2700 ps. E2020 is drawn in yellow color (oxygen atoms are red, and nitrogen atoms are blue). The dashed lines present hydrogen bonds and water bridges. Water molecules are shown as a ball model. Hydrogen atoms are not shown.

TABLE 1: Major Residues of TcAChE Forming Hydrophobic Interactions (HIs) with E2020 and Their Total Atom Pairs of HIs during the SMD Simulation atom pairs of atom pairs of residues hydrophobic interaction residues hydrophobic interaction Trp279 Trp84 Leu282 Tyr334 Tyr121 Phe288 Phe330

10477 7911 5663 4910 3657 2912 2676

Phe284 Ser286 Phe331 Phe290 Tyr70 Asn280 Asp72

2627 2228 2191 2157 2096 1928 1622

HupA from AChE, where the highest force appears when HupA goes through the bottleneck.19 The second force peak appears at about 970 ps (∼440 pN) when the benzene ring of E2020 passes the bottleneck of the gorge (active site). Breaking of the hydrogen bond between Arg289 and E2020 and the hydrophobic interaction between E2020 and TcAChE cause the second force peak (Figure 5B). At about 2150 ps, most of E2020 is at the peripheral anionic site of the gorge. At this position, the benzene ring of E2020 closely contacts with Tyr70, Tyr121, and Trp279. At the same time, O24 of E2020 forms a DHB and a WB with

Ser286 (Figure 5E). Breaking these interactions produces the highest force peck (∼550 pN). At about 2700 ps, although most of the E2020 moiety has moved out of the gorge, the benzene ring of E2020 still contacts closely with Trp279. Meanwhile, a DHB is formed between O24 of E2020 and Phe284 (Figure 5F). These interactions discourage E2020 from leaving and produce the fourth force peak (∼370 pN). It indicates that the peripheral anionic site of the gorge is the main site that prevents E2020 from leaving the gorge. Change of the Bottleneck. Although E2020 orientates along the active site gorge, the steric hindrance of the bottleneck may still affect E2020 dissociation. The crystal structures of AChE revealed that the bottleneck, which mainly is composed of Tyr121 and Phe330, locates at the middle of the binding gorge.4 The minimal distance between Tyr121 and Phe330 (Y-F distance) was monitored during the SMD simulation. Figure 6A displays the fluctuations of the Y-F distance during the SMD simulation. The X-ray crystal structure indicates that the piperidine ring of E2020 inserts into the bottleneck between Tyr121 and Phe330,17 forming a cation-π interaction with Phe330 and a WB interaction with Tyr121. Conformational

E2020 Binding to Acetylcholinesterase

Figure 6. Time dependence of the minimal distance between Tyr121 and Phe330 of TcAChE (A), dihedral angle N14-C17-C18-C23 (B), dihedral angle C13-N14-C17-C18 (C), dihedral angle C8-C10C11-C16 (D), and dihedral angle C7-C8-C10-C11 (E) of E2020. The atomic numbers of E2020 are shown in Figure 1.

J. Phys. Chem. B, Vol. 109, No. 49, 2005 23735 change should occur during E2020 leaving away from the gorge of TcAChE. Accordingly, we monitored the dihedral angles of the four single bonds connecting the piperidine ring to the benzene ring (bonds N14-C17 and C17-C18) and dimethoxyindanone (bonds C11-C10 and C10-C8). Figure 6B-E displays the fluctuations of these dihedral angles during the SMD simulation. At the beginning, the benzene ring of E2020 is under the bottleneck, poised to directly contact with Trp84 at the bottom of the gorge through π-π stacking. As E2020 moves toward the outside of the gorge, the benzene ring of E2020 approaches the benzene ring of Phe330 (from 0 to 350 ps), forming π-π stacking for a while. Then, the benzene ring of E2020 clashes with the side chain of Phe330, discouraging E2020 to go through the bottleneck (Figure 5A). The clash between E2020 and Phe330 enlarges the entrance of the bottleneck, and the Y-F distance increases from 0.4 to 0.75 nm (Figure 6A). After the benzene ring of E2020 crosses the bottleneck, the Y-F distance recovers to its original value and the bottleneck closes. At ∼3500 ps, although there is no direct interaction between E2020 and the bottleneck, the bottleneck opens again. This is in agreement with the conclusion of Tai et al.6,7 and Xu et al.19 that the bottleneck closing and opening is an intrinsic property of AChE. In the same time period (0-350 ps), the dihedral angle of bond C17-C18 (N14-C17-C18-C23) of E2020 flips from -50 to 100° (Figure 6B) and bond C8-C10 (dihedral angle C8C10-C11-C16) rotates from 50 to -50° (Figure 6D); no dramatical changes occur for the other two dihedral angles. This

Figure 7. Profiles for center-center distances of the benzene ring of E2020 to the aromatic rings of the group 1 residues (Trp84 (d1), Phe330 (d2), and Phe331 (d3)) (A), group 2 residues (Phe288 (d4), Phe290 (d5), and Tyr334 (d6)) (B), and group 3 residues (Tyr70 (d7), Tyr121 (d8), and Trp279 (d9)) (C); (D) profiles for the center-center distances between the indene ring of E2020 and aromatic rings of Trp279 (d10), Tyr70 (d11), Tyr121 (d12), and Phe284 (d13).

23736 J. Phys. Chem. B, Vol. 109, No. 49, 2005

Niu et al.

Figure 8. Major hydrophobic interactions between E2020 and TcAChE in the crystal structure (A) and the snapshots at 250 ps (B), 1350 ps (C), and 2000 ps (D) during unbinding processes. E2020 is drawn in yellow color (oxygen atoms are red, nitrogen atoms are blue). The dashed lines are center-center distances between the benzene/indene ring of E2020 and aromatic residues of the gorge of TcAChE.

indicates that both the benzene ring and dimethoxyindanone may adjust their conformations when E2020 passes the bottleneck, avoiding serious conflict with Tyr121 and Phe330. Functions of Aromatic Residues in the Gorge. The above analyses indicate that aromatic residues are major components contributing to hydrophobic interactions for E2020 to TcAChE (Table 1). To further probe the roles that aromatic residues in the gorge play in the E2020 binding, center-center distances of the benzene and indene rings of E2020 to the rings of the aromatic residues listed in Table 1 were monitored during the SMD simulation; the results are shown in Figure 7. According to the center-center distances and E2020-AChE interaction features, aromatic residues in the gorge can be clearly categorized into three groups: group 1 contains Trp84, Phe330, and Phe331; group 2 consists of Phe288, Phe290, and Tyr334; and group 3 includes Tyr70, Tyr121, and Trp279. We define the center-center distances of the benzene ring to the aromatic rings of these residues as d1-d9 and those of the indene ring to the aromatic rings of Trp279, Tyr70, Tyr121, and Phe284 as d10d13. By analyzing the profiles of these distances in detail along the simulation time, an interesting phenomenon has been addressed that each profile contains platform(s). In general, the appearing times and durations for the platforms are also approximately sorted by the groups of the aromatic residues

mentioned above. In the X-ray crystal structure, the benzene ring of E2020 forms π-π stacking interactions with the aromatic rings of Trp84 (Figure 8A). After ∼80 ps of simulation, the aromatic rings of group 1 residues move toward the benzene ring of E2020; then, the profiles of the center-center distances d1 and d3 form the first platforms from 80 to 1000 ps and that of the center-center distance d2 forms a platform from 80 to 1500 ps. During this period, the distances (d1-d3) are relatively stable at ∼0.5 nm. From 1000 to 1500 ps, the profiles of d1 and d3 form the second platforms, fluctuating slightly around ∼0.8 and ∼0.7 nm, respectively (Figure 7A). Each profile of the distances formed between the benzene ring of E2020 and the group 2 residues (d4-d6) also contains platforms. d4 drops dramatically from 1.3 to 0.8 nm in the period from 0 to 80 ps and then forms the first platform at ∼0.8 nm until 900 ps, and d6 forms the first platform at ∼0.7 nm from 80 to 900 ps; then, these two distances form their second platforms at ∼0.5 nm from 900 to 2200 ps. d5 decreases from 1.1 to 0.7 nm during the time period from 0 to 80 ps and then forms a long platform from 80 to 2200 ps at ∼0.6 nm (Figure 7B). d7 and d9 form their first platforms from 400 to 1500 ps at ∼1.0 nm and their second platforms from 1600 to 2600 ps at ∼0.5 nm. The profile of d8 forms platforms at ∼0.5 nm from 400 to 1000 nm and from 1600 to 2600 ps, respectively (Figure 7C). Similarly, the

E2020 Binding to Acetylcholinesterase profiles of center-center distances between the indene ring and the side chains of Trp279 (d10), Tyr70 (d11), Tyr121 (d12), and Phe284 (d13) have platforms from 0 to 1000 ps (Figure 7D). The X-ray crystal structure of the E2020-TcAChE complex indicates that Trp84 is located near the bottom of the gorge; Phe330 and Phe331 are positioned near the bottleneck of the gorge (Figure 8A). After a short-time SMD simulation, the side chains of these two residues move close to the benzene ring of E2020. Phe288, Phe290, and Tyr334 are situated around the middle part of the gorge, and Tyr70, Tyr121, and Trp279 are laid around the lip of the gorge (Figure 8B). Therefore, the characters of the center-center distance profiles indicate that the aromatic residues lining the AChE gorge compose a “conveyer belt” for E2020 entering and leaving. When E2020 leaves the gorge, the aromatic side chains of the group 1 residues follow the benzene ring of E2020 moving toward the entrance of the gorge in the period from 0 to 1500 ps; cooperatively, the aromatic side chains of the group 3 residues move following the indene ring of E2020 from 0 to 1000 ps. It seems that these residues discourage E2020 from leaving. During this period, the aromatic side chains of the group 2 residues wiggle toward the benzene ring of E2020 to receive it. From 900 to 2200 ps, the aromatic side chains of the group 2 residues move along the benzene ring toward the gorge entrance, and the side chains of the group 3 residues switch to pick up the benzene ring (Figure 8C). From 1600 to 2600 ps, the side chains of the group 3 residues send E2020 out following the benzene ring (Figure 8D). Reversibly, when E2020 binds to AChE, the side chains of residues consisting of group 3 at the lip of the gorge first “receive” E2020 through the π-π stacking interaction and then relay the benzene ring of E2020 to the side chains of the group 2 residues, making E2020 move more deep into the gorge. Afterward, the side chains of the group 1 residues attract the benzene ring of E2020 to the bottom of the gorge; at the same time, the aromatic rings of the group 3 residues fasten the indene ring of E2020 via π-π stacking interactions; thus, E2020 is finally fixed inside the gorge. Conclusions The present SMD simulations with different velocities have provided new insight into the unbinding mechanism of E2020 from TcAChE at the atomic level. In the nanosecond simulations, E2020 was pulled out from the long binding gorge of TcAChE with different velocities, and the unbinding force profiles were obtained. The simulation results based on the SMD simulation under a pulling velocity of 0.005 Å/ps have revealed many dynamic features of the inhibitor unbinding process and provided insights into the inhibition mechanism of E2020. The peaks in the force profiles and the SMD trajectory can be assigned to the forming and breaking of specific direct hydrogen bonds, hydrophobic contacts, and water bridges. In particular, our SMD simulations of E2020-TcAChE unbinding allow for the following conclusions: (1) Different from other AChE inhibitors such as HupA,19 the major resistance for TcAChE to prevent E2020 from leaving the gorge is from the peripheral anionic site rather than the bottleneck. When E2020 arrives at this site, the benzene ring of E2020 closely contacts with several aromatic residues (e.g., Tyr70, Tyr121, and Trp279). Meanwhile, O24 of E2020 forms a strong DHB and a WB with Ser286 (Figure 5E). These interactions cause the largest rupture force (∼550 pN) (Figure 4A).

J. Phys. Chem. B, Vol. 109, No. 49, 2005 23737 (2) The rotatable bonds of the piperidine ring to the benzene ring and dimethoxyindanone are essential for E2020 to go through the bottleneck of AChE, facilitating E2020 to adjust its conformation when it passes the bottleneck to avoid serious conflict with Tyr121 and Phe330 (Figure 6). (3) As in the binding and unbinding of other inhibitors to AChE, aromatic residues lining the wall of the gorge play important roles for E2020 binding to TcAChE. However, the roles that these aromatic residues play for E2020 binding are very particular. These aromatic residues are major components contributing to hydrophobic interactions for E2020 to TcAChE (Table 1). Remarkably, the aromatic residues lining the TcAChE gorge compose a conveyer belt for E2020 entering and leaving; three groups of aromatic residues act as “sender” and “receiver” successively while E2020 is entering or leaving the TcAChE gorge (Figures 7 and 8) Acknowledgment. This work was supported by the State Key Program of Basic Research of China (grants 2004CB518901 and 2004CB518905), 863 Hi-Tech Program (grants 2002AA233061, 2002AA104270, 2002AA233011, and 2003AA235030), and Shanghai Science and Technology Commission (grant 03DZ19228). References and Notes (1) Bartus, R. T.; Dean, R. L., III; Beer, B.; Lippa, A. S. Science 1982, 217, 408-414. (2) Kasa, P.; Papp, H.; Kasa, P., Jr.; Torok, I. Neuroscience 2000, 101, 89-100. (3) Luttmann, E.; Linnemann, E.; Fels, G. J. Mol. Model. 2002, 8, 208216. (4) Sussman, J. L.; Harel, M.; Frolow, F.; Oefner, C.; Goldman, A.; Toker, L.; Silman, I. Science 1991, 253, 872-879. (5) Wlodek, S. T.; Clark, T. W.; Scott, L. R.; McCammon, J. A. J. Am. Chem. Soc. 1997, 119, 9513-9522. (6) Tai, K.; Shen, T.; Borjesson, U.; Philippopoulos, M.; McCammon, J. A. Biophys. J. 2001, 81, 715-724. (7) Shen, T.; Tai, K.; Henchman, R. H.; McCammon, J. A. Acc. Chem. Res. 2002, 35, 332-340. (8) Taylor, P.; Radic, Z. Annu. ReV. Pharmacol. Toxicol. 1994, 34, 281-320. (9) Wong, D. M.; Greenblatt, H. M.; Dvir, H.; Carlier, P. R.; Han, Y. F.; Pang, Y. P.; Silman, I.; Sussman, J. L. J. Am. Chem. Soc. 2003, 125, 363-373. (10) Hallak, M.; Giacobini, E. Neuropharmacology 1989, 28, 199-206. (11) Davis, K. L.; Powchik, P. Lancet 1995, 345, 625-630. (12) Sugimoto, H.; Iimura, Y.; Yamanishi, Y.; Yamatsu, K. J. Med. Chem. 1995, 38, 4821-4829. (13) Nightingale, S. L. JAMA 1997, 277, 10. (14) Loewenstein, Y.; Gnatt, A.; Neville, L. F.; Soreq, H. J. Mol. Biol. 1993, 234, 289-296. (15) Giacobini, E.; Zhu, X. D.; Williams, E.; Sherman, K. A. Neuropharmacology 1996, 35, 205-211. (16) Yamanishi, Y.; Araki, A.; Kosasa, T.; Ogura, H.; Yamatsu, K. Soc. Neurosci. Abstr. 1988, 14, 99. (17) Kryger, G.; Silman, I.; Sussman, J. L. Structure 1999, 7, 297307. (18) Pov-Ray-Team. PoV-Ray, version 3; 1999 (www.povray.org). (19) Xu, Y. C.; Shen, J. H.; Luo, X. M.; Silman, I.; Sussman, J. L.; Chen, K. X.; Jiang, H. L. J. Am. Chem. Soc. 2003, 125, 11340-11349. (20) Izrailev, S.; Stepaniants, S.; Isralewitz, B.; Kosztin, D.; Lu, H.; Molnar, F.; Wrigers, W.; Schulten, K. Steered molecular dynamics. In Computational Molecular Dynamics: Challenges, Methods and Ideas, Lecture Notes in Computational Science and Engineering; Denflhard, P., Hermans, J., Leimkuhler, B., Mark, A., Skeel, R., Reich, S., Eds.; SpringerVerlag: Berlin, 1998; pp 39-65. (21) Shen, L.; Shen, J.; Luo, X.; Cheng, F.; Xu, Y.; Chen, K.; Arnold, E.; Ding, J.; Jiang, H. Biophys. J. 2003, 84, 3547-3563. (22) Cheng, F.; Shen, J.; Luo, X.; Jiang, H.; Chen, K. Biophys. J. 2002, 83, 753-762. (23) Berman, H. M.; Westbrook, J.; Feng, Z.; Gilliland, G.; Bhat, T. N.; Weissig, H.; Shindyalov, I. N.; Bourne, P. E. Nucleic Acids Res. 2000, 28, 235-242. (24) Sybyl Molecular Modeling System, version 6.7; Tripos Associates: St. Louis, MO.

23738 J. Phys. Chem. B, Vol. 109, No. 49, 2005 (25) Gilson, M. K.; Straatsma, T. P.; McCammon, J. A.; Ripoll, D. R.; Faerman, C. H.; Axelsen, P. H.; Silman, I.; Sussman, J. L. Science 1994, 263, 1276-1278. (26) Berendsen, H. J. C.; Spoel, D. v. d.; Drunen, R. v. GROMACS: Comput. Phys. Commun. 1995, 95, 43-56. (27) van Gunsteren, W. F.; Berendsen, H. J. Gromos-87 manual.; Biomos BV: Nijenborgh, Groningen, The Netherlands, 1987. (28) van Buuren, A. R.; Marrink, S. J.; Berendsen, H. J. C. J. Phys. Chem. 1993, 97, 9206-9212. (29) van Gunsteren, W. F.; Billeter, S.; Eising, A.; Hunenberger, P. H.; Kruger, P.; Mark, A. E.; Scott, W. R. P.; Tironi, I. G. Biomolecular simulation: The GROMOS96 Manual and User Guide; Biomos BV: Groninger, Zurich, Switzerland, 1996. (30) Van Aalten, D. M. F.; Bywater, R.; Findlay, J. B. C.; Hendlich, M.; Hooft, R. W. W.; Vriend., G. J. Comput.-Aided Mol. Des. 1996, 10, 255-262. (31) Breneman, C. M.; Wiberg, K. B. J. Comput. Chem. 1990, 11, 361397. (32) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Zakrzewski, V. G.; Montgomery, J. A., Jr.; Stratmann, R. E.; Burant, J. C.; Dapprich, S.; Millam, J. M.; Daniels, A. D.; Kudin, K. N.; Strain, M. C.; Farkas, O.; Tomasi, J.; Barone, V.; Cossi, M.; Cammi, R.; Mennucci, B.; Pomelli, C.; Adamo, C.; Clifford, S.; Ochterski, J.; Petersson, G. A.; Ayala, P. Y.; Cui, Q.; Morokuma, K.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Cioslowski, J.; Ortiz, J. V.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Gonzalez, C.; Challacombe, M.; Gill, P. M. W.; Johnson, B. G.; Chen, W.; Wong, M. W.; Andres, J. L.; Head-Gordon,

Niu et al. M.; Replogle, E. S.; Pople, J. A. Gaussian 98, revision A.4; Gaussian, Inc.: Pittsburgh, PA, 1998. (33) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Interaction models for water in relation to protein hydration; Pullman, B., Ed.; D. Reidel: Dordrecht, The Netherlands, 1981; pp 331342. (34) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. J. Comput. Chem. 1997, 18, 1463-1472. (35) Botti, S. A.; Felder, C. E.; Lifson, S.; Sussman, J. L.; Silman, I. I. Biophys. J. 1999, 77, 2430-2450. (36) Park, S.; Khalili-Araghi, F.; Tajkhorshid, E.; Schulten, K. J. Chem. Phys. 2003, 119, 3559-3566. (37) Isralewitz, B.; Gao, M.; Schulten, K. Curr. Opin. Struct. Biol. 2001, 11, 224-230. (38) Wallace, A. C.; Laskowski, R. A.; Thornton, J. M. Protein Eng. 1995, 8, 127-134. (39) Axelsen, P. H.; Harel, M.; Silman, I.; Sussman, J. L. Protein Sci. 1994, 3, 188-197. (40) Harel, M.; Schalk, I.; Ehret-Sabatier, L.; Bouet, F.; Goeldner, M.; Hirth, C.; Axelsen, P. H.; Silman, I.; Sussman, J. L. Proc. Natl. Acad. Sci. U.S.A. 1993, 90, 9031-9035. (41) Raves, M. L.; Harel, M.; Pang, Y. P.; Silman, I.; Kozikowski, A. P.; Sussman, J. L. Nat. Struct. Biol. 1997, 4, 57-63. (42) Harel, M.; Quinn, D. M.; Nair, H. K.; Silman, I.; Sussman, J. L. J. Am. Chem. Soc. 1996, 118, 2340-2346. (43) Koellner, G.; Kryger, G.; Millard, C. B.; Silman, I.; Sussman, J. L.; Steiner, T. J. Mol. Biol. 2000, 296, 713-735. (44) Henchman, R. H.; Tai, K.; Shen, T.; McCammon, J. A. Biophys. J. 2002, 82, 2671-2682.