Visualizing Electron Density

Problem

You want to visualize how a ligand fits into an electron density grid in a 2D molecule diagram (see Figure 1). The electron density grid with various contour levels and the 3D molecule is shown in Table 1.

../_images/1eta_elecdensity2img.png

Figure 1. Example of visualizing ligand electron density fit (1ETS)

Table 1. Example of ligand electron density fit at various contour levels (1ETS)
contour 1.0 contour 1.5 contour 2.0
../_images/1eta_screenshot_contour_1_0_small.png ../_images/1eta_screenshot_contour_1_5_small.png ../_images/1eta_screenshot_contour_2_0_small.png

Ingredients

Difficulty Level

../_images/chilly6.png ../_images/chilly6.png ../_images/chilly6.png

Solution

The SetElectronDensityContourOverlap function tags each atom with a value in the range [0.0, 1.0]. These values represent how much each atom exists in the contour. 0.0 - means no density in the given contour while 1.0 means that the atom is embedded in the grid at the given contour level.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
def SetElectronDensityContourOverlap(mol, skewgrid, contour, contourtag):

    center = OEFloatArray(3)
    extents = OEFloatArray(3)
    oechem.OEGetCenterAndExtents(mol, center, extents)

    subgrid = OEScalarGrid()

    # expand the grid a bit for proper overlaps

    extents[0] += 2.5
    extents[1] += 2.5
    extents[2] += 2.5
    OEMakeRegularSubGrid(subgrid, skewgrid, center, extents, 0.5,
                         skewgrid.GetReentrant() >= 7)

    # clip to the contour

    clipvalue = (lambda x: 0.0 if x < contour else 1.0)
    values = [clipvalue(x) for x in subgrid.GetValues()]

    subgrid.SetValues(values)

    overlap = OEOverlap()
    overlap.SetRefGrid(subgrid)
    result = OEOverlapResults()
    atomscores = OEFloatArray(mol.GetMaxAtomIdx())
    overlap.Overlap(mol, result, atomscores)

    # the sum of all points == 1 in carbon radius
    # should be theortically ~20.5

    for atom in mol.GetAtoms():
        overlap = atomscores[atom.GetIdx()] / 20.5
        if overlap > 0.0:
            atom.SetData(contourtag, overlap)

The SetElectronDensityContourOverlap function shows how to project the electron density fit with various contour levels into a 2D molecule diagram. First the image is divided into two image frames since both a molecule and a color gradient will be depicted. Then the SetElectronDensityContourOverlap function is called that calculates whether the atoms of the molecule are inside the electron density grid at various contour levels (see lines 11-13). The molecule is then prepared for depiction generating its 2D coordinates by calling the OEPrepareDepictionFrom3D function (lines 17-22). After constructing a color gradient that will assign colors to various contour levels, the function loops over the atoms and draws a circle around them if they are embedded into the electron density grid (> 0.2) with a color and radius that corresponds to the given contour level (lines 32-41). Finally, the molecule along with the color gradient is rendered to the image (lines 45-52). You can see the result in Figure 1.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
def DepictElectronDensityFit(image, mol, edgrid, opts):

    # generate image frames

    width, height = image.GetWidth(), image.GetHeight()
    mframe = OEImageFrame(image, width, height * 0.90, OE2DPoint(0.0, 0.0))
    lframe = OEImageFrame(image, width, height * 0.10, OE2DPoint(0.0, height * 0.90))

    # calculate fit to electron density at various contours

    contours = [1.0, 1.5, 2.0]
    for contour in contours:
        SetElectronDensityContourOverlap(mol, edgrid, contour, GetContourTag(contour))

    # prepare molecule for depiction

    width, height = mframe.GetWidth(), mframe.GetHeight()
    opts.SetDimensions(width, height, OEScale_AutoScale)

    OEPrepareDepictionFrom3D(mol)
    opts.SetScale(OEGetMoleculeScale(mol, opts) * 0.95)
    disp = OE2DMolDisplay(mol, opts)

    # create color gradient

    colorg = OELinearColorGradient()
    colorg.AddStop(OEColorStop(min(contours), OEColor(190, 190, 255)))  # light blue
    colorg.AddStop(OEColorStop(max(contours), OEColor(80, 80, 255)))  # medium blue

    # visualize electron density fit

    layer = disp.GetLayer(OELayerPosition_Below)
    for contour in contours:
        contourtag = GetContourTag(contour)
        radius = GetContourRadius(contour, contours, disp)
        color = colorg.GetColorAt(contour)
        pen = OEPen(color, color, OEFill_On, 1.0)
        for atom in mol.GetAtoms():
            if atom.HasData(contourtag) and atom.GetData(contourtag) > 0.2:
                adisp = disp.GetAtomDisplay(atom)
                layer.DrawCircle(adisp.GetCoords(), radius, pen)

    # render molecule

    OERenderMolecule(mframe, disp)

    # draw color gradient

    copts = OEColorGradientDisplayOptions()
    copts.SetColorStopPrecision(1)
    copts.AddMarkedValues(contours)
    OEDrawColorGradient(lframe, colorg, copts)

Download code

elecdensity2img.py and supporting data 1eta.mtz and 1eta_ligand.sdf

Usage:

prompt > python3 elecdensity2img.py -grid 1eta.mtz -ligand 1eta_ligand.sdf -out 1eta.png

Discussion

Visualizing electron density helps to evaluate the quality of protein-ligand structures. The Iridium database divides the protein-ligand structures into three categories:

  • Iridium-NT (not trustworthy)
  • Iridium-MT (moderately trustworthy)
  • Iridium-HT (highly trustworthy)

Not surprisingly the 1ETS complex (Figure 1 is considered “not trustworthy” by the Iridium database. The Table 2 shows examples from the MT and HT categories of the Iridium database

Table 2. Examples from the Iridium database
1COY - Iridium MT 1D3H - Iridium HT
../_images/1coy_elecdensity2img.png ../_images/1d3h_elecdensity2img.png

The visualization also helps to compare the “deposited” starting models of the Iridium database with the models that are refined by AFITT. See example in Table 3.

Table 3. Examples from the Iridium database
deposited starting model model after AFITT refinement
../_images/1cx2_start_elecdensity2img.png ../_images/1cx2_refined_elecdensity2img.png

See Also in OEChem TK Manual

API

Theory

See Also in Shape TK Manual

API

See Also in OEDepict TK Manual

Theory

API

See Also in OEGrapheme TK Manual

API