Drug-target residence time, the length of time for which a small molecule stays bound to its receptor target, has increasingly become a key property for optimization in drug discovery programs. However, its in silico prediction has proven difficult. Here we describe a method, using atomistic ensemble-based steered molecular dynamics (SMD), to observe the dissociation of ligands from their target G protein-coupled receptor in a timescale suitable for drug discovery. These dissociation simulations accurately, precisely, and reproducibly identify ligand-residue interactions and quantify the change in ligand energy values for both protein and water. The method has been applied to 17 ligands of the A2A adenosine receptor, all with published experimental kinetic binding data. The residues that interact with the ligand as it dissociates are known experimentally to have an effect on binding affinities and residence times. There is a good correlation (R2 = 0.79) between the computationally calculated change in water-ligand interaction energy and experimentally determined residence time. Our results indicate that ensemble-based SMD is a rapid, novel and accurate empirical method for the determination of drug-target relative residence time.

The accurate identification of the specific points of interaction between G protein-coupled receptor (GPCR) oligomers is essential for the design of receptor ligands targeting oligomeric receptor targets. A coarse-grained molecular dynamics computer simulation approach would provide a compelling means of identifying these specific protein–protein interactions and could be applied both for known oligomers of interest and as a high-throughput screen to identify novel oligomeric targets. However, to be effective, this in silico modeling must provide accurate, precise, and reproducible information. This has been achieved recently in numerous biological systems using an ensemble-based all-atom molecular dynamics approach. In this study, we describe an equivalent methodology for ensemble-based coarse-grained simulations. We report the performance of this method when applied to four different GPCRs known to oligomerize using error analysis to determine the ensemble size and individual replica simulation time required. Our measurements of distance between residues shown to be involved in oligomerization of the fifth transmembrane domain from the adenosine A_2A receptor are in very good agreement with the existing biophysical data and provide information about the nature of the contact interface that cannot be determined experimentally. Calculations of distance between rhodopsin, CXCR4, and β₁AR transmembrane domains reported to form contact points in homodimers correlate well with the corresponding measurements obtained from experimental structural data, providing an ability to predict contact interfaces computationally. Interestingly, error analysis enables identification of noninteracting regions. Our results confirm that GPCR interactions can be reliably predicted using this novel methodology.

G-protein coupled receptors (GPCRs) are the largest superfamily of membrane proteins, regulating almost every aspect of cellular activity and serving as key targets for drug discovery. We have identified an accurate and reliable computational method to characterise the strength and chemical nature of the inter-helical interactions between the residues of transmembrane (TM) domains during different receptor activation states, something that cannot be characterised solely by visual inspection of structural information. Using the fragment molecular orbital (FMO) quantum mechanics method to analyse 35 crystal structures representing different branches of the class A GPCR family, we have identified 69 topologically-equivalent TM residues that form a consensus network of 51 inter-TM interactions, providing novel results that are consistent with and help to rationalise experimental data. This discovery establishes a comprehensive picture of how defined molecular forces govern specific inter-helical interactions which, in turn, support the structural stability, ligand binding and activation of GPCRs.

We apply the hit-to-lead ESMACS (enhanced sampling of molecular dynamics with approximation of continuum solvent) and lead-optimization TIES (thermodynamic integration with enhanced sampling) methods to compute the binding free energies of a series of ligands at the A1 and A2A adenosine receptors, members of a subclass of the GPCR (G protein-coupled receptor) superfamily. Our predicted binding free energies, calculated using ESMACS, show a good correlation with previously reported experimental values of the ligands studied. Relative binding free energies, calculated using TIES, accurately predict experimentally determined values within a mean absolute error of approximately 1 kcal mol−1. Our methodology may be applied widely within the GPCR superfamily and to other small molecule–receptor protein systems.