The development of novel isoform/class-selective inhibitors is still of great biological and medical significance to conquer the continuously reported side effects for the histone deacetylase (HDAC) drugs. The first potent HDAC allosteric inhibitor was discovered last year, and this allosteric inhibitor design is thought to be a promising strategy to overcome the current challenges in HDAC inhibitor design. However, the detailed allosteric mechanism and its remote regulatory effects on the catalytic/inhibitor activity of HDAC are still unclear. In this work, on the basis of microsecond-time-scale all-atom molecular dynamics (MD) simulations and picosecond-time-scale density functional theory/molecular mechanics MD simulations on HDAC8, we propose that the allostery is achieved by the intrinsic conformational flexibility of the binding rail (constituted by a highly conserved X–D residue dyad), which steers the loop–loop motion and creates the diverse shapes of the allosteric sites in different HDAC isoforms. Additionally, the rotatability of the binding rail is an inherent structural feature that regulates the hydrophobicity of the linker binding channel and thus further affects the HDAC enzyme inhibitory/catalytic activity by utilizing the promiscuity of X–D dyad. Since the plastic X residue is different among class I HDACs, these new findings provide a deeper understanding of the allostery, which is guidable for the design of new allosteric inhibitors toward the allosteric site and structure modifications on the conventional inhibitors binding into the active pocket by exploiting the intrinsic dynamic features of the conserved X–D dyad.

The bromodomain and extra terminal domain (BET) family of bromodomains (BRDs) are well-known drug targets for many human diseases. The active pockets of the two tandem bromodomains BD1/BD2 are highly conserved (sequence similarity is about 95%), thus it is of great medical importance and still a significant challenge to develop BD1/BD2 selective inhibitors. A few BD2 selective inhibitors, such as RVX-208 and RVX-297, have been reported recently. However, their selectivity is insufficient for drug discovery, and the molecular basis of the selective inhibition for BD2 over BD1 remains unknown. In this work, by extensive classical molecular dynamics (MD) simulations and hybrid density functional theory/molecular mechanics (DFT/MM) MD simulations, it is for the first time revealed that the selective inhibitory effect towards BD2 is achieved by the distinctive structural dynamics of the ZA-loop (“in/out” conformations) in BD1 and BD2, which originate from the existence of residue Asp144 in BD1 which is replaced by His433 in BD2. Additionally, the more stable inherent H-bond constructed by a conserved D–Y dyad, as well as the stronger π–π stacking interaction formed between His433 and the ligand, are responsible for the higher inhibitory activity of RVX-297 compared to that of RVX-208 in BD2. All these findings should guide further novel inhibitor design and structural modification of validated BD1/BD2 inhibitors to increase the selectivity for BD1/BD2 among the BET family.

The cleavage of the magnesium-assisted diphosphate group (the PPi group) is one significant and prevalent rate-limiting step triggering the enzyme catalysis synthesis of terpenoid natural products. However, the PPi cleavage procedure has been rarely studied in most theoretical research of the terpenoid biosynthetic mechanism. In this work, QM(DFT)/MM MD simulations were employed to illuminate the detailed PPi cleavage mechanism in three different enzyme systems (ATAS, TEAS, and FPPS). We found that the most rational protonation state of the PPi group is highly dependent on the Mg2+ coordination modes and the enzyme classes. The deprotonation of PPi is favorable for triggering the catalysis reaction in ATAS, while monoprotonation in FPPS and biprotonation in TEAS are advantageous. As a result, similar PPi cleavage occurs by means of nucleophilic substitution reactions in TEAS/FPPS/ATAS but presents an SN1, SN2, and borderline mechanism, respectively. Finally, the alternative functions of PPi protonation and Mg2+ coordination modes are discussed.Keywords: diphosphate (PPi) cleavage; farnesyl diphosphate cyclase (FPPC); farnesyl diphosphate synthase (FPPS); protonation state; QM/MM

The Mg-dependent 5-epi-aristolochene synthase from Nicotiana tabacum (called TEAS) could catalyze the linear farnesyl pyrophosphate (FPP) substrate to form bicyclic hydrocarbon 5-epi-aristolochene. The cyclization reaction mechanism of TEAS was proposed based on static crystal structures and quantum chemistry calculations in a few previous studies, but substrate FPP binding kinetics and protein conformational dynamics responsible for the enzymatic catalysis are still unclear. Herein, by elaborative and extensive molecular dynamics simulations, the loop conformation change and several crucial residues promoting the cyclization reaction in TEAS are elucidated. It is found that the unusual noncatalytic NH2-terminal domain is essential to stabilize Helix-K and the adjoining J-K loop of the catalytic COOH-terminal domain. It is also illuminated that the induce-fit J-K/A-C loop dynamics is triggered by Y527 and the optimum substrate binding mode in a “U-shape” conformation. The U-shaped ligand binding pose is maintained well with the cooperative interaction of the three Mg2+-containing coordination shell and conserved residue W273. Furthermore, the conserved Arg residue pair R264/R266 and aromatic residue pair Y527/W273, whose spatial orientations are also crucial to promote the closure of the active site to a hydrophobic pocket, as well as to form π-stacking interactions with the ligand, would facilitate the carbocation migration and electrophilic attack involving the catalytic reaction. Our investigation more convincingly proves the greater roles of the protein local conformational dynamics than do hints from the static crystal structure observations. Thus, these findings can act as a guide to new protein engineering strategies on diversifying the sesquiterpene products for drug discovery.

Histone Deacetylases (HDACs) are promising anticancer targets and several selective inhibitors have been created based on the architectural differences of foot-pockets among HDACs. However, the “gate-keeper” of foot-pockets is still controversial. Herein, it is for the first time revealed that a conserved R–E salt bridge plays a critical role in keeping foot-pockets closed in class-II HDACs by computational simulations. This finding is further substantiated by our mutagenesis experiments.

Geranyl pyrophosphate synthase (GPPS) is responsible for the formation of geranyl pyrophosphate (GPP), a key intermediate which has the potential to derive numerous functionally and structurally diverse groups of terpenoid natural products via the head-to-tail assembly of two isoprenoid building blocks (dimethylallyl diphosphate, DMAPP; isopentenyl diphosphate, IPP) in the initial step of carbon-chain elongation during isoprenoid biosynthesis. Elucidating the detailed catalytic mechanism in GPPS is of significant interests as it will stimulate the development of new technology in generating novel natural productlike scaffolds. It has been known that the catalytic reaction involves three sequential steps, namely “ionization–condensation–elimination”, but the exact catalytic mechanism has remained controversial since the 1970s. By employing Born–Oppenheimer density functional quantum mechanics (B3LYP/6-31(+)G*)/molecular mechanics dynamics simulations, here we suggest that GPPS adopts a protonation-induced catalytic mechanism, in which there are two key points different from previously hypothesized mechanisms. The first point is the protonation of DMAPP which is essential in the initial “ionization” process but was not considered in previous mechanisms. The second point is the stereoselectivity of proton transfer (HS) from IPP to H76 residue in the final “elimination” step as identified in our simulations, in contrast to the proton transfer from IPP (HR) to DMAPP in previous hypotheses. Moreover, the free energy barrier of the whole assembly reaction is predicted to be 18.8 ± 0.6 kcal/mol, in agreement with the experimental value of 18.0 kcal/mol. Furthermore, the catalytic roles of the two Mg2+ ions at the bottom of the active site are also discussed, and key residues (K44, R47, R94, R95, K180, K235, and E73 around DMAPP and IPP) responsible for the stabilization of transition states, intermediates, and/or product are clarified. These findings can assist site-directed mutagenesis experiments in protein engineering as well as inhibitor designs.Keywords: assembly; catalytic mechanism; GPPS; isoprenoid biosynthesis; protonation

The conversion of 1-deoxy-D-xylulose 5-phosphate (DXP) to 2-C-methyl-D-erythritol 4-phosphate (MEP) catalyzed by DXP reductoisomerase (DXR) is the committing step in the biosynthesis of terpenoids. This MEP pathway is essential for most pathogenic bacteria but absent in human, and thus, it is an attractive target for the development of novel antibiotics. To this end, it is critical to elucidate the conversion mechanism and identify the transition state, as many drugs are transition-state analogues. Here we performed extensive combined quantum mechanical (density functional theory B3LYP/6-31G*) and molecular mechanical molecular dynamics simulations to elucidate the catalytic mechanism. Computations confirmed the transient existence of two metastable fragments of DXP by the heterolytic C3–C4 bond cleavage, namely, 1-propene-1,2-diol and glycoaldehyde phosphate, in accord with the most recent kinetic isotope effect (KIE) experiments. Significantly, the heterolytic C3–C4 bond cleavage and C2–C4 bond formation are accompanied by proton shuttles, which significantly lower their reaction barriers to only 8.2–6.0 kcal/mol, compared with the normal single carbon–carbon bond energy 83 kcal/mol. This mechanism thus opens a novel way for the design of catalysts in the cleavage or formation of aliphatic carbon–carbon bonds.Keywords: 1-deoxy-D-xylulose 5-phosphate; carbon−carbon bond cleavage; combined QM(DFT)/MM; proton shuttle; reductoisomerase

Development of isoform-selective histone deacetylase (HDAC) inhibitors is of great biological and medical interest. Among 11 zinc-dependent HDAC isoforms, it is particularly challenging to achieve isoform inhibition selectivity between HDAC1 and HDAC2 due to their very high structural similarities. In this work, by developing and applying a novel de novo reaction-mechanism-based inhibitor design strategy to exploit the reactivity difference, we have discovered the first HDAC2-selective inhibitor, β-hydroxymethyl chalcone. Our bioassay experiments show that this new compound has a unique time-dependent selective inhibition on HDAC2, leading to about 20-fold isoform-selectivity against HDAC1. Furthermore, our ab initio QM/MM molecular dynamics simulations, a state-of-the-art approach to study reactions in biological systems, have elucidated how the β-hydroxymethyl chalcone can achieve the distinct time-dependent inhibition toward HDAC2.

A full enzymatic catalysis cycle in the inosine–adenosine–guanosine specific nucleoside hydrolase (IAG-NH) was assumed to be comprised of four steps: substrate binding, chemical reaction, base release, and ribose release. Nevertheless, the mechanistic details for the rate-limiting step of the entire enzymatic reaction are still unknown, even though the ribose release was likely to be the most difficult stage. Based on state-of-the-art quantum mechanics and molecular mechanics (QM/MM) molecular dynamics (MD) simulations, the ribose release process can be divided into two steps: “ribose dissociation” and “ribose release”. The “ribose dissociation” includes “cleavage” and “exchange” stages, in which a metastable 6-fold intermediate will recover to an 8-fold coordination shell of Ca2+ as observed in apo- IAG-NH. Extensive random acceleration molecular dynamics and MD simulations have been employed to verify plausible release channels, and the estimated barrier for the rate-determining step of the entire reaction is 13.0 kcal/mol, which is comparable to the experimental value of 16.7 kcal/mol. Moreover, the gating mechanism arising from loop1 and loop2, as well as key residues around the active pocket, has been found to play an important role in manipulating the ribose release.

Bromodomains (BRDs) are protein modules that selectively recognize histones as a “reader” by binding to an acetylated lysine substrate. The human BRD4 has emerged as a promising drug target for a number of disease pathways, and several potent BRD inhibitors have been discovered experimentally recently. However, the detailed inhibition mechanism especially for the inhibitor binding kinetics is not clear. Herein, by employing classical molecular dynamics (MD) and state-of-the-art density functional QM/MM MD simulations, the dynamic characteristics of ZA-loop in BRD4 are revealed. And then the correlation between binding pocket size and ZA-loop motion is elucidated. Moreover, our simulations found that the compound (−)-JQ1 could be accommodated reasonably in thermodynamics whereas it is infeasible in binding kinetics against BRD4. Its racemate (+)-JQ1 proved to be both thermodynamically reasonable and kinetically achievable against BRD4, which could explain the previous experimental results that (+)-JQ1 shows a high inhibitory effect toward BRD4 (IC50 is 77 nM) while (−)-JQ1 is inactive (>10 μM). Furthermore, the L92/L94/Y97 in the ZA-loop and Asn140 in the BC-loop are identified to be critical residues in (+)-JQ1 binding/releasing kinetics. All these findings shed light on further selective inhibitor design toward BRD family, by exploiting the non-negligible ligand binding kinetics features and flexible ZA-loop motions of BRD, instead of only the static ligand–protein binding affinity.

SAHA (vorinostat, Merck) is a famous clinical drug for zinc-containing histone deacetylase (HDAC) targets against cancer and several other human disorders, whose inhibition mechanism (namely the protonation mechanism) upon binding to HDAC has been debated for more than ten years. It is very challenging to verify experimentally and is still controversial theoretically. The popular “Class-dependent” (namely “Tyr-dependent”) hypothesis is that the deprotonation of SAHA is mostly regulated by the conserved Tyr308 in class I HDAC while it is replaced by the His843 in class IIa HDAC. Herein, by elaborate QM(DFT)/MM MD simulations, we exclude the prevalent “Class-dependent” mechanism and advance a novel “Metal-dependent” mechanism, where the remote second metal site (K+ in most HDAC and Ca2+ in HDAC2) determines the protonation of SAHA. This proof-of-principle “Metal-dependent” mechanism opens up a new avenue to utilize the second metal site for isoform-selective inhibitor design.

Phosphodiesterases (PDEs) are the sole enzymes hydrolyzing the important second messengers cGMP and cAMP and have been identified as therapeutic targets for several diseases. The most successful examples are PDE5 inhibitors (i.e., sildenafil and tadalafil), which have been approved for the treatment of male erectile dysfunction and pulmonary hypertension. However, the side effects mostly due to nonselective inhibition toward other PDE isoforms, set back the clinical usage of PDE5 inhibitors. Until now, the exact catalytic mechanism of the substrate cGMP by PDE5 is still unclear. Herein, the first computational study on the catalytic hydrolysis mechanism of cGMP for PDE5 (catalytic domain) is performed by employing the state-of-the-art ab initio quantum mechanics/molecular mechanics (QM/MM) molecular dynamics (MD) simulations. Our simulations show a SN2 type reaction procedure via a highly dissociated transition state with a reaction barrier of 8.88 kcal/mol, which is quite different from the previously suggested hydrolysis mechanism of cAMP for PDE4. Furthermore, the subsequent ligand exchange and the release of the product GMP have also been investigated by binding energy analysis and MD simulations. It is deduced that ligand exchange would be the rate-determining step of the whole reaction, which is consistent with many previous experimental results. The obtained mechanistic insights should be valuable for not only the rational design of more specific inhibitors toward PDE5 but also understanding the general hydrolysis mechanism of cGMP-specific PDEs.

Human oxidosqualene cyclase (OSC) is one key enzyme in the biosynthesis of cholesterol. It can catalyze the linear-chain 2,3-oxidosqualene to form lanosterol, the tetracyclic (6–6–6–5 members for A–B–C–D rings) cholesterol precursor. It also has been treated as a novel antihyperlipidemia target. In addition, the structural diversity of cyclic terpenes in plants originates from the cyclization of 2,3-oxidosqualene. The enzyme catalytic mechanism is considered to be one of the most complicated ones in nature, and there are a lot of controversies about the mechanism in the past half a century. Herein, state-of-the-art ab initio QM/MM MD simulations are employed to investigate the detailed cyclization mechanism of C-ring and D-ring formation. Our study reveals that the C and D rings are formed near-synchronously from a stable “6–6–5” ring intermediate. Interestingly, the transition state of this concerted reaction presents a “6–6-6” structure motif, while this unstable “6–6-6” structure in our simulations is thought to be a stable intermediate state in most previous hypothetical mechanisms. Furthermore, as the tailed side chain of 2,3-oxidosqualene shows a β conformation while it is α conformation in lanosterol, finally, it is observed that the rotatable “tail” chain prefers to transfer β conformation to α conformation at the “6–6–5” intermediate state.

Terpenes (isoprenoids) represent the most functionally and structurally diverse group of natural products. Terpenes are assembled from two building blocks, isopentenyl diphosphate (IPP) and dimethylallyl diphosphate (DMAPP or DPP), by prenyltransferases (PTSs). Geranyl pyrophosphate synthase (GPPS) is the enzyme that assembles DPP and IPP in the first step of chain elongation during isoprenoid biosynthesis. The mechanism by which GPPS assembles the terpene precursor remains unknown; elucidating this mechanism will help in development of new technology to generate novel natural product-like scaffolds. With classic and QM/MM MD simulations, an “open-closed” conformation change of the catalytic pocket was observed in the GPPS active site at its large subunit (LSU), and a critical salt bridge between Asp91(in loop 1) and Lys239(in loop 2) was identified. The salt bridge is responsible for opening or closing the catalytic pocket. Meanwhile, the small subunit (SSU) regulates the size and shape of the hydrophobic pocket to flexibly host substrates with different shapes and sizes (DPP/GPP/FPP, C5/C10/C15). Further QM/MM MD simulations were carried out to explore the binding modes for the different substrates catalyzed by GPPS. Our simulations suggest that the key residues (Asp91, Lys239, and Gln156) are good candidates for site-directed mutagenesis and may help in protein engineering.

Discovery of the isoform-selective histone deacetylases (HDACs) inhibitors is of great medical importance and still a challenge. The comparison studies on the structure–function relationship of the conserved residues, which are located in the linker binding channel among class I HDACs (including 4 isoforms: HDAC1/2/3/8), have been carried out by using ab initio QM/MM MD simulations, a state-of-the-art approach to simulate metallo-enzymes. We found that the conserved tyrosine (Y303/308/286/306 in HDAC1/2/3/8, respectively) could modulate the zinc-inhibitor chelation among all class I HDACs with different regulatory mechanisms. For HDAC1/2/3 selective-inhibitor benzamide, the conserved tyrosine could modulate the coordinative ability of the central atom (Zn2+), while for pan-inhibitor SAHA, the conserved tyrosine could increase the chelating ability of the ligand (SAHA). Moreover, it is first found that the conserved tyrosine is correlated with the intertransformation of π–π stacking styles (parallel shift vs T-shaped) by the aromatic ring in benzamide and the two conserved phenylalanine residues of HDACs. In addition, the catalytic roles of the conserved tyrosine in stabilizing the transition state and intermediate are further revealed. These findings provide useful molecular basis knowledge for further isoform-selective inhibitor design among class I HDACs.

The glucosamine 6-phosphate deaminase (NagB), which catalyzes the conversion of D-glucosamine 6-phosphate (GlcN6P) into D-fructose 6-phosphate (F6P) and ammonia, determines the final metabolic fate of N-acetylglucosamine (GlcNAc). Here using state-of-the-art ab initio QM/MM MD simulations, we have explored the plausible mechanisms for the enzymatic ring-opening of GlcN6P in the basic environment. Two different proton-shuttle mechanisms have been proposed. Calculations show that the protonated state of the amino group in the substrate dominates the concerted and stepwise catalytic pathways and a catalytic triad plays an important role in mediating the proton transfer and the resulting ring-opening process. The free energy barrier for the rate-determining step in the low-energy stepwise reaction is 17.9 kcal mol−1. In acidic solution, the lid motif prefers a closed state while it always stays in the open state in basic solution upon substrate binding, which is basically dominated by the protonated state of the residue His145.

Although various Trypanosoma vivax purine-specific inosine–adenosine–guanosine nucleoside hydrolase (IAG-NH) crystal structures have been determined and the chemical reaction mechanism of substrate hydrolysis has been studied recently, the mechanistic details for the release of base and ribose are still unclear. Herein molecular dynamics (MD) simulations combined with umbrella sampling technique were utilized to explore the regulation mechanisms of key residues and loops 1 and 2 for the base release. Our results have indicated that the base release process is not the rate-limiting step in the entire hydrolysis process, and the very low barrier of ~ 5.6 kcal/mol can be washed out easily by the notable exothermicity from the substrate hydrolysis step. Moreover, the MD simulations have revealed that Glu82/Trp83 in loop 1 and His247/Arg252 in loop 2 are important to modulate the base release. The partial helix-to-coil change of loop 2 along with the base release process has been observed, showing good agreement with the IAG-NH crystal structures. The local binding site around the ribose after the base release is also discussed.Highlights► The inner-work mechanism of base release in IAG-nucleoside hydrolase was illuminated by MD simulations. ► The key residues promoting the base release have been characterized. ► The helix-to-coil change of loop 2 was revealed and it was in agreement with the XRD's structures. ► Our simulations confirmed the key switches to control the loop's open/closed state and found His247 is a novel controller.

The glucosamine 6-phosphate deaminase (NagB), which catalyzes the conversion of D-glucosamine 6-phosphate (GlcN6P) into D-fructose 6-phosphate (F6P) and ammonia, determines the final metabolic fate of N-acetylglucosamine (GlcNAc). Here using state-of-the-art ab initio QM/MM MD simulations, we have explored the plausible mechanisms for the enzymatic ring-opening of GlcN6P in the basic environment. Two different proton-shuttle mechanisms have been proposed. Calculations show that the protonated state of the amino group in the substrate dominates the concerted and stepwise catalytic pathways and a catalytic triad plays an important role in mediating the proton transfer and the resulting ring-opening process. The free energy barrier for the rate-determining step in the low-energy stepwise reaction is 17.9 kcal mol−1. In acidic solution, the lid motif prefers a closed state while it always stays in the open state in basic solution upon substrate binding, which is basically dominated by the protonated state of the residue His145.

Histone Deacetylases (HDACs) are promising anticancer targets and several selective inhibitors have been created based on the architectural differences of foot-pockets among HDACs. However, the “gate-keeper” of foot-pockets is still controversial. Herein, it is for the first time revealed that a conserved R–E salt bridge plays a critical role in keeping foot-pockets closed in class-II HDACs by computational simulations. This finding is further substantiated by our mutagenesis experiments.

SAHA (vorinostat, Merck) is a famous clinical drug for zinc-containing histone deacetylase (HDAC) targets against cancer and several other human disorders, whose inhibition mechanism (namely the protonation mechanism) upon binding to HDAC has been debated for more than ten years. It is very challenging to verify experimentally and is still controversial theoretically. The popular “Class-dependent” (namely “Tyr-dependent”) hypothesis is that the deprotonation of SAHA is mostly regulated by the conserved Tyr308 in class I HDAC while it is replaced by the His843 in class IIa HDAC. Herein, by elaborate QM(DFT)/MM MD simulations, we exclude the prevalent “Class-dependent” mechanism and advance a novel “Metal-dependent” mechanism, where the remote second metal site (K+ in most HDAC and Ca2+ in HDAC2) determines the protonation of SAHA. This proof-of-principle “Metal-dependent” mechanism opens up a new avenue to utilize the second metal site for isoform-selective inhibitor design.