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.

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

 contour 1.0 contour 1.5 contour 2.0

Ingredients¶

 OEChem TK - cheminformatics toolkit Shape TK - molecule shape overlap toolkit OEDepict TK - molecule depiction toolkit Grapheme TK - molecule and property visualization toolkit

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) ```

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

 1COY - Iridium MT 1D3H - Iridium HT

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.

 deposited starting model model after AFITT refinement

API

Theory

API