Interactions of polyelectrolytes (PEs) with proteins play a crucial role in numerous biological processes, such as the internalization of virus particles into host cells. Although docking, machine learning methods, and molecular dynamics (MD) simulations are utilized to estimate binding poses and binding free energies of small-molecule drugs to proteins, quantitative prediction of the binding thermodynamics of PE-based drugs presents a significant obstacle in computer-aided drug design. This is due to the sluggish dynamics of PEs caused by their size and strong charge–charge correlations. In this paper, we introduce advanced sampling methods based on a force-spectroscopy setup and theoretical modeling to overcome this barrier. We exemplify our method with explicit solvent all-atom MD simulations of the interactions between anionic PEs that show antiviral properties, namely heparin and linear polyglycerol sulfate (LPGS), and the SARS-CoV-2 spike protein receptor binding domain (RBD). Our prediction for the binding free-energy of LPGS to the wild-type RBD matches experimentally measured dissociation constants within thermal energy, kBT, and correctly reproduces the experimental PE-length dependence. We find that LPGS binds to the Delta-variant RBD with an additional free-energy gain of 2.4 kBT, compared to the wild-type RBD, due to the additional presence of two mutated cationic residues contributing to the electrostatic energy gain. We show that the LPGS–RBD binding is solvent dominated and enthalpy driven, though with a large entropy–enthalpy compensation. Our method is applicable to general polymer adsorption phenomena and predicts precise binding free energies and reconfigurational friction as needed for drug and drug-delivery design.
