Journal Articles

Ligand Interaction Diagrams

Dr. Alex M. Clark
October 2006
Chemical Computing Group Inc.
 

Introduction

A 2D molecular display method has been added to MOE 2006 (Molecular Operating Environment), which is designed for viewing the active sites of protein:ligand complexes. The ligand is arranged and rendered using an enhanced version of the previously introduced 2D depiction layout algorithm, and the protein residues are arranged around it in order to indicate spatial proximity and hydrogen bonds. Additional interaction styles are used to render coordinated metals, bound solvent molecules, substitution proximity contour and solvent-accessible surface area. While optional parameters can be specified as input, generation of the diagram is entirely automated.

BIOTIN
3D Active Site 2D Diagram
3D Active Site2D Diagram

Methodology

The generation of Ligand Interaction diagrams in MOE differs from preceding efforts such as LIGPLOT [Wallace 1995] in that the diagram is produced with a high regard for aesthetics, and as a primary step the ligand itself is subjected to a layout procedure which is derived from the standalone 2D depiction algorithm, under additional constraints. The interacting residues are drawn as single units for purposes of clarity, and are positioned in order to try to impart as much indication of the relative 3D distances as possible, but with a heavily weighted preference for achieving a clear layout at the expense of strictly reproducing exact distances. When the layout is complete, the large amount of information content found in an active site is clearly and concisely annotated with visual cues.

Residue Interactions

Residues in the active site are identified in two steps: (1) receptor residues, ions and solvent molecules which have a strong interaction with the ligand, such as a favorable hydrogen-bonding interaction, and (2) receptor residues and ions which are close to the ligand, but whose interactions with the ligand are weak or diffuse, such as collective hydrophobic or electrostatic interactions.

Hydrogen bonds and chelated metal interactions are enumerated and scored by a pairwise comparison of heavy atoms, which includes parameters such as:

  • Atom types (element, hybridization, bonding environment)
  • Distance
  • In-plane and out-of-plane angles of substituents

For certain atom pair combinations (e.g. secondary amine nitrogen and carbonyl oxygen), a hydrogen bond is considered possible and a scoring function is defined, having been trained on large quantities of protein data using contact statistics [Labute 2001]. The score is expressed as a percentage probability, where 10% would be considered a significant likelihood of being a strong hydrogen bond, whereas 1% would be considered relatively weak. The default minimum cutoff is set at 10%, but may be modified by the user.

Solvent molecules such as water are often computational or crystallographic artifacts, and are excluded from the diagram unless they are anchored by at least two qualifying hydrogen bonds, one to the ligand and one to the receptor.

Chelated metal ions are identified by a similar technique, with a different set of parameters.

Receptor residues and ions which remain may be included in the diagram as second-tier entities, if they are sufficiently close to the ligand. The method used to include residues is to define a maximum distance between heavy atoms of the ligand and receptor as 4.5 Å. Any residue with a heavy atom this close qualifies. The distance is then extended to 4.6 Å, within which range a residue must have two atoms, and so on, out to 10 atoms at 5.4 Å. This simplistic criterion is sufficient for selecting residues for inclusion in the diagram, prior to a more quantitative definition of diffuse interactions, which will be established subsequently.

Ligand Layout

We have previously described a method for reliably producing aesthetically high quality 2D layouts for small molecules [Clark 2006], which has been available since MOE 2005. In MOE 2006, there is an additional capability for providing a bias function to influence the desirability of each proposed 2D layout. The Ligand Interactions algorithm exploits this feature in order to encourage the selection of interatomic distances such that the 2D layout is similar to the 3D conformation. The magnitude of the scoring function is scaled such that it is sufficiently large to override the scoring method which the depiction layout uses internally for nonbonded repulsion, but not large enough to force the depiction layout algorithm to enter disaster recovery mode.

score = score + Σtorsion w |θ3D - θ2D| + Σi,j w |dist3Dij - dist2Dij|

Therefore, the resulting layout is bound by the constraints of an aesthetic diagram layout, but as a secondary constraint, does not maximize dispersion but instead favors emulating the input conformation, as is shown in the following example:

3D conformation of ligand

Normal depiction

Flattening depiction

Residue Placement

The flattened coordinates of the ligand are used as the starting point for the Ligand Interactions diagram. The receptor residues, solvent molecules and ions are placed around the 2D ligand, beginning with strongly-bound residues.

Residues are positioned by using an initial placement guess, followed by the continuous optimization of a function based on residue positions. The residue layout proceeds in two steps (strongly bound residues, followed by weak). Both groups use the same initial placement method, which consists of creating a grid around the ligand. Upon atoms and ring centers are added a Gaussian distribution of repulsion intensity, as shown in the following illustration:

1. Base avoidance

To find the most sensible starting position for an individual residue, the minimum distance from each point to the ligand atoms with which the residue is in contact is determined, and this value is scaled and overlaid onto the grid. As shown in the following example, the residue with an interaction to the aromatic nitrogen atom (circled) has a most-favorable starting position at the position marked by the X:

2. Specific attraction

Once plausible starting positions have been obtained for the strongly-bound residues, a continuous energy function is composed and optimized, in order to obtain a more finely-tuned placement.

energy = Σij wij (dist2ij - D2)2       (variable-variable desired distances)
+ Σij wij exp [-K dist2ij]       (variable-variable short range repulsion)
+ Σik wik (dist2ik - D2)2       (variable-fixed point desired distances)
+ Σik wik exp [-K dist2ik]       (variable-fixed point short range repulsion)
+ Σθ wθ exp [-½K (1 - cos θ)2]       (variable-fixed point angular repulsion)

The variables are the 2D positions of the strongly-bound receptor residues, and the fixed points are positions of the ligand atoms. The ideal solution to the function is one in which the residue:residue and residue:ligand distances are the same in the 2D layout as they are in the 3D active site, but due to the reduction in dimension and repulsion terms which express the aesthetic preference for a disperse layout, these effects are counter-balanced.

The receptor residues which were not classified by strong specific interactions are placed subsequently, at which time the previously placed residues are considered to be fixed points. The same initial placement method is the same, as is the subsequent refinement, except that the residues are weakly tethered to their initial placement positions, rather than to nearby parts of the ligand.

3D Active Site 2D Placement

Proximity Contour

In order to provide an indication of the amount of space that is available at the possible substitution points on the ligand, a proximity contour around the perimeter is calculated. The contour is calculated by first assigning a magnitude to each atom, which is subsequently used to produce a 2D surface.

The space about each atom is calculated by first obtaining its "substitution vectors", which are geometrical locations that would be appropriate for hydrogen atoms, depending on the atomic hybridization and number of heavy substituents. For terminal non-linear atoms, an umbrella-like sampling is used:

sp3/3 sp3/2 sp3/1
sp2/2 sp2/1 sp/1

Along each of the substitution vectors, a cannonball is fired, considered to have a radius equal to the van der Waals radius of the parent atom. The distance travelled before connecting with the van der Waals surface of any other heavy atom in the system (ligand or receptor) is the extent, which is bounded between 0 and 2 Å. The magnitude of the atom is set to be the average. If more than half of the vectors extend out into space indefinitely, the atom is classified as being entirely exposed.


Proximity contour: 1RO6 Hydrolase

The contour drawn around the ligand is produced by plotting Gaussian spheres based on the magnitudes of each of the constituent atoms, then reducing to a perimeter with a cutoff value that scales the observed distance in a manner that has qualitative significance, e.g. if the distance between a terminal atom and the contour line is approximately that of a single bond length, it is indicative of there being enough room for substitution of an additional heavy atom at that point. The contour line is broken if it is closest to an atom which is fully exposed.

Surface Contacts

Surface contacts are defined by a notion of solvent accessible surface area, and displacement thereof. Each heavy atom is considered to be enclosed by a sphere of van der Waals radius, plus the effective radius of a water molecule (taken to be 1.4 Å):

Whenever the outer radii of two atoms overlap, there is a portion of the surfaces of both atoms where there is no longer sufficient room for a water molecule to roam freely about the surfaces of either:

Therefore the total solvent accessible surface area is correspondingly reduced for both atoms. In practice, because the spheres placed around the atoms are so large, being in proximity to several covalently bound neighbors and a handful of moderately proximate non-bonded atoms quickly reduces the solvent accessible surface area to zero. A residual exposed surface area on a ligand atom is therefore indicative of an empty pocket inside the active site, or that the atom in question projects out of the active site.

In this context, solvent accessible surface area is interesting in two ways: ligand atoms with remaining surface area are either outside of the binding pocket, or they are in an interior region that would ideally be hydrophilic, and in practice filled with water. Hydrophobic atoms with high solvent accessible surface area are likely to be energetically disfavored. From the receptor point of view, the interesting property is the difference in solvent accessible surface area which is induced by the ligand. For hydrobobic residues in particular, the replacement of exposure to water by proximity to non-polar regions of the ligand are typically indicative of hydrophobic interactions which may have a significant role in stabilizing the complex.

Solvent accessible surface area of the ligand is plotted directly onto the atoms in the form of a blue smudge, as is shown in the following diagram:

Difference in solvent accessible surface area for the receptor residues is accumulated over the constituent atoms, and plotted as a turquoise halo, as is shown below:

Rendering

Various options are available with regard to the rendering modes, but the default settings convey specific meaning, as is demonstrated in the following example:

Residues are annotated with their 3-letter amino acid code, and their position classification. If there are multiple chains in the system, the positions are prefixed with letters of the alphabet.

There are 4 default color schemes for residues. Hydrophobic residues are all colored with a green interior, whereas polar residues are colored in light purple. Basic residues are further annotated by a blue interior ring, and acidic residues with a red ring. Water molecules are drawn without interior color, and metals are colored grey.

Hydrogen bonding interactions between the receptor and the ligand are drawn with an arrowhead to denote the direction of the hydrogen bond (i.e. the donor is at the base of the arrow, and the acceptor is at the head). In cases where there could be a hydrogen bond in either direction, the arrow is double headed.

When the hydrogen bond is formed with the residue sidechain, the arrow is drawn in green. Hydrogen bonds to the residue backbone are drawn in blue, with an additional dot drawn at the residue attachment point.

The ligand is always drawn in a conventional 2D schematic form, with most of the hydrogen atoms either implied or subsumed into the atom label. Hydrogen atoms explicitly involved in a hydrogen bond are drawn individually, and oriented toward the residues with which they interact. The annotated proximity contour and solvent surface area exposure properties are drawn in the background of the diagram.

Gallery

1RO6 1ETS
1RO61ETS
1FVT 1G5S
1FVT1G5S
1ADD 1ATP
1ADD1ATP

Future Work

Since the release of MOE 2006, a number of subsequent improvements are in development, some of which are:

  • Residue-residue interactions (such as hydrogen bonds or disulfide bridges), which influence the order and priority of weakly bound residue placement, in order to produce an arrangement which is more readily comparable to the original 3D active site, and may optionally be drawn in the diagram.

  • Arene π-π and π-cation interactions, which are detected and drawn in the final diagram.

  • Ligand series alignment, when multiple similar active sites are available in the input. The orientation of the ligands relative to one another can be reproduced in the series of diagrams.

  • Concurrent residue placement, for a series of similar (but not necessarily identical) active sites, in which the union of residues is arranged for all active sites at the same time, such that the resulting diagrams show an identical binding pocket.

  • Conserved interactions and residues for a series, showing which features are consistent throughout.

References

[Wallace 1995] Wallace, A.C., Laskowski, R.A., Thornton, J.M.; "LIGPLOT: a program to generate schematic diagrams of protein-ligand interactions". Protein Engineering 8 (1995) 127-134.
[Clark 2006] (a) Clark, A.M. ; Labute, P.; Santavy, M.; "2D Structure Depiction", J. Chem. Inf. Model., 46, 1107-1123 (2006); (b) Chemical Computing Group Journal article (2005).
[Labute 2001] Labute, P.; "Probabilistic Receptor Potentials". Chemical Computing Group Journal (2001), http://www.chemcomp.com/journal/cstat.htm.
[FuncRef] SVL function reference, available to MOE users.