Structure-based models have been introduced to study protein folding, based on the hypothesis that the native protein structure should capture the most relevant interactions across the whole free energy landscape of a protein. They have proven useful to study otherwise, inaccessible folding mechanisms at high resolution . Recently, we have introduced multi-eGO, an ensemble-based model rooted in Bayesian statistics, where one or more conformational ensembles of a protein, representing the fluctuations of relevant free energy minima, are weighted with prior information . This approach enables the generation of a model that can simultaneously describe folded and disordered proteins, in addition to self-assembly processes such as peptide amyloidogenic aggregation. Here, we have extended the multi-eGO approach to include small molecules. Our initial results on four case studies allow us on the one hand to propose multi-eGO as a model for benchmarking and developing in silico techniques for ligand binding studies, and on the other hand to suggest that multi-eGO can be used for quantitative estimation of ligand affinities. This is because multi-eGO describes the system at atomic resolution, with binding free energies compatible with those measured in vitro, but with a dramatic acceleration of the kinetics (that we can tentatively quantify between 2 and 4 orders of magnitude). We exploited the LYZ:BNZ system to apply orthogonal free energy calculation techniques (TI with multiple restraining conditions and VMetaD) showing their applicability and robustness on the model. It must be underlined that all simulations are replicated multiple times to obtain statistically meaningful estimates, at minimal computational costs. While Aβ42:10074-G5 represented a case that is inaccessible to current conventional approaches, namely explicit titration as a function of ligand concentration, with direct comparisons with in vitro concentrations. The latter case, on the one hand, demonstrated the efficiency of multi-eGO and, on the other hand, underlined the importance of correcting single-molecule simulations for the finite-size effect, following procedures such as those described in Lopez et al. The use of multi-eGO requires model training, which consists of running training simulations and finding a set of free parameters. Here, we introduced a possible strategy to set our energy scale parameters, excluding the intervention of any possible user bias. In practice, we minimized the RMSE with the protein fluctuation as well as the RMSE for the intermolecular residue-wise contact map, comparing the multi-eGO and training simulations. This approach was first successfully tested on the LYZ:BNZ system and then successfully replicated for the SRC:DAS case. Successfully, because in both cases, knowledge of the experimental binding free energy was not used in the input and instead was accurately obtained by the simulations. In both cases, our training simulations include a ligand in its binding pose as well as several ligands free to sample the surface of the protein. The latter data are used to set the energy scale of the system. For Aβ42:10074-G5, we used a previously generated, large dataset of simulations, from which we inherited an overestimated dissociation constant. This is where the advantage of multi-eGO becomes apparent, as it is sufficient to tune a single parameter to make the model match the experimental value. This allows us to speculate that multi-eGO may be particularly suitable for studying ligand binding processes such as those represented by the LYZ:BNZ, SRC:DAS, and SRC:PP1 cases, where the training data can be generated with very limited computational cost and should represent only the weak (i.e., faster to sample) non-specific interactions of a ligand with the protein surface, given knowledge or a hypothesis about the binding pose. We anticipate that further work will be needed to better understand the strengths and limitations of multi-eGO in ligand binding studies, including cases with more complex binding mechanisms such as induce-fit or conformational selection, as well as to streamline and possibly automatize the simulation setup protocol.