Preparation of a multi-eGO model for lysozyme:benzene: all-atom training. To train a multi-eGO model for LYZ:BNZ binding, we performed multiple apo and holo simulations using the DES-Amber 48 and the GAFF2 force field (see Methods and Table ). With respect to the apo state, the simulations reproduced a breathing motion of the two domains over the microsecond time scale that is well characterized experimentally by residual dipolar couplings by nuclear magnetic resonance (NMR) spectroscopy (Fig. ), as shown by the distribution of the root mean square deviation of the atomic positions (RMSD) calculated with respect to the backbone of a reference structure (PDB 1L84). The simulations of the holo state were performed with one BNZ molecule in the binding cavity and four additional molecules in solution. We found that BNZ is more likely to be found in regions surrounding the binding cavity, suggesting a role of surface interactions in directing BNZ to its binding site (Fig. ). !C254D785EF8<4 of 0.53 kJ/mol. The resulting multi-eGO model reproduced the intermolecular contact probabilities between LYZ and BNZ (Fig. ), demonstrating that multi-eGO can effectively learn heterogeneous interactions. As a first test, we performed 40 unbiased binding simulations, starting from random configurations of 4 unbound BNZ molecules at the same concentration as the training simulation (see SI Movie 1 for a representative trajectory). Simulations were stopped when the minimum distance between BNZ atoms and the C⍺ atom of A99 in the binding pocket fell below 0.4 nm. The extracted binding times were then fitted to a Poisson distribution, yielding a mean binding time 𝜏 of 963±20 ns with a p-value of 0.988 obtained from a Kolmogorov-Smirnov (KS) test , corresponding to a kon of 7.6x10 7 M -1 s -1 (see Fig. ). This rate should be considered as nominal due to the simplified nature of the model, but when compared with the experimental value 67 of 8x10 5 -10 6 M -1 s -1 it gives an indication of the intrinsic acceleration of multi-eGO due to removal of degrees of freedom. Thermodynamic integration and volume-based metadynamics on lysozyme:benzene multi-eGO model. Given the efficiency of multi-eGO, its accuracy in estimating the LYZ:BNZ binding free energy was tested, in comparison with alternative methods. We focused on Thermodynamic Integration (TI) and Volume Based Metadynamics (VMetaD), the former being the industry standard for such calculations and the latter, a potentially more informative but also more computationally expensive alternative . TI of LYZ:BNZ was performed in 6 replicates using three alternative restraints (CL1, CL2, and CL3) to keep the ligand in the binding pose and to speed up convergence (see Methods and Fig. ). The unbinding free energies, calculated as the mean and the standard deviation of the mean, are shown in Fig. . The three estimates of 21.41±0.04 kJ/mol, 21.45±0.05 kJ/mol, and 21.46±0.07 kJ/mol from the three restraints, respectively, show a high statistical precision and a remarkable accuracy when compared to the experimental unbinding free energy of 21.7±0.8 kJ/mol (Fig. ). An additional unbinding free energy estimate was obtained using VMetaD (see Methods) by running four independent replicates of 1.5 µs each, reaching convergence after approximately 200 ns (Fig. ). A free energy surface (FES) representing the unbinding process was obtained by reweighting VMetaD simulations as a function of the distance of the BNZ center of mass from the center of the sphere centered on the binding site, which was used to define the sample volume, and the coordination number of BNZ with the residues of LYZ (see Fig. and Methods section). From the FES, we estimated a ΔG of 21.4±0.2 kJ/mol, obtained as the mean and standard deviation of the mean over the 4 replicates (Fig. ), which is within the statistical precision of the TI estimates, demonstrating the reliability of the two different approaches. The use of VMetaD also allowed us to compare our sampled binding/unbinding pathways with those already reported in the literature (Fig. and Fig. ). Using a dynamic timewarping clustering algorithm , we were able to observe most, if not all (depending on the classification approach), of the previously observed pathways , further supporting the ability of multi-eGO to describe correctly the binding and unbinding mechanism of a ligand. and we can define all the six main unbinding pathways. 1) BNZ passes behind D and under C; 2) BNZ passes over C and behind F; 3) BNZ passes over C and under F; 4) BNZ passes in between the top part of H and G; 5) BNZ passes behind H; 6) BNZ passes under the C-terminal domain, between helices G and D. These pathways qualitatively match the ones identified with explicit water and atomistic potentials in ref. .