Article Information

• PDF (1808 kb)
• Full Text with Large Figures
• Supporting Material
• Cited by in Scopus (1)

PubMed

• Articles by Alemayehu A. Gorfe
• Articles by J. Andrew McCammon

Related Articles

• Structure of acetylcholinesterase complexed with E2020 (Aric...

  Structure of acetylcholinesterase complexed with E2020 (Aricept): implications for the design of new anti-Alzheimer drugs Structure, Volume 7, Issue 3, 15 March 1999, Pages 297-307 Gitay Kryger, Israel Silman and Joel L Sussman Summary Our study shows, , that the design of E2020 took advantage of several important features of the active-site gorge of AChE to produce a drug with both high affinity for AChE and a high degree of selectivity for AChE versus butyrylcholinesterase (BChE). It also delineates voids within the gorge that are not occupied by E2020 and could provide sites for potential modification of E2020 to produce drugs with improved pharmacological profiles. Summary | Full Text | PDF (239 kb)

• Flexibility of Aromatic Residues in the Active-Site Gorge of...

  Flexibility of Aromatic Residues in the Active-Site Gorge of Acetylcholinesterase: X-ray versus Molecular Dynamics Biophysical Journal, Volume 95, Issue 5, 1 September 2008, Pages 2500-2511 Yechun Xu, Jacques-Philippe Colletier, Martin Weik, Hualiang Jiang, John Moult, Israel Silman and Joel L. Sussman Abstract The high aromatic content of the deep and narrow active-site gorge of acetylcholinesterase (AChE) is a remarkable feature of this enzyme. Here, we analyze conformational flexibility of the side chains of the 14 conserved aromatic residues in the active-site gorge of AChE based on the 47 three-dimensional crystal structures available for the native enzyme, and for its complexes and conjugates, and on a 20-ns molecular dynamics (MD) trajectory of the native enzyme. The degree of flexibility of these 14 aromatic side chains is diverse. Although the side-chain conformations of F330 and W279 are both very flexible, the side-chain conformations of F120, W233, W432, Y70, Y121, F288, F290 and F331 appear to be fixed. Residues located on, or adjacent to, the Ω-loop (C67–C94), namely W84, Y130, Y442, and Y334, display different flexibilities in the MD simulations and in the crystal structures. An important outcome of our study is that the majority of the side-chain conformations observed in the 47 AChE crystal structures are faithfully reproduced by the MD simulation on the native enzyme. Thus, the protein can assume these conformations even in the absence of the ligand that permitted their experimental detection. These observations are pertinent to structure-based drug design. Abstract | Full Text | PDF (1458 kb)

• Bacterial Division: Another Way to Box in the Ring

  Bacterial Division: Another Way to Box in the Ring Current Biology, Volume 16, Issue 20, 24 October 2006, Pages R881-R884 William Margolin Summary Proper placement of the cell division site in some rod-shaped bacteria requires two different negative regulatory systems, nucleoid occlusion and the Min proteins. lacks these systems, but recent work has uncovered a novel regulator that achieves the same goals. Summary | Full Text | PDF (144 kb)

• …more

Copyright © 2008 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 94, Issue 4, 1144-1154, 15 February 2008

doi:10.1529/biophysj.107.117879

Biophysical Theory and Modeling

Dynamics of the Acetylcholinesterase Tetramer

Alemayehu A. Gorfe^{*, †, 1, ,}  , Chia-en A. Chang^{*, †, 1, ,}  , Ivaylo Ivanov^*, † and J. Andrew McCammon^{*, †, ‡, §}

^* Department of Chemistry and Biochemistry, University of California at San Diego, La Jolla, California 92093-0365
^† Center for Theoretical Biological Physics, University of California at San Diego, La Jolla, California 92093-0365
^‡ Howard Hughes Medical Institute, University of California at San Diego, La Jolla, California 92093-0365
^§ Department of Pharmacology, University of California at San Diego, La Jolla, California 92093-0365

Address reprint requests to Alemayehu A. Gorfe, Tel.: 858-822-0255.

Address reprint requests to Chia-en A. Chang, Tel.: 858-822-1469.

¹ Alemayehu A. Gorfe and Chia-en A. Chang contributed equally to this work.

Termination of nerve impulses in cholinergic synapses requires very rapid hydrolysis of the neurotransmitter acetylcholine (ACh). This is carried out by one of the most efficient enzymes known—acetylcholinesterase (AChE) 1. AChE is found in all vertebrates and invertebrates so far investigated, as well as in other species 2. Structurally, AChE belongs to the α/β hydrolase protein fold and has been purified and crystallized in the monomeric, dimeric, and tetrameric forms 3,4.

The tetramer is the most important functional form of the enzyme under physiological conditions. The crystal structure of recombinant mouse tetramer (mAChE) was solved at a resolution of 2.9Å 5. The structure was characterized by an antiparallel alignment of two canonical homodimers with a compact, pseudosquare planar arrangement of the subunits. Furthermore, the soluble, trypsin-released tetrameric form of AChE from Electrophorus electricus (eAChE) was crystallized in two distinct crystal forms 6. The two low-resolution eAChE tetrameric structures, resolved to 4.5 and 4.2Å, significantly differ in the arrangement of the catalytic domains: one is compact square nonplanar (Protein Data Bank (pdb) code 1C2O), whereas the other is loose, pseudosquare planar (pdb code 1C2B). It has therefore been suggested that the tetramer is flexible and could exist in multiple conformations separated by low barriers 5,6. Computations on the diffusion of ACh to the active sites of the tetrameric enzyme showed that the reaction rates differ for individual active sites in the compact structure 7,8. The rates are similar for individual active sites in the loose structure and an intermediate structure obtained by morphing between the two. The key findings of these studies were i), active sites are partially occluded in the compact structure with a competition for substrate between individual active sites, and ii), electrostatic forces substantially enhance the rate in the tetramer relative to the monomer.

The C-terminal ∼40-residue t-peptide 9,10, also called the tryptophan amphiphilic tetramerization (WAT) domain, interacts with the proline-rich attachment domain (PRAD) of the collagen-like Q subunit (ColQ). ColQ is a structural protein that anchors AChE to the synaptic basal lamina. However, the mAChE tetramer entirely lacked the WAT domain, which was flexible and not observed in either of the eAChE structures. The t-peptide has a series of aromatic amino acids, including three equally spaced tryptophan residues (WWW). The crystal structure of the PRAD/WAT complex was solved at a resolution of 2.35Å and showed a tight association through hydrogen bonding and hydrophobic interactions between the WWW motifs in WAT and consecutive proline residues in PRAD 11. Based on the crystal structures of the compact eAChE (pdb 1C2O) and the PRAD/WAT complex (pdb 1VZJ), Zhang and McCammon constructed a complete, roughly symmetric model of the complex of AChE tetramer (AChEt, Fig. 1) with the anchoring protein tail 12. Using block normal mode (BNM) analysis, they found several low-frequency and low-barrier normal modes corresponding to intersubunit motions due to the presence of flexible hinges. This flexibility may play a role in the near diffusion controlled breakdown of ACh 13,14.

Display large version of this figure

Figure 1

Model-built AChEt structure (left) and its schematic representation (right) showing subunits A (black), B (red), C (green), and D (blue). The ColQ onto which the WAT domains (residues 544–583) wrap around is in purple. Note the nearly symmetric arrangement of the monomers in this model. The three PAS residues (Tyr⁷², Trp²⁸⁶, and Tyr³⁴¹) are shown in yellow vdW spheres. The distances (d) between centers of subunits, the angle (θ) formed by three centers of subunits, and interfacial contact (χ) are illustrated. Note that χ represents the intersection (yellow sphere) of two selections: the numbers of atoms in one subunit that are within 10Å of any atom in a complementary subunit and vice versa (see green and lavender arcs).

To what extent does intersubunit fluctuation modulate catalytic efficiency? The answer to this question lies in being able to exhaustively sample the conformational space accessible to AChEt. Once the dominant configurations, i.e., the most populated “low energy” states, are identified, it will be possible to compute the probability of active site occlusion due to a complementary subunit and to correlate these structural features to function by, for example, calculating the reaction rates based on structures representing various conformational states. In principle, rigorous conformational sampling can be performed by atomistic molecular dynamics (MD) simulations. However, for a system with the complexity and size of AChEt, exhaustive sampling by atomistic, explicit water simulations is impractical; the low frequency and large amplitude intersubunit fluctuations cannot be sufficiently sampled in the few tens of nanoseconds generally afforded by fine-grained MD. Moreover, the high-frequency motions, which take up the major portion of atomistic simulation time, contribute very little to the large-scale intersubunit fluctuations. Therefore, less detailed but fast and inexpensive approaches are needed. On the other hand, access to the active center can also be modulated by local conformational changes in a complementary subunit. These motions may not be captured by the coarse-grained (CG) model. There is also the need to validate and parameterize the CG model. Thus, a multiscale simulation approach was adopted. In multiscale simulations, fine-grained models, such as atomistic MD, are combined with CG ones, such as Cα-based bead models. CG models, both alone and as part of a multiscale simulation scheme, are becoming increasingly useful to investigate large-scale molecular motions 15,16,17,18,19,20,21, such as those of AChEt.

In this work, we present results from multiscale simulations designed to explore the intersubunit motions. We combine all-atom MD (AA-MD) and CG Brownian dynamics (CG-BD) 22 simulations to characterize conformational states that AChEt may adopt at physiological conditions. The results show that there is indeed a substantial intersubunit movement (translational and rotational) that allows dynamic exposure of the active site in each subunit, albeit to various degrees. We also performed BNM analysis on selected snapshots from the AA-MD simulation and found three groups of principal modes, each composed of three modes representing closely related motions. Similar results were obtained from a principal component analysis (PCA) carried out on the CG-BD simulated structures.

The model-built tetrameric AChE structure (AChEt, Fig. 1) was used as a starting point for the simulations after the following modifications. Since they do not contribute to structural stability but would increase system size, the first nine residues of ColQ were omitted. Note that ColQ has 47 residues, whereas each of the four AChE monomers is 583 amino acids long, with residues 544–583 constituting the WAT domain. We introduced physiologically existing disulfide bridges that connect i), the WAT of subunits A and B at residue 580, ii), ColQ and subunit C (residues 15 and 580), and iii), ColQ and subunit D (residues 16 and 580).

A 20ns explicit water AA-MD simulation was performed. The total number of atoms in the simulation was ∼256,000. Standard simulation procedures (e.g., 23) were followed using the program NAMD 24 and the CHARMM22 force field 25. Briefly, after preparation of the system by sequential steps of energy minimization and equilibration, the production phase was carried out at constant temperature (310K) and pressure (1atm) conditions. Periodic boundary conditions with full particle mesh Ewald electrostatics, a 10Å cutoff for the van der Waals (vdW) interactions, a 12Å cutoff for nonbonded list update, the SHAKE algorithm, and RESPA multiple time stepping (with 2, 1, 2 fs) were used.

CG model was also used to represent AChEt, where each residue is represented by a single interaction center (bead) placed on the α carbon. A charged residue was assigned ±1e charge, and the beads were stepped through time by Brownian dynamics (BD) simulations. The force field parameters were based on those developed by Tozinni and McCammon 22,26, which were implemented in the UHBD simulation package 27,28. The force field has the functional form:

(1)

In our CG simulations, the degree of complexity and some unique features presented by the AChEt necessitated a number of modifications to be made. For example, as the model lacks the screening effect of explicit water molecules, the electrostatic interactions between two charged amino acids on the protein surface, Asp³¹⁰ of the catalytic domain and Arg⁵⁵¹ of the WAT, occasionally become too strong. We thus introduced a weak Morse potential for this pair to avoid unrealistically close distances between them (see Supplemental Material for parameters). A cutoff of 15Å was set for both nonbonded terms. Furthermore, the WAT helices possess a special WWW motif making strong hydrophobic stacking with the PRAD. There are also disulfide bonds within a monomer and between monomers and ColQ. To capture these special interactions which are much stronger than typical vdW interactions, a harmonic form, k_r(r–r₀)², is applied to replace the Morse potential in Eq. (1), where r is the distance between two interaction beads. Probability distribution functions for r were obtained from the 20ns AA-MD simulation, which provided force constants, k_r, of, e.g., ∼41kcal/mol Ų for a Pro-Trp stacking pair (Supplementary Material ).

Three snapshots from the AA-MD simulation collected at 1.16, 7.72, and 14.72ns were used to initiate BD simulation runs I, II, and III, respectively. The use of several snapshots ensures that, first, the model-built (AChEt) structure has become sufficiently relaxed that any remaining bad contacts are removed. Second, as some of the parameters for the CG-BD runs are derived from the initial structure, alternative starting conformations reduce undue dependence on the initial conditions. Furthermore, different random number seeds were assigned to an initial conformation to increase sampling. The BD simulations were run for 10μs or 15μs. The time step was set to 0.1ps, and the viscosity of water was set to 1 cp, which corresponds to a water temperature of 293K. The BD algorithm used here has been reported in several original publications 32,33, and all the parameters and equations can be found in a previous publication 29.

Deviation from the initial model and the structural integrity of each individual subunit during the simulations was monitored by the root mean-square deviations (RMSD). The fluctuations within a subunit were studied using the root mean-square fluctuations (RMSF).

The method used to compute BNM, where each residue is treated as a block, has been published previously 34. The method has also been applied on the AChEt model 12. The BNMs were calculated for four selected snapshots from the AA-MD simulation. We chose the last snapshot and the three snapshots used for starting independent CG-BD runs. PCA was carried out for the CG-BD trajectories using the structural analysis program Bio3D and its default parameters 35.

To investigate the large-scale intersubunit dynamics, three variables were defined (Fig. 1): the distance (d) between the center of mass of subunits in contact with one another (cf. A-B, B-C, C-D, and D-A); the angle (θ) subtended by vectors connecting the center of a subunit to the centers of the other two subunits it is in contact with (cf., and ); and interfacial contact (χ), defined as the average of the number of atoms in subunit i within 10Å of any atom in complementary subunit j, and the number of atoms in j that are within 10Å of any atom in i. As illustrated in Fig. 1 (right), χ is the intersection between two selections and represents contacts of surface atoms only. We obtained the same results by using the interfacial surface area of a subunit occluded from solvent by an opposing subunit (see Supplementary Material ). To calculate the change in solvent-accessible surface area (ΔSASA) from the CG-BD trajectories, the CG bead radii (which account for the size of each residue) were introduced into the CHARMM parameter file. Then, the solvent-accessible surface for individual subunits and for pairs of subunits (dimer) was computed with the CHARMM program 36, using the new radii and a regular probe radius of 1.4Å. ΔSASA of subunit A due to subunit B is then obtained by subtracting the SASA of the dimer A-B from the sum of SASAs of subunits A and B in isolation.

As a qualitative measure of occlusion of the active site of a subunit by a complementary subunit, we checked for the number of atoms in subunit j that are within 16Å (approximately twice the average length of a residue plus the radius of ACh) of the peripheral anionic site (PAS) residues Tyr⁷², Trp²⁸⁶, and Tyr³⁴¹ residing in subunit i. These residues were selected because they can interact with a flexible loop in a complementary subunit in the tetramer 5. We would like to stress that although the three residues Tyr⁷², Trp²⁸⁶, and Tyr³⁴¹ at the PAS are preselected, the occlusion as defined here is a dynamic measure that for each frame counts all atoms in an opposing subunit that make contact (cutoff 16Å) with the preselected fixed residues. This can occur through translation and/or rotation of subunits or through local movements of surface loops or through a combination of all three.

The intersubunit motions may serve as a gate to control binding site accessibility. In this study, we assume that the active sites have two states, open and closed, which may be described as follows:

(2)

(3)

(4)

(5)

We first present the intrasubunit structural properties during the simulations. Then, we present the large-scale motions as suggested by normal mode analysis and PCA, followed by the analysis of the intersubunit rotational and translational motions together with interfacial properties. We will then show how the intersubunit fluctuations affect access to active centers. As we are interested mainly in the large-scale intersubunit motions that would affect the accessibility of the active sites to substrates, in the following we focus on the four catalytic domains (residues 1–543) only, unless specifically mentioned. Note that many CG-BD simulations were carried out, but only one representative run (run II) is shown in Figure 2 and Figure 3 and Figure 4 and Figure 6.

Examples of the RMSD and RMSF derived from AA and CG simulations are plotted in Fig. 2 for each individual subunit. We can see that the subunits are stable in both the AA and CG simulations (Figure 2A). The RMSD rose to 4Å in the 20ns production run of the AA simulation, whereas the CG simulation fluctuated between 2 and 5.5Å. Although the values appear somewhat high, they are reasonable for a system of this size. Furthermore, some reorganization is expected upon the relaxation of the initial model. As expected, however, the AA-MD run is not yet fully equilibrated, although it is sufficiently stabilized for the purposes of this work. Note that the CG simulations, including the one displayed in Fig. 2, were started from a snapshot of the AA simulation; thus the starting conformation had already moved away from the initial model.

Display large version of this figure

Figure 2

Time evolution of AChEt catalytic subunits’ (residues 5–543, excluding the WAT domain) structural properties during the 20ns AA-MD (left) and one of the (CG-BD, right) simulations. (A) RMSD and (B) RMSF in Å calculated after individually fitting the Cα atoms of each catalytic subunit. Color code: A (black), B (red), C (green), and D (blue).

The calculated RMSFs show that both simulations suggest substantial internal mobility of the subunits (Figure 2B). For example, consistent with previous studies 5,39,40, the Ω-loop (residues 257–272) 39 exhibits significant mobility in all cases. Although the RMSF of the Cα atoms derived from the AA simulations is higher for a number of residues, the similar pattern between the AA and CG simulations suggests that the extremely fast CG model qualitatively captures the most relevant motions suggested by the AA-MD. The CG model was expected primarily to sample the long timescale low frequency motions; differences in, for example, time steps (2 fs vs. 100 fs) between AA-MD and CG-BD imply that fluctuations would not be identical.

The initial structure used in this study was a model built from two x-ray structures of different resolutions (1C2O and 1VZJ). BNM analysis performed on the model suggested large-scale motions. To investigate whether those motions are preserved when the structure evolves during the MD simulation, we recomputed the BNM on four selected AA-MD snapshots. We found nine low-frequency normal modes (excluding translation and rotation) that can be grouped into three sets (or types of motions), each of which contain three modes (see Supplementary Fig. S2 ). On the whole, the first and second sets represent out-of-plane and in-plane motion of the catalytic domains, respectively, whereas the third set accounts for a “swinging” motion of the WAT domain (Fig. 3). The three modes within a set represent essentially the same motion in orthogonal directions with barely separated eigenvalues. We found the same type of motions in all four snapshots except for switches in the order of modes within the same set. This is to be expected because modes within a set are almost degenerate, with the three types of large-scale motions remaining as the most relevant fluctuations. These motions are similar to those reported earlier that were based on BNM analysis of the model-built structure 12. Most interestingly, the first three principal components computed from the CG-BD trajectories closely match the results from the BNM analysis (Fig. 3). These results clearly show that there are few low-frequency motions involving translation, rotation, and bending deformations that account for the major conformational changes in AChEt. Furthermore, similarity of the BNM results among different AA-MD snapshots and between the AA-MD relaxed structures and the initial model indicates that the motions are robust with regard to fine structural details and thus unaffected by deficiencies potentially present in the initial model of the AChEt. This is further reinforced by the similarity between the BNM and the PCA results, which also validate the capability of our CG model to capture the relevant large-scale motions.

Display large version of this figure

Figure 3

BNM (left) and PCA (right). The BNMs were calculated for all three snapshots used for starting the CG-BD runs as well as for the last snapshot; a representative result from a snapshot at 7.72ns, which was used to start run II, is shown here. Similarly, the PCA shown here is only for run II. The rest of the BNM and PCA calculations are very similar. (A–C) Representative modes from the first, second, and third sets of lowest frequency BNMs (see text) and the first, second, and third PCs. The blue arrows, drawn every four residue for clarity, indicate the direction of motion. For each mode (i.e., for those in A, B, and C), the same orientation is used and for different modes (i.e., A versus B, etc.), structures are rotated for the best view.

The time evolutions of the intersubunit distance (d), angle (θ), and interfacial contact (χ) in the course of the AA and CG simulations are shown in Figure 4AC. The plots show that first, the subunit arrangement is asymmetric and there is indeed a substantial intersubunit movement. Second, both the AA and CG simulations sample a similar range of distances (Figure 4A) and angles (Figure 4B). Third, compared with what was accessible by the short AA-MD run, the fast and highly simplified CG runs sampled various arrangements of subunits, achieving equilibrium fluctuations in which a conformational space is reversibly revisited. For example, d between subunits A and B (A-B) was >65Å in the first 10ns of AA-MD but subsequently dropped to <55Å, whereas subunits B and C kept moving apart. In contrast, the CG-BD run sampled the whole range of distances, roughly between 52 and 65Å for both A-B and B-C. An even larger conformational space was sampled by running several CG-BD simulations. Similar information is obtained using the variable θ, where the asymmetry and dynamics of the subunit arrangement is clearly seen in Figure 4B. The subunits rotate substantially relative to one another, such that a significant proportion of the conformations adopt a nonplanar subunit arrangement. As mentioned before, both (pseudosquare) planar and nonplanar subunit arrangements have been observed by x-ray crystallography 6. It is therefore instructive to quantify the extent of nonplanarity in the simulated structures. To characterize deviation from planarity, the difference between the sum of all four angles in AChEt and 360° has been calculated. Fig. 5 displays histograms of the deviation from planarity. For comparison, the nonplanarity calculated for the nonplanar crystallographic structure (1C2O) is 3°, which is shown in the figure as a magenta dashed line. We find that 27%, 15%, and 26% of all the collected snapshots from the CG-BD runs I, II, and II, respectively, are even more nonplanar than 1C2O. Using a more relaxed cutoff of 1°, over 50% (i.e., 70%, 51%, and 60%) of the structures are nonplanar. These data clearly show that AChEt fluctuates between planar and nonplanar subunit arrangements and that x-ray structures 1C2O and 12CB represent individual snapshots of the available conformations in solution.

Display large version of this figure

Figure 4

Distance (d), angle (θ), number of interfacial contacts (χ), and occlusion involving the catalytic domains (residues 5–543). (A) d between the center-of-mass of two subunits: A-B (black), B-C (red), C-D (green), and D-A (blue). (B) θ subtended by the centers of three subunits: (black), (red), (green), and (blue). (C) χ defined as the mean of the numbers of atoms in subunit i within 10Å of any atom of subunit j and vice versa: A-B (black), B-C (red), C-D (green), and D-A (blue). (D) Occlusion defined as the number of atoms in a complementary subunit j within 16Å of the PAS residues 72, 286, and 341 of subunit i: A (black), B (red), C (green), and D (blue).

Display large version of this figure

Figure 5

Histograms of the deviations from planarity, defined as , where θ is as defined in the legend of Fig. 4 and k is an index for the angles. The magenta dashed line represents the nonplanarity of the x-ray structure 1C2O.

In addition to studying the fluctuations of the catalytic domains in terms of a rigid-body motion, it is important to investigate how the local contacts at the subunit interfaces (χ) behave. There are a variety of ways to quantify the interfacial contacts. We used the variable χ as defined in Methods. Furthermore, since the value of χ and its ability to exclude buried atoms depends on the arbitrary cutoff used, we checked the change in SASA upon the association of two subunits. Essentially the same results were obtained (see Supplementary Fig. S2 ). Therefore, with the current cutoff (10Å), χ is a good measure of interfacial contact; its time evolution is shown in Figure 4C. The plot shows substantial fluctuation in both the AA and CG simulations. Occasionally, the contact between subunits B and C is lost. In fact, several CG runs indicate that any two pairs of subunits tend to have tighter interaction (higher χ) than the remaining two pairs. Furthermore, as in the case of the distances and the angles, the number of contacts between subunits varies from pair to pair. For example, all three variables (d, θ, χ) display higher fluctuation when subunit B is involved in the measurement. Taken together, these data demonstrate that AChEt is significantly dynamic due to the weak interaction between the catalytic subunits that allows translational and rotational motions of the subunits relative to one another.

The major goal of this work is to elucidate the nature and extent of active site occlusion by neighboring subunits. To monitor active site occlusion due to a complementary subunit, it is important to identify which of the gorge residues are affected upon tetramer formation. Buried residues, such as those lining the gorge, are inaccessible to residues in an opposing subunit. Although these residues experience high-frequency internal fluctuations that may play a role in accommodating a bound substrate, such fluctuations would not affect occlusion as defined here. Among the surface residues that play a critical role in substrate capture and binding, Tyr⁷², Trp²⁸⁶, and Tyr³⁴¹ at the PAS were found to consistently become blocked by residues from an opposing subunit 41. Therefore, contacts to these three residues by any residue of an opposing subunit can be a good estimate of the degree of occlusion. We therefore defined occlusion as the number of atoms in an opposing subunit that are within 16Å of the three PAS residues of a particular subunit (see Methods). The results are shown in Figure 4D. In these plots, the higher the number, the more occluded the active site would be. We can see that the peripheral sites of subunits A and D become frequently occluded by subunits D and C, respectively (Figure 4D). Further, although the fractions of occluded PAS vary among the subunits, the active sites of subunits B and C are accessible more frequently than those of subunits A and D. Similar results were obtained in all three CG-BD runs.

Overall, the presented data clearly show that the intersubunit fluctuations are intrinsic properties of the tetramer. Such a dynamic assembly of the subunits permits individual subunits to experience events of occlusion and opening the extent and rate of which varies among the subunits. These issues will be discussed further in the following sections.

The large-scale intersubunit motions cannot be fully captured by nanoseconds-long AA-MD simulations; the size of the system (∼2400 residues) prohibits sufficient sampling. We therefore decided to approach the problem by combining AA-MD and Cα-based CG-BD simulations. There are two advantages to this approach: we can i), improve the parameterization of the CG model using data from the AA simulation, and ii), exploit the extensive sampling power of the CG model for investigating the long timescale dynamics of AChEt. To be able to handle the very large and complicated structure of the AChEt, we made a number of modifications on the CG model with the help of results from the AA-MD (see Methods). First, the AA-MD provided the parameters for the harmonic constraints between the unusually tightly bound PRAD-WAT helices. This modification helped to create a stable structure of the PRAD-WAT region without affecting the structure and dynamics of the globular catalytic subunits. Second, due to the lack of electrostatic screening by the high dielectric water, Arg⁵⁵¹ of the WAT domain becomes locked in interaction with a surface residue Asp³¹⁰ in the catalytic domain. This was prevented by applying a Morse potential between the two residues based on their average separation in the AA simulation. The computational efficiency of the CG model allowed several CG-BD runs to be performed using several AA-MD-relaxed structures as starting conformations, thus reducing dependency on initial conditions.

After the required parameterizations, the CG-BD simulations both captured the large-scale motions (Figure 3 and Figure 4 and Figure 5) and, as illustrated in Figure 4D, allowed an extensive sampling of conformational space that was inaccessible to the AA-MD; important regions were visited repeatedly. For example, the AA simulation achieved only one and two events of PAS occlusion at subunits B and D, respectively, but none at subunits A and C. In the BD simulations, many events of occlusion and opening, most prominently in subunit D, have been achieved. The latter enables analysis of some of the equilibrium and kinetic properties of the tetramer (discussed below). Therefore, the multiscale simulations provided information that went beyond the capabilities of either the AA or CG model alone. Yet, despite their profoundly different spatial and temporal scales, the two models complement one another, as evidenced by the above examples.

Furthermore, we calculated the normal modes for selected snapshots from the AA-MD simulation. We found nine low-frequency modes (Fig. 3) which represent different flavors of three major motions. The first three represent rotation-dominated “shearing” of the catalytic domains relative to one another. The second set of modes represents a largely translational in-plane motion of the catalytic domains. The third set accounts for the “swinging” motion of the WAT helices. To investigate whether the major movements are altered during the simulation, we compared the BNM obtained from several snapshots. We found that except for small exchanges in the eigenvalues within a given set, there are no significant changes in the type of motion. To further check this issue, we computed the principal components from the CG-BD simulations (Fig. 3). As expected, the low-resolution CG model was not able to resolve the small differences observed within a given set. Remarkably, however, the three major motions seen in the BNM are faithfully reproduced. Therefore, the large-scale motions involving subunits of the tetramer, which were implied by AA-MD and captured by both the BNM and CG-BD simulations, are intrinsic properties of the enzyme.

Having discussed the reliability of the CG model and the initial structure, as well as the conformational sampling afforded by the multiscale simulations, we would like to investigate the structures to determine which of the motions, localized high-frequency interfacial fluctuations or large-scale domain motions, are more important for the function of the tetrameric enzyme. Inhibitors that bind to the PAS at the entrance to the active site gorge 42 influence catalysis by blocking entrance of ligands and removal of products 1,5,43,44,45,46,47. In the tetramer, a similar blockade could come from a nearby subunit. For example, in the mouse tetramer, which is compact and pseudosquare planar 5, the Gly²⁶⁰–Gly²⁶⁴ cluster of hydrophobic residues at the tip of the short Ω-loop 39 pack against the PAS residues Tyr⁷², Trp²⁸⁶, and Tyr³⁴¹, sterically occluding access to the active center gorge. Similarly, in the compact eAChE tetramer (1C2O), two of the gorges from diagonally opposed subunits face the tetramer internal space, with their PAS occluded. In contrast, the loose eAChE tetramer (1C2B), despite its resemblance to the mAChE tetramer in terms of the pseudosquare planar arrangement of the subunits, has all four active centers fully accessible to solvent. Therefore, occlusion of the PAS is dependent on the translations and rotations of the subunits relative to one another but can also be modulated by local conformations such as those of the Ω-loop or other surface loops. Therefore, to determine whether or not the global intersubunit fluctuations are more important for the occlusion of the PAS than the local interfacial motions, we plotted the occlusion against the intersubunit distance (Fig. 6, main plot) and the interfacial contacts (Fig. 6, inset). There is a reasonably strong anticorrelation between occlusion and distance (correlation coefficient, R=−0.78) but none between occlusion and interfacial contact. There was no correlation with the interfacial surface area either (ΔSASA, see Supplementary Fig. S3 ). Therefore, the simulations suggest that occlusion is a function of intersubunit fluctuations and that the rate of blockade and opening of the route to the active center is regulated mostly by the motions of the subunits relative to one another.

Display large version of this figure

Figure 6

Plots of occlusion against the intersubunit distance (main plot) and interfacial contact (inset) for the A-D pair. The correlation coefficients are indicated. A similar behavior was obtained for the other subunit pairs.

The tetramer is the major functional form of AChE in cholinergic synapses, though tetramerization may affect ligand accessibility to the peripheral site 2,48,49. Our simulations show that the subunits are highly dynamic and that alternative opening and occlusion of an active site is allowable. This would imply that instead of being trapped in an occluded active site, for which ligand binding is limited, the subunits fluctuate relative to each other to expose all four active sites.

The accessibility of an active site may be divided into three states based on the probability distribution of the active site occlusion (Figure 7A). The first of these, the open state, has an occlusion number <2. The second, which has an occlusion number between 2 and 5, is considered partially open. The third one, the occluded state, covers an occlusion number >5. For example, in the crystal structure 1C2O, two of the four active sites are in the open and two are in the occluded state (occlusion numbers are both 9) whereas in 1C2B, all four active sites are open (occlusion numbers are 0). This definition is reasonable, as demonstrated by ACh-protein docking. Fig. 8 illustrates the docking results using the Vdock package 50,51. ACh was always docked near the entrance to the active site when AChE was in the open state. In the partially open state, ACh was placed mostly near the entrance, but the molecule could also be trapped by the complementary subunit. ACh can still reach the PAS when an active site is occluded, but it is very likely to interact and stay with the complementary subunit (Fig. 8). As a result, the association may be reduced by occlusion. Although it may be more efficient for catalysis to keep the four active sites always in the open state, the (weak) intersubunit interactions and thermal fluctuations bring two subunits close enough to hinder ligand binding, as shown in Fig. 8.

Display large version of this figure

Figure 7

Statistical analysis of the PAS occlusion. (A) Probability distribution of the occlusion in one of the CG-simulations. Vertical dashed lines indicate the cutoffs used to define three states represented by the snapshots having an occlusion less than or equal to 2 (state 1), between 2 and 5 (state 2), or greater than 5 (state 3). (B) The fraction of structures in state 1 (square) and state 2 (circle) for subunits A, B, C, and D derived from each of the four CG-BD runs.

Display large version of this figure

Figure 8

Docking of ACh into the PAS of AChE. (A) Docking results. Two hundred ACh structures were docked into three of the four PASs representing the open (occlusion=0), partially occluded (occlusion=2), and occluded (occlusion=6) states. See text and legends of Figure 1 and Figure 5 for the definitions of PAS and occlusion. (B) A zoom-in on the PAS residues and the ACh structures docked onto PAS.

Figure 7B shows partially open (i.e., partially occluded) and fully open fractions computed from four BD runs. The fractions of fully open conformations vary between subunits, and the average value is ∼52%. This suggests that, on average, two open states may be observed in the tetramer. However, in contrast to the crystal structures where only subunits B and D are in the open state, our BD simulations reveal that all four active sites can adopt three different conformational states, albeit with varying populations. Such dynamics may provide more symmetric binding sites for the tetramer, as compared with the axially asymmetric binding sites in 1C2O. It is worth noting that most conformations generated by the simulations are not very close to structures 1C2O or 1C2B. Zhang et al. generated an intermediate structure INT by morphing between the two crystal structures. Their calculation suggested that reaction rates for INT and 1C2B are similar, whereas the more compact 1C2O has slower rates 7,8. As a result, although all four active sites may not be in their open state simultaneously, the conformational dynamics may enable the tetramer to retain rates similar to those of 1C2B.

Finally, it is of interest to investigate if the binding of ACh is “gated” by intersubunit motions. Eq. (2) is applied here; thus we classify the partially and fully open states as an open gate, whereas we consider the occluded state a closed gate. The calculated average dwell times of the open gate and closed gate are ∼300ns and ∼50ns, respectively. The opening and closing rate constants k_o and k_c of this system are thus 1/50 and 1/300ns⁻¹. The computed (k_o+k_c)⁻¹ is a few orders of magnitude larger than the diffusional relaxation time (see Eq. (5)) if we assume that the protein and ACh have approximate radii of 65Å and 3Å, respectively. Thus, according to Eq. (3), the intersubunit motions may be viewed as “slow gating”, and the gated association rate constant may therefore be estimated by the ungated rate constant multiplied by the opening fraction. The approximate gated rate constant here is ∼85% of that for an ungated active site. We may also view the number as though one of the four active sites is closed, on average, during the approach of ACh to the tetramer. This result is consistent with the experimental finding that the binding of a substrate to the PAS is only somewhat hindered by the association of the subunits 49, resulting in a slightly higher substrate inhibition constant 48 in the tetramer than in the monomer 49. Note that our method provides reasonable approximations in subunit motions. However, since the CG models do not have explicit water molecules, the timescales exhibited here might not correspond exactly to those in cells. It is also important to underline that the arrangement of the subunits in the tetramer is asymmetric, as illustrated by a comparison between the initial nearly symmetric model (Figure 1A) and the MD relaxed structure (Figure 8A). Note also that the effect of electrostatics, which was shown to enhance the rate of binding 7, is not considered here. Only the geometry of the tetramer is discussed as having a gating effect. Taken together, the catalytic activity of each subunit might not be very adversely affected by the tetramer formation, as has also been suggested by experiments 49. Moreover, the tetramer provides more active sites that, for an approaching ACh, present at least three times greater a chance for reaction than an isolated monomer.

In this work, we presented results from multiscale simulations that allowed us to directly connect long timescale domain motions with the active site substrate accessibility of the enzyme AChE. The in-plane translation and out-of-plane rotational motions of the four catalytic domains lead to dynamic occlusions of the peripheral cationic site. These motions gate substrate-protein association and allow the active center to be accessible for more than 80% of the time on average. A combination of methods involving atomistic and CG simulations as well as normal mode analyses was used to establish the link between domain motion and substrate accessibility. This approach is likely to be useful for future investigations of biomolecular systems. The results also lay the foundation for further study on the role of tetramerization in the catalytic efficiency of AChE, which may include rate calculations on several structures representing different active site occlusions.