The study of ligand binding plays a pivotal role in the understanding of biochemical pathways and in drug design. Ligands, which can be small molecules, peptides, or other macromolecules, interact with specific targets, modulating their conformations, interactions and biological activities . Characterizing a ligand-receptor interaction is a three-level problem (structure, thermodynamics and kinetics) that involves the characterization of the bound configuration (the pose), of the stability of the interaction (the binding free energy), and of the kinetics of the process (the on and off rate constants) including the detailed binding mechanism. Although a wide range of experimental and computational methods have been developed to address such aspects , they must be scaled up to reduce time and associated costs, in particular when used for drug design. A common computational approach for drug design is based on the use of molecular docking methods , including recent machine learning approaches , to search and score the binding pose for large chemical libraries. The binding mechanism of top-scoring ligands can then be studied by a variety of methods , including endpoint approaches, such as the well-established thermodynamic integration (TI) , free energy perturbation (FEP) , and more recent single-simulation based or non-equilibrium approaches ; pathway-dependent calculations, such as long time-scale MD simulations , Markov state models ; and reaction coordinate based techniques like umbrella sampling and metadynamics (MetaD) . The intensive computational cost of these methods can be reduced by lowering the resolution of the structural representation of the ligand-receptor system (i.e. by coarse graining), but this usually leads to reduced accuracy . There is, therefore, a need to develop methods that achieve a more advantageous trade-off between speed and accuracy. Recently, we developed multi-eGO to simulate protein selfassembly processes . Multi-eGO simplifies the description of protein interactions by deriving effective Lennard-Jones (LJ) parameters (C6, C12 or ε, 𝜎 equivalently), by combining training conventional MD simulations of a system with prior models representing baseline interactions. We showed that multi-eGO provides accurate, atomic-level data on protein folding, intrinsically disordered proteins (IDPs) conformational dynamics and peptide aggregation, with highly reduced computational costs . Here we extended the applicability of multi-eGO to ligand binding simulations (Figure ). The training data consists of simulating a protein in its apo and holo states. The holo training is performed with additional free ligands in the box that allows simultaneous sampling of the bound pose and estimation of non-specific interactions (both protein-ligand and ligand-ligand). The apo training intramolecular contact probabilities are reweighted, following Bayes, by those obtained from a prior simulation, which represents a protein as a self-avoiding polymer. To reweight the protein-ligand and ligand-ligand contact probabilities, we introduce, as prior, a simulation performed at the same ligand concentration of the training data but with ligands and protein interacting only by hard-core repulsions (i.e. LJ C12 interactions). This prior simulation captures the rotational and translation entropy associated to the free ligand at the training set concentration. We expect this approach to have a twofold advantage: (i) the cost of training simulations, while not negligible, is often in part offset before free energy calculations are performed, to generate initial configurations and to test the stability of a binding pose; (ii) once the training cost is paid, the resulting multi-eGO model can be used to perform multiple simulations, using multiple free energy techniques as well as unbiased simulations, in a fraction of the time with respect to conventional force-fields (cf. Table ).