[ Plotting vw faq | REBEL faq | PK shift | Binding energy | Ensemble average | Helicity | Stacks merge ]
How to plot the distance dependence of a van der Waals interaction |
The following script will plot three energy-interatomic distance plots for three possible
van der Waals terms defined by the
vwMethod
preference ("exact"-black,"soft"-blue and "old soft"-red). Simply mark and paste the following
lines into your ICM session:
build string "o;o" # two oxygens
set term "vw" only
set a_2//o Sum(Xyz(a_1//o ))+{.0 .0001 .0}
n=200
a=Rarray(n)
b=a
c=a
r=Rarray( n .03 3.)
vwSoftMaxEnergy = 14.5
for i=1,n
translate a_2 add {0.03 0. 0.}
#
vwMethod="exact"
r[i]=Distance( a_2//o a_1//o ) # use Sum(Distance(..)) for the old version
show ey mute
a[i]=Energy("vw")
#
vwMethod="old soft"
show ey mute
b[i]=Energy("vw")
vwMethod="soft"
show ey mute
c[i]=Energy("vw")
endfor
s=Sarray(3*n)
s[n+1]="_red line"
s[2*n+1]="_blue line"
plot ds r//r//r a//b//c s {0. 6.01 1. 1. ,-1. 17. 1. 1.}
How to calculate the electrostatic free energy by the REBEL-method |
This short script solves the Poisson equation by the "Rapid-Exact-
Boundary ELement (REBEL) method for crambin.
Examples:
electroMethod="boundary element"
read object s_icmhome+"crn"
delete a_w* # get rid of water molecules
show energy "el"
show Energy("el")-r_out, r_out # Coulomb and solvation components
To extract the surface polarization charge per atom use the Rarray (a_//*) function, e.g.
electroMethod = 4 # REBEL
show energy "el"
Rarray( a_//* ) # returns polarization charges
# if you want to correct the partial charges induced polarization
set charge a_//* Charge(a_//*) - Rarray( a_//* )
How to evaluate the pK shift |
Suppose we mutate a surface Asp into Ala and want to evaluate how the pK of
the neighboring His is changed. The pK shift may be evaluated as the
difference of potential at Nd1 His and Ne2 His nitrogens due to the
mutation. Since Ala may be considered as uncharged, the shift is simply
the potential at the nitrogens due to the Asp charge.
Example (pKshift of His34 from Asp22):
read object s_icmhome+"rinsr" # load insulin
# we assume that the positive charge is
# equally distributed between the two nitrogens
make boundary
pKshift=Sum({0.5, 0.5}*Potential(a_/his/nd1,ne2 a_/22/*))/ \
(Boltzmann*300.*Log(10.))
How to evaluate the binding energy |
There are many different approaches to the evaluation of binding energy.
One of the reasonable approximations has the following features:
- van der Waals/hydrogen bonding interaction is excluded since
it has close magnitudes for protein-protein and for protein-solvent
interactions;
- electrostatic free energy change is calculated by the
REBEL
method (see also the section
"How to calculate the electrostatic free energy ... ")
above);
- side-chain entropy change is calculated by standard ICM
entropic term
based on exposed surface area of flexible side-chains;
- hydrophobic energy change is calculated using
surface term
with constant surface tension of 20. cal/Angstrom.
Example:
electroMethod="boundary element"
surfaceMethod="constant tension"
surfaceTension=0.020
dielConst = 12.7
set terms "sf,el,en"
read object s_icmhome+"2ptc"
show energy a_1 a_1 mute
e1 =Energy("el,sf,en")
show energy a_2 a_2 mute
e2 =Energy("el,sf,en")
show energy mute
e12 =Energy("el,sf,en")
print "Binding energy = ", e12 - e1 - e2
How to calculate an ensemble average |
The following is an example of calculating the average of an inter-atomic
distance over a set of conformations collected in the
conformational stack.
This calculation is written as a macro.
Feel free to change it. You may also use trajectory file and
load frame instead of stack and load conf,
respectively.
# first, define the macro
macro ensembleAverage r_temperature
l_commands = no
l_info = no
load conf 0 # extract the lowest energy
e0 = Energy("func")
ansAver = 0. # the statistical sum initialization
statsum = 0.
r_temperature = r_temperature * Boltzmann
for i = 1,Nof(conf) # loop through all the stack
# conformations
load conf i
prob = Exp((e0-Energy("func"))/r_temperature)
# averaging distance between two ca
# is just an example
ansAver = ansAver + Distance(a_/2/ca a_/4/ca)*prob
statsum = statsum+ prob
endfor
r_out = ansAver/statsum
print " Ensemble average is: ", r_out
endmacro
# Now you can calculate your average
read object "my_peptide"
read stack # stack file is assumed to have
# the same name
ensembleAverage 600. # sometimes you may use the elevated
# temperature to account for relaxation
How to evaluate helicity of a peptide from the BPMC simulation |
-
Run the
_folding script first.
Make sure the procedure converges by running several simulations
(say _f1 _f2 _f3) from different random starting conformations. E.g.
cp $ICMHOME/_folding _f1 # adjust the script
icm _f1 > f1.ou &
cp _f1 _f2
icm _f2 > f2.ou &
cp _f2 _f3
icm _f3 > f3.ou &
-
You can evaluate helicity for each simulation. If they converge the result
will be about the same.
-
Helicity is just the ensemble average of the parameter which can be
calculated as the relative number of the helical residues. Therefore you need
to assign secondary structure for a particular trajectory frame or stack conformation
and count number of helical residues. See macro _helicity averaging helicity
over the trajectory frames.
macro helicity s_trjFileName r_temperature
# attention: 'temperature' is extremely important.
# You may use elevated temperature to account for relaxation.
l_commands=no
l_info=no
read conf 0 s_trjFileName
e0=Energy("func") # the lowest energy
av=0.
ssum = 0.
r_temperature = r_temperature * Boltzmann
res = Real(Nof(a_/*))
read trajectory s_trjFileName
for i=1,Nof(frame)
load frame i s_trjFileName
assign sstructure
prob = Exp((e0-Energy("func"))/r_temperature)
av = av + prob*Nof(Sstructure(a_1/*),"H")/res
ssum= ssum+prob
endfor
print " The best E=", e0, " Helicity= " av*100./ssum
endmacro
How to merge and compress several conformational stacks |
You may run several
montecarlo simulations and accumulate several
conformational stacks (
*.cnf files).
To unite them it is essential that they have been created with the same energy function,
because the compression algorithm takes the energy into account to decide which structure
is more valuable. If it is not the case, you can always recalculate energies for the
stack conformations by the following procedure:
read object "f"
read stack "f1"
read stack "f2" append # the second stack will be appended
for i=1,Nof(conf)
load conf i
show energy s_correctTerms # say, "vw,14,to,el,sf,en"
store conf i
endfor
Now, to unite all the stacks and compress them you may do the following (just the idea):
read object "f"
delete stack
read stack "f1" # first simulation
read stack "f2" append # second simulation
read stack "f3" append # third simulation
show stack # look at what you have now
compare v_//phi,psi # use the comparison criterion from the simulation script
# compose a new one, see also the compare command
vicinity = 40. # you may redefine the vicinity parameter
compress stack
show stack # look at the compressed stack
write stack "f4" # save the result