- Research article
- Open Access
Role of magnesium ions in the reaction mechanism at the interface between Tm1631 protein and its DNA ligand
Chemistry Central Journalvolume 10, Article number: 41 (2016)
A protein, Tm1631 from the hyperthermophilic organism Thermotoga maritima belongs to a domain of unknown function protein family. It was predicted that Tm1631 binds with the DNA and that the Tm1631–DNA complex is an endonuclease repair system with a DNA repair function (Konc et al. PLoS Comput Biol 9(11): e1003341, 2013). We observed that the severely bent, strained DNA binds to the protein for the entire 90 ns of classical molecular dynamics (MD) performed; we could observe no significant changes in the most distorted region of the DNA, where the cleavage of phosphodiester bond occurs. In this article, we modeled the reaction mechanism at the interface between Tm1631 and its proposed ligand, the DNA molecule, focusing on cleavage of the phosphodiester bond. After addition of two Mg2+ ions to the reaction center and extension of classical MD by 50 ns (totaling 140 ns), the DNA ligand stayed bolted to the protein. Results from density functional theory quantum mechanics/molecular mechanics (QM/MM) calculations suggest that the reaction is analogous to known endonuclease mechanisms: an enzyme reaction mechanism with two Mg2+ ions in the reaction center and a pentacovalent intermediate. The minimum energy pathway profile shows that the phosphodiester bond cleavage step of the reaction is kinetically controlled and not thermodynamically because of a lack of any energy barrier above the accuracy of the energy profile calculation. The role of ions is shown by comparing the results with the reaction mechanisms in the absence of the Mg2+ ions where there is a significantly higher reaction barrier than in the presence of the Mg2+ ions.
Tm1631 is a member of the domain of unknown function 72 (DUF) family in the Protein family (Pfam) database, a protein domain that has no characterized function . Protein function can only be unambiguously determined experimentally, but in case of a new protein with no computationally predicted putative function it is difficult to choose the correct experiment. New procedures are necessary to facilitate this research and improve determination of function for all the DUF proteins. A new approach to this problem has been developed by our group and was described in a previous report to predict the function of the Tm1631 protein .
In an earlier paper , we used the binding site comparison capability of the ProBiS algorithm [3, 4] to predict the binding site in the Tm1631 protein and to speculate on the nature of the DNA ligand that could bind to the Tm1631 protein. The Tm1631 protein was predicted to be analogous to endonuclease IV despite sharing <10 % sequence identity, and the proposed Tm1631–DNA complex was subjected to 90 ns long classical MD, using the CHARMM simulation package . This resulted in structures of the Tm1631–DNA complex that are used in the QM/MM study reported in this paper, where we attempted to investigate the catalytic mechanism of the reaction between Tm1631 and DNA. It was of especial interest in this connection to determine how Mg2+ ions act in this DNA binding site (Fig. 1). Endonuclease repair mechanism is an important mechanism that allows organisms to escape DNA damage and plays a major role in the prevention of cancer in higher organisms . Endonucleases cleave phosphodiester bond of DNA at the damaged site and create a nick in the phosphodiester backbone that is recognized by further repair enzymes in the base excision repair pathway. The endonuclease catalytic mechanism is thought to involve a hydroxide ion derived from water, which forms a bond with phosphorus of the DNA and induces cleavage of the phosphodiester bond.
Endonucleases require one, two or three divalent metal ions, such as Mn2+ or Mg2+ in the catalytic site [7–13]. It is believed that in thermal environment, such as the one Thermotoga maritima lives in at temperatures around 80 °C, the most suitable metal ion for this kind of system is Mg2+ . However, Tm1631 crystal structure (PDB: 1VPQ) does not contain any metal ions in the predicted binding site. On the other hand, there are relatively few crystallographically characterised magnesium binding sites . For some binding sites, a metal ion is observed in the binding site but which metal ion it is and how many of them are needed for the catalytic activity remains undetermined [16–19].
Here, we report the results of theoretical QM/MM studies of the reaction mechanism for the Tm1631 protein. We postulated a catalytic mechanism with two Mg2+ ions that resembles the one of the apurinic/apyrimidinic (AP) endonuclease enzyme presented by Mol et al. . We tested the system with and without the ions, and found that the energetically most favourable pathway of the phosphodiester bond cleavage catalysed by Tm1631 requires presence of Mg2+ ions. In the proposed catalytic mechanism (Fig. 1b), Lys73 of the Tm1631 should be in a deprotonated state before water ionizes. One of the Mg2+ ions attacks the water molecule, and subsequently water ionizes, the proton forms a bond with Lys73 and the OH group forms a bond with the phosphorus atom. The transition state then collapses leading to cleavage of the phosphodiester P-O3′ bond, while the O3′ leaving group is stabilized by the second Mg2+ ion. Our findings suggest that the catalytic mechanism of Tm1631 requires Mg2+ ions and is similar to known Steitz’s mechanism .
Experimental and computational studies [24, 25] point out that the cleavage of the phosphodiester bond occurs via SN2 nucleophilic substitution in three steps (Fig. 1b) explained in work by Sgrignani et al.  and Yang et al. .
Calculations were based on the crystal structure of protein Tm1631, (PDB: 1VPQ) and the DNA chains from another crystal structure (PDB: 2NQJ).
System setup and MD
The structures used in this work were previously equilibrated by classical MD simulation. Since there is no crystal structure of Tm1631 with DNA available in the Protein Data Bank (PDB), we used as the starting structure for this study the predicted Tm1631–DNA complex after 90 ns of classical MD, which was performed in our previous study . To validate this starting complex structure, we plotted its all-atom, protein Tm1631 and DNA root-mean-square deviations (RMSDs) (compared to the first snapshot at 0 ns of classical MD) dependence against the simulation time, which showed that in the last 20 ns of simulation the RMSDs have reached a plateau, suggesting that the starting structure is well equilibrated (Additional file 1: Figure S1).
In order to obtain good reactant and product structures we probed different Mg2+ ion and water molecules positions. In this search we required that Mg2+ initially coordinates with six oxygen atoms (Additional file 1: Table S1; Distances to Mg1 2+ and Mg2 2+), and with the distant environment atoms in the same positions for both, the reactant and product conformations. The final energy minimization procedures were performed without any constraints or restraints and we were able to obtain suitable starting positions for the reaction mechanism studies with the above mentioned properties.
Model building and QM/MM simulation
The starting structure of our simulation was the Tm1631–DNA complex after 90 ns of classical MD simulation. We replaced two water molecules with two Mg2+ ions and then minimized the system. A minimized structure with Mg2+ ions was further optimized at quantum mechanics/molecular mechanics (QM/MM) level using the CHARMM software package. The CHARMM force field parameters were used to describe the molecular mechanics (MM) part, while the quantum mechanics (QM) region (42 atoms in total, including both Mg2+ ions (Fig. 2a, b) was treated at the density functional theory (DFT) level using the B3LYP functional and the 6–31G* basis set. We used the Replica path (RPATh) method to divide the system into 16 structures equidistantly apart in the RMSD space between the reactants and products and minimized each obtained structure using 3000 steps of adopted basis Newton–Raphson (ABNR) minimization. After minimization we checked the distances between Mg1 2+ and Mg2 2+ ions, which were both around 4 Å. Both Mg2+ ions were coordinated with 6 oxygen atoms (reactants, Fig. 2c; product, Fig. 2d), and the distances between coordinated atoms and both Mg2+ ions were approximately 2 Å.
Ab initio QM/MM  methods were performed with the CHARMM biomolecular program . In the molecular mechanics part of the calculations, we used CHARMM force field version 36. The ab initio DFT calculations were performed using the general atomic and molecular electronic structure system (GAMESS)  software package, interfaced to the CHARMM program. We used the B3LYP/6–31G* level of theory, which is implemented in the GAMESS program. The ABNR minimization algorithm was used for energy minimizations. Molecular mechanics calculations were performed with a constant dielectric of ε = 1 using a classical force shift method and a cut-off distance of 12 Å. Molecular graphic images were produced using the VMD software package . All ab initio QM/MM and QM calculations were carried out on the CROW clusters at the National Institute of Chemistry in Ljubljana .
Chain-of-replica methods involve discretising any reaction pathway by defining replicated conformations along the path between the reactants and products . In order to keep the pathway along the reaction smooth a penalty term is included into the potential function which keeps the RMSD values equidistant between all the neighbouring conformations . This modified potential function is then used in the geometry minimization procedure to obtain the minimum energy pathway between the reactant and product structures. Since the penalty is in the RMSD space there is no preference in the reaction coordinate to the individual distances among the atoms and is thus an efficient tool to investigate the order of bond breaking and bond making during the reaction process. The RPATh method can be used for both minimizations and the MD simulations. In order to investigate the basic steps of the reaction process, one first explores the potential energy surface by minimum energy pathway from reactants to products by a minimization procedure. In the RPATh setup this means that all the replicas are minimised along with the restraint which forces them to be equidistant in RMSD space. To completely understand the energetics of the reaction process one must calculate the free energy along the reaction pathways. Unfortunately, it is not possible to perform such calculations with the satisfactory accuracy using the currently available computational methods, although much effort is invested to develop practically feasible QM/MM methods to calculate the free energy of a reaction processes [34, 35].
Study of concurrent reaction mechanisms
The following steps define the procedure for the reaction with two Mg2+ ions in the binding pocket:
For coordination and structural files, we read the last coordinate frame from the file produced by the 90 ns classical MD simulation as explained earlier .
To perform QM/MM calculations, the QM region and if necessary also the boundary between the classical potential and the quantum potential involving link atoms must be defined. We assigned the QM region as follows: Lys73 side chain, two water molecules close to reaction center: water number 1863 (2.710 Å; Fig. 2 and Additional file 1: Table S1) is involved in the reaction and nearby water 6134, phosphate group with adjacent C5′ atom from abasic dideoxyribose on place 7 in the DNA (residue 3DR7), and the 5-member ring from Cyt6 residue on the DNA. We used 4 link atoms, because we divided QM and MM region in the middle of four covalent bonds: one in Lys73, one in 3DR7, and two in Cyt6 residues. The total number of QM atoms in the system including two Mg2+ ions was 42 consisting of 170 electrons described by 343 basis functions. The total number of atoms in the system was 13,932, including 2860 water molecules.
A reactant was constructed for which the QM/MM minimization procedure was initiated by using RESD  restraints for bonded atoms P (3DR7) and O3′ (CYT6) so that they remained separated at 1.6 Å and for OH− (water 1863) and P (3DR7) which remained 3.5 Å apart. On both bivalent metal ions, the Mg2+ harmonical restraint (the “cons harm” command in CHARMM) was used to preserve the coordination of ions. At this stage 100 steps of ABNR minimization were performed. Then we removed CONS HARM restraint from the minimization run of 300 steps with only RESD restraints for bonds left. Subsequently we removed all the restraints and ran 1500 steps of minimization.
From the structure of the reactant we made a product, breaking the bond between P (3DR7) and O3′(CYT6) and created a bond between OH− (water 1863) and P (3DR7) 1.6 Å using restraints. We also restrained the distance between O5′ (3DR7) and P (3DR7) at 1.6 Å. After 300 steps with RESD, we removed all the restraints and ran 1500 steps of a geometry minimization procedure.
We replicated the whole system 16 times and produced initial replica conformations by linear interpolation of the coordinates between the reactant and the product. The first and last replica represented reactant and product conformations, respectively, and were fixed throughout the pathway minimizations. With this setup a 3000 steps ABNR minimization was performed.
Procedure steps for reaction without Mg2+ ions in the binding pocket:
We started with the minimized reactant from the simulation with two Mg2+ ions, then removed them from the structure and held the distance between OH− (water 1863) and P (3DR7) at 4 Å for 300 steps. Afterwards we removed all RESD restraints and ran 3000 steps of ABNR minimization.
We took the product from the previous simulation with two Mg2+ ions and removed them; subsequently running 3000 steps of an ABNR minimization procedure.
We performed an RPATh calculation, with 3000 steps of ABNR minimization procedure using identical number of replicas as in the calculations with ions present in the system.
To check the influence of the protein environment on the lowering of the energy barrier, we performed pure QM calculations with the GAMESS program as well as the reaction path calculations with the RPATh method. B3LYP/6–31G* level of theory was also used for all calculations in the pure QM case.
Procedure steps for reaction without Mg2+ ions and pure QM calculations:
First we made the reactant state from the already QM/MM minimised structure without Mg2+ ions. Subsequently we removed the whole protein, DNA structure and water molecules that were not in the QM region. Finally we ran 3000 steps of ABNR geometry minimization with restrains, holding all three protons attached to the Lys73.
We made a product from the QM/MM calculations without Mg2+ ions in the binding pocket. Then we ran 3000 steps of ABNR minimization without restrains.
RPATh calculation was finally performed with 3000 steps of ABNR minimization, using 16 replicas.
Results and discussion
The Tm1631 protein with yet unknown function from the organism T. maritima has a similar binding site to that of DNA repair proteins, as established earlier . By exploring the possible reaction pathways using QM/MM methods we tried to gain insight into the catalytic mechanism of the Tm1631–DNA complex. The similarity of the energetically most favourable pathway of the Tm1631–DNA complex with that of Mol et al.  strongly suggests that the mechanism is the same as in other endonucleases. The mechanism for Mg2+ ions catalysis that we propose is most likely the so-called Steitz’s mechanism [23, 37].
A variety of conformations within the active site were energetically evaluated and compared and the following systems were studied:
QM/MM calculation for the reaction with 2 × Mg2+ ions in the binding pocket (Additional file 2: Video S1).
QM/MM calculation for the reaction without Mg2+ ions in the binding pocket.
QM calculation for the reaction without Mg2+ ions in the binding pocket and environment with no protein or solvent molecules.
Quantum mechanics/molecular mechanics and QM calculations were performed using the B3LYP/6–31G* DFT method. As expected, obtained results suggest that the protein has an impact on lowering the reaction barrier and also establish that metal ions are required in the binding pocket. Our interest in this study is focused on the core DNA repair function, leaving the deprotonation of Lys73 for future investigations, however there is evidence that such a state exists . As explained in the previous section we made sure that the ions are coordinated as expected . We are aware of the possibility that one, two or three ions may be involved in the DNA repair reaction mechanisms. However we present in this paper the results for the two ion systems and only the core part of the repair mechanism where the hydroxide ion attacks the phosphate and the P-O3′ bond gets cleaved.
The hydroxide ion involved in nucleophilic attack on the phosphodiester bond P-O3′ was derived from water, ionization of which is accomplished with the help of an Mg2+ ion. Then proton from the same water molecule forms a bond with the nitrogen atom of the side chain of the Lys73, which in its deprotonated state can act as proton acceptor in the enzyme’s active site , and is part of the nucleic acid repair mechanism . We calculated the pKa of the Lys73 to be 9.55 using the DEPTH tool , which is the lowest of all lysine residues of the Tm1631 protein. We also calculated the average pKa of the Lys73′s surrounding residues (10.5), which suggested that Lys73 is most likely in a deprotonated state before reaction occurs and can accept a proton from water molecule (Additional file 1: Figure S2). Literature reports that Mg2+ plays a functional role in the catalytic mechanism and the stability of protein-DNA complex. Metal ions also lower the local pKa, and this, considering the harsh environment that the organism experiences, is in a good agreement with our study [7, 8, 41–43].
Magnesium ions coordination is essential for most phosphoryl transfer enzymes . The common catalytic mechanism was proposed previously . This mechanism works on the same principle as our proposed mechanism: Mg1 2+ coordinates the nucleophile and Mg2 2+ coordinates the leaving oxygen atom (O3′). Many similar systems with two Mg2+ ions in the binding pocket have been studied [12, 37]. Distance between both metal ions should be ~4 Å and in our case it is 3.85 Å (reactant) and 4.18 Å (product). Our system is also coordinated in the octahedral shape and most of the angles are 90° between O–Mg–O. Also seen from other enzymatic literature enzymatic phosphate hydrolysis proceeds as SN2-like nucleophilic attack on the scissile phosphate performed by an hydroxide ion, which is formed upon water activation [24–26, 45]. The important role of Mg2+ ions is lowering the pKa of its ligands and also for the presence of second metal ion (which coordinates nucleophilic water or hydroxide in the binding pocket). This leads to a mechanism with early proton transfer [46, 47], preceding the cleavage of phosphodiester bond in case on RNase catalytic system. Mg2+ ions have an essential contribution for the specific catalytic reactions by lowering the pKa of the leaving group and can impose specific geometry for the triphosphate chain—pentacovalent transition state and products .
We observe six critical points in the energy profile for two Mg2+ ions shown in Fig. 3. The minimum with the lowest energy is the frame 5, which is very close to a symmetric structure, the so-called pentacovalent intermediate, with P atom in the center, three oxygen atoms (O1P, O2P, O5′) in a planar arrangement, and the two reacting oxygen atoms positioned almost equidistant on the opposite sides of the plane. Usually such a structure would suggest a transition state in the reaction pathway but the Mg2+ ions effectively stabilize the energy of pentacovalent intermediate to make it a stable minimum. In order to check this structure and its energy we performed a separate calculation. We minimized the geometry of the single number five replica structure. After 10,000 steps of ABNR minimization the geometry had no visible changes to the one in the chain of replicated structures and the energy was lowered by 2.5 kcal/mol due to the removal of the restraint which slightly distorts the distances and makes the energy higher. The effect of the distortion could be made smaller by increasing the number of replicas, however this would not change the overall properties of the reaction mechanism. At the current level of accuracy one can estimate that energy noise level in Fig. 3 is <3 kcal/mol. This makes the reaction kinetically controlled exothermic and not thermodynamically controlled because there are no energy barriers between the reactants and the products. Our calculated energy change (−14.0 kcal/mol) agrees with results of another study , in which energy difference of phosphodiester bong cleavage starting from OH− was found to be −18.1 kcal/mol. The distances to important amino acids are reported in Additional file 1: Table S1.
Next, we present the results of QM/MM calculations without Mg2+ ions (Fig. 4), which supports two observations:
There is a high 17 kcal/mol barrier in the middle of the reaction path between structures five and nine. This structure has features analogous to those of structure six from Fig. 3, but in this case it is a transition state structure suggesting that the role of the Mg2+ ions is to transform this configuration into a stable and low energy pentacovalent intermediate.
Many chemical reactions have multiple reaction channels which depend mostly on the positions of the species entering the reaction. In the case where the two Mg2+ ions are not present the conformational space of the entering species is larger than the one with the two ions present. This makes the reaction energetically less favourable because the system may explore more pathways which are of higher energy than the ones with the ions present. This works in addition to the fact that the environment atoms in the enzyme systems usually lower the energy barriers of any transition state structure. At this point we can add that a possible solution to choose the most favourable reaction channel would be the use of QM/MM molecular dynamics. But there is a sampling problem to be emphasized which, to be resolved, would require tens of nanoseconds of simulation time what translates into tens of millions steps of complete QM/MM calculations. In the present studies less than 100 thousands steps of QM/MM calculations were performed and still it took a few months of CPU time using between 32 and 128 processors depending on the task.
To show the impact of the protein, DNA, solvent and ion environment on the studied reaction, we also studied the system in vacuum without Mg2+ ions. The vacuum calculation was set up keeping just the QM region in the reaction path calculations. From the Fig. 5 one can observe similar behaviour as in the QM/MM calculation without ions, but the barrier is extremely high which suggests that the protein environment indeed significantly contributes to the energy stabilisation of the pentacovalent transition state structure.
In order to verify the stability of the initial QM/MM setup we performed two additional classical MD 50 ns simulations starting with the reactant and product structures that we used for the minimum energy pathway calculations. The ions in these simulations kept the hexacoordinated structures throughout the simulations. In the case of reactant simulations none of the waters were exchanged and all 12 atoms around the two ions were identical. This means that the atom positions of reactant are in the stable and favorable positions to enter the reaction.
Protein Tm1631 from the organism T. maritima was predicted to be an endonuclease-like DNA binding protein, and consequently we investigated its function focusing specifically on the role of Mg2+ ions in its binding pocket. We performed a QM/MM study of Tm1631 in a complex with damaged DNA. We found that Mg2+ ions are required in the binding pocket in order that the reaction occurs. This allows us to conclude that Tm1631 is indeed an endonuclease binding protein with a reaction mechanism similar to that of other endonucleases. Some reconciliation is still needed regarding the number of metal ions, e.g. it is possible that only one ion suffices for the reaction to take place. The present paper is one of the few theoretical insights in the available literature to study a series of the reactions that play a role in the complex of the endonuclease repair process. Future work should be aimed at determination of the precise number and type of ions that are needed for the reaction to occur. Another interesting study would be to explain the formation of the hydroxide ion in connection with the protonation of Lys73 and the role of ions in such mechanisms. It would also be interesting to compare the present results with the results obtained by similar studies in different proteins (Additional file 3).
Bateman A, Coin L, Durbin R, Finn RD, Hollich V, Griffiths-Jones S, Khanna A, Marshall M, Moxon S, Sonnhammer EL et al (2004) The Pfam protein families database. Nucleic Acids Res 32(suppl 1):D138–D141
Konc J, Hodošček M, Ogrizek M, Konc JT, Janežič D (2013) Structure-based function prediction of uncharacterized protein using binding sites comparison. PLoS Comput Biol 9:e1003341
Konc J, Janežič D (2010) ProBiS algorithm for detection of structurally similar protein binding sites by local structural alignment. Bioinformatics 26:1160–1168
Konc J, Janežič D (2012) ProBiS-2012: web server and web services for detection of structurally similar binding sites in proteins. Nucleic Acids Res 40:W214–W221
Brooks BR, Brooks CL, MacKerell AD, Nilsson L, Petrella RJ, Roux B, Won Y, Archontis G, Bartels C, Boresch S (2009) CHARMM: the biomolecular simulation program. J Comput Chem 30:1545–1614
Johnson G (2010) Unearthing Prehistoric Tumors, and Debate. NY Times 61:233–246
Zahran M, Berezniak T, Imhof P, Smith JC (2011) Role of magnesium ions in DNA recognition by the EcoRV restriction endonuclease. FEBS Lett 585:2739–2743
Imhof P, Fischer S, Smith JC (2009) Catalytic mechanism of DNA backbone cleavage by the restriction enzyme EcoRV: a quantum mechanical/molecular mechanical analysis. Biochemistry (Mosc) 48:9061–9075
Kluge S, Weston J (2005) Can a hydroxide ligand trigger a change in the coordination number of magnesium ions in biological systems? Biochemistry (Mosc) 44:4877–4885
Paspaleva K, Thomassen E, Pannu NS, Iwai S, Moolenaar GF, Goosen N, Abrahams JP (2007) Crystal structure of the DNA repair enzyme ultraviolet damage endonuclease. Structure 15:1316–1324
Paspaleva K, Natuurwetenschappen, F. der W. en. In: Structure and function of the UVDE repair protein https://openaccess.leidenuniv.nl/handle/1887/13711 Accessed 18 Apr 2015
Dupureur CM (2010) One is enough: insights into the two-metal ion nuclease mechanism from global analysis and computational studies. Metallomics 2:609–620
Martínek V, Bren U, Goodman MF, Warshel A, Florián J (2007) DNA polymerase β catalytic efficiency mirrors the Asn279–dCTP H-bonding strength. FEBS Lett 581:775–780
Dupureur CM (2008) Roles of metal ions in nucleases. Curr Opin Chem Biol 12:250–255
Cowan JA (1997) Inorganic biochemistry: an introduction. John Wiley & Sons, New York
Cowan JA (1998) Metal activation of enzymes in nucleic acid biochemistry. Chem Rev 98:1067–1088
Cowan JA (2002) Structural and catalytic chemistry of magnesium-dependent enzymes. Biometals 15:225–235
Horton NC, Perona JJ (2004) DNA Cleavage by Eco RV endonuclease: two metal ions in three metal ion binding sites. Biochemistry (Mosc) 43:6841–6857
Palermo G, Stenta M, Cavalli A, Dal Peraro M, De Vivo M (2013) Molecular simulations highlight the role of metals in catalysis and inhibition of type II topoisomerase. J Chem Theory Comput 9:857–862
Lahiri SD, Zhang G, Dunaway-Mariano D, Allen KN (2003) The Pentacovalent Phosphorus Intermediate of a Phosphoryl Transfer Reaction. Science 299:2067–2071
Mol CD, Izumi T, Mitra S, Tainer JA (2000) DNA-bound structures and mutants reveal abasic DNA binding by APE1 DNA repair and coordination. Nature 403:451–456
Zhang H, Bren U, Kozekov ID, Rizzo CJ, Stec DF, Guengerich FP (2009) Steric and electrostatic effects at the C2 atom substituent influence replication and miscoding of the DNA deamination product deoxyxanthosine and analogs by DNA polymerases. J Mol Biol 392:251–269
Steitz TA, Steitz JA (1993) A general two-metal-ion mechanism for catalytic RNA. Proc Natl Acad Sci 90:6498–6502
Rosta E, Nowotny M, Yang W, Hummer G (2011) Catalytic mechanism of RNA backbone cleavage by ribonuclease H from quantum mechanics/molecular mechanics simulations. J Am Chem Soc 133:8934–8941
De Vivo M, Dal Peraro M, Klein ML (2008) Phosphodiester cleavage in ribonuclease H occurs via an associative two-metal-aided catalytic mechanism. J Am Chem Soc 130:10955–10962
Sgrignani J, Magistrato A (2015) QM/MM MD simulations on the enzymatic pathway of the human flap endonuclease (hFEN1) elucidating common cleavage pathways to RNase H enzymes. ACS Catal 5:3864–3875
Yang W (2011) Nucleases: diversity of structure, function and mechanism. Q Rev Biophys 44:1–93
Field MJ, Bash PA, Karplus M (1990) A combined quantum mechanical and molecular mechanical potential for molecular dynamics simulations. J Comput Chem 11:700–733
Eurenius KP, Chatfield DC, Brooks BR, Hodoscek M (1996) Enzyme mechanisms with hybrid quantum and molecular mechanical potentials. I. Theoretical considerations. Int J Quantum Chem 60:1189–1200
Humphrey W, Dalke A, Schulten K (1996) VMD: visual molecular dynamics. J Mol Graph 14:33–38
Hodoscek M, Borstnik U, Janezic D (2001) CROW for large scale macromolecular simulations. Cell Mol Biol Lett 7:118–119
Woodcock HL, Hodošček M, Brooks BR (2007) Exploring SCC-DFTB paths for mapping QM/MM reaction mechanisms. J Phys Chem A 111:5720–5728
Tao P, Hodošček M, Larkin JD, Shao Y, Brooks BR (2012) Comparison of three chain-of-states methods: nudged elastic band and replica path with restraints or constraints. J Chem Theory Comput 8:5035–5051
Hudson PS, White JK, Kearns FL, Hodoscek M, Boresch S, Woodcock HL (2015) Efficiently computing pathway free energies: New approaches based on chain-of-replica and non-boltzmann bennett reweighting schemes. Biochim Biophys Acta 1850:944–953
König G, Brooks BR (2015) Correcting for the free energy costs of bond or angle constraints in molecular dynamics simulations. Biochim Biophys Acta 1850:932–943
Lee YS, Hodoscek M, Brooks BR, Kador PF (1998) Catalytic mechanism of aldose reductase studied by the combined potentials of quantum mechanics and molecular mechanics. Biophys Chem 70:203–216
Palermo G, Cavalli A, Klein ML, Alfonso-Prieto M, Dal Peraro M, De Vivo M (2015) Catalytic metal ions and enzymatic processing of DNA and RNA. Acc Chem Res 48:220–228
Weber J, Hammond ST, Wilke-Mounts S, Senior AE (1998) Mg2+ Coordination in Catalytic Sites of F1-ATPase. Biochemistry (Mosc) 37:608–614
Teles JH, Brode S, Berkessel A (1998) Hydrogenation without a metal catalyst: an ab initio study on the mechanism of the metal-free hydrogenase from Methanobacterium thermoautotrophicum. J Am Chem Soc 120:1345–1346
Tan KP, Nguyen TB, Patel S, Varadarajan R, Madhusudhan MS (2014) Depth: a web server to compute depth, cavity sizes, detect potential small-molecule ligand-binding cavities and predict the pKa of ionizable residues in proteins. Nucleic Acids Res 41:W314–W321
Sissi C, Palumbo M (2009) Effects of magnesium and related divalent metal ions in topoisomerase structure and function. Nucleic Acids Res 37:702–711
Mones L, Kulhánek P, Florián J, Simon I, Fuxreiter M (2007) Probing the two-metal ion mechanism in the restriction endonuclease BamHI. Biochemistry (Mosc) 46:14514–14523
Akabayov B, Richardson CC (2011) Binding of Mn-deoxyribonucleoside triphosphates to the active site of the DNA polymerase of bacteriophage T7. Powder Diffr 26:159–162
Lopata A, Jambrina PG, Sharma PK, Brooks BR, Toth J, Vertessy BG, Rosta E (2015) Mutations decouple proton transfer from phosphate cleavage in the dUTPase catalytic reaction. ACS Catal 5:3225–3237
Cassano AG, Anderson VE, Harris ME (2004) Understanding the transition states of phosphodiester bond cleavage: insights from heavy atom isotope effects. Biopolymers 73:110–129
Florián J, Goodman MF, Warshel A (2003) Computer simulation of the chemical catalysis of DNA polymerases: discriminating between alternative nucleotide insertion mechanisms for T7 DNA polymerase. J Am Chem Soc 125:8163–8177
Mones L, Tang W-J, Florián J (2013) Empirical valence bond simulations of the chemical mechanism of ATP to cAMP conversion by anthrax edema factor. Biochemistry (Mosc) 52:2672–2682
The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript. All authors read and approved the final manuscript.
We sincerely thank Dr. Walter E. Knapp for hosting Mitja Ogrizek at FU Berlin and Dr. Petra Imhof for fruitful discussion.
The authors declare that they have no competing interests.
Financial support was provided through Grants P1-0002, J1-6736 and J1-6743 the Slovenian Research Agency. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.