Conventional MD simulations were performed using the AMBER 2020 32 simulation package, along with the ff14SB force field and the TIP3P water model for the parameters of protein atoms and water molecules, respectively. For each of the eight simulation systems, the protonation state of each residue at pH 9.0 was estimated with the H++ server, in order to match the experimental conditions in enzymatic degradation assays using either PET films or granules. The consistency of the estimated protonation states for the H/F and S/I variants of a given enzyme were manually checked. Then, all systems were solvated in truncated octahedral water boxes with 15 Å of padding and neutralized with Na + or C l-counterions. A standard minimization, heating, and 10 ns equilibration protocol was performed, followed by NPT production at 1 bar and the experimentally determined optimal temperature of catalytic activity for each system (i.e. 303 K for IsPETase, 308 K for BurPL, 323 K for LCC, and 333 K for TfCut). Constant temperature and pressure conditions for each system were controlled using a Langevin thermostat 36 with a collision frequency of 1 ps -1 , and a Berendsen barostat 37 with a pressure relaxation time of 1 ps. MD simulations were run using a time step of 2 fs along with the SHAKE 38 algorithm to constrain hydrogen-containing bonds, a 12 Å cut-off for nonbonded interactions, and the particle-mesh Ewald 39 method for long-range electrostatics. Five 500-ns-long replicas were produced for each of the four H/F and four S/I systems, for a total of 20 𝜇s of MD simulations. Convergence of the simulations is shown in Figure .