After proper setup, the system was minimized for 20000 steps using the steepest descent algorithm followed by 2000 steps using the conjugate gradient algorithm. Then, the system was heated from 0 to 300 K in 50 ps under NVT ensemble with a weak restraint of 25 kcal/mol/Å 2 on protein. To achieve a uniform density after heating, 1 ns of density equilibration was performed under the NPT ensemble, where the target temperature of 300 K was maintained by using the Langevin thermostat with a collision frequency of 2 ps, while the target pressure of 1.0 atm was controlled by using the Berendsen barostat 60 with a pressure relaxation time of 1 ps. Afterwards, a production MD simulation was propagated for 50 ns under the NPT ensemble. During the equilibration and production runs, a weak restraint of 25 kcal/mol/Å 2 has been applied to the substrate, because the substrate binding is quite flexible at the enzyme-water interface. A 2 fs time-step was used for all MD runs. The SHAKE algorithm was used to constrain the covalent bonds involving hydrogen, while the "three-point" SHAKE algorithm was used for the water molecules. All MD simulations were performed by using the pmemd.cuda program in the Amber 18 software package. QM/MM MD simulations. A representative snapshot between 40-50 ns extracted from the classical MD trajectory was used for the subsequent QM/MM MD simulations, which was converged (Figures ) and close to the most populated structures from the cluster analysis (Figures ). All QM/MM Born Oppenheimer MD simulations were performed by using the CP2K 5.0 package, 62 along with the QM program QUICKSTEPS and the MM driver FIST. In this code, a real space multigrid technique is used to compute the electrostatic coupling between the QM and MM regions. The QM region was modeled at the DFT(UB3LYP-D3) level, employing the dual basis set formalism of Gaussian and plane-waves (GPW), whereas the remaining part of the system was modeled at the MM level using the same parameters as in the classical MD simulations. The Gaussian double-ζ valence polarized (DZVP) basis set was used to expand the wave function, while the auxiliary plane-wave basis set with a density cutoff of 360 Ry and GTH pseudopotentials was utilized to converge the electron density. The dangling bonds between the QM and MM regions were capped with hydrogen atoms. The QM region consists of the electron donor lignin, the electron acceptor, and the surrounding key residues, which has been specified for each step (see Figure ). To speed up the calculation of Hartree-Fock exchange within the B3LYP hybrid functional, the auxiliary density matrix method (ADMM) was used. All QM/MM MD simulations were performed under the NVT ensemble using a time-step of 0.5 fs. The system was equilibrated without any constraints for 2 ps. For the second ET step without tartaric acid (see Figure ), the QM/MM metadynamics calculations were performed to sample the free energy profile, in which the width of the Gaussian hills was set to 0.2 Å, the Gaussian height was set to 0.6 kcal/mol, while the time deposition interval between two consecutive Gaussians was set to 10 fs.