AMBER Archive (2009)

Subject: [AMBER] glycopeptide Parametrization

From: guardiani_at_fi.infn.it
Date: Wed Jun 10 2009 - 10:41:14 CDT

Dear Amber Experts,

I am currently trying to perform simulations of a 21-residues long
glycopeptide. The NMR structure of this peptide includes a beta-hairpin
(residues 1-14) followed by a short alfa-helix. The orientation of the
alfa-helix is rather variable, but in general it points away from the
beta turn. On top of the beta turn there is an Asparagine linked to
a beta-glucose unit through an N-glycosidic bond.

Before describing my calculation procedure, I anticipate a bizarre
result. The hairpin region produced by the REMD simulation in
implicit solvent (igb=7), superposes reasonably well to the corrisponding
NMR region. Conversely the region that was expected to adopt an helical
structure, runs parallel to the C-terminal beta-strand. I noticed that
in almost all conformations, this structure is stabilized by several
H-bonds between the oxydrils of glucose and the C-terminal COO- group
of the peptide. I therefore suspect that this kind of interactions may
be overstabilized by Glycam. Do you have any experience of this ?

First of all, here I show how I attained the topology of the glycopeptide.

STEP 1
-----------

I built with xleap a new unit called ASG comprizing a capped asparagine linked
to a beta-D_glucose molecule.

xleap -s -f $AMBERHOME/dat/leap/cmd/leaprc.ff99

source leaprc.glycam04
set ACE tail ACE.1.5
set ASN head ASN.1.1
set ASN tail ASN.1.13
set NME head NME.1.1
x = sequence { ACE ASN NME }
remove x x.2.11
set x tail x.2.10
set 1GB head 1GB.1.1
y = sequence { x 1GB }
savepdb y Asn_Glc.pdb
quit

STEP 2
------------

As in the file Asn_Glc.pdb the two hydrogens bound to atom C6 of glucose are
almost superposed, I manually corrected they coordinates. First of all
I measured the correct bond distances C6-H16 and C6-H26 and the bond angle
H16-C6-H26 in beta-D-glucose.

xleap -s -f $AMBERHOME/dat/leap/cmd/leaprc.ff99

source leaprc.glycam04
x = copy 1GB
savepdb x 1GB.pdb
quit

Using VMD I measured the parameters of interest:

vmd 1GB.pdb

measure bond {5 6} molid 0
measure bond {5 7} molid 0
measure angle {5 6 7} molid 0

If we now define the coordinates of the two hydrogens as (x1,y1,z1) and
(x2,y2,z2) and if we call theta the bond angle H16-C6-H26, and b the
C-H bond length, it is possible to establish the correct relative position
of the two hydrogens:

x2 = x1 + 2*b*sin(theta/2)
y2 = y1
z2 = z1

The edited ASG unit is saved in the file Asn_Glc_edit.pdb.

STEP 3
-----------

Using VMD, the Asn_Glc_edit.pdb file is saved in mol2 format:
Asn_Glc_edit.mol2

STEP 4
------------

Preparation of input file foe Gaussian run:

antechamber -fi mol2 -fo gzmat -i Asn_Glc_edit.mol2 -o Asn_Glc_edit.pdb

STEP 5
------------

Gaussian run to optimize the structure and to compute the ESP charges.
The input file for Gaussian has the form:

--Link1--
%Nproc=2
%chk=molecule
#HF/6-31G* SCF=Tight Test Pop=MK
iop(6/33=2) iop(6/42=6) opt
optimization of Asn_Glc

HERE I PUT THE Z-MATRIX OF ASN_GLC

The Gaussian job is started with the command:

bsub -W 24:00 -n 2 -i Asn_Glc_edit.com -o Asn_Glc_edit.log
-e Asn_Glc_edit.err g03

STEP 6
--------------

Here we generate an .ac file where the ESP charges are converted to RESP
charges, but the atomtypes are still those of the GAFF force field.

antechamber -fi gout -fo ac -i Asn_Glc_edit.log
-o Asn_Glc_edit.ac -c resp

STEP 7
---------------

Creation of a new .ac file where the GAFF atomtypes are converted to
the AMBER atomtype:

atomtype -i Asn_Glc_edit.ac -o Asn_Glc_edit_amber.ac -p amber

STEP 8
---------------

Preparation of a file mainchain.AsnGlc where in which we identify the atoms
of the capping groups as atoms to be removed.

HEAD_NAME N1
TAIL_NAME C6
MAIN_CHAIN C3
OMIT_NAME C7
OMIT_NAME H10
OMIT_NAME H11
OMIT_NAME H12
OMIT_NAME N3
OMIT_NAME H9
OMIT_NAME C2
OMIT_NAME O1
OMIT_NAME C1
OMIT_NAME H1
OMIT_NAME H2
OMIT_NAME H3
PRE_HEAD_TYPE C
POST_TAIL_TYPE N
CHARGE 0.0

STEP 9
---------------

Preparation of prepi file.

prepgen -i Asn_Glc_edit_amber.ac -o Asn_Glc_edit_amber.prepi
         -f prepi -m mainchain.AsnGlc -rn ASG -rf ASG.res

STEP 10
---------------

Preparation of frcmod file

parmchk -i Asn_Glc_edit_amber.prepi -o Asn_Glc_edit_amber.frcmod -f prepi

STEP 11
---------------

Creation of topology file for implicit-solvent (igb=7) REMD simulation.

xleap -s -f $AMBERHOME/dat/leap/cmd/leaprc.ff99SB

set default PBradii bondi
loadamberprep Asn_Glc_edit_amber.prepi
loadamberparams Asn_Glc_edit_amber.frcmod
x = sequence { NTHR PRO ARG VAL GLU ARG ASG GLY HIS SER VAL PHE LEU ALA
                PRO TYR GLY TRP MET VAL CLYS}
saveamberparm x csf114_glc_ext.top csf114_glc_ext.crd
quit

=========================================================================
                          REMD SIMULATION
=========================================================================

STEP 1: MINIMIZATION
----------------------

Minimization
  &cntrl
   imin = 1,
   maxcyc = 400,
   ncyc = 200,
   ntb = 0,
   igb = 7,
   gbsa = 1,
   saltcon = 0.2,
   cut = 999.0,
   rgbmax=999.0,
   ntpr = 10,
   ntx = 1,
  /

mpirun $sander -O -i minim.in -o minim.out -c csf114_glc_ext.crd -p
csf114_glc_ext.top -r minim.rst > out

STEP 2: EQUILIBRATION
------------------------

We consider 20 different temperatures chosen in geometric progression
T = 310, 324, 338, 353, 368, 384, 401, 418, 437, 456, 476, 497, 519,
541, 565, 590, 616, 642, 671, 700 K. Each replica is gradually heated
to the corresponding temperature.

Equilibration
  &cntrl
   imin = 0,
   cut = 999.0,
   rgbmax = 999.0,
   ig= 909,
   igb = 7,
   gbsa = 1,
   saltcon = 0.2,
   nstlim = 5000000,
   dt = 0.001,
   ntt = 3, gamma_ln = 1.0,
   tempi = 0.0, temp0 = 310.0,
   ntx = 1, irest = 0, ntb = 0,
   nscm = 50,
   ntpr = 100, ntwr = 100, ntwx = 100
  /

STEP 3: REMD RUN
--------------------

We performed REMD 250 ns of REMD run proceding through restarts of 10 ns.

MD run
  &cntrl
   imin = 0,
   cut = 999.0,
   rgbmax = 999.0,
   ig = 11377,
   igb = 7,
   gbsa = 1,
   saltcon = 0.2,
   nstlim = 2500,
   numexchg = 2000,
   dt = 0.001,
   ntt = 3, gamma_ln = 1.0,
   tempi = 310.0, temp0 = 310.0,
   ntx = 5, irest = 1, ntb = 0,
   nscm = 50,
   repcrd = 0,
   ntpr = 500, ntwr = 500, ntwx = 500
  /

What did I do wrong ? Is there anything that I may try to
improve the parametrization of the glycopeptide ? Is there
any known overstabilization of the polar group-charged group
interactions in GLYCAM ?

Thank you very much for your help.

Best regards,

Carlo Guardiani

----------------------------------------------------------------
This message was sent using IMP, the Internet Messaging Program.

_______________________________________________
AMBER mailing list
AMBER_at_ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber