Molecular dynamics, crystallography and mutagenesis studies on the substrate gating mechanism of prolyl oligopeptidase

Altered prolyl oligopeptidase (PREP) activity is found in many common neurological and other genetic disorders, and in some cases PREP inhibition may be a promising treatment. The active site of PREP resides in an internal cavity; in addition to the direct interaction between active site and substrate or inhibitor, the pathway to reach the active site (the gating mechanism) must be understood for more rational inhibitor design and understanding PREP function. The gating mechanism of PREP has been investigated through molecular dynamics (MD) simulation combined with crystallographic and mutagenesis studies. The MD results indicate the inter-domain loop structure, comprised of 3 loops at residues, 189-209 (loop A), 577-608 (loop B), and 636-646 (loop C) (porcine PREP numbering), are important components of the gating mechanism. The results from enzyme kinetics of PREP variants also support this hypothesis: When loop A is (1) locked to loop B through a disulphide bridge, all enzyme activity is halted, (2) nicked, enzyme activity is increased, and (3) removed, enzyme activity is only reduced. Limited proteolysis study also supports the hypothesis of a loop A driven gating mechanism. The MD results show a stable network of H-bonds that hold the two protein domains together. Crystallographic study indicates that a set of known PREP inhibitors inhabit a common binding conformation, and this H-bond network is not significantly altered. Thus the domain separation, seen to occur in lower taxa, is not involved in the gating mechanism for mammalian PREP. In two of the MD simulations we observed a conformational change that involved the breaking of the H-bond network holding loops A and B together. We also found that this network was more stable when the active site was occupied, thus decreasing the likelihood of this transition.     

Computational study of IAG-nucleoside hydrolase: determination of the preferred ground state conformation and the role of active site residues

The mechanism of action of inosine-adenosine-guanosine nucleoside hydrolase (IAG-NH) has been investigated by long-term molecular dynamics (MD) simulation in TIP3P water using stochastic boundary conditions. Special attention has been given to the role of leaving group pocket residues and conformation of the bound substrate at the active site of IAG-NH. We also describe the positioning of the residues of an important flexible loop at the active site, which was previously unobservable by X-ray crystallography due to high B-factors. Five MD simulations have been performed with the Enzyme x Substrate complexes: Enzyme x anti-Adenosine with Asp40-COOH [E(40H) x Ade(a)], Enzyme x anti-Adenosine with Asp40-COO- [E(40) x Ade(a)], Enzyme x syn-Adenosine with Asp40-COOH [E(40H) x Ade(s)], Enzyme x syn-Adenosine with Asp40-COO- [E(40) x Ade(s)], and Enzyme x anti-Inosine with Asp40-COO- [E(40) x Ino(a)]. Overall, the structures generated from the MD simulation of E(40H) x Ade(s) preserve the catalytically important hydrogen bonds as well as electrostatic and hydrophobic interactions to provide a plausible catalytic structure. When deprotonated Asp40 (Asp4-COO-) is present, the active site is open to water solvent which interferes with the base stacking between Trp83 and nucleobase. A calculation using Poisson-Boltzmann equation module supports that Asp40 indeed has an elevated pK(app). Solvent accessible surface area (SASA) calculations on all the five MD structures shows that systems with protonated Asp40, namely, E(40H) x Ade(a) and E(40H) x Ade(s), have zero SASA. It is found that a water molecule is hydrogen-bonded to the N7 of the nucleobase and is probably the essential general acid to protonate the departing nucleobase anion. The N7-bonded water is in turn hydrogen-bonded to waters in a channel, held in place by the residues of the flexible loop, Tyr257, His247, and Cys245. Using normal-mode analysis with elastic network model, we find that the flexible loop explores a conformational space much larger than in the MD trajectory, leading to a "gating"-like motion with respect to the active site.     

Molecular dynamics simulations of arachidonic acid complexes with COX-1 and COX-2: insights into equilibrium behavior

The cyclooxygenase (COX) enzymes are responsible for the committed step in prostaglandin biosynthesis, the generation of prostaglandin H(2). As a result, these enzymes are pharmacologically important targets for nonsteroidal antiinflammatory drugs, such as aspirin and newer COX-2 selective inhibitors. The cyclooxygenases are functional homodimers, and each subunit contains both a cyclooxygenase and a peroxidase active site. These enzymes are quite interesting mechanistically, as the conversion of arachidonic acid to prostaglandin H(2) requires two oxygenation and two cyclization reactions, resulting in the formation of five new chiral centers with nearly absolute regio- and stereochemical fidelity. We have used molecular dynamics (MD) simulations to investigate the equilibrium behavior of both COX-1 and COX-2 enzyme isoforms with bound arachidonate. These simulations were compared with reference simulations of arachidonate in solution to explore the effect of enzyme on substrate conformation and positioning in the active site. The simulations suggest that the substrate has greater conformational freedom in the COX-2 active site, consistent with the larger COX-2 active site volume observed in X-ray crystal structures. The simulations reveal different conformational behavior for arachidonate in each subunit over the course of extended equilibrium MD simulations. The simulations also provide detailed information for several protein channels that might be important for oxygen and water transport to or from active sites or for intermediate trafficking between the cyclooxygenase and peroxidase active sites. The detailed comparisons for COX-1 versus COX-2 active site structural fluctuations may also provide useful information for design of new isozyme-selective inhibitors.     

Carbon monoxide binding to the heme group at the dimeric interface modulates structure and copper accessibility in the Cu,Zn superoxide dismutase from Haemophilus ducreyi: in silico and in vitro evidences

X-ray absorption near-edge structure (XANES) spectroscopy and molecular dynamics (MD) simulations have been jointly applied to the study of the Cu,Zn superoxide dismutase from Haemophilus ducreyi (HdSOD) in interaction with the carbon monoxide molecule. The configurational flexibility of the Fe(II)-heme group, intercalated between the two subunits, has been sampled by MD simulations and included in the XANES data analysis without optimization in the structural parameter space. Our results provide an interpretation of the observed discrepancy in the Fe-heme distances as detected by extended X-ray absorption fine structure (EXAFS) spectroscopy and the classical XANES analysis, in which the structural parameters are optimized in a unique structure. Moreover, binding of the CO molecule to the heme induces a long range effect on the Cu,Zn active site, as evidenced by both MD simulations and in vitro experiments. MD simulation of the CO bound system, in fact, highlighted a structural rearrangement of the protein-protein hydrogen bond network in the region of the Cu,Zn active site, correlated with an increase in water accessibility at short distance from the copper atom. In line, in vitro experiments evidenced an increase of copper accessibility to a chelating agent when the CO molecule binds to the heme group, as compared to a heme deprived HdSOD. Altogether, our results support the hypothesis that the HdSOD is a heme-sensor protein, in which binding to small gaseous molecules modulates the enzyme superoxide activity as an adaptive response to the bacterial environment.     

MD simulations reveal alternate conformations of the oxyanion hole in the Zika virus NS2B/NS3 protease

Recent crystallography studies have shown that the binding site oxyanion hole plays an important role in inhibitor binding, but can exist in two conformations (active/inactive). We have undertaken molecular dynamics (MD) calculations to better understand oxyanion hole dynamics and thermodynamics. We find that the Zika virus (ZIKV) NS2B/NS3 protease maintains a stable closed conformation over multiple 100-ns conventional MD simulations in both the presence and absence of inhibitors. The S1, S2, and S3 pockets are stable as well. However, in two of eight simulations, the A132-G133 peptide bond in the binding pocket of S1' spontaneously flips to form a 3₁₀-helix that corresponds to the inactive conformation of the oxyanion hole, and then maintains this conformation until the end of the 100-ns conventional MD simulations without inversion of the flip. This conformational change affects the S1' pocket in ZIKV NS2B/NS3 protease active site, which is important for small molecule binding. The simulation results provide evidence at the atomic level that the inactive conformation of the oxyanion hole is more favored energetically when no specific interactions are formed between substrate/inhibitor and oxyanion hole residues. Interestingly, however, transition between the active and inactive conformation of the oxyanion hole can be observed by boosting the valley potential in accelerated MD simulations. This supports a proposed induced-fit mechanism of ZIKV NS2B/NS3 protease from computational methods and provides useful direction to enhance inhibitor binding predictions in structure-based drug design.     

On the Permeation by Dioxygen of Urate Oxidase from Aspergillus flavus in Complex with Xanthine Anion: Dioxygen Pathways and a Portrait of the Enzyme Cavities from Molecular Dynamics Simulations in Water Solution

This work describes molecular dynamics (MD) simulations in aqueous media for the complex of the homotetrameric urate oxidase (UOX) from Aspergillus flavus with xanthine anion ( 5 ) in the presence of dioxygen (O₂). After 196.6 ns of trajectory from unrestrained MD, a O₂ molecule was observed leaving the bulk solvent to penetrate the enzyme between two subunits, A/C. From here, the same O₂ molecule was observed migrating, across subunit C, to the hydrophobic cavity that shares residue V227 with the active site. The latter was finally attained, after 378.3 ns of trajectory, with O₂ at a bonding distance from 5 . The reverse same O₂ pathway, from 5 to the bulk solvent, was observed as preferred pathway under random acceleration MD (RAMD), where an external, randomly oriented force was acting on O₂. Both MD and RAMD simulations revealed several cavities populated by O₂ during its migration from the bulk solvent to the active site or backwards. Paying attention to the last hydrophobic cavity that apparently serves as O₂ reservoir for the active site, it was noticed that its volume undergoes ample fluctuations during the MD simulation, as expected from the thermal motion of a flexible protein, independently from the particular subunit and no matter whether the cavity is filled or not by O₂.     

Ion permeation through a narrow channel: using gramicidin to ascertain all-atom molecular dynamics potential of mean force methodology and biomolecular force fields

We investigate methods for extracting the potential of mean force (PMF) governing ion permeation from molecular dynamics simulations (MD) using gramicidin A as a prototypical narrow ion channel. It is possible to obtain well-converged meaningful PMFs using all-atom MD, which predict experimental observables within order-of-magnitude agreement with experimental results. This was possible by careful attention to issues of statistical convergence of the PMF, finite size effects, and lipid hydrocarbon chain polarizability. When comparing the modern all-atom force fields of CHARMM27 and AMBER94, we found that a fairly consistent picture emerges, and that both AMBER94 and CHARMM27 predict observables that are in semiquantitative agreement with both the experimental conductance and dissociation coefficient. Even small changes in the force field, however, result in significant changes in permeation energetics. Furthermore, the full two-dimensional free-energy surface describing permeation reveals the location and magnitude of the central barrier and the location of two binding sites for K(+) ion permeation near the channel entrance--i.e., an inner site on-axis and an outer site off-axis. We conclude that the MD-PMF approach is a powerful tool for understanding and predicting the function of narrow ion channels in a manner that is consistent with the atomic and thermally fluctuating nature of proteins.     

Spectral signatures of heterogeneous protein ensembles revealed by MD Simulations of 2DIR spectra

A model for the calculation of amide I FTIR and 2DIR spectra taking into account fluctuations in hydrogen bonding and structure from molecular dynamics (MD) simulations is tested on three systems. It is found that although the homogeneous lineshape approximation yields satisfactory FTIR spectra, 2DIR spectra are sensitive to the inhomogeneity naturally present in any solvated protein and the common approximations of a static structure and averaged-effect solvent are invalid. By building on the local amide Hamiltonian and incorporating site energy variation with electrostatic-based models and disorder from MD trajectories, good agreement is obtained between calculated and measured 2DIR spectra. The largest contribution to the observed inhomogeneity is found to be the fluctuating site energies, which in turn are most sensitive to the water solvent. With the ability to accurately predict 2DIR spectra from atomistic simulations, new opportunities to test force fields and mechanistic predictions from MD are revealed.     

Superoxide dismutase: fluctuations in the structure and solvation of the active site channel studied by molecular dynamics simulation

The molecular dynamics (MD) simulation of superoxide dismutase (SOD) in water is carried out for a total of 23 ps. The simulation system is a 26 A sphere centered at the active site of SOD, including 1602 atoms from SOD and 1761 water molecules. There is no gross deviation from the x-ray structure for the average MD structure. The structure and potential fluctuations around the active site are examined. The results provide new insight to the interactions between SOD and its substrate superoxide.     

Molecular dynamic simulations of the N-terminal receiver domain of NtrC reveal intrinsic conformational flexibility in the inactive state

The N-terminal receiver domain of NtrC is the molecular switch in the two-component signal transduction. It is the first protein where structures of both the active (phosphyroylated) and inactive (unphosphyroylated) states are determined experimentally. Phosphorylation of the NtrC at the active site induces large structural change. NMR experiments suggested that the wild type unphosphorylated NtrC adopts both the active and the inactive conformations and the phosphorylation stabilizes the active conformations. We applied free (unconstrained) molecular dynamic (MD) simulation to examine the intrinsic flexibilities and stabilities of the NtrC receiver domain in both the active and inactive conformations. Molecular dynamic simulations showed that the inactive state of NtrC receiver domain is more flexible than the active state. There were large movements in helix 4 and loop beta3-alpha3 which coincide with major structural differences between the inactive and active states. We observed large root-mean-square deviations from the initial starting structure and the large root-mean-square fluctuations during MD simulation for the inactive state. We then investigated the activation pathway with Targeted MD simulation. We show that the intrinsic flexibility in the loop beta3-alpha3 plays an important role in triggering the conformational change. Phosphorylation at the active site may serve to stabilize the conformational change. These results together suggest that the unphosphorylated NtrC receiver domain could be involved in a conformational equilibrium between two different states.     

Modeling a new water channel that allows SET9 to dimethylate p53

SET9, a protein lysine methyltransferase, has been thought to be capable of transferring only one methyl group to target lysine residues. However, some reports have pointed out that SET9 can dimethylate Lys372 of p53 (p53-K372) and Lys4 of histone H3 (H3-K4). In order to understand how p53 can be dimethylated by SET9, we measured the radius of the channel that surrounds p53-K372, first on the basis of the crystal structure of SET9, and we show that the channel is not suitable for water movement. Second, molecular dynamic (MD) simulations were carried out for 204 ns on the crystal structure of SET9. The results show that water leaves the active site of SET9 through a new channel, which is made of G292, A295, Y305 and Y335. In addition, the results of molecular docking and MD simulations indicate that the new water channel continues to remain open when S-adenosyl-L-methionine (AdoMet) or S-adenosyl-L-homocysteine (AdoHcy) is bound to SET9. The changes in the radii of these two channels were measured in the equilibrium phase at the constant temperature of 300 K. The results indicate that the first channel still does not allow water to get into or out of the active site, but the new channel is large enough to allow this water to circulate. Our results indicate that water can be removed from the active site, an essential process for allowing the dimethylation reaction to occur.     

Ligand diffusion in the catalase from Proteus mirabilis: a molecular dynamics study

The role of the channels and cavities present in the catalase from Proteus mirabilis (PMC) was investigated using molecular dynamics (MD) simulations. The reactant and products of the reaction, H(2)O(2) -->1/2 O(2) + H(2)O, catalyzed by the enzyme were allowed to diffuse to and from the active site. Dynamic fluctuations in the structure are found necessary for the opening of the major channel, identified in the X-ray model, which allows access to the active site. This channel is the only pathway to the active site observed during the dynamics, and both the products and reactant use it. H(2)O and O(2) are also detected in a cavity defined by the heme and Ser196, which could play an important role during the reaction. Free energy profiles of the ligands diffusing through the major channel indicate that the barriers to ligand diffusion are less than 20 kJ mol(-1) for each of the species. It is not clear from our study that minor channels play a role for access to the protein active site or to the protein surface.     

Dynamics of the native and the ligand-bound structures of eosinophil cationic protein: network of hydrogen bonds at the catalytic site

Eosinophil Cationic Protein (ECP) is sequentially and structurally similar to ribonuclease A (RNase A). It belongs to the RNase A family of proteins and the RNA catalysis is essential to its biological function. In the present study, we have generated the dinucleotide-bound structures of ECP by docking the dinucleotides to a number of molecular dynamics (MD) generated ECP structures. The stability of the docked enzyme-ligand complexes was ascertained by extensive MD simulations. The modes of ligand binding are explored by essential dynamics studies. The role of water molecules in the stability of the complex and in the catalysis was investigated. The active site residues form a complex network of connections with the ligand and with a water molecule. The catalytic mechanism of the RNA cleavage is examined on the basis of the active site geometry obtained by the simulations.     

Acetylcholinesterase: mechanisms of covalent inhibition of H447I mutant determined by computational analyses

The reaction mechanisms of two inhibitor TFK(+) and TFK(0) binding to H447I mutant mouse acetylcholinesterase (mAChE) have been investigated by using a combined ab initio quantum mechanical/molecular mechanical (QM/MM) approach and classical molecular dynamics (MD) simulations. TFK(+) binding to the H447I mutant may proceed with a different reaction mechanism from the wild-type. A water molecule takes over the role of His447 and participates in the bond breaking and forming as a "charge relayer". Unlike in the wild-type mAChE case, Glu334, a conserved residue from the catalytic triad, acts as a catalytic base in the reaction. The calculated energy barrier for this reaction is about 8kcal/mol. These predictions await experimental verification. In the case of the neutral ligand TFK(0), however, multiple MD simulations on the TFK(0)/H447I complex reveal that none of the water molecules can be retained in the active site as a "catalytic" water. Taken together our computational studies confirm that TFK(0) is almost inactive in the H447I mutant, and also provide detailed mechanistic insights into the experimental observations.     

Gramicidin A channel as a test ground for molecular dynamics force fields

We use the well-known structural and functional properties of the gramicidin A channel to test the appropriateness of force fields commonly used in molecular dynamics (MD) simulations of ion channels. For this purpose, the high-resolution structure of the gramicidin A dimer is embedded in a dimyristoylphosphatidylcholine bilayer, and the potential of mean force of a K(+) ion is calculated along the channel axis using the umbrella sampling method. Calculations are performed using two of the most common force fields in MD simulations: CHARMM and GROMACS. Both force fields lead to large central barriers for K(+) ion permeation, that are substantially higher than those deduced from the physiological data by inverse methods. In long MD simulations lasting over 60 ns, several ions are observed to enter the binding site but none of them crossed the channel despite the presence of a large driving field. The present results, taken together with many earlier studies, highlights the shortcomings of the standard force fields used in MD simulations of ion channels and calls for construction of more appropriate force fields for this purpose.     

Atomic crystal and molecular dynamics simulation structures of human carbonic anhydrase II: insights into the proton transfer mechanism

Human carbonic anhydrase II (HCA II) is a zinc-metalloenzyme that catalyzes the reversible interconversion of CO2 and HCO3-. The rate-limiting step of this catalysis is the transfer of a proton between the Zn-bound solvent molecule and residue His64. In order to fully characterize the active site structural features implicated in the proton transfer mechanism, the refined X-ray crystal structure of uncomplexed wild type HCA II to 1.05 A resolution with an Rcryst value of 12.0% and an Rfree value of 15.1% has been elucidated. This structure provides strong clues as to the pathway of the intramolecular proton transfer between the Zn-bound solvent and His64. The structure emphasizes the role of the solvent network, the unique positioning of solvent molecule W2, and the significance of the dual conformation of His64 in the active site. The structure is compared with molecular dynamics (MD) simulation calculations of the Zn-bound hydroxyl/His64+ (charged) and the Zn-bound water/His64 (uncharged) HCA II states. A comparison of the crystallographic anisotropic atomic thermal parameters and MD simulation root-mean-square fluctuation values show excellent agreement in the atomic motion observed between the two methods. It is also interesting that the observed active site solvent positions in the crystal structure are also the most probable positions of the solvent during the MD simulations. On the basis of the comparative study of the MD simulation results, the HCA II crystal structure observed is most likely in the Zn-bound water/His64 state. This conclusion is based on the following observations: His64 is mainly (80%) orientated in an inward conformation; electron density omit maps infer that His64 is not charged in an either inward or outward conformation; and the Zn-bound solvent is most likely a water molecule.     

Molecular Dynamic Simulations of the N-Terminal Receiver Domain of NtrC Reveal Intrinsic Conformational Flexibility in the Inactive State

Abstract

The N-terminal receiver domain of NtrC is the molecular switch in the two-component signal transduction. It is the first protein where structures of both the active (phosphyroylated) and inactive (unphosphyroylated) states are determined experimentally. Phosphorylation of the NtrC at the active site induces large structural change. NMR experiments suggested that the wild type unphosphorylated NtrC adopts both the active and the inactive conformations and the phosphorylation stabilizes the active conformations. We applied free (unconstrained) molecular dynamic (MD) simulation to examine the intrinsic flexibilities and stabilities of the NtrC receiver domain in both the active and inactive conformations. Molecular dynamic simulations showed that the inactive state of NtrC receiver domain is more flexible than the active state. There were large movements in helix 4 and loop β3-α3 which coincide with major structural differences between the inactive and active states. We observed large root-mean-square deviations from the initial starting structure and the large root-mean-square fluctuations during MD simulation for the inactive state. We then investigated the activation pathway with Targeted MD simulation. We show that the intrinsic flexibility in the loop β3-α3 plays an important role in triggering the conformational change. Phosphorylation at the active site may serve to stabilize the conformational change. These results together suggest that the unphosphorylated NtrC receiver domain could be involved in a conformational equilibrium between two different states.     

Are time-averaged restraints necessary for nuclear magnetic resonance refinement? A model study for DNA.

A recently suggested method for refinement of structural data obtained from two-dimensional nuclear magnetic resonance experiments using molecular dynamics (MD) is explored. In this method, the time-averaged values of the appropriate internal co-ordinates of the molecule, calculated from the MD trajectory, are driven by restraints towards the experimental target values. This contrasts with most refinement procedures currently in use, where restraints are applied based on the instantaneous values of the appropriate co-ordinates. Both refinement methods are applied to the EcoRI restriction site DNA hexamer d(GAATTC)2, using target nuclear Overhauser enhancement distances derived from a one nanosecond unrestrained MD simulation of this structure. The resulting refined structures are compared to the results of the unrestrained MD trajectory, which serves as our "experimental" data. We show that although both methods can yield an average structure with the correct gross morphology, the new method allows both a much more realistic picture of inherent flexibility, and reproduces fine conformational detail better, such as sequence dependency. We also analyze the very long MD trajectory generated here (longer than any previously reported for a DNA oligomer), and find that significantly shorter simulations, typical of those frequently performed, may not yield acceptably reliable values for certain structural parameters.     

Outside‐binding site mutations modify the active site's shapes in neuraminidase from influenza A H1N1

The recent occurrence of 2009 influenza A (H1N1) pandemic as well as others has raised concern of a far more dangerous outcome should this virus becomes resistant to current drug therapies. The number of clinical cases that are resistant to oseltamivir (Tamiflu®) is larger than the limited number of neuraminidase (NA) mutations (H275Y, N295S, and I223R) that have been identified at the active site and that are associated to oseltamivir resistance. In this study, we have performed a comparative analysis between a set of NAs that have the most representative mutations located outside the active site. The recently crystallized NA‐oseltamivir complex (PDB ID: 3NSS) was used as a wild‐type structure. After selecting the target NA sequences, their three‐dimensional (3D) structure was built using 3NSS as a template by homology modeling. The 3D NA models were refined by molecular dynamics (MD) simulations. The refined models were used to perform a docking study, using oseltamivir as a ligand. Furthermore, the docking results were refined by free‐energy analysis using the MM‐PBSA method. The analysis of the MD simulation results showed that the NA models reached convergence during the first 10 ns. Visual inspection and structural measures showed that the mutated NA active sites show structural variations. The docking and MM‐PBSA results from the complexes showed different binding modes and free energy values. These results suggest that distant mutations located outside the active site of NA affect its structure and could be considered to be a new source of resistance to oseltamivir, which agrees with reports in the clinical literature. © 2012 Wiley Periodicals, Inc.