The understanding of binding interactions between any protein and a small molecule plays a key role in the rationalization of affinity and selectivity and is essential for an efficient structure-based drug discovery (SBDD) process. Clearly, to begin SBDD, a structure is needed, and although there has been fantastic progress in solving G-protein-coupled receptor (GPCR) crystal structures, the process remains quite slow and is not currently feasible for every GPCR or GPCR–ligand complex. This situation significantly limits the ability of X-ray crystallography to impact the drug discovery process for GPCR targets in ‘real-time’ and hence there is still a need for other practical and cost-efficient alternatives. We present here an approach that integrates our previously described hierarchical GPCR modelling protocol (HGMP) and the fragment molecular orbital (FMO) quantum mechanics (QM) method to explore the interactions and selectivity of the human orexin-2 receptor (OX2R) and its recently discovered nonpeptidic agonists. HGMP generates a 3D model of GPCR structures and its complexes with small molecules by applying a set of computational methods. FMO allows ab initio approaches to be applied to systems that conventional QM methods would find challenging. The key advantage of FMO is that it can reveal information on the individual contribution and chemical nature of each residue and water molecule to the ligand binding that normally would be difficult to detect without QM. We illustrate how the combination of both techniques provides a practical and efficient approach that can be used to analyse the existing structure–function relationships (SAR) and to drive forward SBDD in a real-world example for which there is no crystal structure of the complex available.
G-protein-coupled receptor (GPCR)–ligand interactions are fundamental to almost all processes occurring in living organisms, and as such it is perhaps unsurprising that they are the targets of about 40% of all prescribed drugs [1–3]. What is surprising is that these drugs only target approximately 50 of the 800 known GPCRs . Thus there is huge potential in terms of the number of targets for new therapies to be designed against . Further progress of drug discovery for GPCRs is highly dependent on the understanding of structure–function relationships (SAR) and the interactions between the receptor and the small molecule (drug candidate) [4,6–8]. Currently there are no experimental or computational tools available that can provide an accurate list of interactions between the receptor and the ligand  and their complex chemical nature .
The efficiency and cost-effectiveness of any drug-discovery process is highly dependent on the availability of structural data regarding the target receptor and on the reliability of the data mining tools [6–8]. The recent breakthroughs in structural biology of GPCRs have resulted in the solving of over 100 structures of GPCR–ligand complexes representing over 30 unique GPCRs [7,11]. Nevertheless, the interpretation of interactions observed in the atomic-resolution structure is usually performed by ‘visual inspection’ or alternatively with a simple molecular mechanics (MM) model that cannot explain the full complexity of the molecular interactions . In many cases, the affinity and reasons why a particular ligand interacts the way it does with its receptor remain unclear.
Recently several notable reports have been published [10,12–14] that emphasized the pivotal role of a large number of ‘non-obvious’, hidden-from-eye interactions such as CH/π [15,16], halogen/π , cation/π  and non-classical H-bonds  in receptor–ligand binding that are not properly parameterized in currently available force fields (FF) . This problem can be handled with quantum mechanics (QM) methods, which have always been considered a reliable approach for the exploration of receptor–ligand interactions [20,21]. However, in spite of their many advantages, traditional QM methods are generally not feasible for large biological systems, due to their high computational cost .
The fragment molecular orbital (FMO) method [16,21,23] offers a considerable computational speed-up over traditional QM methods . One of the key features of the FMO approach is that it can provide a list of the interactions formed between the ligand and the receptor and a chemically intuitive breakdown of these interactions . Such information is essential for medicinal chemists to be able to rationally approach modification of lead compounds in order to increase favourable interactions. It works by dividing the system into smaller pieces called fragments (Figure 1). For example, in proteins, each residue can be represented as a fragment. Similarly, the ligand can be represented by single or multiple fragments as necessary. By performing QM calculations on fragments, one can achieve a high level of accuracy with very high efficiency.
Workflow for PIEDA calculations and details on each of PIE terms that are computed 
The pair interaction energy (PIE) between any two fragments calculated by FMO is a sum of four energy terms: electrostatics, exchange-repulsion, charge transfer and dispersion, provided by pair interaction energy decomposition analysis (PIEDA) –see Figure 1. The electrostatic and charge transfer terms are important in salt-bridges, hydrogen bonds and polar interactions, whereas the dispersion term can be thought of as hydrophobic in nature. The role of hydrophobic interactions is vital for biomolecular recognition but there is still no reliable predictive method for its quantification . The exchange-repulsion term describes the steric repulsion between electrons  that prevents atoms from collapsing into each other.
The key difference between FMO and MM methods derives from the fact that FMO takes into account polarization and charge transfer [16,26]. The description of electrostatics in most FF is based on static charges that neglect polarization and in polar systems such as proteins they are an approximation to the actual state. The van-der-Waals forces, despite being generally well parameterized on average, are not capable of detecting the directional nature of the dispersion terms involving halogens . Reported examples  comparing FMO and MM, have shown that the FMO method clearly outperformed FF-based scoring functions and demonstrated a high correlation with experimentally measured values of protein–ligand affinity [28,29]. In our recent report  we described how FMO has been applied for the analysis of 18 GPCR–ligand crystal structures representing different branches of the GPCR genome. This work revealed key and consensus interactions that are involved in receptor–ligand binding and were previously omitted from structure-based descriptions, including hydrophobic interactions, non-classical hydrogen bonds and the involvement of backbone atoms.
There is no need to compromise today in performing detailed analysis of protein–ligand structures using MM/FF whereas a similar analysis can be done with FMO that is reasonably quick. A typical FMO calculation on a ligand–receptor complex takes approximately 4 h on a 36 CPU cores to complete, significantly faster than weeks to a month (or more) for traditional QM approaches that have been used for estimating binding free energies.
Although X-ray crystallography is the preferred start point for structure-based drug-discovery (SBDD), in the context of a drug-discovery programme it is often the case that time prohibits its use for a series of complexes. The hierarchical GPCR modelling protocol (HGMP) has been developed to support GPCR SBDD programmes. HGMP has been successfully applied in GPCR drug discovery projects such as MCH-1R for obesity treatment , the orexin-1 and -2 receptors (OX1R and OX2R) for insomnia [31,32], the 5-HT2C for the treatment of metabolic disorders [33,34] and in other confidential drug discovery programmes. Additionally, the HGMP technology was used in the solving of the two H1R crystals structures  bound to the second and third generation antihistamines: cetirizine and fexofenadine.
In this work, we review how the integrated HGMP–FMO approach was applied to investigate the binding and selectivity of a recently reported nonpeptidic OX2R agonist (compound 26, see Figure 2) and its closely related analogues . We demonstrate that the FMO approach can be successfully applied to explore interactions in GPCR–ligand models rather than just crystallographic structures as previously reported [15,16,28,29,36,37].
First nonpeptidic OX2R compound 26 discovered by Nagahara et al. .
Hierarchical GPCR modelling protocol
The HGMP was developed by Evotec Ltd in conjunction with the University of Oxford to support SBDD programmes [30,31]. The HGMP generates a 3D model of GPCR structures and its complexes with small molecules by applying a set of computational methods. The models produced by HGMP are then used in drug discovery. The HGMP involves homology modelling followed by optimization protocols and flexible ensemble docking to predict binding poses and function of ligands bound to GPCRs. The HGMP includes a large set of integrated protocols like MD simulation , LowModeMD [30,38] and others to refine the GPCR models and exclusive scoring functions like the GPCR-likeness assessment score (GLAS) to evaluate model quality . The HGMP has been successfully applied in a large number GPCR drug discovery projects and also to support crystallography.
FMO is a general QM method that can be applied to any set of atoms, no matter if it is a soluble or membrane protein as described previously for exploration of the role of none-classical CH/π hydrogen bond in ligand recognition of β2-adrenergic GPCR receptor , for analysis of the GPCR–ligand crystal structures , for potency calculations on a novel Hsp90 fragment-linked inhibitor  and others. The FMO method can give a breakdown of the interaction energy between any two fragments and is consists of a sum of four pairwise interaction energy (PIE) terms: electrostatics, exchange-repulsion, charge transfer and dispersion (Figure 1). It is important to emphasize that the PIE is not a difference between energies of ‘free’ and ‘bound’ ligand but it rather represents the ‘strength’ of the interaction between the ligand and protein residues in the complex. Residues within a radius of ≤4.5 Å (1 Å=0.1 nm) around the ligand atoms often included in the FMO calculations. Based on the previous reports  we considered any interaction with an absolute PIE greater than or equal to 3.0 kcal/mol to be significant.
Medical significance of nonpeptidic OXR agonists
Class A GPCRs; OX1R and OX2R, are located mostly in the brain and associated with a range of different physiological functions, including the control of feeding, energy metabolism, modulation of neuro-endocrine function and the regulation of the sleep−wake cycle. Two non-selective neuropeptides orexin-A (OxA) and orexin-B (OxB) are natural agonists of OX1R and OX2R, which have dual activity, at both receptors. The dual activity aspect of these peptides has limited the usefulness of these natural agonists as probe compounds to dissect out the precise role of each receptor, in several conditions related to OXR activation. The peptides are also inefficient for in vivo studies due to lack of ability to penetrate the blood–brain barrier (BBB). Small-molecule agonists of OXR are important for both research and medicine as having the potential to address both these problems of selectivity and BBB penetration. A wealth of data so far suggests that OXR agonists could be used for the treatment of sleep disorders, narcolepsy, cataplexy, obesity, hypophagia, as well as attention deficit hyperactivity, depression and related bipolar disorders [35,40–43]. Furthermore, it was demonstrated that OX1R agonists might be promising candidates for colon cancer therapy . Activation of OX1R can drive apoptosis in human colon cancer cells and even reverse the development of established tumours.
However, in spite of their medical importance, the design of small-molecule agonists (rather than antagonists of peptide-activated GPCRs), is considered as one of the big challenges in drug discovery . This is because for agonists, there is the added requirement that it must not only bind the receptor but also activate it. Peptide-activated GPCRs like OX1R and OX2R, are considered especially challenging due to the large number of specific and non-specific interactions that are usually involved in peptide binding and activation.
FMO study of OX2R–agonist (compound 26) complex
Through an extensive synthesis and screening programme, Nagahara et al.  recently reported the discovery of the first selective nonpeptidic OX2R agonists culminating in compound 26 (Figure 2). This new chemical screening information along with the recently solved OX2R crystal structure  (PDB entry 4S0V) provides a new opportunity to develop drugs against this important target. As we do not yet have a crystal structure for the OX2R in complex with compound 26, the application of protocols such as the HGMP–FMO becomes the method of choice to advance the discovery of new ligands via the generation of plausible binding hypotheses that can be experimentally tested.
In previous site-directed mutagenesis (SDM) studies it was shown that alanine mutations of T1112.61, Q1343.32, D211ECL2, W214ECL2, Y2235.38, F2275.42, F3467.35 and H3507.39 caused a large (>50-fold) decrease in the potency of endogenous agonist without affecting the efficacy compared with WT. The mutations Y232A5.47 and Y317A6.48 resulted in a reduction of both EC50 (by 28.4- and 17.7-fold respectively) and Emax of 44.9% and 49.6%. These mutations caused a moderate decrease in potency of endogenous agonist (by 22.3-fold) without affecting its efficacy. These SDM data suggest that there is no clear correlation between the importance of residues for potency and for efficacy.
We recently proposed  two potential binding modes of compound 26 with OX2R produced by the HGMP: (1) ‘L’ shape docking pose (Figure 3) and (2) ‘U’ shape (Figure 4) as the antagonist Suvorexant adopts according to the recently solved crystal structure of the complex with OX2R (PDB entry 4S0V ).
‘L’-shaped docking pose of compound 26
In the first ‘L’ pose FMO detected 16 key interaction between compound 26 and OX2R residues: T1192.69, M1914.64, C210ECL2, D211ECL2, H2245.39, F2275.42, F2285.43, Y2325.47, Y3176.48, I3206.51, N3246.55, K3276.58, Arg3286.59, F3467.35, H3507.39 and with one water molecule HOH4025. Furthermore, these modelling observations are directly supported by the published SDM data [31,47]: the role of D211ECL2, F2275.42 F3467.35, H3507.39 and particularly of the residue Y232A5.47 in the potency and efficacy of OX2R endogenous agonist suggest their involvement in 26 agonist activity. However the relatively highly exchange repulsion term of D211ECL2 can be slightly artificial due to the fact that this residue is located on the loop which is more challenging to model. The exact role of F/Y5.47 (Y2325.47 in OX2R) in class A GPCRs activation is not clear but it is frequently engaged in interaction with agonists .
In the ‘U’ shape (Figure 4) pose, FMO detected 13 residues: T1112.61, P1313.29, Q1343.32, T1353.33, E212ECL2, H2245.39, F2275.42, Y3176.48, I320, V3256.56, R328, H3507.39, Y3547.43 and two water molecules: HOH4021 and HOH4025, that are involved in 26 binding. This pose is also supported by SDM data with residues T1112.61, Q1343.32, F2275.42 and H3507.39 and particular the interactions with Y3176.48. The interactions with the toggle switch residue Y3176.48 and with the aromatic cluster residue F2275.42, support the GPCR activation switch mechanism that allows compound 26 to have OX2R agonist function. This agonist-bound switch was proposed to be part of a larger ‘transmission switch’ that accounts for the relocation of conserved residues W6.48 (Y3176.48 in OX2R) and F6.44 towards P5.50 . A potential explanation for OX1R selectivity arises from potential interactions with the non-conserved residues T1112.61 (S1022.61 in OX1R) and T1353.33 (A1353.33 in OX1R).
‘U’-shaped docking pose of compound 26
The total PIE energy in docking pose ‘L’ is −110.17 kcal/mol and in ‘U’ is −143.18 kcal/mol. This suggests that in binding mode ‘U’, compound 26 forms a more stable complex with OX2R. In both docking poses the relative contributions of electrostatic and hydrophobic (dispersion) interactions are equal. However, we should add a note of caution here in that these are docked poses of agonists based on an antagonist-bound crystal structure template. There is a possibility that the agonist-bound form of the receptor is one that is quite different from the antagonist-bound form.
Another interesting aspect that we evaluated here is the correlation between experimentally measured EC50s and PIE as calculated by FMO for analogues of compound 26 as published by Nagahara et al. . It is known that ligand receptor affinity and efficacy might be driven by different energy terms including direct enthalpic contributions, entropy, solvation and the ‘strain energy’ of the ligand's bioactive conformation . Despite the fact that not all these factors are accounted for we observed significant correlation (r2=0.872, Figure 4D) between experimental values of EC50, measured for 16 analogues of compound 26 and the PIE calculated for their ‘U’ poses (no significant correlation was observed for the ‘L’ pose). These results are in agreement with our previously published report  where we demonstrated significant correlation between PIEs and experimentally measured affinities of OX2R Suvorexant-based antagonists. These observations increase our confidence that FMO can be used to provide additional insight into SBDD against GPCR targets even for modelled GPCR–ligand complexes.
FMO can be a highly useful tool for rational SBDD [16,39,49], as it provides an accurate and comprehensive list of strong, weak or repulsive interactions between the ligand and its surrounding residues. Such information is highly useful to guide modifications, substitution, scaffold hoping, linking or extension of chemical moieties to form stronger or new interactions with the protein or alternatively to remove repulsions. FMO can also be helpful in analysis of the ligand–water–protein network, to distinguish between energetically favourable and unfavourable water molecules and to design ligands that can interact or displace certain waters. As previously demonstrated , significant correlation between protein–ligand affinity and FMO energy terms  indicates that they can be efficiently used as descriptors in QSAR modelling to predict the binding affinities of new molecules. FMO has been successfully applied in the discovery of a novel Hsp90 inhibitors  and in many other our confidential drug discovery programmes. We learned that application of the FMO in hit-to-lead and lead optimization stages of drug discovery is a highly efficient for the design, evaluation and filtering of targets for synthesis that significantly decreases the effort, time and cost of chemical synthesis.
Here we present a new approach that opens an alternative avenue for the structural exploration of GPCRs and structure-based GPCR drug discovery. We have applied this approach to explore the binding and selectivity of an OX2R agonist. The outcome of this study is currently being applied in the generation of new OX2R agonists and for the discovery of first in class OX1R agonists. To our knowledge this is the first time that FMO calculations have been applied to docking poses of ligands in a GPCR. Applying the FMO analysis to this kind of problem results in two distinct benefits: (a) complex QM theories are condensed into four simple and intuitive quantities and (b) calculations become much faster than traditional QM approaches. This knowledge can be used to understand the chemical nature of existing receptor–ligand complexes, make suggestions for mutations and more importantly can suggest rational ligand optimization routes. The HGMP–FMO approach creates a cost-efficient new avenue for a SBDD against GPCR targets.
This work was supported by the Biotechnology and Biological Sciences Research Council [grant number BB/L026287/1 (to P.C.B. and A.H.); the Next Generation Super Computing Project, Nanoscience Program (MEXT, Japan) and Computational Materials Science Initiative (CMSI, Japan) (to D.G.F.); and the EPSRC and Evotec via the Systems Approaches to Biomedical Sciences Doctoral Training Centre [grant number EP/G037280/1 (to M.A.)].
fragment molecular orbital method
hierarchical GPCR modelling protocol
pair interaction energy
pair interaction energy decomposition analysis
structure-based drug discovery
GPCRs: Beyond Structure Towards Therapy: Held at Monash University, Prato Centre, Prato, Italy, 16–18 September 2015