We subject a series of five protein–ligand systems
which
contain important SARS-CoV-2 targets, 3-chymotrypsin-like protease
(3CLPro), papain-like protease, and adenosine ribose phosphatase,
to long time scale and adaptive sampling molecular dynamics simulations.
By performing ensembles of ten or twelve 10 μs simulations for
each system, we accurately and reproducibly determine ligand binding
sites, both crystallographically resolved and otherwise, thereby discovering
binding sites that can be exploited for drug discovery. We also report
robust, ensemble-based observation of conformational changes that
occur at the main binding site of 3CLPro due to the presence of another
ligand at an allosteric binding site explaining the underlying cascade
of events responsible for its inhibitory effect. Using our simulations,
we have discovered a novel allosteric mechanism of inhibition for
a ligand known to bind only at the substrate binding site. Due to
the chaotic nature of molecular dynamics trajectories, regardless
of their temporal duration individual trajectories do not allow for
accurate or reproducible elucidation of macroscopic expectation values.
Unprecedentedly at this time scale, we compare the statistical distribution
of protein–ligand contact frequencies for these ten/twelve
10 μs trajectories and find that over 90% of trajectories have
significantly different contact frequency distributions. Furthermore,
using a direct binding free energy calculation protocol, we determine
the ligand binding free energies for each of the identified sites
using long time scale simulations. The free energies differ by 0.77
to 7.26 kcal/mol across individual trajectories depending on the binding
site and the system. We show that, although this is the standard way
such quantities are currently reported at long time scale, individual
simulations do not yield reliable free energies. Ensembles of independent
trajectories are necessary to overcome the aleatoric uncertainty in
order to obtain statistically meaningful and reproducible results.
Finally,
we compare the application of different free energy methods to these
systems and discuss their advantages and disadvantages. Our findings
here are generally applicable to all molecular dynamics based applications
and not confined to the free energy methods used in this study.