Harmonic Linear Discriminant Analysis (HLDA) systematically constructs the CV for the process of interest as a linear combination of an arbitrary number of descriptors of the system. This method uses as an input the equilibrium distributions of a set of descriptors at the endpoints of each transition, i.e., in this case, the reactants and products of the acylation and deacylation steps) and looks for the direction that maximizes the distance between the two states' distributions (or classes) by using a modified version of Fisher's linear discriminant analysis (LDA) with harmonic weighting. The method is described in detail in Section S1. The solution of the eigenvalues problem reported in eq. (S1) provides the coefficients of the linear combination of the starting descriptors (d1-11, c1-3 in Figure ). This defines a CV along which the projections of the reactant and product distributions are well separated, enabling highly efficient and accurate sampling that better distinguishes the two states in the FE space. In principle, each degree of freedom of the system can be included in the definition of the reaction coordinate, so that any a priori knowledge on the process is required. Indeed, the coefficients provided by HLDA describe the relative weight of the descriptors, indicating the relevance of each feature in driving the system from the reactants to the products. The optimization of the HLDA variables from the descriptor distributions was performed with the python script deposited on the PLUMED-NEST repository (www.plumed-nest.org), with plumID:19.042. Descriptors for the acylation and deacylation steps.