Molecular Heatmaps via weighted Gaussian & rational-blob kernel summation over per-atom descriptors

[pleasureoffiguring (J. M.)]

Open in Wolfram Cloud ↗

DOWNLOAD-DESKTOP-NOTEBOOK.nb (5.79 MB)

[wolframcommunity (EDITORIAL BOARD)]

– you have earned Featured Contributor Badge
Your exceptional post has been selected for our editorial column Staff Picks http://wolfr.am/StaffPicks and Your Profile is now distinguished by a Featured Contributor Badge and is displayed on the Featured Contributor Board. Thank you!

[Vitaliyk (Vitaliy Kaurov)]

I added a small animation to show how the heatmap is built rather than just showing the final picture. The atoms are added in order of decreasing Topological Steric Effect Index, so the most sterically crowded parts of idarubicin contribute first. This makes the plot readable in a different way: one can see which atomic neighborhoods dominate.

ClearAll["Global`*"];

mol2 = Molecule[Entity["Chemical", "Idarubicin"]];

coords2da =
  MoleculeValue[
   mol2,
   "AtomDiagramCoordinates",
   IncludeHydrogens -> None
   ];

coords2daAll =
  MoleculeValue[
   mol2,
   "AtomDiagramCoordinates",
   IncludeHydrogens -> All
   ];

rot = RotationTransform[-Pi/2, Mean[coords2da]];

coords2daV = rot /@ coords2da;
coords2daAllV = rot /@ coords2daAll;

mol2V = Molecule[
   mol2,
   AtomDiagramCoordinates -> coords2daAllV
   ];

topid =
  MoleculeValue[
   mol2,
   "TopologicalStericEffectIndex",
   IncludeHydrogens -> None
   ];

order = Ordering[topid, All, Greater];

{{xm, xp}, {ym, yp}} =
  MinMax[#, Scaled[0.1]] & /@ Transpose[coords2daV];

gauss[u_?NumericQ] := If[u > 40, 0., Exp[-u]];

hmfun2[n_Integer][x_?NumericQ, y_?NumericQ] :=
 With[
  {sc = 0.3, ii = order[[;; n]]},
  Total[
   MapThread[
    #2 gauss[
       SquaredEuclideanDistance[{x, y}, #1]/(2 sc^2)
       ] &,
    {coords2daV[[ii]], topid[[ii]]}
    ],
   Method -> "CompensatedSummation"
   ]
  ];

Manipulate[
 Show[
  DensityPlot[
   hmfun2[n][x, y],
   {x, xm, xp},
   {y, ym, yp},
   PerformanceGoal->"Quality",
   AspectRatio -> Automatic,
   ColorFunction->"ValentineTones",
   Frame -> None,
   PlotPoints -> 100,
   PlotRange -> All,
   PlotRangePadding->0,
   ImageSize->550
   ],
  MoleculePlot[mol2V,ColorRules -> {_ -> Directive[White,Thickness[.005]]}]
  ],
 {{n, Length[topid], "atoms included"}, 0, Length[topid], 1,
  Appearance -> "Labeled"}
 ]

frnsTM=ParallelTable[
 Show[
  DensityPlot[
   hmfun2[n][x, y],
   {x, xm, xp},
   {y, ym, yp},
   PerformanceGoal->"Quality",
   AspectRatio -> Automatic,
   ColorFunction->"LightTemperatureMap",
   Frame -> None,
   PlotPoints -> 100,
   PlotRange -> All,
   PlotRangePadding->0,
   ImageSize->550
   ],
  MoleculePlot[mol2V,ColorRules -> {_ -> Directive[Black,Thickness[.005]]}]
  ],
 {n, 0, Length[topid], 1}
 ];//AbsoluteTiming

frnsSC=ParallelTable[
 Show[
  DensityPlot[
   hmfun2[n][x, y],
   {x, xm, xp},
   {y, ym, yp},
   PerformanceGoal->"Quality",
   AspectRatio -> Automatic,
   ColorFunction->"SolarColors",
   Frame -> None,
   PlotPoints -> 100,
   PlotRange -> All,
   PlotRangePadding->0,
   ImageSize->550
   ],
  MoleculePlot[mol2V,ColorRules -> {_ -> Directive[White,Thickness[.005]]}]
  ],
 {n, 0, Length[topid], 1}
 ];//AbsoluteTiming

Export["mollights.gif", ImageAssemble /@ Transpose[{frnsTM, frnsSC}], 
 ImageSize -> 900]

[Vitaliyk (Vitaliy Kaurov)]

A slightly different versions where widths variable is made wider for higher sorting values. I.e. the goal: keep the animation code essentially the same, still adding atoms in descending TopologicalStericEffectIndex , but let that same descriptor also make high-value atoms slightly wider. So topid now controls both the Gaussian height and, modestly, the Gaussian width.

ClearAll["Global`*"];

mol2 = Molecule[Entity["Chemical", "Idarubicin"]];

coords2da =
  MoleculeValue[
   mol2,
   "AtomDiagramCoordinates",
   IncludeHydrogens -> None
   ];

coords2daAll =
  MoleculeValue[
   mol2,
   "AtomDiagramCoordinates",
   IncludeHydrogens -> All
   ];

rot = RotationTransform[-Pi/2, Mean[coords2da]];

coords2daV = rot /@ coords2da;
coords2daAllV = rot /@ coords2daAll;

mol2V = Molecule[
   mol2,
   AtomDiagramCoordinates -> coords2daAllV
   ];

topid =
  MoleculeValue[
   mol2,
   "TopologicalStericEffectIndex",
   IncludeHydrogens -> None
   ];

order = Ordering[topid, All, Greater];

widths = Rescale[topid, MinMax[topid], {0.24, 0.45}];

{{xm, xp}, {ym, yp}} =
  MinMax[#, Scaled[0.1]] & /@ Transpose[coords2daV];

gauss[u_?NumericQ] := If[u > 40, 0., Exp[-u]];

hmfun2[n_Integer][x_?NumericQ, y_?NumericQ] :=
 With[
  {ii = order[[;; n]]},
  Total[
   MapThread[
    #2 gauss[
       SquaredEuclideanDistance[{x, y}, #1]/(2 #3^2)
       ] &,
    {coords2daV[[ii]], topid[[ii]], widths[[ii]]}
    ],
   Method -> "CompensatedSummation"
   ]
  ];

Manipulate[
 Show[
  DensityPlot[
   hmfun2[n][x, y],
   {x, xm, xp},
   {y, ym, yp},
   PerformanceGoal -> "Quality",
   AspectRatio -> Automatic,
   ColorFunction->"LightTemperatureMap",
   Frame -> None,
   PlotPoints -> 100,
   PlotRange -> All,
   PlotRangePadding -> 0,
   ImageSize -> 550
   ],
  MoleculePlot[
   mol2V,
   ColorRules -> {_ -> Directive[Black, Thickness[.005]]}
   ]
  ],
 {{n, Length[topid], "atoms included"}, 0, Length[topid], 1,
  Appearance -> "Labeled"}
 ]

frnsTM =
  ParallelTable[
   Show[
    DensityPlot[
     hmfun2[n][x, y],
     {x, xm, xp},
     {y, ym, yp},
     PerformanceGoal -> "Quality",
     AspectRatio -> Automatic,
     ColorFunction -> "LightTemperatureMap",
     Frame -> None,
     PlotPoints -> 100,
     PlotRange -> All,
     PlotRangePadding -> 0,
     ImageSize -> 550
     ],
    MoleculePlot[
     mol2V,
     ColorRules -> {_ -> Directive[Black, Thickness[.005]]}
     ]
    ],
   {n, 0, Length[topid], 1}
   ]; // AbsoluteTiming

SetDirectory[NotebookDirectory[]]

Export["mollights-widths.gif", frnsTM, ImageSize -> 550,"DisplayDurations"->Append[Table[.1,36],1]]