is estimated, where i,j are atom indices, 𝜀 ' is the energy scale parameter, and 𝑃 234563789 %& is the minimum probability used for regularization. A detailed description of the model can be found in ref. and the associated code and parameters are available on GitHub. For each apo system, we performed a prior simulation consisting of a self-avoiding chain, obtained with bonded interactions and C12 repulsion from the multi-eGO transferable potential. For multi-domain proteins like LYZ, a second prior simulation was introduced. This consists of a multi-eGO simulation in which the inter-domain attractions, namely between residues 1-71 and 72-162 are turned off, while preserving intra-domain and local interactions, in the alpha-helix connecting the two domains, to maintain its integrity. The two sets of prior probabilities allowed interand intra-domain interactions to be correctly decoupled and to assign them specific energy scale (ε0). For intermolecular interactions, the reference prior simulation corresponds to protein:ligand complex at the same concentration of the training simulation (considering only unbound ligands) with only hard-core repulsions between protein and ligands. This allows the concentration-dependent, rototranslational entropy of the ligands to be accounted for by removing the effect of box size on the contact probability estimation. Multi-eGO ligands were parametrized starting from the training topology, removing hydrogens and optimizing the bonded parameters to reproduce the local geometries of the training simulations. All multi-eGO MD simulations were performed using stochastic dynamics integration at a 5 fs timestep and 25 ps relaxation time. The cutoff for the LJ interactions was set to 2.5σmax, corresponding to 1.44 nm. A 10% larger radius was used for the neighbor lists, which were updated every 20 steps. Different values of ε0 were tested to maximize the overall agreement between training and multi-eGO simulations. All simulations were run using GROMACS 2022 . Thermodynamic Integration. For the multi-eGO unbinding free energy estimation, we used TI and the Bennett's accepting ratio (BAR) algorithm . For each system a set threecross linked restraints, namely CL1, CL2, and CL3, to keep the ligand in the binding pose during the decoupling. The effect of this restraining potential can be removed analytically a posteriori 39 . The thermodynamic cycle was computed through a set of λ values scaling the restraint and non-bonded interactions (Table ). For LYZ:BNZ, a steepest descent energy minimization was performed at each λ, followed by a 500 ps relaxation and a 1 ns run at 300K in the NVT ensemble. For SRC:DAS and SRC:PP1, a steepest descent energy minimization was performed at each λ, followed by a 2 ns relaxation and a 1 ns run at 310K in the NVT ensemble. All simulations were run using GROMACS 2022 . Volume-based Metadynamics. VMetaD 38 simulations LYZ:BNZ were run using well-tempered metadynamics 60 considering 3 CVs, namely the relative position in spherical coordinates (r, q, j) of the benzene center of mass, with respect to the center of mass of the lysozyme binding domain (residues 72-162). We set an initial height of the Gaussians of 1.2 kJ/mol, the widths at 1 Å, π/16 rad, and π/8 rad for r, q, and j, respectively and a bias factor of 20, depositing a Gaussian every 10 ps. 4 different 1.5 µs-long replicas of the simulations were run, employing the same MD parameters used in the unbiased multi-eGO simulations. To avoid an extensive sampling time of the unbound state, a repulsive spherical potential was set at a distance of 𝜌 6:3 (here set at 30 Å) from the center of mass of LYZ: