ICM Manual v.3.8
by Ruben Abagyan,Eugene Raush and Max Totrov
Copyright © 2018, Molsoft LLC
Oct 18 2018

Contents
 
Introduction
Reference Guide
Command Line User's Guide
 ICM-shell
 ICM graphics
 Str.Analysis
 Sequence
 Molcart
 Pharmacophores
 Energy
 Molecules
 Animation
 Symmetry
 X-ray
 Plotting
  Plot simple
  Plot histogram
  Plot 3D 2Dfunction
  Plot 3D shape
  Peptide docking
 Docking/VLS
 Examples
 _chemSuper
 Chemical Conformation Generator
References
Glossary
 
Index
PrevICM Language Reference
How to plot
Next

[ Plot simple | Plot histogram | Plot 3D 2Dfunction | Plot 3D shape | Peptide docking ]

How to make a simple plot y=F( x)


If your function is determined on an uneven set of abscissa values, stored in array x, then use:
 
 plot display x y 
If your function is determined on a regular 1D-grid of values varying from x0 to x1, and is stored in the array y, use:
 
 plot display Rarray(Nof(y), x0, x1) y 
Use on-line help:
 
 help plot 
or to see the concise list of arguments
 
 help command plot 

Another way to make a file with a plot is to use the make plot command, e.g.


add column t Random(100,0.,1.) Random(100,0.,1.) Random(100,0.,1.)
make plot t "x=A;y=B;color=C;size=6;style=dots;;" output="aaa.png" size={500,500} # or pdf,eps,jpg
The type of file will be determined by the extension.

How to plot a histogram


Function Histogram returns the 2*N matrix which may be used on the fly with a command like
 
   plot display BAR red Histogram(a 50) 
More examples:
 
 a=Random(0. 100. 10000) 
 b=.04*(Count(1 50)*Count(1 50)) 
 annot = {"Histogram with quadratic ruler" "Random value" "N"} 
 plot display BAR red Histogram(a b) annot 
 
 b=Sqrt(100.*Count(1 100)) 
 annot = {"Histogram with square root ruler" "Random value" "N"} 
 plot display green BAR Histogram(a b) annot 


How to make a 3D-surface plot of a 2D-function


In the most common case you may have a matrix N*M my_matrix representing a two-dimensional function of x and y real arrays. First create a graphics object, then display it and save as a picture if you wish:
 
 read matrix "ram" 
                             # let's try a "wire" representation first 
 x=Rarray(Nof(ram),-180.,180.)  #this scale may be nonlinear 
 y=x 
 make grob "matwire" ram x y # scales for x,y. Def. z-scale=1. 
 display matwire 
                             # now let's have a look at a solid surface 
 make grob solid "matsolid" ram x y 0.9 # scales for x,y and z 
 display matsolid solid 
                             # adjust the view by rotating/translating 
#   connect matsolid         # you may move them with respect to each other 
                             # and save image 
 write image rgb "img_matsolid" 


How to create a new graphics object of a specific shape


First, use one of the provided graphics objects, stored in *.gro files to try these kind of ICM non-molecular objects. You can scale graphics objects by multiplying them by a real number and shifting (e.g. g1=10.*g + {1.,0.,0.}) and/or adjust their positions either by translate or rotate commands or by connecting to the selected grob and dragging/rotating them with the mouse in the graphics window, for example:
 
 read grob "icos" 
 my_icos = g_icos*0.2 
 display my_icos g_icos 
 translate add my_icos {5. 15. 0.} 
You may also create a new or edit an existing graphics objects (see provided *.gro files to learn about format).

Flexible peptide docking


Let us consider a problem of flexible peptide docking to a receptor. The peptide may be completely free or constrained (e.g. have a helical fragment in the middle), and the receptor will be assumed to have a fixed conformation. We will be going through a series of interactive or non-interactive steps. Feel free to put the commands below into several scripts and run them sequentially.
We will follow a particular example through all the steps. In this example we will be docking a constrained 11-residue alpha-helical peptide to a Bcl-XL molecule. It is a good idea to create a separate directory for the project, say, projectX , to keep all the files in one place.
Preparing the receptor binding surface for peptide docking.
First, we need to get the receptor model as an ICM object. It may come from several sources. The most common one is a pdb-structure converted with the convertObject macro, or the build model command which builds a model by homology, followed by regularization and refinement with the refineModel macro.
We will create two objects:
a_1.|the original pdb object with receptor with some other peptide
a_2.|converted receptor (without the peptide)

Example:
 
% cd projectX   
% icmgl -G 
 nice "1g5j"  # an open form of the receptor 
 delete a_w*  # delete waters but keep bound peptide 
 copy a_ "rec" 
 delete a_rec.2  # keep only the main chain 
 undisplay window 
 convertObject a_rec. yes yes yes no # also optimizes hydrogens 
# 
 display box Box( a_1.2 4. )  # display a box around the peptide 
Now you may rotate the scene and adjust the box size by pulling a box corner with the leftMouseClick. The box needs to cover the receptor surface you need to include in the simulation, but needs to be minimized to reduce the computational burden and memory use.
Sometimes the orientation of the box with respect to the receptor surface is not favorable. In that case, use connect a_*. and rotate the objects with the respect to the origin ( use axisLength = 15. and display origin to to see the axis in addition to seeing the box ). If you change the absolute position of the receptor model, you may want to save the new orientation of the pdb object (e.g. write pdb a_1. "template" )
The result Now we have two objects: (i) a pdb object in which we deleted the unnecessary molecules, but kept the peptide for comparison and convenience in the box definition; and (ii) an ICM-object with a receptor for which we will be calculating the grid maps.
Create grid potential maps
Now calculate five grid potentials with the make map potential command:
 
 make map potential "gc,gh,ge,gs,gb" Box() 0.5 
This command uses all molecules in the current object as a source for the map calculations. The Box() function will return the box you see on your screen. Alternatively you can just provide the 6 numbers defining that box.
0.5 is the grid size. Even though the ICM grid potential is extremely fast compared to other programs, for large boxes you may need to increase the grid cell size. We do not advise to go beyond 0.7A, however.
The result of this step is five grid maps around the docking surface of interest.
Building the peptide
The peptide of interest can be build with the following simple command:
 
 build string "ECLKRIGDELD"   # peptide sequence from Bax 
 rename a_ "pep" 
 set a_/* "HHHHHHHHHHH" 
 superimpose a_1.2/310:320/ca a_/1:11/ca align  

Saving the results
Let us now save both the objects and the maps for future use.
 
 write m_gc m_gh m_ge m_gs m_gb # gc.map .. files are created 
 write object a_*. "all"  # write 3 objects in one multi-object file 

Running a simple simulation from one starting conformation
The actual simulation does not need to be run interactively. An ICM script file _dockpep is described below. Run the script by
 
% _dockpep >! f1.out  # the output file 
The file contains the following commands:
 
 #!/usr/local/bin/icmng  
 read libraries  
 read object "all"  # contains a_ 
 set object a_pep.  # the peptide is the current object 
 fix v_/3:10/phi,psi,omg       # fix some backbone angles  
 nvar = Nof( v_//phi,psi,H,P ) # number of essential variables 
   # let us choose iteration limits depending on nvar  
 mncallsMC    = 10000 + Integer(0.008*nvar*nvar*nvar*nvar*nvar) 
 mncalls      = 170+nvar*3   
 print mncallsMC mncalls 
   # make sure that these two numbers are not insane  
   # or set smaller mncallsMC and mncalls (e.g. 10000 and 100) as a test 
 temperature  = 600   # optimal temperature for the simulation 
 tolGrad      = 0.001 # exit minimization when gradient is < 0.001 
 mnhighEnergy = 15 + nvar/2   # 
 set vrestraints a_/* # load preferred backbone and side-chain angle zones 
 # 
 vicinity = 2.5 
 compare a_pep. static  # use sRmsd of the peptide 
 set terms "vw,14,hb,el,to,en,gc,gh,ge,gb,gs"  
 minimizeMethod = "conjugate" 
 # 
 # maps 
 read map "gc,gh,ge,gb,gs" 
 GRID.margin = Distance( a_pep./1/ca a_pep./11/ca ) 
 # 
 # run montecarlo from a single start  
 # the reverse option will make more reasonables moves 
 # 
 # impose tethers ? 
 montecarlo reverse trajectory   
  
 # run montecarlo from multiple starts and merge the stacks 
  for i=1,Nof(<starting conditions>) 
   montecarlo reverse trajectory append 
  endfor 

This simulation created several files:
f1.cnfa stack file with mnconf (or less) best non-redundant conformations
f1.trja trajectory file with all accepted conformations (if you used the trajectory option)
f1.ou|the file with the text output of the script
The output file f1.ou contains information about every accepted conformation of the simulation, e.g.
 
... 
 __ __      300   5   asp BPMC   db   da   da    -32.42    -23.66  100   1.29 97810 
 DY Visi    300   1   phe BPMC   fb   fb   fb    -32.42    -32.42  101   0.57 97911 
 __ __      300   4   ser BPMC   sa   sg   sd    -32.42    -24.45   88   0.06 97999 
 DY Visi    300   5   asp BPMC  dmn  dmn  dmn    -32.42    -32.43   49   0.08 98048 
... 

Analyzing the results of a single simulation
Energy profile: Is the simulation long enough
The energies can be plotted to analyze the progression of the global energy optimization.
 
plotEnergy "f1.ou" 50.     
plotBestEnergy "f1.ou" 50. "display append"              
The first macro plots energies of each conformation generated and accepted in the course of the ICM stochastic global energy optimization. Due to the nature of the procedure, the energy may go up and down. The second macro plots (and appends the plot to the previously generated def.eps file) only the lowest energy achieved by a given iteration. This plot shows if and when a better conformation is found. Naturally, the energy can not be improved after the global minimum is found.
However, having a long flat plateau does not guarantee that the global minimum is found since a new drop of energy may come after a very long simulation. The time required for finding the global minimum in general depends on the number of variables, number of atoms (the cost of a single energy calculation) and complexity of the potential energy surface.
The best empirical way to make sure that the minimum is found is to compare the lowest energies found from completely different starting conformations (read about multiple start simulations below).
We use a window of 50 kcals/mole from the lowest energy. It helps to avoid the dependence of the scale on the initial energy which may be quite high.
To analyze the improvements only use plotBestEnergy .
Graphical analysis of the results
You can generate all stack conformations as independent ICM objects with the mkStackConf macro,e.g.
 
 read object "all"  # three objects, the peptide is the current object 
 mkStackConf 1 10   # make ICM objects from 10 best energy structures 
The objects will be called a_pep1. , a_pep2. , etc. Now these objects can be displayed simultaneously or separately and further analyzed, e.g.
 
 display a_pep*. 
 color a_pep?*. Energy( stack ) # in case you extracted all confs 
You can loop through the conformations interactively either with the original stack, e.g.
 
# create some view you like 
for i=1,10 # or Nof(conf) 
  load conf i 
  pause   # rotate and zoom, hit Return  
endfor 
, or use the objects you created with the mkStackConf macro.
Running a multiple start simulation
You can figure it out yourself :-) .
Convergence analysis
The key question we would like to ask is if at least two independent simulations from random starting conformations identified a nearly identical conformation with a similar energy as the lowest energy conformation. A low-tech example with just two simulations f1 and f2 :
 
 read object "all"  
 load conf 0 "f1"  # 0 is the lowest energy conformaiton.  
 load conf 0 "f2" 
You can also plot the best energy (see above) and compare the plots, as well as perform graphical analysis to see if the best conformations are similar.
Analyzing the results of a multiple start simulation, merging stacks
Use the read stack append command to merge all stack conformations together. Now you can redefine the compare command and the vicinity parameter depending on how you want to further compare and filter out the accumulated conformations. The compress command performs the compression.

Prev
Density correlation
Home
Up
Next
Docking/VLS

Copyright© 1989-2018, Molsoft,LLC - All Rights Reserved. Copyright© 1989-2018, Molsoft,LLC - All Rights Reserved. This document contains proprietary and confidential information of Molsoft, LLC. The content of this document may not be disclosed to third parties, copied or duplicated in any form, in whole or in part, without the prior written permission from Molsoft, LLC.