|Home | About | Journals | Submit | Contact Us | Français|
The different coordination modes and fast ligand exchange of zinc coordination has been suggested to be one key catalytic feature of the zinc ion which makes it an invaluable metal in biological catalysis. However, partly due to the well known difficulties for zinc to be characterized by spectroscopy methods, evidence for dynamic nature of the catalytic zinc coordination has so far mainly been indirect. In this work, Born-Oppenheimer ab initio QM/MM molecular dynamics simulation has been employed, which allows for a first-principle description of the dynamics of the metal active site while properly including effects of the heterogeneous and fluctuating protein environment. Our simulations have provided direct evidence regarding inherent flexibility of the catalytic zinc coordination shell in Thermolysin (TLN) and Histone Deacetylase 8 (HDAC8). We have observed different coordination modes and fast ligand exchange during the picosecond's time-scale. For TLN, the coordination of the carboxylate group of Glu166 to Zinc is found to continuously change between monodentate and bidentate manner dynamically; while for HDAC8, the flexibility mainly comes from the coordination to a non-amino-acid ligand. Such distinct dynamics in the zinc coordination shell between two enzymes suggests that the catalytic role of Zinc in TLN and HDAC8 is likely to be different in spite of the fact that both catalyze the hydrolysis of amide bond. Meanwhile, considering that such Born-Oppenheimer ab initio QM/MM MD simulations are very much desired but are widely considered to be too computationally expensive to be feasible, our current study demonstrates the viability and powerfulness of this state-of-the-art approach in simulating metalloenzymes.
Zinc is relatively abundant in biological materials. Approximately 10% of the total human proteome have been identified to bind with zinc in vivo from a bioinformatics investigation1 and they play very crucial roles in all forms of life2-6. For mononuclear zinc enzymes, a typical metal coordination environment contains three amino acid side chain ligands (His, Glu, Asp and Cys) and one/two small molecule(s).3,7,8 The flexibility of zinc coordination, which allows different coordination modes and fast ligand exchange, has been suggested to be one key catalytic feature of the zinc ion which makes it an invaluable metal in biological catalysis.9 However, partly due to the well known difficulties for zinc to be characterized by spectroscopy methods10,11, evidence for dynamic nature of the catalytic zinc coordination has so far mainly been indirect, coming from determined X-ray crystal structures 5,12 and geometry optimizations of model complexes with electronic structure methods13,14.
A classical example of the flexible zinc coordination is thermolysin (TLN), one of the best experimentally characterized zinc proteases with one Glu and two His amino-acid-side-chains as ligands.15 Extensive structural studies have indicated that besides the simple tetrahedron coordination, its catalytic zinc ion can also be five-, or six-coordinated,16-19 especially for the carboxylate of Glu/Asp whose coordination is highly flexible. The different coordination mode of carboxylate to zinc with a continuous range between monodentate and bidentate manner has been observed in model complexes as well as other zinc enzymes with one Glu or Asp as ligands.12 It has been suggested that the coordination of carboxylate to zinc could be dynamic in the catalytic process, known as carboxylate shift, which has been suggested to be important in the function of zinc enzymes.13,14,20-22 Intriguingly, for zinc-dependent histone deacetylases (HDACs), in which the catalytic zinc ion is coordinated to one His and two Asp residues, only the monodentate mode has been observed for the carboxylate coordination in several crystal structures.23-26 It should be noted that HDACs, which catalyze the removal of acetyl group from histone tails, have been emerged to be critical in gene regulation and are among the most promising targets for the development of anti-tumor therapeutics. For example, the recently FDA approved anti-cancer drug SAHA is a HDAC inhibitor which is directly coordinated to the zinc27,28.
In order to provide deep insights into the dynamics and flexibility of the zinc catalytic site, which would be essential in characterizing their catalytic mechanisms and rational design of novel inhibitors for zinc enzymes, we have carried out DFT QM/MM Born-Oppenheimer molecular dynamics (BOMD) simulations on TLN and HDAC8. Although semi-empirical QM/MM BOMD simulations of some zinc-dependent enzymes have been carried out29-32, one main concern is the accuracy and transferability of the semi-empirical QM Hamiltonian in describing the Zinc coordination shell and the enzyme reaction profile. Often some parameters of the semi-empirical Hamiltonian need to be re-optimized for specific systems and reactions to obtain reasonable results30,32-35, which significantly reduce its transferability and predictability. Thus, in our current DFT QM/MM molecular dynamics simulations, the zinc ion and its ligands are treated by B3LYP hybrid functional with a Stuttgart ECP/Basis set36 for the zinc atom and a 6-31G* basis set for all other QM atoms. Our employed theoretical approach, which has recently been demonstrated to be feasible and powerful in several enzymatic studies37-40, allows for a first-principle description of the dynamics of the metal active site while properly including effects of the heterogeneous and fluctuating protein environment.
As illustrated in Fig. 1 and listed in Table 1 and and2,2, four TLN and two HDAC8 simulation systems have been prepared base on their respective experiment crystal structures: TLNa, a substrate-free structure of TLN (pdb code: 1LNF16); TLNb: an TLN-inhibitor complex which mimics the reactant state (pdb code: 1ZDP17); TLNc: an TLN-inhibitor complex which mimic the transition state (pdb code: 2TMN18); TLNd: an TLN-product complex (pdb code: 3TMN19); model 1: a complex of the Y306F mutant of HDAC8 and its substrate (pdb code: 2V5W25); model 2: an HDAC8-SAHA complex (pdb code: 1T6924). Representative active site structures of TLN and HDAC8 enzymes were illustrated in Figure 2.
For each simulation system, its initial structure was prepared based on their respective crystal structure with molecular modeling and dynamics simulations. First, missing residues were added with the Swiss-PdbViewer41. Protonation states of charged residues in the enzyme complex were determined via H++ program42 and by carefully examining local hydrogen bond networks. In all TLN models, His142 and His146 was identified as singly protonated, whereas His231 was doubly protonated, which are the same as previous studies43,44. For HDAC8 models, His180 was determined as singly protonated on δ site, while His142 and His143 were determined as singly protonated on ε site. Then each prepared system was solvated into a rectangular box with a 10 Å buffer distance between the solvent box wall and the nearest solute atom. Finally, in order to neutralize each simulation system, one to six sodium ions were added at the protein surface respectively by employing the Amber Tool. All sodium ions were located more than 17 Å away from the Zinc active site. The resulted simulation system was about 45000 atoms for each model. The Zn2+ ion was modeled with the Stote's scheme45. Considering that zinc coordination shell is very difficult to be well described by a molecular mechanical force field46-49, ~1500 kcal/mol harmonic restraint was placed on ~15 atoms in the zinc bind site to retain the zinc coordination structure during the MM equilibrations. The rest protein and solvent molecules were first minimized and then more than 3 ns MD simulations were carried out for each system with periodic boundary condition. A time step of 2 fs was used. Berendsen thermostat method50 was used to control the system temperature at 300K. All MD simulations were performed with AMBER10 molecular dynamics package51. The Amber99SB52-54 force field for the protein and TIP3P model55 for water molecules were employed, and the force field parameters for substrate in HDAC8 and inhibitors in these simulation models were generated from AMBER GAFF force field56 via AMBER tools. The SHAKE algorithm57 was applied to constrain all bonds involving hydrogen atoms with tolerance of 10-5 Å, and 12 Å cutoff were used for both van der Waals and electronic interactions.
Considering that the trajectory was very stable after 2 ns for all these classical MD simulations, the resulted snapshot after 3ns MD simulation was employed as the initial structure for the preparation of ab initio QM/MM MD simulations. Each QM/MM model was prepared by deleting the solvent molecules beyond 30 Å from the zinc atom. The resulting QM/MM system had totally ~13000 atoms. The detailed QM/MM partition for all these TLN and HDAC8 models were presented in Figure 1S. The QM subsystem, including the zinc and its coordinating ligands, was treated with B3LYP functional using Stuttgart ECP/basis set36 for the zinc atom and 6-31G(d) basis set for all other atoms, which has been previously tested and employed successfully to describe zinc coordination shell 14,21,58-60. The QM/MM boundary were described by improved pseudobond approach61-64. All other atoms were described by the same molecular mechanical force field used in classical MD simulations. The spherical boundary condition had been applied so that atoms beyond 22 Å from the zinc atom were fixed. The 18 and 12 Å cutoffs were employed for electrostatic and van der Waals interactions, respectively. There was no cutoff for electrostatic interactions between QM and MM regions. The prepared system was first minimized by QM/MM calculations. Finally, more than 20 ps ab initio QM/MM MD simulations were carried out with 1 fs as the time step, and the Beeman algorithm65 was used to integrate the Newton equations of motion. In order to further check the convergence of our simulations we have also extended the ab initio QM/MM MD simulations to 40 ps for the TLNa system. As shown in Table 1S shown, the difference of average distances and fluctuation for the zinc coordination shell in TLNa are very similar for different time periods. The Berendsen thermostat method50 was used to control the system temperature at 300 K. The last 20 ps trajectory which has achieved the equilibrated temperature (300K) is adopted as data analysis. In order to make sure that the fluctuation of Zinc coordination shell is not due to the higher temperature of the QM sub-system, we have monitored it along our QM/MM MD simulation trajectories. We found that most of the time the temperature of the QM subsystem is not higher than 300 K, and the fluctuation is reasonable. All ab initio QM/MM calculations were performed with modified Q-Chem66 and Tinker67 programs.
For TLN, we simulated four enzyme complexes of TLN with different coordinating ligands starting from their respective crystal structures, as illustrated in Figure 1(a). In the model TLNa, in which no substrate is bound to the active site, the zinc ion is six-coordinated with a bidentate carboxyl group and two water molecules in the crystal structure. From Table 1, we can see that our MD simulation results not only reproduced the coordination shell very well but also clearly indicated the flexibility and dynamics of zinc catalytic site. The coordination between Zn2+ and Glu166:OE1 is demonstrated to be most dynamic, and the zinc-carboxylate coordination fluctuates between mono- and bi-dentate manners as shown in Figure 3 and and4.4. Interestingly, in the model TLNb, in which the inhibitor is bound to Zn2+ with a thiolate, its zinc coordination is much less flexible which is manifested by the small standard deviation (SD), and is kept to be tetrahedral as in the crystal structure most of the time (See Figure 3). Meanwhile, in models TLNc and TLNd, which have been suggested to mimic the reaction transition state and product respectively, their zinc active sites are found to be even more flexible than that in model TLNa. For both models, the zinc-carboxylate coordination is the most flexible, and continuously changes between mono-dentate and bi-dentate manner dynamically. We can also see that the ligand exchange in zinc coordination is very fast, and can occur at the picoseconds scale, as illustrated in Figure. 4.
For HDAC8, a key distinct feature of its active site is that there are two carboxylate groups coordinated to zinc instead of only one carboxylate in TLN. From structural studies, only mono-dentate mode of carboxylate coordination has been observed.23-26 Considering the medical and biological significance of HDACs, it would be of particular interest to probe the dynamics of its zinc active site. Here we have simulated two HDAC8 complexes starting from their respective crystal structures: model 1 is a Y306F mutant in complex with its substrate acetyl-lysine,25 and model 2 is the wild-type HDAC8 binding with its superstar inhibitor SAHA,24 as shown in Figure 1(b). From Table 2, we can see that the average distances are consistent with the experimental data, and fluctuations indicate the flexibility of the zinc coordination sphere. Although the zinc active sites in both models are still quite dynamic, the carboxylate coordination is very stable, which is manifested by the small fluctuation, and always kept in the mono-dentate mode. As shown in Figure 3 and and4,4, in the mutant-substrate complex (model 1), its zinc coordination number fluctuates between 4 and 5 and its main flexibility comes from the coordination of the water molecule. In the HDAC8-SAHA complex, its five-fold zinc coordination is considerably more stable. SAHA is coordinated to the zinc in the bidentate manner most of the time, but it remains to be more flexible than other amino acid ligands, which is indicated by the relative large standard deviation in Table 2. Thus our simulations clearly indicate that the dynamics of zinc active site of HDAC8 is quite different from thermolysin, and suggest that HDAC8 is not likely to employ the carboxylate-shift mechanism in its catalytic cycle.
As noted in Figure 1, the first coordination shell in TLN is “2His+Glu”, for HDAC8 is “1His+2Asp”, so the total charge for the coordination shell is different. This would be an important factor which leads to their different flexibility of zinc-binding sites. Such distinct dynamics in the zinc coordination shell also suggests that the catalytic role of zinc in TLN and HDAC8 is likely to be different in spite of the fact that both catalyze the hydrolysis of amide bond.
In summary, with ab initio QM/MM MD simulations, we have provided direct evidence regarding the inherent flexibility of the catalytic zinc coordination in both TLN and HDAC8, and have observed different coordination modes and fast ligand exchange. For TLN, the coordination of the carboxylate group of Glu166 is found to be most flexible, which can continuously change between monodentate and bidentate manner dynamically. While for HDAC8, its coordination to all three amino acid ligands, including two carboxylate groups of Asp, is very stable. Its flexibility mainly comes from a non-amino-acid ligand. Such distinct dynamics in their zinc coordination shell suggest that the catalytic role of zinc in TLN and HDAC8 is likely to be different in spite of the fact that both catalyze the hydrolysis of amide bond. Furthermore, our work demonstrates the feasibility and applicability of Born-Oppenheimer ab initio QM/MM MD simulations in simulating metalloenzymes, and set the stage for more detailed understanding of catalysis and inhibition in TLN and HDAC8.
This work was supported by from NIH (R01-GM079223), NSF (CHE-CAREER-0448156) and the China Scholarship Council. We thank NYU-ITS and NCSA for providing computational resources.