Why is it so difficult to converge MM/GBSA free energies of binding for cyclodextrin inclusion complexes?

The molecular mechanics/generalized Born surface area (MM/GBSA) method is a popular choice to compute free energies of binding, ΔG_bind, where a small ligand (L) binds to a larger host (H) to form an inclusion complex LH:

L+H⇌LH

While most of applications of MM/GBSA have historically been centered around protein-ligand complexes, there are also examples in the literature of cyclodextrin (CD) inclusion complexes being studied with this method.

Typical simulation protocols used for CD inclusion complexes were, however, originally designed and used for protein-ligand complexes. When applied to CD inclusion complexes, the resulting free energies of binding are pretty much hit-and-miss in terms how well the simulations compare with experimental free energies of binding. The origin of this problem cannot easily be traced to one consistent failure of the MM/GBSA model but is more a combination of multiple issues with the simulation protocol: type of molecular dynamics simulation, equilibration length (t_eq), length of production run (t_prod) and sampling time (t_s) all of which have been studied extensively for protein-ligand complexes.

In this work, I present insight into why it can be quite demanding to obtain free energies of binding with the MM/GBSA method that are statistically converged for CD inclusion complexes. The presented results will also be accompanied by recommendations on what a consistent simulation protocol at minimum should include for MM/GBSA when applied to CD inclusion complexes.