REAL META-PROGRAM FOR
PROTEIN 3D STRUCTURE MODELLING
#
CovBond := proc(d,i,j) k50/DIST(i,j)^12 + k51/DIST(i,j)^6 end:
##############
# k50 .. k59 #
##############
#
# k50 k51
# --------- + --------
# 12 6
# (a - b) (a - b)
# Handbook of Chemistry and Physics F-217, F-224
eqns := { subs( dd = 1887/1000, diff( subs(DIST(i,j)=dd,CovBond(dd,i,j)), dd)),
subs( DIST(i,j) = 1887/1000, CovBond(dd,i,j)) = -10158/100 }:
assign(solve( eqns, {k50,k51} ));
#
# Extra dimensions squeezing
##############
# k64 .. k69 #
##############
# DIM
# -----
# \ 2
# ) k[60 + i] x[i]
# /
# -----
# i = 4
for i from 4 to DIM do k6.i := 1/1000 od:
i := 'i':
#
# Interior/Surface interactions
##############
# k70 .. k79 #
##############
#
# Interior-Interior attraction
# k70 k71
# -------- + --------
# 4 2
# (a - b) (a - b)
InterInter := proc(d,i,j) k70/DIST(i,j)^4 + k71/DIST(i,j)^2 end:
subs(DIST(i,j)=dd,InterInter(dd,i,j)):
eqns := subs( dd=7, { %=5, diff(%,dd)=0 } ):
assign(solve( eqns, {k70,k71} ));
#
# Interior-Surface repulsion
# k72
# --------
# 2
# (a - b)
InterSurf := proc(d,i,j) k72/DIST(i,j)^2 end:
assign( solve( { subs(DIST(i,j)=RadiusR[i]+RadiusR[j],InterSurf(dd,i,j))=5 },
k72));
###########################################################
printf( `Model constants:\n` );
printf( `Distance penalty, k10=%g, k11=%g\n`, k10, k11 );
printf( `van der Waals constants: k20=%g, k21=%g\n`, k20, k21 );
printf( `Charge potential: k30=%g, Charge H+=%g, Charge O-=%g\n`,
k30, k31, k32 );
printf( `alignment = exp( %g*(AlignNHO-1) )\n`, k37 );
printf( `O-H (same aa) factor=%g\n`, k38 );
printf( `OH - Charged residues factor=%g\n`, k39 );
printf( `Covalent bond (van der Waals): k50=%g, k51=%g\n`, k50, k51 );
printf( `Extra dimensions squeezing:` );
for i from 4 to DIM do printf( ` k6%d=%g`, i, k6.i ) od;
printf( `\n` );
printf( `Interior-Interior attraction: K70=%g, k71=%g\n`, k70, k71 );
. . . . .
############################
# Models of peptide chains #
#########################################################################
# #
# ModelII has the following characteristics: #
# #
# 3 blobs per amino acid #
# CA(i) = xyz[3*i-2] #
# R(i) = xyz[3*i-1] (collection of atoms of the side chain) #
# O(i) = xyz[3*i] #
# #
# The C', N and H atoms are in the plane defined by CA(i), O(i) #
# and CA(i+1) and are computed based on these points #
# #
# The plus charge of the NH is assumed to be +k31 and at H(i) #
# The negative charge of the C=0 is assumed to be k32 at O(i) #
# #
# Radii of side chains are taken from the volumes (as spheres) #
# Distance from centre of side chains to CA taken from a table #
# #
# Charged amino acids, +- 1 taken at the centre of the side chain #
# Disulphide bonds, formed from the centre of the side chains #
# #
#########################################################################
# The number of blobs includes a fictitious CA at the end of the chain
NPEPT := nops(Seq):
NBLOB := 3*NPEPT+1:
#
# Planarity of the CA, C', N and CA atoms
#
# except for Proline, the following distances should hold:
#
#
# N O
# \ //
# \ //
# CA ------ C' .
# / \ .
# / \.
# / (1) .\ (2)
# CB . \
# / . N ------ CA
# / /
# / /
# R /
# H
#
Digits := 14:
CA1x := 0:
CA1y := 0:
C1x := 1.53:
C1y := 0:
# CA - C' - N angle is 115 deg
N2x := evalf( C1x + 1.325*cos( (115-180)*2*Pi/360 )):
N2y := evalf( 1.325*sin( (115-180)*2*Pi/360 )):
# distance from N to CA
eq1 := (N2x-CA2x)^2 + (N2y-CA2y)^2 = 1.453^2:
# distance from C' to CA(2)
eq2 := (C1x-CA2x)^2 + (C1y-CA2y)^2 = 2.427^2:
fsolve( {eq1,eq2}, {CA2x,CA2y} ):
assign(%);
# CA - C' - O angle is 120.5
O1x := evalf( C1x + 1.230*cos( (180-120.5)/360*2*Pi )):
O1y := evalf( 1.230*sin( (180-120.5)/360*2*Pi )):
# CA - N - H angle is 115
aH := arctan( (CA2y-N2y) / (CA2x-N2x) ) - evalf( 115/360*2*Pi ):
H2x := N2x + 1.0 * cos(aH):
H2y := N2y + 1.0 * sin(aH):
# position of C' in terms of CA1, O, CA2
eq1 := a1*CA1x + a2*O1x + (1-a1-a2)*CA2x = C1x:
eq2 := a1*CA1y + a2*O1y + (1-a1-a2)*CA2y = C1y:
convert( fsolve( {eq1,eq2}, {a1,a2} ), rational, exact):
C := subs( %, proc(i) a1*CA(i) + a2*O(i) + (1-a1-a2)*CA(i+1) end ):
# position of N in terms of CA1, O, CA2
eq1 := a1*CA1x + a2*O1x + (1-a1-a2)*CA2x = N2x:
eq2 := a1*CA1y + a2*O1y + (1-a1-a2)*CA2y = N2y:
convert( fsolve( {eq1,eq2}, {a1,a2} ), rational, exact):
N := subs( %, proc(i) a1*CA(i-1) + a2*O(i-1) + (1-a1-a2)*CA(i) end ):
# position of H in terms of CA1, O, CA2
eq1 := a1*CA1x + a2*O1x + (1-a1-a2)*CA2x = H2x:
eq2 := a1*CA1y + a2*O1y + (1-a1-a2)*CA2y = H2y:
convert( fsolve( {eq1,eq2}, {a1,a2} ), rational, exact):
H := subs( %, proc(i) a1*CA(i-1) + a2*O(i-1) + (1-a1-a2)*CA(i) end ):
# position of CA2 in terms of CA1, O, C'
eq1 := a1*CA1x + a2*O1x + (1-a1-a2)*C1x = CA2x:
eq2 := a1*CA1y + a2*O1y + (1-a1-a2)*C1y = CA2y:
convert( fsolve( {eq1,eq2}, {a1,a2} ), rational, exact):
# position of CB1 in terms of CA1 and R
CB := proc(i) CA(i) + 153/100/DistCA_SC[i]*(R(i)-CA(i)) end:
# Position of CB(1) in CB(1), CA(1), C'(1) plane
CB1x := 153/100 * cos(-AngleCpr_Cbeta[aanum]*2*Pi/360):
CB1y := 153/100 * sin(-AngleCpr_Cbeta[aanum]*2*Pi/360):
# Position of N(1) in CB(1), CA(1), C'(1) plane
N1x := 1453/1000*cos(1/180*AngleN_Cpr[aanum]*Pi):
N1y := 1453/1000*(cos(1/180*AngleCpr_Cbeta[aanum]*Pi)*
cos(1/180*AngleN_Cpr[aanum]*Pi)-cos(1/180*AngleN_Cbeta[aanum]*Pi))/
sin(1/180*AngleCpr_Cbeta[aanum]*Pi):
N1z := -1453/1000*((cos(1/180*AngleN_Cpr[aanum]*Pi)^2-2*cos(1/180*
AngleCpr_Cbeta[aanum]*Pi)*cos(1/180*AngleN_Cpr[aanum]*Pi)*cos(1/180*
AngleN_Cbeta[aanum]*Pi)+cos(1/180*AngleN_Cbeta[aanum]*Pi)^2-1+cos(1/180*
AngleCpr_Cbeta[aanum]*Pi)^2)/(cos(1/180*AngleCpr_Cbeta[aanum]*Pi)^2-1))^(1/2):
# Position of N1 in terms of CA, C' and CB (R)
# N-CA = aN * Perpen( C'-CA, CB-CA ) + bN * (C'-CA) + cN * (CB-CA)
N1 := evalm( aN * Perpen( [C1x,0,0], [CB1x,CB1y,0] ) +
bN * array([C1x,0,0]) + cN * array([CB1x,CB1y,0]) -
array([N1x,N1y,N1z]) ):
N1inCACR := solve( convert({N1[1],N1[2],N1[3]},rational,exact), {aN,bN,cN} ):
# distance equations (for each possible i)
d := DistCA_SC[i]:
distances :=
# distances to side chain
PotDistance( d, R(i), CA(i) ) +
# DistCpr_R[i] = sqrt( 1.53^2 + d^2 -\
2*1.53*d*cos( AngleCpr_Cbeta[Seq[i]]*2*Pi/360 ));
PotDistance( DistCpr_R[i], R(i), C(i) ) +
# DistN_R[i] = sqrt( 1.453^2 + d^2 -\
2*1.453*d*cos( AngleN_Cbeta[Seq[i]]*2*Pi/360 ));
PotDistance( DistN_R[i], R(i), N(i) ) +
# distance to main plane (CA,O,CA)
PotDistance( sqrt( (CA1x-O1x)^2 + (CA1y-O1y)^2 ), CA(i), O(i) ) +
PotDistance( sqrt( (CA1x-CA2x)^2 + (CA1y-CA2y)^2 ), CA(i), CA(i+1) ) +
PotDistance( sqrt( (O1x-CA2x)^2 + (O1y-CA2y)^2 ), O(i), CA(i+1) ) +
# other distances
# DistN_Cpr[i] = sqrt( 1.453^2 + 1.53^2 - 2 * 1.453 * 1.53 *\
cos( AngleN_Cpr[Seq[i]]*2*Pi/360 ));
PotDistance( DistN_Cpr[i], N(i), C(i) ):
distances := convert(distances,rational,exact):
# No van der Waals interaction within an amino acid
# is considered significant
# van der Waals interaction (for each amino acid i against j)
vdWij := 0:
ai := [ [H,12/10], [N,15/10], [CA,2], [C,2], [O,14/10], [R,RadiusR[i]]]:
aj := [ [H,12/10], [N,15/10], [CA,2], [C,2], [O,14/10], [R,RadiusR[j]]]:
for a1 in ai do for a2 in aj do
vdWij := vdWij + VanDerWaals( a1[2]+a2[2], a1[1](i), a2[1](j) );
od od:
# Coulumb's law within one amino acid
Couli := PotCharge( ChargeR[i]*k31, R(i), H(i) ) +
k39*PotCharge( ChargeR[i]*k32, R(i), O(i) ) +
k38*PotCharge( k31*k32, H(i), O(i) ): # these will not form a H-bond
# Coulumb's law within two amino acids
Coulij := k39*PotCharge( ChargeR[i]*k31, R(i), H(j) ) +
k39*PotCharge( ChargeR[i]*k32, R(i), O(j) ) +
PotCharge( ChargeR[i]*ChargeR[j], R(i), R(j) ) +
PotCharge( k31*k31, H(i), H(j) ) +
PotCharge( k31*k32, H(i), O(j) ) * AlignNHO(i,j) +
k39*PotCharge( k31*ChargeR[j], H(i), R(j) ) +
PotCharge( k32*k31, O(i), H(j) ) * AlignNHO(j,i) +
PotCharge( k32*k32, O(i), O(j) ) +
k39*PotCharge( k32*ChargeR[j], O(i), R(j) ):
# the Coulumb's potential between amino acids i and i+1
# must exclude the potential between O(i) and H(i+1)
Couli1 := subs( j=i+1, Coulij) - PotCharge( k32*k31, O(i), H(i+1) ) *
AlignNHO(i+1,i):
# Disulphide covalent bond
CovBondij := CovBond(dd,R(i),R(j)) - VanDerWaals(RadiusR[i]+RadiusR[j],
R(i), R(j) ):
# Dimension squeezing for amino acid i
DimSqueeze := sum( 'INDEX( k6.k*(CA(i)^2+O(i)^2+R(i)^2),k)', k=4..DIM ):
# Position of N(i) (wrt CA, CB and C')
PositionN := PotDistance( 0, N(i)-CA(i), aN[i]*Perpen3( C(i)-CA(i),
CB(i)-CA(i) ) + bN[i]*(C(i)-CA(i)) + cN[i]*(CB(i)-CA(i)) ) / k10 * k11: