Xiong B et al / Acta Pharmacol Sin 2004 Jun; 25 (6): 705-713
Bing XIONG, Xiao-qin HUANG, Ling-ling SHEN, Jian-hua SHEN 2, Xiao-min LUO, Xu SHEN, Hua-liang JIANG2, Kai-xian CHEN
Drug Discovery and Design Center, State Key Laboratory of Drug Research, Shanghai Institute of Materia Medica, Shanghai Institutes for Biological Sciences, Chinese Academy of Sciences, Shanghai 201203, China
1 Project supported by the Shanghai Key Program of Basic Research (No 02DJ14070), the National Natural Science Foundation of China (No 20072042), and the State Key Program of Basic Research of China (No 2002CB512802)
2 Correspondence to Prof Jian-hua SHEN, Prof Hua-liang JIANG. Phn 86-21-5080-7188. Fax 86-21-5080-7088. E-mail jhshen@mail.shcnc.ac.cn, jiang@iris3.simm.ac.cn
Received 2003-07-04 Accepted 2004-01-19
KEY WORDS molecular dynamics; secretase; essential dynamics; computer simulation; Alzheimer disease
ABSTRACT
AIM: Based on the structural analysis to reveal the mechanism of ligand binding to b-secretase and the specificity of each binding sub-site. METHODS: Molecular dynamics was used to simulate on the ligand free b-secretase and ligand bound b-secretase. The trajectories were analyzed using the essential dynamics, and the significant conformational change was illustrated employing the DynDom program. RESULTS: The essential dynamics and DynDom analyses clearly showed that the b-secretase experienced a large conformational change upon the substrate or inhibitor binding. The flap structure adopted a swing motion, gradually covering the active site to facilitate the ligand binding process. Residues Ser86 and Ile87 served as the hinge point. Inhibitor-enzyme interaction analysis revealed that residues at P2, P1, and P1' positions of the inhibitor were very important for the binding, and residues at P2' and P3' positions may be modified to improve the binding specificity. S3 subsite of the enzyme still had space to modify the inhibitors in increasing the binding affinity. CONCLUSION: The information presented here is valuable and could be used to identify small molecular inhibitors of b-secretase.
INTRODUCTION
Alzheimer disease (AD) is a progressive neuro-degenerative disorder characterized by the presence of amyloid plaques composed mainly of the 39-42 amino acid amyloid b (Ab) peptide produced by the b-secretase and g-secretase proteolytically cleaving the amyloid precursor protein (APP)[1-4]. The peptide then forms insoluble aggregates, which ultimately leads to neuron loss and dementia. Several studies have found that b-secretase firstly cleaves APP to generate the N-terminus of Ab and then the b-secretase cleaves at the C terminus, leading to the release of Ab, in which the g-secretase is a rate-limiting enzyme in this key event[5-7]. Recently, this enzyme was independently cloned by five groups at the same period and was proved to be the predominant activity of proteolytic cleavage of APP at the b-site in human brain tissue[8-11]. Experimental results of overexpression and knockout of this enzyme further demonstrate that the g-secretase is the most potential target compared with the a- and gsecretase for therapeutic intervention in AD[12,13].
It has also been generally recognized that the b-secretase belongs to the superfamily of aspartic protease and of the class I trans-membrane protein, and more like the pepsin subfamily based on the structural analysis and the similarity of amino acid sequences[14]. Although the b-secretase is fully composed of the signal peptide, pro-peptide, extra-cellular protease domain, trans-membrane region, and a cytoplasmic domain, the protease domain actually plays the most important role for proteolytic cleavage of APP, and the pro-peptide has less effect on the active site of b-secretase. Dynamically, substrate and inhibitors bind b-secretase with a two-step process, when the "tightening up" of the initial encounter complex occurs[15]. The binding process induces structural rearrangement (most likely the "flap" closing) of the initial enzyme-ligand encounter, leading to a binary complex of higher affinity. Several motion forms, such as tip curling motion and hinge-bending motion, have been suggested for the systems of HIV-1 protease and pepsin-like protease in dynamic and energetic aspects of the substrate or inhibitor binding mode[16]. Although sharing common conserved catalytic triad Asp-Thr(Ser)-Gly with pepsin-like protease, the b-secretase has its own structural specificity as it shows almost no observable binding affinities with most of the available protease inhibitors. Therefore, whether the conformational change resulted from ligand binding is similar to or different from the proposed motion forms becomes to be our great interest on the conformational flexibility of the b-secretase.
To investigate the dynamic properties of the b-secretase and find out the structural factors of different subsites for the selectivity of inhibitor binding to this specific secretase, long-time molecular dynamics (MD) simulations[17] have been performed for both the ligand-free (LF) and the ligand-bound (LB) states of b-secretase. By iteratively tracking the trajectory of conformational changes and applying the essential dynamics (ED) analysis technique[18], the particular motion mode has been identified, and the specific motion form of the flap region involving in the ligand binding process was revealed. Integrating simulation results of the present paper with the binding properties of the most potent ligand, OM00-3[19], the interaction mechanism and structural selectivity were reasonably demonstrated at the level of atomic resolution. The whole work here has provided valuable guidelines for further studying the functions b-secretase and for discovering more potent inhibitors.
MATERIALS AND METHODS
Molecular coordinates The crystal structure OM99-2-b-secretase complex at 1.90 Å resolution isolated from the PDB (entry code 1FKN) [14] was used in the three-dimensional (3D) model construction of OM00-3-b-secretase complex. Computational mutageneses of Val/Leu, Gln/Asp, and Ala/Val respectively at P3, P2, and P2' positions were performed on the crystal structure of OM99-2-b-secretase complex, producing a primary 3D model of OM00-3-b-secretase complex[20]. According to the experimental and theoretic studies on the aspartic protease family, the Asp228 of b-secretase catalytic triad was protonated. Atomic charges of OM00-3 at subsites P1 (Val variant) and P1' (Ala variant) were assigned according to the results of ESP calculation. The parameters for all the bonds and torsion angles of OM00-3 were set as the CHARMM19 force field[20]. The coordinates of uncomplexed state for the b-secretase were prepared by deleting OM99-2 from the 1FKN structure, and Asp32 and Asp228 at the catalytic triad were protonated to mimic the active state of the enzyme at pH=4.50. The two structures of LB and LF states for b-secretase were energetically optimized using molecular mechanics method with Amber force field[21] and Kollman-all-atom charges for 2000 steps in order to relax the steric crashes.
Molecular dynamics simulations The MD simulations were run on the parallel computer using the program EGO-VIII[22], and the CHARMM19 force field[21] was used. The TIP3P[23,24] water model was used to solve the structure models. The most popular Verlet algorithm[25] was adopted to integrate the equation of motion, and the FAMUSAMM algorithm[26] was applied to rapidly evaluate electrostatic interactions, while the lengths of bonds involving hydrogen atoms were fixed with the SHAKE algorithm[27].
In the simulation models of b-secretase, the histidine residues on the protein surface were set as dual-protonated form (Nd and Ne) according to the pH value of the bulk solution (pH=4.5). The ionization states of other acidic and basic residues were determined through pKa calculations. The solute was then solvated by a drop of water sphere by using the SOLVATE1.0 program (http://www.mpibpc.gwdg.de.abteilungen/071/solvate/node4.html), which ensures the whole simulation system to be covered by water molecules at least 10 Å thick. Thus total number of atoms for the LF and LB models of were 34 913 and 35 144, respectively (Fig 1).
Fig 1. Schematic representation of the MD simulation system. The b-secretase is represented as ribbon and the peptide inhibitor is shown as CPK model. The water molecules are illustrated as wireframe style.
The solved simulation models were minimized at constant volume. To remove poor steric contacts all the water molecules adopted a criterion that the maximum force of the whole system was set less than 10.00 kcal (mol·Å)-1. Afterwards, the whole systems were energetically optimized. The friction factor t was set as 0.1 at the end of each integration step during the minimization process, and then changed to 1.0 for the MD equilibrium simulations. As the minimization converged, the whole system was directly subjected to a slow heating procedure for about 100 ps in a heat reservoir of 300 K. After that a 4.5-ns (4.5×10-9 s) and a 1.2-ns (1.2×10-9 s) MD simulations were respectively performed for the LF state LB states of b-secretase. In the MD simulations, the time step was set as 1 fs (1×10-15 s), and the frequency for analyzing the MD output was set as 1 ps (1×10-12 s) for the LF simulation and 0.5 ps for the LB simulation, respectively.
Essential dynamics analysis Besides the routine energy and root mean square deviation (rmsd) analyses to monitor the b-secretase motions in the MD simulations, the essential dynamics (also termed as principal component (PCA)) method[18] was used to identify collective motion in the protein. The basic essential dynamics theory was outlined briefly as follows.
A 3N×3N covariance matrix can be built from N carbon-a atoms by using a method of Amadei et al[18]:
C=<(x-<x>)(x-<x>)T> (1)
Here the <> means time-averaged.
Then the QR method may be used to diagonalize this covariance matrix to compute the eigenvector and eigenvalue:
¦«=TTCT (2)
The ith column of the matrix ¦« is the ith eigenvector corresponding to the ith eigenvalue. Sorted the eigenvalue of the matrix C, the eigenvectors related to the largest eigenvalues are the important vectors that indicate the low frequency motion of the protein.
Finally we projected the trajectory to the largest 20 eigenvectors to identify the detailed structural information. DynDom[28] were adopted to analyze the domain motion of the enzyme.
RESULTS AND DISCUSSION
Structural dynamics A 4.5-MD simulation was preformed on the LF model of b-secretase, and only a 1.1 ns MD simulation was carried out on the LB model because the apo-enzyme had been equilibrious at 1.1 ns as illustrated in Fig 2. The energy trajectories and the root mean square deviations (RMSD) with respect to the crystal structure for the two simulations were shown in Fig 2.
Fig 2. The total energy and RMSD with respect to the crystal structure of the two simulation models versus the simulation time. (A) Simulation result of the LB model. (B) Simulation result of the LF model. (C) RMSD comparison between the LB model and LF model for the first 1.1 ns MD simulations.
As the energy curves illustrated, the total energy seemed to be stable after 1000 ps and 300 ps for the LF and LB models, respectively, which indicates that the LB state of b-secretase reaches its equilibration much quickly than the LF state does. The total energy of the LF system had only a trivial change (about 0.6 %), indicating that there is no high energy barriers in the energy potential surface of b-secretase, conformational changes can easily occur.
After 500 ps relaxation, the RMSD of the LF model was obviously larger than of the LB model (Fig 2C). This is reasonable because the MD simulations of both LF and LB models were started from the peptide-binding X-ray structure of b-secretase, and without restriction of ligand, the LF state of b-secretase should process a larger movement to recover its close state. After a slow but continuous closing motion (around 3.5 ns), the RMSD of the LF simulation reached a platform of about 3.5 Å.
The scaffold of the b-secretase is the b-sheet dominant structure consisting of two domains: the left lobe is composed of ten b strands and one a helix while the right lobe has sixteen b strands and several short a helices (Fig 1). DSSP[29] analysis indicated that most secondary structures were stable except the 5 short helices (each consisting of about 3 residues), which quickly unfolded at the first 1000 ps for the LF simulation. These helices included the segments of 53-56, 113-116, 171-173, 277-279, and 379-382. The very weak a helix feature of these segments may contribute to this unfolding event.
Although some conformational changes of the b-secretase were detected (see discussion below for details), the scaffold of b-secretase keeps well during the LF and LB simulations due to the large number of b sheets. The radius of gyration of the LF simulation also showed that the whole protein packing was stable and the average number was 21.70 while the fluctuation was about 0.27.
Collective motion of different Domains
Overall of the ED analysis The trajectories of MD simulations were analyzed using the essential dynamics (ED) method. It should be noted that inhibitor atoms were not considered in the ED analyses because we only focused on the collective motion of the b-secretase. The original 3N-dimensional configuration space formed with the N Ca atom coordinates could be mainly represented by the important subspaces spanned by a set of eigenvectors. Motions along the eigenvectors corresponding to the large eigenvalues were mainly low-frequency and large-amplitude fluctuations of the protein. The ligand peptide in the LB simulation was not involved in the ED analysis because we concerned mainly the remarkable motion of the enzyme. As shown in Fig 4, about 90 % of the total positional fluctuations could be described by the first 10 eigenvectors in the ligand-free simulation. However in the LB simulation, the motion was very complex and the first 100 eigenvectors must be used to describe 90 % fluctuations. The projection over the trajectories of LB and LF simulations along first three eigenvectors were plotted in the Fig 5. The LB model had lower eigenvalues in comparison with the LF model, because ligand binding could reduce the anharmonic fluctuation of the b-secretase.
Fig 3. Eigenvalue accumulations of the LF and LB simulations.
Fig 4. The projections along the first six eigenvectors over the whole trajectory of LF simulation.
The collective motions of the enzyme could be easily identified by the eigenvector projections over the trajectory and the DynDom analysis[28]. The sub-spaces of the first three eigenvectors could be contributed to the insert parts, domain-domain linking segment (residues 180-186) and flap structure.
Insert parts The main differences between the b-secretase and other pepsin sub-family structures are identified to be six segments termed as insert parts (insertion A: residues 158-167; insertion B: residues 218-221; insertion C: residues 251-258; insertion D: residues 270-273; insertion E: residues 290-295; insertion F: residues 311-317)[14]. These insert parts might execute special function. Detailed analysis of the projection data along the eigenvectors 3-5 revealed that the structural fluctuation originates from the insert parts. The most significant deviation from the crystal structure was inserts A and F (the maximum RMSD is about 9.0 Å). The root mean square displacement (data not shown) demonstrated that the inserts A, C, and F process large-amplitude conformational changed and the large motions happened at about 0.5 ns. Inserts A and F moved down, making the active site of b-secretase flatter and wider. The disulfide bond between residues 269 and 319 was close to inserts D and F, which may constrain the motion of the two inserts in a limited space. The motion of inserts D and F implies that they have the intrinsic flexibility to allow the access of different substrate to b-seretase[28].
Domain linking segment The linking segment (residue 175-189) contributes lots of fluctuations to the first six eigenvectors, especially to eigenvector 2. According to the dynamic structural information it was found that this segment deviated obviously from the crystal structure (about 11 Å). Although this segment had a large conformational change, no related motion between the two lobes was found. The b-sheets below the active sites contribute the stable force to maintain the shape of the active site.
Specific motion of the flap The first eigenvector projection data indicates that the most significant feature of conformational motion is clearly periodic (the periodicity is about 100 ps) in the LF trajectory. The snapshot structures revealed that the largest collective motion along the first eigenvector was the flap movement. In order to illustrate the flap motion mode, a vector between the centroid of six flap residues 70-75 and total protein centroid was defined and the distance and angle changes were calculated over the whole simulation time. The largest distance and angle deviations from the crystal structure are about 4.5 Å and 26 degree, respectively. The angle deviation also demonstrated that the angle of the flap structure changed at very short time intervals, which indicates the flap adopts a swing motion, therefore, it withdraws gradually from the binding site and widens the binding site. We also adopted the DynDom program[30] to analyze the domain with two snapshot structures having the minimum and maximum fluctuation in the projection along the first eigenvector. The result confirmed the hinge-bending motion of the flap, and Ser86 and Ile87 were found to act as a mechanical hinge. From Fig 5, it can be seen that the flap structure has a positive correlation with the top level of b-sheet (residues 104-106 and 43-47) in the conformational movement. The correlation coefficients of RMSDs of segments 104-106 and 43-47 with the flap structure were 0.93 and 0.78, respectively. The ground level of b-sheet had a negative correlation with the flap, which indicates the motion between them is in the opposite direction. This further supports the bending-motion of the flap. Accordingly, we assume that the flap b-sheet, residues 104-106, and residues 43-47 twisted around the mechanic hinge are involved in the collective motion. This can be easily seen from the superimposed snapshot structures that these residues consist of the top level of b-sheets in the left lobe while the hinge is in the middle level of the b-sheets. This indicates that the b-sheets are the origin of the intrinsic flexibility.
Fig 5. The flap bending motion map produced by DynDom program (left). The hinge-bending angle is 18.36 degree. The yellow residues are the moving sub-domain, and the blue residues represente the fixed sub-domain. The mechanic hinge is shown in spacefill form. The snapshot (every 500 ps) superimpose of flap structure in the LF simulation (right).
Fig 6. The important water molecules mediated hydrogen-bonding network.
Inhibitor specificity and binding mechanism of b-secretase The crystal structure of inhibitor bound b-secretase shows that the active site of the enzyme interlaces with hydrophobic and hydrophilic sub-pockets. Detailed analysis of the interaction of each amino acid of the inhibitor showed that P4, P2, and P3' positions had several hydrogen bonding interaction with the enzyme while the P3, P2, P1, P1', and P2' positions had relatively stronger hydrophobic interaction (Tab 1). The number of hydrogen bonds between P4 (Glu) and the enzyme fluctuates slightly along the trajectory. In comparison with the crystal structure, the final structure of P4 only lost the weak hydrogen bonding interaction with Gly11 and Tyr71, while the salt bridge interactions of P4 with the Arg307 and Arg235 were kept in the simulation. This is in good agreement with the result of combinatorial inhibitor library probing experiment[19]. In their work, substitution of other residues for P4 will decrease the catalytic efficiency since loss of the electrostatic interaction with the Arg235 and Arg307. The hydrogen bonding and hydrophobic interaction of the P2 (Asp) were very stable over the trajectory. Another hydrogen bonding interactions of the inhibitor with Asn233 and Gln73 were detected successively in simulation. Although the P3' (Glu) has a large RMSF value in comparison with the other residues of the inhibitor, it forms moderate hydrogen bonding interactions and hydrophobic interactions with the enzyme. The transition-state analogous residues at P1 and P1' positions interact with the b-secretase mainly through the hydrophobic interaction. Residues at P3 and P2' positions form hydrophobic interaction to b-secretase.
Tab 1. Interactions (the number of hydrogen bonds formed and the number of hydrophobic contacts) between the b-secretase and the peptide inhibitor averaged over the whole MD simulations
| Subsites |
No HB* |
No HC** |
RMSF (Å) |
| P4 |
11 |
4 |
0.84 |
| P3 |
5 |
9 |
0.60 |
| P2 |
11 |
9 |
0.53 |
| P1 |
3 |
15 |
0.50 |
| P |
2 |
13 |
0.51 |
| P |
3 |
11 |
0.64 |
| P |
8 |
8 |
1.41 |
|
P |
1 |
5 |
1.77 |
* hydrogen bonds; ** hydrophobic contacts within distance of 4 Å.
Recently, Tung et al investigated the significance of the each amino acid of the peptide inhibitor[19], which provides a valuable data for modification of the inhibitor. They reported that truncation of the P4' (Phe) did not significantly reduce the potency, which was in agreement with our simulation result: residue at P4' position form weak hydrogen bond and hydrophobic interactions to the enzyme (Tab 1). Residue at P4 position could be modified to acetyl group, which did not dramatically decrease the potency[31]. This indicated that the carbonyl could form a hydrogen bond to the NH group of Asn at P2 position, maintaining the conformation of the inhibitor. In the present work, this hydrogen bond was often detected over the whole LB simulation.
Water molecules play an important role in the aspartic proteases as has been demonstrated in HIV-1 proteases. These enzymes transfer proton through water molecules. Moreover, the water molecules mediate the hydrogen bonding interactions between the enzyme and substrate. Crystal structures of aspartic proteases indicated that one water molecule was important in stabilizing the flap and active site structures. To find similar important water molecule(s) in the b-secretase, we calculated the translational diffusion coefficient from the root mean square displacement:
<(r-r0)2>=Dt (3)
The diffusion coefficient of the crystal water in the LB simulation was about 5.01×10-5 cm2/s, which was consistent with the TIP3P diffusion coefficient. However, this value was larger than the diffusion coefficient of the water molecules from the shell of 4 Å thickness around the peptide inhibitor (1.77×10-5 cm2/s). Two water molecules were found to have very low diffusion coefficients, about 1.98×10-6 cm2/s and 2.50×10-6 cm2/s, respectively. Further analysis for the two water molecules over the MD trajectory showed that they positioned in the cavity below the flap, mediating the hydrogen bonding interactions of Tyr71 and Ser35 with Asn37, which could stabilize the flap structure and maintain the binding pocket.
Based on the detailed structural dynamics analysis, together with the discovery of two important water molecules, the binding mechanism can be assumed: when substrate or inhibitor approach b-secretase, the flap structure gradually moves with a swing and hinge bending motion to cover up the active sites; the two water molecules in the cavity below the flap serve as the hydrogen bond mediator for maintaining the shape of flap and catalytic binding sites through a hydrogen bonding network.
CONCLUSION
In the present work, we performed MD simulations on two models of b-secretase, ligand-free and ligand-bound models. The results produced by essential dynamics analysis clearly showed that the b-secretase experienced a large conformational change upon the substrate or inhibitor binding. Due to the intrinsic b-sheet flexibility the flap structure could adopt a swing motion to gradually cover the active site, which facilitated the ligand binding. The force inducing flap swing may originate from the bending hinge composed by Ser86 and Ile87. Inhibitor-enzyme interaction analysis indicated that Tyr71 played an important role in stabilizing the flap structure. It was also identified that residues of the inhibitor at P2, P1, and P1' positions were important for the binding, and the residues at P2' and P3' positions could be modified to improve the inhibitor specificity. Residues at P3 position have more space to be modified to improve the binding affinity. These results presented here are informative and could be used for further inhibitor design.