Computational Chemistry

This page describes the current research in our lab to better understand the molecular mechanisms involved in receptor-ligand interactions.

For this purpose we use

• classical medicinal chemistry,
• computational medicinal chemistry, and
• molecular biology tools.

This particular page is an excerpt from a publication:

• van Rhee, A.M.; Fischer, B.; Jacobson, K.A.
  Modelling the P_2Y purinoceptor using rhodopsin as template.
  Drug Design Discovery, 1995, 13: 133-154.
  

which describes how we approached the problem of visualizing the molecular architecture of the G protein-coupled P2Y₁ purinoceptor.

It is the intention to add more features (both classical medicinal chemistry and computational chemistry) as we go along, but time is the limiting factor. You could make a bookmark now, and check in every now and then.

A very characteristic feature of G protein-coupled receptors (GPCRs) is their overall topology:

They all have similar hydrophobicity profiles with seven stretches of increased hydrophobicity, most probably corresponding to alpha-helical regions spanning the cell membrane.

Unfortunately, isolation, purification and crystallization of P₂ purinoceptors, as for other GPCRs, remains to be demonstrated.

Not only is a crystallographic structure as yet unobtainable, but the very limited quantities of receptor protein present in the membrane, and isolation problems, prohibit the use of modern techniques, such as nuclear magnetic resonance spectroscopy, towards structure elucidation.

Chemical modification ²⁶ of the ligand binding site and site-directed mutagenesis have been shown to be useful techniques in understanding receptor-ligand interactions, but is particularly in the case of P₂ purinoceptors restricted by the lack of suitable radioligands.

A possible further approach is to construct models of the receptors and their ligand complexes using molecular modelling, as was first proposed by Hibert et al.

On these pages we present models of the chick P_2Y purinoceptor, and the human A_2A receptor based on the electron density map of rhodopsin ²⁸ , rather than on the bacteriorhodopsin template ²⁹ used by Hibert et al.

The models are verified by incorporation of site-directed mutagenesis ^{41, 50, 51} and ligand affinity data. ^{14, 15}

• Quanta,
  
• Insight II/Discover, and
• MOPAC molecular modelling packages.
  

Model building.
Kyte-Doolittle hydrophobicity and Emini surface probability parameters were calculated using a 7 amino acid window. With these parameters, the putative transmembrane domains (TMs) were identified.

Sequences of the putative transmembrane domains were aligned manually using the alignment procedure described on our "GPCR page".

Helices were built from the full length of the TM, using the Biopolymer module in the Insight II program. The N terminus was capped with an acetamido group, and the C terminus was capped with a carboxamido group. The secondary structure was assumed to be a right-handed alpha helix (default phi and psi angles from the Biopolymer module). After generation of the Amber atom type parameters and template charges, the helices were energy minimized as described in the Energy Minimization section.

The energy minimized helices were converted to Protein Data Bank (PDB) format files and imported in the Quanta program. Coils (an alpha carbon backbone representation of the helix) of the minimized helices were drawn, that resemble the helical wheel representation to aid in orienting the individual helices towards each other.

The following constraints were used constructing the helical bundle:

• 1. The helical axes of the first, third, fifth and seventh helices were quasi-antiparallel to those of the second, fourth and sixth helices.
  
• 2. The hydrophobic side of each helix was facing the lipid phase and the hydrophilic side of each helix was facing either another helix or the pore formed by the helical bundle.
  
• 3. Conserved residues, either identical throughout a subset or highly homologous, determined the orientation of each helix relative to the other helices.
  
• 4. Distances in the electron density map of rhodopsin ²⁸ correlated to atomic coordinates in the model. This was obtained through triangulation of the electron density map and conversion of the coordinates to a grid on the screen.
  
• 5. The assembly of helices maintained a clockwise order, when seen from the intracellular side, as argued by Baldwin. ³¹
  
• 6. None of the helices were intersecting.

Energy Minimization.
The individual helices generated with Insight II's Biopolymer module were energy minimized in Discover using a stepwise process. Initially, 200 steps of Steepest Descent were performed, followed by minimization using the Conjugate Gradient (CG) method until the gradient reached a value below 0.1 kcal/mol/Å. Use of the Amber forcefield in Discover required that in all calculations the 1-4 nonbonded interactions were scaled by a factor 0.5, and the dielectric constant was assumed to be distance independent with a magnitude of 4. The bond, theta, phi, out-of-plane, nonbonded and Coulombic interactions were all used to obtain the final energy, but were not expressly scaled.

The helical bundles generated with Quanta were energy minimized in Discover in a stepwise process. Initially, 500 steps CG were performed with the backbone of the helices tethered with a force constant of 100 kcal/Å. In consecutive runs (500 steps CG each), the force constant was reduced to 50 kcal/Å, then 25 kcal/Å and, finally, to 0 (no tethering).

Ligand Docking.
Selective and/or specific ligands for the receptors were fully minimized using the MNDO hamiltonian of MOPAC. All ligands were rigidly docked into the helical bundle using graphical manipulation coupled to continuous energy monitoring, i.e. the ligand was manually docked into the binding site without relaxing the atomic coordinates of either ligand or protein while continuously calculating the energy of the whole. When a final position was reached, consistent with a low local energy and known pharmacological data, the complex of receptor and ligand was subjected to a minimization run of 4000 steps CG (or until the gradient was < 0.1 kcal/mol/Å). Charges for the ligands were imported from the MOPAC output files.

We will now briefly discuss the results obtained through this approach, but the information specific to either the P_2Y purinoceptor or the adenosine A_2A receptor is kept on separate pages.

Alignment of G protein-coupled receptors and dendrogram.
In TM1, the motif GXXXN (P_2Y) or GXXGN (P_2U) occurs in P₂ purinoceptors, rather than the GN motif present in biogenic amine or adenosine receptors. The alternate motifs, however, are not exclusively used by P₂ receptors, but are shared with, e.g., the PGE₃-II (GXXGN) and the PAF (GXXXN) receptors. In all sequences the last 5 C-terminal residues of this helix are frequently occupied by lysine or arginine residues, indicating the end of the transmembrane domain. Such basic residues may serve as 'membrane anchors', ^{34, 42} and are useful in determining the position of the helix in the lipid bilayer.

In TM2, the residue preceding the conserved leucine in the LXXXD motif is a conserved serine for the biogenic amine receptors, but is either an asparagine or a histidine in the P_2Y and P_2U receptors, respectively. The asparagine is conserved among the P_2Y, AT_1A, CKR1, IL-8A, PAF, NK2 receptors, the opsins and rhodopsin. The histidine is shared between P_2U and thrombin receptors.

Another difference between biogenic amine receptors and P₂ receptors is the position of the conserved proline relative to the conserved aspartate.

Although absent in the muscarinic receptor, the proline is consistently spaced by 8 residues from the aspartate in the other biogenic amine receptors, whereas in many other receptors, the P₂ receptors included, only 7 residues separate these two residues.

This last difference is possibly significant since the proline is located near the luminal side of the receptor protein, which is thought to be important for ligand binding.

The only exceptions, observed so far, to the conservation of the aspartate occur in the substance P receptor where it is substituted with a glutamate, and in the gonadotropin-releasing hormone receptor where an asparagine replaces this residue.

The conserved aspartate in TM2 was shown by Horstman et al to be a sodium-dependent allosteric regulatory site in the alpha₂-adrenergic receptor.

The DRY motif, characteristic of the third transmembrane domain in biogenic amine receptors, is replaced by a HRY motif in the P2Y₁ receptor (HRC in the P_2U receptor). This motif has not been found in other receptors. The only other substitution for the conserved aspartate, seems to be the D to E in, e.g., the TXA₂, PGE₃-II, opsin and rhodopsin sequences. The significance of this deviation is not clear, but is supposedly important for coupling of the P₂ receptors to G proteins and not for ligand binding, the subject of this study.

One or more of the first 6 positions in the N-terminal sequence of TM4 are generally (37 out of 40 aligned sequences of family A GPCRs ) occupied by lysine or arginine residues. Again these residues could well serve as 'membrane anchors'. They seem to be occurring more frequently towards the cytosolic side of helices than on the luminal side, thus reflecting the polarity of the membrane. These 'membrane anchors' also occur much less frequently in the C-terminal sequences of helices, than in the N-terminal sequences, as is demonstrated by the alignment for TM5. The reason for this disparity is believed to be the direction in which the side-chains are pointing, i.e. opposite to the propagation direction of the helix, allowing more rotational freedom for lysine and arginine residues near the cytosolic N-terminus than for the cytosolic C-terminus.
There are no other marked differences between the various receptor subfamilies for either TM4 or TM5.

The CXXP motif used for the alignment of TM6 could be substituted by WXP for most receptors, but the P₂ receptors deviate at this position. The P_2Y receptor has a tyrosine at the position of the conserved tryptophan, and this seems to be rather unique. The phenylalanine in the P_2U receptor at the same position occurs more frequently, such as in the thrombin, and various orphan receptors. Characteristic of the P₂ receptors is the presence of a lysine (P_2Y: K269) or an arginine (P_2U: R265) at an otherwise non-conserved position. It shares this feature only with the endothelin receptors (both ET-A and ET-B subtypes) and the orphan receptor RSC338. The proposed 'membrane anchors' occur quite frequently at the first position in the alignment of this transmembrane domain.

TM7 is best aligned by means of the NPXXY motif. However, the P₂ receptors, the gonadotropin-releasing hormone receptor, the TXA₂, PGE₃-II, and several orphan receptors constitute an exception to this rule. The conserved asparagine is replaced by an aspartate residue in the latter cases. A non-conserved aspartate (D352) in the 5HT_1B receptor aligns perfectly with aspartates in the IL-8A receptor (D288), and in the D_1B receptor (D342), a glutamine in the P_2Y receptor (Q296) and a lysine in the P_2U receptor (K289). Arginines 299 (P_2Y) and 292 (P_2U) align with Y530 in the m3, E291 in the IL-8A and E287 in the CKR1 receptor, respectively. Both positions are near the luminal side of the receptor and are likely involved in ligand binding. The K and R residues occurring near the C-terminus are probably better characterized as 'membrane anchors'.

If you are interested in the alignment of over 150 GPCRs, or in mutagenesis studies of GPCRs, just follow these links. You might want to read our introductory page first, though.

Bacteriorhodopsin, included in this study to facilitate comparison with other modelling studies, ^{[e.g., 27, 35]} is clearly not related to any of the GPCR subfamilies. The degree of relatedness between bacteriorhodopsin and GPCRs shown in the dendrogram is probably an overestimate of the actual distance, caused by the residue alignment procedure of the GCG program.

Kyte-Doolittle hydrophobicity and Emini surface probability analysis.
At the portion of the sequence where proteins are supposed to cross the lipid bilayer membrane, the hydrophobicity of that segment is usually increased relative to the cytosolic and luminal portions. The procedure developed by Kyte and Doolittle uses this phenomenon to identify transmembrane domains in protein sequences with unknown structure. In the case of GPCRs, these transmembrane domains are thought to be alpha helical, and more importantly amphiphilic. The amphiphilicity of the TMs supposedly reflects the way a GPCR is built out of seven of these helical transmembrane domains. Consequentially, the hydrophobicity profile of these GPCRs is not as clear as one would wish.

An example of this effect is the poor separation of the sixth (VI) and seventh (VII) TM in both the P2Y₁ receptor profile (below) and the bacteriorhodopsin profile (not shown).

To facilitate identification of TMs, the Kyte-Doolittle hydrophobicity method is often supplemented with methods describing other sequence dependent vectorial parameters, such as the 'conservation moment' and 'hydrophilic moment', ³⁶ for each helix or other computational methods ⁴³ . Use of the 'conservation moment' method requires a particularly well defined alignment of highly homologous sequences, and can not be applied to the P₂ GPCRs because of the low sequence similarity of these receptors with any other GPCR subfamily. The 'hydrophilic moment' method relies on a database of partial hydrophilic factors for any given amino acid and is dependent on the method with which these factors were measured. Furthermore, the high incidence of basic residues in the P_2Y receptor sequence greatly influences the results obtained with this method.

In contrast, the Emini surface probability can be calculated from a given sequence, indicating the propensity of a stretch of amino acids (the window) to be at the surface of a protein. The combination of an increase in the Kyte-Doolittle hydrophobicity index and a decrease in the Emini surface probability index was successfully applied to predict the TMs in the P2Y₁ receptor (see figure above) and in the reference protein bacteriorhodopsin. The start and end residue of each helix was usually predicted within 3 residues of the TMs used by Henderson et al. The only major deviation in the P2Y₁ receptor prediction occurred at the N-terminus of TM7 (VII). This particular sequence contains 3 hydrophilic residues (K, Q and R) that are possibly involved in ligand binding, thus delaying the start of the predicted TM.

Building the receptor model.
Modelling of G protein-coupled receptors has become an important tool in understanding drug-receptor interactions and in the development of new ligands for these receptors.

The first widely accepted method was the homology modeling method by Hibert et al. This method involved the alignment of the receptor sequence with the sequence of bacteriorhodopsin, and the subsequent mapping of the sequence onto the structure of bacteriorhodopsin that was determined by Henderson et al.

Bacteriorhodopsin is a proton pump in the outer membrane of Halobacterium halobium, and lacks any functional or sequence homology with GPCRs. Nevertheless, the procedure was based on the assumption that there would be considerable structural homology. This structural homology was inferred by the extraordinary similarity in the hydrophobicity plots, or Kyte-Doolittle plots, of the biogenic amine subfamily of GPCRs and bacteriorhodopsin.

Recently, a low resolution electron density map of rhodopsin, a true member of the GPCR superfamily, was published. ²⁸ The low sequence homology with bacteriorhodopsin, the structural differences that must arise from the different placement of proline residues (causing bends in helices) in bacteriorhodopsin- and GPCR-sequences, and the availability of an electron density map of a true member of the GPCR superfamily prompted us to adapt a new method to build models of GPCRs. ³⁴ This novel method is based on a computational approach rather than strict compliance with the atomic coordinates of a distantly related protein, albeit with higher resolution.

To ascertain the viability of modelling transmembrane proteins (based on the structure of rhodopsin ³⁴ ) by the methods described in the paper, we built models of bacteriorhodopsin (data not shown) based on the electron density map of bacteriorhodopsin, ²⁹ and rhodopsin (data not shown) based on the electron density map of bovine rhodopsin as published by Schertler et al.

Root mean square (r.m.s.) distance calculations on superimposed structures were performed to establish how well the models fitted the experimental data. The bacteriorhodopsin model that was constructed, compared with the one deposited in the PDB, ²⁹ had an r.m.s. value of 2.15 Å when measured on all backbone atoms and 2.06 Å on all C^alpha atoms. Both values are lower than the resolution of the model, i.e. 3.5 Å. When the rhodopsin model was compared with the bacteriorhodopsin model by Henderson et al, however, the r.m.s. value increased to 16.96 Å when measured on all backbone atoms and 16.86 Å on all C^alpha atoms. These rather high values indicate that there are more structural differences between bacteriorhodopsin and rhodopsin than have been assumed thus far. ⁴⁴

In our opinion, the position of the proline residues in the helices is a major determinant in this matter. As can be seen from the alignment for the helices, proline residues occur less frequently and at different positions in the sequence of bacteriorhodopsin than in the sequence of rhodopsin. The influence of proline residues in alpha-helices was extensively studied by Sankararamakrishnan and Vishveshwara and von Heijne.

We have found that applying the Amber forcefield to our calculations yields results in agreement with Sankararamakrishnan and Vishveshwara and von Heijne, and that those results are consistent with data obtained from a crystallographic study of the transmembrane protein photosynthetic reaction center (PDB:1prc), and globular proteins such as phosphoglycerate kinase (PDB:3pgk), lysozyme (PDB:1l27) and alcohol dehydrogenase (PDB:5adh). The differences between the structures derived from the electron density maps of bacteriorhodopsin and rhodopsin also illustrate why helical wheel models, widely used by molecular biologists, are highly imprecise when applied to GPCRs. ³¹

The general picture looks like this:

For a detailed view of either our

• A_2A adenosine receptor model, or our
  
• P2Y₁ purinoceptor model,
  just follow the appropriate links.