Setup of Protein (Amber equilibration)

Here focus is on the maestro program (and how maestro can be double checked by Ulf Rydes programs). A complete list of the setup can be found under "Amber equilibrium" at http://signe.teokem.lu.se/~ulf/Methods/

Characterization of the pdb
Check (use open-office chart):
 * 1) Is the protein a monomer, dimer ? or n-mer?
 * 2) Are there Cys-Cys- crosslinks
 * 3) Are there missing parts of the structure (the first amino acids are often not resolved, check REMARK 465)
 * 4) Are the non-standard residues?
 * 5) At what pH was the structure obtained?

Protein setup with Maestro

 * 1) Make log-file (write down important information concerning  flipped residues etc)
 * 2) Open Prep Wiz and
 * 3) Import the pdb file
 * 4) Uncheck the default "Delete water beyond 5 Å" and hit preprocess
 * 5) Review "Problems => Current  overlapping atoms" If all close contacts are hydrogens, these will be optimized later and are likely not a problem in the end
 * 6) Review the "Problems => alternative positions" (This can be checked with Ulfs Progtrams, see XXX)
 * 7) Go to "Review and Modify" and then "Analyze Workspace"  and Delete undesired (co-crystallized solvent moelcules that are not water or other artefacts).
 * 8) Save a pdb before optimization
 * 9) Optimize H-positions (can take some minutes; usually 5-10 but might take up to several hours for big systems).
 * 10) Run Propka and check his residues (see also below). Check if Asp and Glu, Lys, Arg, and Tyr have pKa below 7
 * 11) Check histidines with changepdb (see below)
 * 12) Check burried charges through visualization and with changepdb (there should be none that are not involved in ionic pairs - excpetions for metal ligands)
 * 13) Decide on protonation state for Asp, Glu, Lys, Arg, Tyr and His
 * 14) Use "Interactive optimizer => Analyze network"
 * 15) Note (and check) which Asn and Gln groups that are flipped.
 * 16) Check His protonation states (HIE, HID or HIP). Does it fit with previous analysis from changepdb?. If not, then change them by hitting the arrow to the left side.
 * 17) Hit "optimize" (lock the histidines by checking the "lock" box. Protonation states of other residues that should remain as they are should also be checked as "locked")
 * 18) Save a pdbfile (HID, HIE, HIP must be specified afterwards - this can be done with changepdb command "hi")

Changepdb
Some important analysis can be done with changepdb (should be compared with maestro).
 * 1) Histidine analysis
 * 2) run changepdb
 * 3) give the input pdb
 * 4) command: hh (creates output to screen and files "hisX" with X = 1,2,3...).
 * 5) copy the output to screen into a file (e.g. named "hishb")
 * 6) grep after the relevant information (grep HHb hishb
 * 7) Inspect the histidine files and decide on protonation state (HIE, HID, HIP) - for highly negatively charged proteins solvent exposed histidines are generally HIP (unless there are good arguments against it). Their PKa values can also be calculated with PROPKA
 * 8) Burried charges
 * 9) Run changepdb
 * 10) Give input pdb
 * 11) Command bc
 * 12) Metal center (if present)
 * 13) Run changepdb
 * 14) Give input pdb
 * 15) command: des (gives file *.actX (X = 1,2,..) that contains an active site cut-out
 * 16) Change histidines
 * 17) Run changepdb
 * 18) Give input pdb
 * 19) run command "hi"
 * 20) type "D" for HID, "E" for HIE, "P" for HIP

Resp charges
Non-standard residues requires determination of charges and sometimes a "dummy" force field. More advanced force fields can be constructed if necessary (see PON parametrization)

1. Run a geometry optimization on the non-standard residues [typically B3LYP/def2-SV(P) or TPSS/def2-SV(P)]. Make sure that the structure is not too different from the crystal structure or other trustworthy experimental structures (if there are large deviations, consider using a single point calculation or only optimized hydrogens to generate charges)

2. With Turbomole (must be version > 7.0) insert "$esp_fit resp in the "control" file 'and run turbomole (Note: It can be convenient to ensure that the order of atoms is the same as in the pdb file)

ridft -resp -proper > logm 

The use of "-proper" indicates that no wave function optimization is required. Look for "GRID" in logm to check whether the Electrostatic Potential (which will be used to fit the resp charge) is generated.

3. Use now the script "changepot" (generates resp.pot, resp.in1, resp.in and qin). Note: changepot might give an error due to a missing line in logm. Insert this line with:

 awk '/charge for group #   1/{x++} x==2{sub(/charge for group #    1/,"writing ESP data for Amber to file esp.dat \n&")}1' logm > logm_2 && mv logm_2 logm 

changepot The files resp.in and resp.in1 are input files for the fitting of resp charges. qin are initial charges. Fitting of resp charges are usually done in two steps, first a fit involving all atoms (check this) with loose convergence. This is done with resp.in. In a next step, all heavy atoms are fixed while hydrogens in methylene (CH2) and methyl (CH3) groups are constrained to be identical and the fit is done with more tight convergence (see below).
 * 1) Select prgoram
 * 2) Press enter (new)
 * 3) name of turbomole output (logm)
 * 4) Do not use Boltzman weights (press enter)
 * 5) Press enter (special requirements if # grid pints > 99999)
 * 6) Press enter (output: resp.pot)
 * 7) Press enter (resp.in / resp.in1)
 * 8) Press enter (qin)

4. Check if there are any polar hydrogens bound to the same atom (e.g. water). These should be constrained to have the same charge (this is done in resp.in, see below for the format):

Format of resp.in files

&cntrl  ihfree=0, qwt=0.0005, iqopt=2 # qwt = optimization parameter

&end

1.0 Mol1

2  63 # charge and tot. # of atoms 7   0 # Atom nr. (here 7 = N) and command (-1 = freeze, 0 = optimized freely, n = fix to the same as atom #n e.g. 43 means fix to same value as atom number 43) 6    0</tt> 6    0</tt>

5. Run resp

resp -O -i resp.in -o resp.out -q qin -e resp.pot</tt>

6. Open resp.in1: Fix all charges except carbon atoms with more than one hydrogen (usually CH3 or CH2 groups). Fix the the hydrogens to have the same charge). qwt should be 0.001. Save qout in qin1:

cp qout qin1</tt>

and run

resp -O -i resp.in1 -o resp.out1 -q qin1 -e resp.pot</tt>

The final charges are now in resp.out1

Residues that should be renamed

 * 1) Cys => Cym if Cys is a metal ligand
 * 2) Cys => CYX if Cys is a Cys-Cys cross link
 * 3) HIS => HID, HIE or HIP if Maestro was used for the setup

tleap - add non-polar hydrogens and setup amber files
Data required for solvating the protein can be obtained through changepdb:
 * 1) Solvent cap

changepdb</tt>

 c # cap </tt>

 m # move </tt>

 w # write (overwrite) </tt>

2. Run tleap

tleap -f leap.in</tt>

A sample leap.in file is provided below:

source leaprc.ff14SB # Amber protein force field</tt>

loadAmberPrep hic.in # special residue prep.in file </tt>

loadAmberPrep cu2.in # special residue prep.in file </tt>

loadAmberParams hic.dat # special residue data file </tt>

loadAmberParams cu2.dat # special residue data file </tt>

x=loadpdb lpmo_setup_part_3_NH2corrected.pdb # pdb input file </tt>

<tt>bond x.56.SG x.178.SG # Cys-Cys cross-link </tt>

<tt>solvatecap x TIP3PBOX {0.0 0.0 0.0} 40 # solvation (see above) </tt>

<tt>saveamberparm x prmtop prmcr</tt>d

<tt>quit</tt>

There will likely be errors from tleap. A few examples on how to fix these are given below:

End-groups (OXT and AMT)
If one wishes to have the end amino group not charged, then delete the lines in addPdbResMap (copy file leaprc.ff03 to the local directory; "cp $AMBERHOME/dat/leap/cmd/leaprc.ff03 ."

Prep.in files
For special residues, a prep.in file must be constructed for the atom-types. A pdb file with the residue can be cut out from the pdbfile. The prep file can be constructed by (amber)

<tt> antechamber -i hic.pdb -fi pdb -o hic.in -fo prepi -nc 0 -c bcc -rn HIC -at amber # AMBER force field </tt>

<tt> antechamber -i hic.pdb -fi pdb -o hic.in -fo prepi -nc 0 -c bcc -rn HIC -at gaff # GAFF force field </tt>

In some cases, neither the regular amber protein force field nor the GAFF force field have parameters. In this case, the parameters must be constructed.

Check force field by energy calculation (Changeparm)
It is a good idea to check the resulting energies with

<tt> changeparm </tt>

<tt> prmtop </tt>

<tt> m (mdrest) </tt>

<tt> prmcrd </tt>

Energies above 1000 kcal/mol should be checked

Run simulated annealing
A simulated annealing run with shake

<tt>Title</tt>

<tt> &cntrl</tt>

<tt>  irest=0,ntx=1,</tt>

<tt>  nstlim=1200000,dt=0.002,</tt>

<tt>  temp0=300.0,ntt=1,tautp=0.2,</tt>

<tt>  ntc=2,ntf=2,</tt>

<tt>  nsnb=15,cut=15.0,dielc=1.0,</tt>

<tt>  ntpr=1000,ntwx=0,ntwv=0,ntwe=0,</tt>

<tt>  ntb=0,ntp=0,taup=0.2,</tt>

<tt>  ipol=0,igb=0,</tt>

<tt>  vlimit=20.0,nmropt=1,</tt>

<tt>  ibelly=1,bellymask=":WAT | @H="</tt>

<tt> &end</tt>

<tt> &wt type='TEMP0',istep1=     0,istep2= 400000,value1=370.00,value2=370.00  &end</tt>

<tt> &wt type='TEMP0',istep1= 400001,istep2=1200000,value1=370.00,value2=  0.00  &end</tt>

<tt> &wt type='TAUTP',istep1=     0,istep2= 400000,value1=  0.20,value2=  0.20  &end</tt>

<tt> &wt type='TAUTP',istep1= 400001,istep2= 800000,value1=  1.00,value2=  1.00  &end</tt>

<tt> &wt type='TAUTP',istep1= 800001,istep2=1000000,value1=  0.50,value2=  0.50  &end</tt>

<tt> &wt type='TAUTP',istep1=1000001,istep2=1200000,value1=  0.05,value2=  0.05  &end</tt>

<tt> &wt type='END' &end</tt>

<tt> &rst iat=0 &end</tt>

After this, run a minimization with 10000 steps

<tt>Title

<tt> &cntrl </tt>

<tt>  irest=0,ntx=1, </tt>

<tt>  nstlim=0,dt=0.002, </tt>

<tt>  imin=1,maxcyc=10000,drms=0.001, </tt>

<tt>  temp0=300.0,ntt=1,tautp=0.2, <tt>  ntc=2,ntf=2, </tt>

<tt>  nsnb=15,cut=15.0,dielc=1.0, </tt>

<tt>  ntpr=1000,ntwx=0,ntwv=0,ntwe=0, </tt>

<tt>  ntb=0,ntp=0,taup=0.2 </tt>

<tt>  ipol=0,igb=0, </tt>

<tt>  ncyc=10,ntmin=1,dx0=0.01, </tt>

<tt>  ibelly=1,bellymask=":WAT | @H=" </tt>

<tt> &end </tt>

Known problems

 * 1) Problems with SHAKE: In the case that SHAKE complains, it can often be relieved by running a short minimization followed and short MD without SHAKE (see input "sander.in00" and sander.in0) and then continue the annealing with SHAKE. If the problem persists, then continue without SHAKE

A simulated annealing run without shake (sander.in1_noshake)

<tt>Title</tt>

<tt> &cntrl</tt>

<tt>  irest=0,ntx=1, # The time-step is lowered</tt>

<tt>  nstlim=1200000,dt=0.0005,</tt>

<tt>  temp0=300.0,ntt=1,tautp=0.2,</tt>

<tt>  ntc=1,ntf=1,# changed to run without SHAKE</tt>

<tt>  nsnb=15,cut=15.0,dielc=1.0,</tt>

<tt>  ntpr=1000,ntwx=0,ntwv=0,ntwe=0,</tt>

<tt>  ntb=0,ntp=0,taup=0.2,</tt>

<tt>  ipol=0,igb=0,</tt>

<tt>  vlimit=20.0,nmropt=1,</tt>

<tt>  ibelly=1,bellymask=":WAT | @H="</tt>

<tt> &end</tt>

<tt> &wt type='TEMP0',istep1=     0,istep2= 400000,value1=370.00,value2=370.00  &end</tt>

<tt> &wt type='TEMP0',istep1= 400001,istep2=1200000,value1=370.00,value2=  0.00  &end</tt>

<tt> &wt type='TAUTP',istep1=     0,istep2= 400000,value1=  0.20,value2=  0.20  &end</tt>

<tt> &wt type='TAUTP',istep1= 400001,istep2= 800000,value1=  1.00,value2=  1.00  &end</tt>

<tt> &wt type='TAUTP',istep1= 800001,istep2=1000000,value1=  0.50,value2=  0.50  &end</tt>

<tt> &wt type='TAUTP',istep1=1000001,istep2=1200000,value1=  0.05,value2=  0.05  &end</tt>

<tt> &wt type='END' &end</tt>

<tt> &rst iat=0 &end</tt>