Parameters for substrates 1 and azide were generated with the antechamber and parmchk2 modules of AMBER20 66 using the 2nd generation of the general amber force-field (GAFF2). Partial charges (RESP model) were set to fit the electrostatic potential generated at the HF/6-31G(d) level of theory. The charges were calculated according to the Merz-Singh-Kollman 70,71 scheme using the Gaussian16 software package. The protonation states were predicted using PROPKA. The enzyme structure was obtained from the PDB with the code 5o30 and cleaned from other non-peptidic molecules to obtain the wild-type system in a tetrameric oligomerization state. The single mutation M45F was introduced using the Pymol mutagenesis tool. Proteins were solvated in a pre-equilibrated truncated octahedral box of 12 Å edge distance using the OPC water model, resulting in the addition of ca. 21.300 water molecules, and neutralized by the addition of explicit counterions (i.e., Na + ) using the AMBER20 leap module. All MD simulations were performed using the amber19 force field (ff19SB) in our in-house GPU cluster, GALATEA. The Pmemd.cuda program from Amber20 was used to perform a two-stage geometry optimization. In the first stage, solvent molecules and ions were minimized, while solute molecules were restrained using 500 kcal•mol -1 •Å-2 harmonic positional restraints. In the second stage, an unrestrained minimization was performed. The systems were then gradually heated by increasing the temperature by 50 K during six 20 ps sequential MD simulations (0-300 K) under constant volume. Harmonic restraints of 10 kcal•mol -1 •Å-2 were applied to the solute, and the Langevin equilibration scheme was used to control and equalize the temperature. The time step was kept at one fs during the heating stages to allow potential inhomogeneities to self-adjust. Each system was then equilibrated without restraints for 2 ns at a constant pressure of 1 atm and temperature of 300 K using a 2 fs time step in the isothermal-isobaric ensemble (NPT). After equilibration, five replicas of 250 ns were run for each system (i.e., 1.25 μs per system and 5 μs in total simulated time) in the canonical ensemble (NVT). MD simulations were analyzed by monomers to make it easier to study, multiplying the simulated time by four. All analysis was done using available Python libraries (pyemma , mdtraj , and mdanalysis ) in a jupyter lab environment.