|
Jul 1 2004
|
| Prev |
3.5 Energetics and electrostatics | Next |
[ vwplot_ht | rebel_ht | ps2 | be | ea | helicity | sm ]
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:
buildpep "one;one" # 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.}
This short script solves the Poisson equation by the "Rapid-Exact-
Boundary ELement (REBEL) method for crambin.
Examples:
electroMethod="boundary element"
read object "crn"
delete a_w* # get rid of water molecules
show energy "el"
show Energy("el")-r_out, r_out # Coulomb and solvation components
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 "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.))
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
The following is an example of calculating the average of an interatomic
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 movie 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
-
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 movie frame or stack conformation
and count number of helical residues. See macro _helicity averaging helicity
over the movie frames.
macro helicity s_movieName 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_movieName
e0=Energy("func") # the lowest energy
av=0.
ssum = 0.
r_temperature = r_temperature * Boltzmann
res = Real(Nof(a_/*))
read movie s_movieName
for i=1,Nof(frame)
load frame i s_movieName
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
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
|