In a recent work of some of us, the structural properties and binding of PETase with tetrameric chains of PET (PET4, Figure ) and PEF (PEF4, Figure ) were investigated through an extensive application of molecular modeling techniques. Specifically, molecular docking, molecular dynamics (MD) simulations and unbinding free energy (FE) calculations indicated that PET4 is stabilized within the binding pocket of the enzymes thanks to the aromatic interactions established by the TPA moieties of PET4 with the W185 and Y87 sidechains. Conversely, PEF4 formed weak interactions with the enzymes, quickly dissociating from the PETase binding pocket during MD simulations. However, that preliminary assessment did not examine the reaction mechanisms of the enzyme with these substrates. The hydrolysis of PET by PETase has been undertaken by various studies that used multiscale quantum mechanics / molecular mechanics (QM/MM) MD simulations to investigate the elementary steps of the catalytic cycle. The process is initiated by the nucleophilic attack of S160 to the hydrolysable ester C-atom of PET (acylation). This stage is activated by the abstraction of a proton from the S160 hydroxyl side chain by H237 and leads to the formation of an acyl-enzyme intermediate. Then, a water molecule enters the active site and attacks the acyl-enzyme intermediate, completing the hydrolysis of PET and restoring the enzyme to its initial state (deacylation). These steps proceed through the formation of tetrahedral transition states where the carbonyl O-atom of the substrate is stabilized by the backbone of the M161 and Y87 (oxyanion hole) (Scheme 1). While some evidences from these studies pinpoint common features of the reaction mechanism, other results remain debated. In particular, the two-steps versus four-steps nature of the process, as well as identification of the rate-limiting step, are still open questions, thus highlighting the necessity of further investigation to deeper understand the underlying hydrolysis mechanism at the molecular level. Furthermore, none of the studies published so far delved into the assessment of the reaction mechanism of PEF depolymerization. Building on these findings, in this study we aim to address these gaps of knowledge by studying the catalysis of PET4 and PEF4 by PETase using QM/MM MD simulations. A crucial step in the study of chemical reactions is the definition of the reaction coordinate, which is defined by collective variables (CVs) capturing the slowest motions of the system, i.e. those functions of the atomic coordinates describing the slow degrees of freedom associated to the activated process of interest. At first glance, a straightforward strategy to select the CVs may seem using chemical intuition for the system under study. For example, the transition can be followed along specific bonds that are broken or formed during the reaction. Although this approach may provide a set of CVs that are physically and chemically meaningful, it inherently assumes that the reaction mechanism is known a priori, overlooking any other degree of freedom that might undergo changes during the transition. In the context of PETase reaction mechanism, some prior works estimated the FE landscapes using forming or breaking bonds as CVs, whereas others adopted more sophisticated strategies. García-Meseguer et al. utilized an adaptive string method to define the minimum FE path connecting reactants to the product, while Burgin and collaborators 30 used transition path sampling and inertial likelihood maximization to provide optimal mechanistic description of the acylation and deacylation steps, and validated the resulting reaction coordinates via committor probability test. In this work, we defined the reaction coordinate by using harmonic linear discriminant analysis (HLDA) , a supervised learning classification technique. Given a limited information about the two metastable states of interest in their local equilibrium states, namely reactants and products, HLDA allows designing an efficient one-dimensional CV in the form of a linear combination of a set of simple descriptors, such as bond distances or coordination numbers. A well-tempered metadynamics (WTMD) bias has been applied along these CVs to help the system to overcome the FE barriers and favor transitions from reactants to products. To gather exhaustive statistics at an affordable computational cost, the QM region has been treated with the third-order density functional tight binding (DFTB3) Hamiltonian. This level of theory found successful application in many contexts, but it does not guarantee the same accuracy as density functional theory (DFT) or that of even higher level techniques such as post-Hartree-Fock (HF) methods. Unfortunately, despite the increasing power of computational infrastructures, running long simulations with such high-level methods is still unaffordable, forcing one to find compromises among the accuracy, size of the QM region and timescale of the simulations.