Introducing Free Energy Perturbation: An effective method for calculating and predicting ligand binding affinities
GPCRs (G-protein-coupled receptors) represent the largest family of membrane receptors currently targeted by approved drugs. It is estimated that approximately 700 approved drugs target GPCRs representing about 35% of all approved drugs.1 Free Energy Perturbation (FEP) calculations, a highly effective computational method, can accurately predict ligand binding affinities to GPCR targets. The binding affinity predictions (ΔG) can be within 1 kcal/mol from the corresponding experimental measurements. Relative FEP is particularly suitable for accurately calculating ΔG of binding for a congeneric series. The method computes the free energy difference of transforming one molecule to another. FEP methods can help prioritize compounds for synthesis and greatly decrease the time and cost involved in identifying an active lead compound, thus expediting, and increasing the efficiency of designing therapeutics for these targets.
A GPCR case study: P2Y1
The P2Y1 purinergic receptor, a class A 7-transmembrane GPCR, is the specific example GPCR presented and discussed here. P2Y1 is of significant interest as a drug target for diabetes, as it is involved in the regulatory mechanism for insulin secretion2 and for hemostasis and thrombosis as it plays a role in mediating platelet activation.3 Drugs can bind to different regions of P2Y1, as in addition to an extracellular binding site, P2Y1 also has a transmembrane exposed site on the extra-helical bundle (as shown in Figure 1). This second ligand binding site, which lies between the membrane and the protein provides a new challenge and opportunity for computational ligand design.
Figure 1: PDB: 4XNV- crystal structure of GPCR P2Y1 with ligand BPTU shown in a POPC membrane, minimized after molecular dynamics (MD). The zoomed in box shows the interactions that the ligand forms with the protein (hydrophobic interactions: grey, π-π staking: purple).
Performing the FEP calculations
In this computational FEP study, we have illustrated the power of this method, by calculating the binding affinity of 30 P2Y1 ligands for which experimental data is available, to make an effective assessment of the accuracy of the method. This computational study used Cresset’s Flare™ FEP software.4 Taking reference from the published work of Dickson5 and Ross,6 it involved predicting the activities of a 30-compound series of congeneric antagonists with a common substructure (shown in Figure 2). This antagonist series was designed starting from BPTU which is a potent (Ki = 16 nM, DGbinding = -11.2 kcal/mol) P2Y1 antagonist, whose crystal structure has been solved with P2Y1 (PDB: 4XNV: Figure 1). This ligand-GPCR crystal structure was used as the reference for the alignment of the congeneric series of ligands.
Figure 2: (A) The P2Y1 common antagonist substructure and (B) the crystallographic ligand of P2Y1 PDB:4XNV complex, BPTU, which was used as a reference to align 29 ligands in Flare.
Careful system preparation is required for FEP studies and there are several steps involved in order to obtain accurate results. MD and water analysis studies (using a technique called 3D-RISM7 to identify waters which play a role in the ligand binding site) were performed first, to appropriately prepare the P2Y1 crystal structure system for FEP studies. Placement of waters is challenging in membrane proteins which are often thought to be largely dry: however, key water molecules may still play a part, and a missing or incorrectly placed water can cause loss of accuracy. MD studies were performed to ensure the system was stabilized, essential protein-ligand interactions were conserved during dynamics, and the ligand bioactive binding mode was captured. In this case the P2Y1 membrane crystal structure with the co-crystalized ligand were placed in a POPC lipid bilayer and were solvated with TIP3P water molecules. The complete system on which 20 ns MD is performed, contained 93,255 atoms. The OpenFF 2.0 was used to model the small ligand molecule, and the AMBER ff14SB was used to model the protein and POPC lipid membrane bilayer system. The calculations were performed using an NPT ensemble with a 4 fs timestep at 298 K. The MD performed did not show much movement in the ligand coordinates from the original crystal structure indicating the starting structure for our FEP studies was stable.
Predicting binding affinities: FEP results and conclusions
The 30-compound antagonist dataset, for which we want to predict activities, was aligned to the crystal structure ligand, BPTU. The ligands were aligned based on their maximum common substructure and Cresset electrostatic ligand field points.8 This results in a highly compact alignment for each of the ligands in the dataset. Chiral, rotamer and tautomer states were verified: the overlayed ligand positions were used as the starting point for FEP ligand binding calculations. For the FEP experiments, the perturbation network was created with dual way perturbations (Figure 3) and adaptive Lambda scheduling was used. The adaptive scheduling has been uniquely developed and implemented by Flare FEP to automatically and efficiently optimize the number of Lambda windows required for each perturbation.
Figure 3: The FEP graph links pairs of ligands and defines the transformations. If a transformation is predicted to be computationally difficult, Flare FEP offers the user an option to automatically insert an ‘intermediate’ molecule which will be highlighted in the graph in grey.
The results of the FEP calculations showed an R2 = 0.57, MUE = 1.18, a Kendall’s Tau = 0.56, and RMSDpw = 1.67 (Figure 4). In essence this means that most of the activities were predicted accurately, to about 1.1 kcal/mol from experiment. At the same time, it informs us that in a prospective experiment with new molecular designs, a similar level of accuracy and correlation between the predicted and experimental data can be expected. However, it is essential to consider the specific characteristics of the new designs and the quality of the data used for calibration. The FEP method's performance may vary depending on the complexity of the transformations and the similarity between the new designs and the molecules used for calibration.
Figure 4: Plot of the predicted activities versus the experimental activities for the 30 P2Y1 ligand set.
Flare FEP provided good results with ligand binding affinities calculated by FEP in agreement with those reported in the Dickson paper.5 Flare FEP can relatively quickly, easily, and accurately perform binding affinity predictions for membrane proteins, even in potentially difficult cases such as this allosteric, lipid exposed binding example P2Y1. Now that this method has been verified, it can be used to predict the affinity of newly designed ligands. Such a workflow would result in time, effort, and energy saved, since only those novel compounds predicted to have high binding affinity will need to be synthesized and experimentally tested, overall avoiding time wasted on compounds with less potential.
References
- Sriram, K.; Insel, P. A. G Protein-Coupled Receptors as Targets for Approved Drugs: How Many Targets and How Many Drugs? Molecular Pharmacology 2018, 93 (4), 251–258. https://doi.org/10.1124/mol.117.111062.
- Marcheva, B.; Weidemann, B. J.; Taguchi, A.; Perelis, M.; Ramsey, K. M.; Newman, M. V.; Kobayashi, Y.; Omura, C.; Manning Fox, J. E.; Lin, H.; Macdonald, P. E.; Bass, J. P2Y1 Purinergic Receptor Identified as a Diabetes Target in a Small-Molecule Screen to Reverse Circadian β-Cell Failure. eLife 2022, 11. https://doi.org/10.7554/elife.75132.
- Agnès Ribes; Garcia, C.; Marie-Pierre Gratacap; Evi Kostenis; Martinez, L. O.; Payrastre, B.; Jean-Michel Sénard; Céline Galès; Pons, V. Platelet P2Y1 Receptor Exhibits Constitutive G Protein Signaling and β-Arrestin 2 Recruitment. BMC Biology 2023, 21 (1). https://doi.org/10.1186/s12915-023-01528-y.
- Flare™, version 7, Cresset®, Litlington, Cambridgeshire, UK; http://www.cresset-group.com/flare/; Cheeseright T., Mackey M., Rose S., Vinter, A.; Molecular Field Extrema as Descriptors of Biological Activity: Definition and Validation J. Chem. Inf. Model. 2006, 46 (2), 665-676; Bauer M. R., Mackey M. D.; Electrostatic Complementarity as a Fast and Effective Tool to Optimize Binding and Selectivity of Protein–Ligand Complexes J. Med. Chem. 2019, 62, 6, 3036-3050; Maximilian Kuhn, Stuart Firth-Clark, Paolo Tosco, Antonia S. J. S. Mey, Mark Mackey and Julien Michel Assessment of Binding Affinity via Alchemical Free-Energy Calculations J. Chem. Inf. Model. 2020, 60, 6, 3120–3130
- Dickson, C. J.; Viktor Horn̆ák; Duca, J. S. Relative Binding Free-Energy Calculations at Lipid-Exposed Sites: Deciphering Hot Spots. Journal of Chemical Information and Modeling 2021, 61 (12), 5923–5930. https://doi.org/10.1021/acs.jcim.1c01147.
- Ross, G. A.; Chao Lü; Guido Scarabelli; Albanese, S. K.; Houang, E. M.; Abel, R.; Harder, E.; Wang, L. The Maximal and Current Accuracy of Rigorous Protein-Ligand Binding Free Energy Calculations. Research Square (Research Square) 2022. https://doi.org/10.21203/rs.3.rs-2179899/v1.
- Skyner, R. E.; McDonagh, J. L.; Groom, C. R.; van Mourik, T.; Mitchell, J. B. O. A Review of Methods for the Calculation of Solution Free Energies and the Modelling of Systems in Solution. Physical Chemistry Chemical Physics 2015, 17 (9), 6174–6191. https://doi.org/10.1039/c5cp00288e.
- Cheeseright, T.; Mackey, M.; Rose, S.; Vinter, A. Molecular Field Extrema as Descriptors of Biological Activity: Definition and Validation. Journal of Chemical Information and Modeling 2006, 46 (2), 665–676. https://doi.org/10.1021/ci050357s.