Molecular Dynamics (MD) simulation and binding energy calculation

The Molecular Dynamics (MD) simulations, including energy minimization, system equilibration and production protocols, were performed with the GROMACS 4.5.5 package as described previously^{40, 41}. Each enzyme-substrate complex was placed in a cubic box after adding GROMOS 9643a1 position, and filled with atomistic TIP3P water. A two-step energy minimization process was performed after system equilibration. The Molecular Mechanics Poisson–Boltzmann surface area (MM-PBSA) method, which has been widely used to predict the binding affinities for a variety of macromolecular complexes and protein-ligand complexes^41–43, was employed to calculate the binding energy of different pol-DNA template complexes (DNA template in PDB: 4NLG is 5′-GAAGCCGGCGGAA-3′). In this study, the free energy of each molecule is defined as follows: ΔG _binding = G _complex − (G _protein + G _ligand). Here, G _complex, G _protein, and G _ligand are the free energy of the polymerase-DNA complex, the free energy of the polymerase, and the free energy of DNA, respectively. The free energy G can be calculated by the following scheme, based on the MM-PBSA method: G = E _MM  TS _MM + G _solv, where E _MM is comprised of the intramolecular electrostatic (E _elec) and van der Waals (E _vdW) interaction energies. The free energy of solvation, G _solv, was approximated as the sum of electrostatic and nonpolar contributions in liquid phase. TS _MM was ignored as it did not contribute significantly to the binding energy in these conditions. All binding energies were calculated by Gromacs 4.5.5 combined with the g_mmpbsa program developed by Kumari et al.⁴⁴.