Amber Workshop - tutorial 4

(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 2006

Simulating a Solvated Protein that Contains Non-Standard Residues
(More Advanced Version)

In tutorials 1 to 3 all the systems we studied contained only standard amino acid or nucleic acid residues. As such we did not need to create units for non-standard residues and provide our own parameters. In this tutorial we will look at one method for creating a non-standard residue.

In this tutorial we will be aiming to setup a protein simulation of Plastocyanin in explicit solvent. In order to do this there are a number of things we will need to do:

This will take a bit of work but it is significantly easier to do in AMBER 8 than it used to be.

You should always read all of the header information in a pdb file since it often contains important information about disordered pairs etc. Often 1 pdb file can also contain a series of different structures / inhibitors etc for the same protein system. In this pdb we see that waters #187 and #183 form a disordered pair and so both should not be present.

Thus I will arbitrarily remove water 187 from the pdb file and keep 183. I will also remove all of the connectivity data at the end of the PDB as this will not be used by XLeap. Each of our water molecules here is an individual resiude. Normally we need to add TER cards between residues that are not bound together in a chain. Fortunately we don't need to do that for our crystallographic waters since XLeap is smart enough to know that each of our waters is a separate molecule. This is defined in the WAT unit.

Since pdb files do not distinguish between cysteine residues that are involved in bonds or other things (and hence which have no proton on the sulphur atom) we need to edit the relevant cysteine residues to correct for this. The residue name used by Leap for regular protonated cysteine residues is CYS, for deprotonated (-ve charge) it is CYM and for cysteine residues involved in disulphide bridges and other bonds it is CYX. Since cysteine 84 in plastcyanin is bonded to the copper ion we need to change the residue name of residue 84 from CYS to CYX. The same is true for histidine residues which can be protonated in the 'delta' position (HID), the epsilon position (HIE) or at both (HIP). Fortunately in plastocyanin this is fairly easy since there are only two histidine residues (37 and 87) both of which are bonded to the copper via the 'delta' nitrogen. Thus they must both be protonated on the epsilon nitrogen. We are making our own unit (HIC) for residue 37 but we initially want it protonated on the epsilon nitrogen so we change it from HIS to HIE here. We also need to change the other histidine residue 87 from HIS to HIE.

The next stage before we can use the PDB is to sort out the non-standard hydrogen names. One option would be to simply strip out the protons and have XLeap add them back in at standard position as we did in the previous two tutorials. In this example though we shall try to keep the positions specified in the PDB file. We can use an AMBER tool for this called protonate. Protonate is typically used to protonate pdb files that do not have protons present. It can also be used, however, to correct pdb files that don't conform to IUPAC proton naming conventions. We run it as follows:

$AMBERHOME/exe/protonate -k -d $AMBERHOME/dat/PROTON_INFO -i 1PLC_mod.pdb -o 1PLC_mod2.pdb -l protonate.out

The -k option changes the names but "keeps" the positions as specified in the pdb file. The -d option specified the PROTON_INFO file that is in the dat directory of the AMBER 8 installation. This specifies information on where to place the protons. Protonate also strips all of the unnecessary header info for us.

The line above should produce the following files: 1PLC_mod2.pdb, protonate.out.

We don't need the protonate.out file but it is a good idea to take a look at it to check for any errors. In my case it was empty and so there were no errors or mysterious protons.

The next stage is to work out how we are going to treat the copper atom. In this example we shall make the copper atom part of residue 37. Residue HIE 37 will then become a non-standard residue, since it is now a histidine atom with a copper atom attached. Hence we will edit the pdb and move the copper atom from the end of the protein and place it as the last atom of residue 37. We also need to change its residue name to HIC and its residue number to 37. We also need to change all of the other atoms in residue 37 from HIE to HIC.

ATOM    527  N   HIC    37       4.103  34.339  15.185
ATOM    528  H   HIC    37       3.437  33.878  14.758
ATOM    529  CA  HIC    37       4.992  33.545  15.990
ATOM    530  HA  HIC    37       5.838  33.996  16.025
ATOM    531  C   HIC    37       5.176  32.169  15.301
ATOM    532  O   HIC    37       4.427  31.772  14.434
ATOM    533  CB  HIC    37       4.395  33.294  17.344
ATOM    534 3HB  HIC    37       4.944  32.577  17.764
ATOM    535 2HB  HIC    37       3.447  32.920  17.176
ATOM    536  CG  HIC    37       4.163  34.407  18.291
ATOM    537  ND1 HIC    37       5.162  35.218  18.775
ATOM    538  CD2 HIC    37       2.952  34.795  18.848
ATOM    539  HD2 HIC    37       2.063  34.452  18.734
ATOM    540  CE1 HIC    37       4.572  36.081  19.598
ATOM    541  HE1 HIC    37       5.044  36.764  20.072
ATOM    542  NE2 HIC    37       3.236  35.857  19.674
ATOM    543  HE2 HIC    37       2.630  36.345  20.187
ATOM   1444  CU  HIC    37       7.050  34.960  18.716

For the methionine that is bound to the copper via a long bond (MET 92) we will create a new residue called MEM (a made up residue name that is not currently in use) that will have a sulphur of type SM rather than S. SM, is a type I just made up that we will use to create special force field parameters for MET 92's S-Cu bond. Thus we need to edit the pdb and change all of residue names for residue 92 from MET to MEM.

The original pdb file also contained several "alternate" configurations for some residues. E.g. for LYS 30 we had

ATOM    426  N   LYS    30      -0.930  27.774  20.957  1.00  8.07      1PLC 549
ATOM    427  CA  LYS    30      -2.028  28.602  20.421  1.00 10.86      1PLC 550
ATOM    428  C   LYS    30      -1.629  30.029  20.254  1.00  9.16      1PLC 551
ATOM    429  O   LYS    30      -1.226  30.665  21.243  1.00  7.63      1PLC 552
ATOM    430  CB ALYS    30      -3.201  28.489  21.415  0.50 13.41      1PLC 553
ATOM    431  CB BLYS    30      -3.250  28.517  21.354  0.50 15.09      1PLC 554
ATOM    432  CG ALYS    30      -4.397  29.366  21.249  0.50 16.84      1PLC 555
ATOM    433  CG BLYS    30      -4.600  28.495  20.646  0.50 21.50      1PLC 556
ATOM    434  CD ALYS    30      -5.681  28.891  21.893  0.50 20.64      1PLC 557
ATOM    435  CD BLYS    30      -5.745  28.171  21.589  0.50 24.43      1PLC 558
ATOM    436  CE ALYS    30      -5.527  28.212  23.225  0.50 23.18      1PLC 559
ATOM    437  CE BLYS    30      -5.585  26.973  22.460  0.50 24.88      1PLC 560
ATOM    438  NZ ALYS    30      -6.825  28.052  23.929  0.50 20.02      1PLC 561
ATOM    439  NZ BLYS    30      -5.971  25.681  21.860  0.50 26.52      1PLC 562
ATOM    440  H   LYS    30      -0.661  27.945  21.802  1.00  8.96      1PLC 563
ATOM    441  HA  LYS    30      -2.302  28.229  19.582  1.00 10.61      1PLC 564
ATOM    442 1HB ALYS    30      -3.380  27.572  21.665  0.50 14.09      1PLC 565
ATOM    443 1HB BLYS    30      -3.134  27.712  21.919  0.50 15.93      1PLC 566
ATOM    444 2HB ALYS    30      -2.700  28.880  22.297  0.50 13.96      1PLC 567
ATOM    445 2HB BLYS    30      -3.188  29.326  21.939  0.50 15.16      1PLC 568
ATOM    446 1HG ALYS    30      -4.210  30.316  21.598  0.50 18.44      1PLC 569
ATOM    447 1HG BLYS    30      -4.784  29.428  20.234  0.50 20.74      1PLC 570
ATOM    448 2HG ALYS    30      -4.609  29.538  20.265  0.50 17.98      1PLC 571
ATOM    449 2HG BLYS    30      -4.629  27.910  19.855  0.50 20.04      1PLC 572
ATOM    450 1HD ALYS    30      -6.278  29.712  22.058  0.50 21.47      1PLC 573
ATOM    451 1HD BLYS    30      -5.929  28.976  22.170  0.50 24.46      1PLC 574
ATOM    452 2HD ALYS    30      -6.224  28.346  21.263  0.50 21.94      1PLC 575
ATOM    453 2HD BLYS    30      -6.606  28.094  21.037  0.50 24.54      1PLC 576
ATOM    454 1HE ALYS    30      -5.138  27.298  23.130  0.50 22.69      1PLC 577
ATOM    455 1HE BLYS    30      -4.709  26.867  22.883  0.50 25.84      1PLC 578
ATOM    456 2HE ALYS    30      -4.956  28.735  23.845  0.50 23.28      1PLC 579
ATOM    457 2HE BLYS    30      -6.257  27.063  23.262  0.50 25.75      1PLC 580
ATOM    458 1HZ ALYS    30      -7.360  28.765  23.779  0.50 21.25      1PLC 581
ATOM    459 1HZ BLYS    30      -6.721  25.740  21.353  0.50 26.03      1PLC 582
ATOM    460 2HZ ALYS    30      -7.206  27.240  23.677  0.50 21.83      1PLC 583
ATOM    461 2HZ BLYS    30      -5.279  25.211  21.533  0.50 25.17      1PLC 584
ATOM    462 3HZ ALYS    30      -6.682  27.986  24.852  0.50 21.60      1PLC 585
ATOM    463 3HZ BLYS    30      -6.289  25.095  22.599  0.50 25.92      1PLC 586

By default Leap and protonate will use the A conformation and disregard the others. This is fine for our purposes here. If we specifically wanted to start out in one of the other confirmations we would need to remove the A conformation from the file.

The final thing we need to do is add a TER card between our last protein atom and our first crystallographic water atom.

If simply loaded our edited pdb file above into XLeap at this point the vast majority of the file would load ok. There would be problems, however, with our two non-standard residues HIC and MEM. Before we can successfully load our 1PLC_modified_final.pdb file into XLeap we need to tell it what our two non-standard units are.

We have several options on how to deal with this. The easiest option for a simple molecule would be to use Antechamber as you will do in tutorial 5. Antechamber provides an automated method for creating non-standard units. However, it only works with complete molecules, not fragments as we have here.

A second option would be to simply edit the residues within xleap this once and then disregard them when we quit. This is the quickest method but doesn't allow for any re-use of our residues. E.g if we wanted to create a second very similar protein we would have to go through and repeat the process of editing the residues in xleap all over again.

The third option, and the one we will use here, is to create new library files for each of our non-standard units. In this way we can re-use those units, by simply loading the library files into xleap before we load our protein pdb, for similar proteins.

So, if you haven't already, fire up xleap and load our 1PLC_mod_final.pdb file into a new unit called PLC.

$AMBERHOME/exe/xleap -s -f $AMBERHOME/dat/leap/cmd/leaprc.ff99
> 1PLC=loadpdb 1PLC_mod_final.pdb

We should check the messages xleap gives us to ensure it has done what we expected. It should report only two unknown residues, our HIC and MEM. It should also have added any missing protons. Since we already ran our pdb through protonate the only missing atoms should be the protons on the our 109 water molecules. Hence we expect xleap to add a total of 218 H atoms. Verify that this is indeed the case. If it isn't then you missed something during the pdb edit. There should also have been 35 atoms that were not in residue templates, these are 18 in our HIC residue and 17 in our MEM residue. This is just a test to check everything is okay. We can quit xleap now.

We now need to create library files for our MEM and HIC units. The quickest way to do this, since we need initial structures for these units, is to simply cut them out of the 1PLC_mod_final.pdb file and save them as their own pdb files.

Now we can go ahead and load xleap again (are you getting good at this yet? :-)) and load the two pdb files into their own units:

$AMBERHOME/exe/xleap -s -f $AMBERHOME/dat/leap/cmd/leaprc.ff99
> HIC = loadpdb hic.pdb
> MEM = loadpdb mem.pdb

Now we need to tell xleap about our non-standard residues. We need to tell it what is bonded to what, what the atom types are and what the point charges on the atoms are. Lets start with the modified methionine residue.

This should bring up our MEM residue in the editing window, you can turn on atom names with the Display menu.

As you can see xleap currently knows nothing about the bonding within this residue. We can correct that quickly by simply drawing in the bonds. On the manipulation bar near the top of the edit window select draw and then connect the atoms as they should be bonded within methionine (remember, hold middle button to rotate, and middle + right to zoom in and out).

H->N->CA->C->O
CA->HA
CA->CB->HB3
CB->HB2
CB->CG->SD
CG->HG2
CG->HG3
SD->CE->HE1
CE->HE2
CE->HE3

Note, once you are proficient with AMBER you can copy prep files from the default units and modify these and therefore save yourself a lot of time. For the purposes of this tutorial it easier to understand what is going on if we do it from scratch. Here is the finished product:

At this point we should save a library file of what we have done so far. In this way if we make a mistake we can simply quit, load xleap again, load the library file and carry on from where we were. Go back to the command window and type:

Now go back to the edit window. The next step is to specify the atom types and charges.

Typically you would need to calculate atomic charges for the copper atom and all of the other atoms in our HIC residue and also the MEM residue. This is necessary because the presence of the copper atom will alter the electron distribution in these surrounding residues. For AMBER this is done using the Restrained Electrostatic Potential Method (RESP). Details of calculating RESP charges can be found on the AMBER website. I will not cover this step in this tutorial, you will just have to trust me on the charge assignments. Note, the RESP fitting procedure has now been automated with a free program called R.E.D. Details of how to obtain and use R.E.D are available at http://www.u-picardie.fr/labo/lbpd/RED/.

For the purposes of this tutorial we will use the charges shown in the screen shot below. In order to specify the charges and atom types we need to select the entire unit. Hit the select button in the Manipulation bar and then rubberband the entire molecule. It should all change colour. Then go to the Edit menu:

You should now go through and assign an atom type and charge for all of the atoms. For the charges you can use the ones shown below that I have calculated for you. For the atom types I have simply used all of the same atom types as defined in a regular MET unit (try editing MET if you want to see them) with the exception of the sulphur atom SD which, as mentioned earlier, I have made type SM. So edit the table so it looks like this (or if you don't have time just read on and I will give you the library file at the end).

Once this is done we can select Table->Save and Quit. We can then close the edit window we should just be left with the command window. The last thing we need to do before we can save a library file for our completed unit is to tell xleap what is the head atom for this unit and what is the tail atom. This information is used in order to connect the protein back bone. E.g. if we type desc MEM we will see that there is no head or tail defined:

However, if we look at MET we shall see that the backbone N is defined as the Head and the backbone C is defined as the Tail:

> desc MET
UNIT name: MET
Head atom: .R<MET 1>.A<N 1>
Tail atom: .R<MET 1>.A<C 16>
.
.
.

We should now repeat the entire process for our HIC residue. To save you time though I have done it for you already. In creating the bonds I ensured that I bound the copper atom CU to the delta nitrogen ND1. The atoms types and charges I assigned were:

We are now almost done. We have our two non-standard residues defined we now just need to let xleap know about any missing parameters. Before we do this lets set up our protein in xleap with all the necessary bonds, in this way we can have xleap tell us what parameters are missing.

In order for xleap to know about our new residues when we load the 1PLC_mod_final.pdb file we need to make sure we load our two new library files first. These will define the HIC and MEM units so that when xleap encounters them in the pdb file it will know the topology, charges and atom types:

There should now be no error messages. Now before we go any further we have to ensure that all the bonds to our copper atom are defined. At the moment only the 'delta' nitrogen of the HIC residue is bonded to the copper since this is the only residue we built with this bond defined. We also need to add bonds between the copper and the cysteine (84) sulphur atom, between the copper and the MEM (92) sulphur atom and between the copper and the other histidine (87) delta nitrogen. We can do this with the bond command:

> bond 1PLC.87.ND1 1PLC.37.CU
> bond 1PLC.84.SG 1PLC.37.CU
> bond 1PLC.92.SD 1PLC.37.CU

This will create the four missing bonds. We can now go ahead and solvate our system in a truncated octahedral box of TIP3P water.

This will add a 12 angstrom buffer of water around our system. In addition to the already defined crystallographic waters. Next we charge neutralise our system since it currently has a charge of -8.0. We shall therefore add a total of 8 Na+ ions.

Now we can try and save our prmtop and inpcrd files. At this point we should find that xleap cannot find the type CU or SM.

This is to be expected since CU and SM were atoms types that we made up that currently don't exist in the standard PARM99 force field. If you check our 1PLC unit you should also find, along with a large number of close contacts due to the solvation, a number of missing parameters all corresponding to the new CU and SM atom types.

So, before we can go any further we need to add all of these parameters to the AMBER force field. Before we quit xleap however we shall save a 1PLC library file so that we don't have to repeat all of the steps we did above. We can just reload the library file.

In order to introduce the new parameters we have two options. We can modify the main force field files or we can create an frcmod file with the changes specific to this project. The second option is a much better choice since modifying the main files could lead to clashes with other people who use the same installation. Creating a set of parameters is a bit of an art, since one is often faced with unknown parameters (such as force constants for copper to its ligands as in this example). The purpose of this tutorial is simply to cover the mechanics of running AMBER and so I provide you all of the parameters below. Note, these are for tutorial purposes only, I make no claims about the validity or appropriateness of these parameters. There are a number of articles in the literature about parameter estimation, and users are encouraged to consult them when faced with unusual chemical environments. A starting point is section 12.1 of the AMBER 8 manual entitled Parameter Development.

plc.frcmod

#  modifications to force field for poplar plastocyanin  

MASS
SM 32.06
CU 65.36

BOND 
NB-CU  70.000   2.05000 #kludge by JRS    
CU-S  70.000   2.10000 #kludge by JRS    
CU-SM  70.000   2.90000 #for pcy    
CT-SM 222.000   1.81000 #met(aa)   

ANGLE   
CU-NB-CV   50.000     126.700  #JRS estimate    
CU-NB-CR   50.000     126.700  #JRS estimate   
CU-NB-CP   50.000     126.700  #JRS estimate   
CU-NB-CC   50.000     126.700  #JRS estimate   
CU-SM-CT   50.000     120.000  #JRS estimate   
CU-S -CT   50.000     120.000  #JRS estimate    
CU-S -C2   50.000     120.000  #JRS estimate    
CU-S -C3   50.000     120.000  #JRS estimate    
NB-CU-NB   10.000     110.000  #dac estimate    
NB-CU-SM   10.000     110.000  #dac estimate    
NB-CU-S    10.000     110.000  #dac estimate   
SM-CU-S    10.000     110.000  #dac estimate  
CU-SM-CT   50.000     120.000  #JRS estimate 
CT-CT-SM   50.000     114.700    #met(aa)  
HC-CT-SM   35.000     109.500   
H1-CT-SM   35.000     109.500  
CT-SM-CT   62.000     98.900  #MET(OL)  

DIHE 
X -NB-CU-X    1       0.000 180.000       3.000    
X -CU-SM-X    1       0.000 180.000       3.000   
X -CU-S -X    1       0.000 180.000       3.000   
X -CT-SM-X    3       1.000   0.000       3.000   

NONBON 
CU    2.20      0.200   
SM    2.00      0.200

Anything with a # in front of it is a comment. As you can see we specify the mass, missing bonds, angle and dihedrals as well as VDW parameters.

We can now load this file into xleap and it will add all of these parameters to the PARM99 force field we selected.

$AMBERHOME/exe/xleap -s -f $AMBERHOME/dat/leap/cmd/leaprc.ff99
> loadamberparams plc.frcmod
> loadoff 1PLC.lib

If you wish to look at the starting structure in vmd we can use ambpdb to create a pdb file for us:

We could now use these files to run plastocyanin simulations. You can try this yourself if you want. Just remember that this is in explicit solvent and a peridiodic box so you use periodic boundary conditions. You will also need to minimise the system initially to remove bad contacts. I would then heat it over 20ps from 0 to 300K with constant volume periodic boundaries before moving to a long equilibration at 300 K with constant pressure.

ATOM	TYPE	CHARGE
N	N	-0.3479
H	H	0.2747
CA	CT	-0.1354
HA	H1	0.1212
C	C	0.7341
O	O	-0.5894
CB	CT	-0.0414
HB3	HC	0.0810
HB2	HC	0.0810
CG	CC	-0.0012
ND1	NB	-0.1513
CD2	CW	-0.1141
HD2	H4	0.2317
CE1	CR	-0.0170
HE1	H5	0.2681
NE2	NA	-0.1718
HE2	H	0.3911
CU	CU	0.3866