Amber Workshop - tutorial 3

(Note: These tutorials are meant to provide illustrative examples of how to use the AMBER software suite to carry out simulations that can be run on a simple workstation in a reasonable period of time. They do not necessarily provide the optimal choice of parameters or methods for the particular application area.)
Copyright Ross Walker 2004

In order to obtain our reference states for glutamic acid and aspartic acid in solution we will need to run 4 minimisations. Once for GLU, once for GLH, once for ASP and once for ASH.

Now, running these residues on their own in solution is not really correct so what we will do is cap each one with ACE and NME residues. E.g our GLU will look as follows:

So, the question that arises is how can we quickly create such starting structures for our minimisations. Well, there are a huge number of ways in which we can do this but in this example we will stick to using the tools that come with Amber. One very neat trick that xleap has is that it can automatically add atoms that are missing from residues that it has templates for. This we will use to our advantage to quickly create the starting structures we need. We will start by cutting out the GLU and ASP residues from our FBP1_noH_min.pdb. This will give us reasonable starting structures for glutamic acid and aspartic acid. We will create 4 new pdb files in total, one for GLU, one for GLH, one for ASP and one for ASH.

So, for FBP1_noH_min.pdb we open it in an editor and cut out all of residue 7 (GLU).

ATOM     71  N   GLU     7      -3.473  -8.991  -6.525
ATOM     72  H   GLU     7      -3.580  -9.974  -6.761
ATOM     73  CA  GLU     7      -2.188  -8.352  -6.892
ATOM     74  HA  GLU     7      -1.575  -9.134  -7.340
ATOM     75  CB  GLU     7      -2.380  -7.273  -7.977
ATOM     76 2HB  GLU     7      -3.083  -6.518  -7.625
ATOM     77 3HB  GLU     7      -1.418  -6.790  -8.149
ATOM     78  CG  GLU     7      -2.867  -7.854  -9.314
ATOM     79 2HG  GLU     7      -2.122  -8.568  -9.673
ATOM     80 3HG  GLU     7      -3.798  -8.400  -9.142
ATOM     81  CD  GLU     7      -3.106  -6.785 -10.403
ATOM     82  OE1 GLU     7      -2.607  -5.639 -10.295
ATOM     83  OE2 GLU     7      -3.798  -7.101 -11.403
ATOM     84  C   GLU     7      -1.321  -7.810  -5.738
ATOM     85  O   GLU     7      -0.100  -7.730  -5.893

Now, we want to cap this with ACE at the nitrogen end and NME at the carbon end. However, we don't know the 3D coordinates for our ACE and NME caps. No worry we will get xleap to add these for us. What we can do therefore is just add the C of ACE and N of NME to the file. Leap will add the rest. Thus our new pdb file looks like this:

ATOM     70  C   ACE     6
ATOM     71  N   GLU     7      -3.473  -8.991  -6.525
ATOM     72  H   GLU     7      -3.580  -9.974  -6.761
ATOM     73  CA  GLU     7      -2.188  -8.352  -6.892
ATOM     74  HA  GLU     7      -1.575  -9.134  -7.340
ATOM     75  CB  GLU     7      -2.380  -7.273  -7.977
ATOM     76 2HB  GLU     7      -3.083  -6.518  -7.625
ATOM     77 3HB  GLU     7      -1.418  -6.790  -8.149
ATOM     78  CG  GLU     7      -2.867  -7.854  -9.314
ATOM     79 2HG  GLU     7      -2.122  -8.568  -9.673
ATOM     80 3HG  GLU     7      -3.798  -8.400  -9.142
ATOM     81  CD  GLU     7      -3.106  -6.785 -10.403
ATOM     82  OE1 GLU     7      -2.607  -5.639 -10.295
ATOM     83  OE2 GLU     7      -3.798  -7.101 -11.403
ATOM     84  C   GLU     7      -1.321  -7.810  -5.738
ATOM     85  O   GLU     7      -0.100  -7.730  -5.893
ATOM     86  N   NME     8

So, what I have done here is added a new residue, residue 6, called ACE before the glutamic acid and a new residue, residue 8, called NME after the glutamic acid. I have not specified any coordinates, however, leap will add these for us. We save this file as GLU_capped.pdb. We now still need GLH, ASP and ASH. So, for the GLH we simply change all the GLU's above to read GLH. Leap will again add the missing ACE and NME atoms as well as the missing GLH proton. Repeat the above process for the ASP and ASH residue (15) from the FBP1_noH_min.pdb pdb file.

Ok, now we need to create our topology and structure files. Let's start with GLU_capped.pdb and load this into xleap.

$AMBERHOME/exe/xleap -s -f $AMBERHOME/dat/leap/cmd/leaprc.ff99
GLU_capped=loadpdb GLU_capped.pdb

So a total of 3 missing heavy atoms were added and 7 missing hydrogens. This is exactly as we expect. However, we can't just save our topology and coordinate files at this point since xleap will have simply positioned the missing atoms based on it's internal templates, regardless of the surrounding residues. If we edit the structure we will see the problem:

You can zoom in and out by holding the middle and right buttons down at the same time. To rotate hold down the middle button.

Oh dear, this does not look good. All is not lost, however, since all we need is an approximate starting geometry since we will be minimising our system in sander. So, how can fix this geometry? Well, xleap contains an internal minimiser that we can use to clean up our structure. This minimiser will very quickly clean a structure. However, it only deals with internal degrees of freedom, it knows nothing of electrostatics or VDW forces and so we will still need to carefully minimise in sander afterwards. It should get us close enough, however, that our initial forces don't break the sander minimiser. Here's how we do it. Ensure that the current manipulation mode is set to "Select". Then by holding down the left mouse button "rubber band the entire structure". Every time you release the mouse button anything inside the band will be added to the current selection and change colour. To remove atoms from a selection hold down the shift key while rubber banding. So, go ahead and select everything, your screen should end up looking like this:

Then we go to the Edit menu and we select "Relax Selection". At this point the structure should move and you should get something similar to what is shown below. If you rotate the structure you should see that we still have some fairly bad van der waals contacts. If we were to just start doing MD with this structure the huge initial forces would, unless we used an exceptionally short time step, rapidly cause our system to blow up. Here, however, we will be minimising our system and so such large forces will rapidly be removed in the first few minimisation steps.

Let's, go ahead and create our topology and coordinate files. Close the editor by selecting "Unit"->"Close". then save the prmtop and inpcrd files:

Now do exactly the same thing with the GLH_capped, ASP_capped and ASH_capped files.

Now we need to minimise the structures in implicit solvent and extract the final energy. Since the structures we have created are quite far from the minima we will run a lot more steps of minimisation that we did before. The input file is below. I have chosen to do a total of 2,500 minimisation steps. 500 of steepest descent followed by 2,000 of conjugate gradient. This is probably much more than we need but since these are small molecules they will be very very quick to run and it won't do us any harm to do more steps than we need.

Rather than run each of these jobs by hand, you can if you want, I am going to create a short shell script that will run them all for us. It will grep out the final energy and print it to the screen. It will also create pdb's of the final structures which we should view in VMD to see if they look reasonable.

#!/bin/bash

#minimise
echo -n "GLU Solution: "
$AMBERHOME/exe/sander -O -i min_capped.in -o GLU_capped_min.out \
                         -p GLU_capped.prmtop \
                         -c GLU_capped.inpcrd \
                         -r GLU_capped_min.crd
#Obtain the energy
grep -C 7 "FINAL RESULTS" GLU_capped_min.out | grep "  2500 " | awk '{print$2}'

#create a pdb of the final structure
$AMBERHOME/exe/ambpdb -p GLU_capped.prmtop <GLU_capped_min.crd \
                         > GLU_capped_min.pdb 2>/dev/null

#now move on to the next structure
echo -n "GLH Solution: "
$AMBERHOME/exe/sander -O -i min_capped.in -o GLH_capped_min.out \
                         -p GLH_capped.prmtop \
                         -c GLH_capped.inpcrd \
                         -r GLH_capped_min.crd
#Obtain the energy
grep -C 7 "FINAL RESULTS" GLH_capped_min.out | grep "  2500 " | awk '{print$2}'

#create a pdb of the final structure
$AMBERHOME/exe/ambpdb -p GLH_capped.prmtop <GLH_capped_min.crd \
                         > GLH_capped_min.pdb 2>/dev/null

#now move on to the next structure
echo -n "ASP Solution: "
$AMBERHOME/exe/sander -O -i min_capped.in -o ASP_capped_min.out \
                         -p ASP_capped.prmtop \
                         -c ASP_capped.inpcrd \
                         -r ASP_capped_min.crd
#Obtain the energy
grep -C 7 "FINAL RESULTS" ASP_capped_min.out | grep "  2500 " | awk '{print$2}'

#create a pdb of the final structure
$AMBERHOME/exe/ambpdb -p ASP_capped.prmtop <ASP_capped_min.crd \
                         > ASP_capped_min.pdb 2>/dev/null

#now move on to the next structure
echo -n "ASH Solution: "
$AMBERHOME/exe/sander -O -i min_capped.in -o ASH_capped_min.out \
                         -p ASH_capped.prmtop \
                         -c ASH_capped.inpcrd \
                         -r ASH_capped_min.crd
#Obtain the energy
grep -C 7 "FINAL RESULTS" ASH_capped_min.out | grep "  2500 " | awk '{print$2}'

#create a pdb of the final structure
$AMBERHOME/exe/ambpdb -p ASH_capped.prmtop <ASH_capped_min.crd \
                         > ASH_capped_min.pdb 2>/dev/null

#All done
echo "ALL DONE"
echo ""

run_capped_min.x

Before we can execute this script we need to make it executable, then we can execute it with ./run_capped_min.x:

GLU Solution: -1.0188E+02
GLH Solution: -8.0227E+01
ASP Solution: -1.1390E+02
ASH Solution: -9.2899E+01

We can now complete our table of energies. But before we do we should quickly take a look at the final structures in VMD to check they look reasonable. Load VMD and then open each of the pdb files above in turn and make sure they look sensible.

Structure 1	PROTEIN NO H (A) -1529.5 KCal/mol	PROTEIN GLU 7 H (C) -1508.1 KCal/mol	PROTEIN GLU 10 H (E) -1506.8 KCal/mol	PROTEIN ASP 15 H (G) -1507.5 KCal/mol
Structure 2	PROTEIN NO H (B) -1516.4 KCal/mol	PROTEIN GLU 7 H (D) -1493.3 KCal/mol	PROTEIN GLU 10 H (F) -1494.6 KCal/mol	PROTEIN ASP 15 H (H) -1494.1 KCal/mol
		SOLUTION GLU NO H -101.88 KCal/mol	SOLUTION GLU H -80.227 KCal/mol
		SOLUTION ASP NO H -113.90 KCal/mol	SOLUTION ASP H -92.899 KCal/mol

Now we have all our energies we can compute the pKa shifts using the formula stated at the beginning:

D_pKa = -(D_G(X^- + H⁺ -> XH)_protein - D_G(X^- + H⁺-> XH)_solution) / k_BT*ln(10)

Bare in mind our estimates will not be accurate since we skipped minimising the proton that we added in the protein calculations. We have also not taken into account the lone H⁺ in each case. We also, in order to make the calculation time tractable to this tutorial, did not minimise our protein systems sufficiently. This tutorial should, however, have familiarised you with the type of steps you would go through to do such calculations.

	GLU 7	GLU 10	ASP 15
Structure 1	-0.18	0.76	0.73
Structure 2	1.05	0.11	0.94