Following the completion of the MD production run for the original model, well-tempered metadynamics was applied, using PLUMED, to eject two or three ligands according to the scheme of Figure : this procedure produced the two-ligand and three-ligand models, respectively. The distance between the center of mass (COM) of the binding residues and that of the ligand, and the solvation were employed as collective variables (as described in Ref. ), with a bias factor of 15 and a deposition rate of one Gaussian per ps to gradually eject the ligand; the Gaussians had a width of 0.15 Å and 0.03 respectively and an initial height of 1.2 kJ/mol. The two-and three-ligand models were subjected to a 100 ns MD simulation to equilibrate the structure after ligand ejection, followed by a 1000 ns MD production run for analysis.