The main advantage of using HLDA is that the definition of the reaction coordinate is not driven by any a priori information on the mechanism under investigation, apart from the initial and final states of the process. Indeed, the coefficients of the linear combination indicate the relative weight of each descriptor (feature importance), suggesting how changes in that specific degree of freedom will drive the system toward one of the two endpoints of the transition. The coefficients obtained by solving eq. (S1) depend on the scatter matrices 𝑺 : and 𝑺 ; of the distributions collected during the unbiased QM/MM MD simulation, which are shown in Figure . The larger the separation between the distributions obtained for the two endpoints, the greater is the relevance of the descriptor in separating the two states of the system. Conversely, highly overlapping distributions indicate that the descriptor is not significantly different in the two classes and its relative weight in the linear combination will be negligible. Concerning the acylation of the PETase-PET4 system (Figure ), it can be observed that the profiles for distances d1, d2, d3, and d6 are neatly separated. For each of these descriptors, one of the two distributions of either the reactant (orange) or the product (blue) is sharp and characterized by a small variance because the atoms considered for these distances are involved in a covalent bond. Conversely, larger fluctuations are observed in the other state, in which these bonds are broken. From the coefficients reported in Table , the highest weighted feature is d2 (-0.817), describing the ester bond of PET that is cleaved upon the nucleophilic attack by S160. The second highest-weighted descriptor is d3 (-0.494), that accounts for the distance between the Og atom of S160 and the associated Hg atom (Figure ). Both these bonds are broken during acylation. Indeed, d2 and d3 have the same sign, meaning that and elongation or contraction of these distances drive the system in the same direction along the reaction coordinate. On the other hand, d1 and d6 correspond to bonds that are formed during the acylation step (Scheme 1A). The sign of the coefficients is positive indicating that changes in these distances guide the system towards the opposite endpoint than d2 and d3. From the coefficients obtained for the deacylation of the PETase-PET4 acyl-enzyme intermediate (Table ), it can be observed that d11 is the main driver of the process with a weight of -0.922. This distance accounts for the formation of the new covalent bond between the O atom of water and the C5 atom of PET4. This process is accompanied by a deprotonation of water, whose coordination number (c1) is expected to decrease consistently with d11, thus displaying coefficients with the same sign. The other relevant descriptor of the process is d1 (0.353), that corresponds to the distance between the C5 atom of PET and the Og atom of S160 (Figure ). These atoms form a covalent bond in the acyl-enzyme intermediate, that is broken by the nucleophilic attack of water during the deacylation, with S160 intercepting one proton from water (Scheme 1C). Consequently, the coefficients of d1 and c3 vary in a consistent manner relative to each other and drive the reaction in an opposite direction to d11 and c1. The other descriptors are not relevant in separating the two FE basins of the deacylation step, as the relative distributions overlaps (Figure ) and their coefficients in the HLDA CV are close to zero. As expected, the results obtained for the PETase-PEF4 system (Figure ) resemble those obtained for PET4. In the acylation step, the highest-weighted descriptors are d2 and d3, while d1 and d6 exhibit opposite signs. Moreover, the contributions of d4, d5, and d8 are more relevant in this system compared to PETase-PET4, although they remain secondary relative to d2 and d3. Finally, as with PET4, the deacylation is predominantly governed by d11, bearing a coefficient of -0.955, followed by d1, which is weighted 0.291.