AMBER Archive (2009)

Subject: Re: [AMBER] glycopeptide Parametrization

From: guardiani_at_fi.infn.it
Date: Fri Jun 12 2009 - 09:06:24 CDT

Dear Amber users,

I would like to thank you all for your help and in particular
I would like to thank Francois Dupradeau. Your idea to simplify
the parametrization using only the Amber (ff99SB) force field
and the GLYCAM force field is quite convincing.

Now my parametrization is as simple as that (I am using the bondi
radii because I want to try an implicit-solvent simulation with
igb=7, but of course I will also run the explicit-solvent one):

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

set default PBradii bondi
source leaprc.glycam04
x = sequence { NVAL GLU ARG NLN GLY HIS CSER }
set x tail x.4.11
set 1GB head 1GB.1.1
y = sequence { x 1GB }
savepdb y y.pdb
saveamberparm y y.top y.crd
quit

Please let me know if this is what you meant. Which charges I will be
using with this parametrization ? Charges based on the Connolly surface
for the amino-acid residues and charges based on the CHELPG algorithm
for the sugar ?

As for the 1-4 scaling factors, I performed my previous simulations
using the default values, e.g.

SCNB = 2.0
SCEF = 1.2

Are these values appropriate for a glycopeptide simulation ?
If I explicitly set these values in the .in input files of
my REMD simulation, will the apply to both the peptide and
the sugar parts ?

Thank you very much again for your help.

Best regards,

               Carlo

Quoting FyD <fyd_at_q4md-forcefieldtools.org>:

> Hi Carlo,
>
> You provide a lot of information...
>
> Here you mix a peptide part + a sugar part + am1-cc charges + GB + ....
> These means you have an Amber force field for your protein (which one
> ?) + a force field for the sugar (glycam) + am1-bcc charges (with gaff
> ?)...
>
> Moreover, which 1-4 scaling factors did you use for the non-bonding
> interactions ? Amber FF and GLYCAM FF use different ones. Finally,
> charges in Amber FF are based on the Connolly surface while GLYCAM FF
> used the CHELPG algorithm.
>
> I am not sure what to answer here. In these conditions, I would try to
> "simplify" your computational conditions; i.e. using a minimum number
> of different force fields; mixing ONLY Amber + Glycam FFs using correct
> 1-4 scaling factors for the 1-4 interactions; This means avoiding
> am1-bcc charges and using RESP charges based on conformations you fully
> control. I would also use in a 1st step explicite solvent.
>
> regards, Francois
>
>
>> 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
>>
>
>
>
>
> _______________________________________________
> AMBER mailing list
> AMBER_at_ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber

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

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