pKa Calculations using Thermodynamic Integration

(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

pKa Calculations using Thermodynamic Integration - SECTION 3

pKa Calculations Using Thermodynamic Integration

By Ross Walker & Mike Crowley

3) Setup the prmtop and inpcrd files for the AspH and Asp- thioredoxin protein.

Now that we have calculated the delta-G value for our model compound we have to repeat the process for the thioredoxin protein. Here is the pdb file for the protein.

thioredoxin_chaina.pdb

The original crystal structure (2TRX) for this protein contained more than 1 chain and also a number of alternate conformations for some of the residues. In the pdb file given above I have already extracted just chain A and have also selected conformation A for all duplicated residues. There is still some initial preparation we need to do however.

Firstly we need to check for any disulphide bonds. We can do this visually in VMD. Open the pdb file and have it display all cysteine residues in ball and stick:

As you can see there are only two cysteine residues in this protein and there is obviously a disulphide bond between them so we need to edit the pdb file and change both of these to CYX. We will also have to manually add this bond in Leap when we create the prmtop and inpcrd files.

Next we need to specify the protonation state of any histidine residues. As you can see below there is only one histidine residue in this protein and it is exposed to solvent and a long way from the ASP26 residue who's pKa we will be calculating.

Thus the protonation state of this histidine is unlikely to have much influence on the results we get from our calculation. For the purposes of this tutorial we will follow what they state in the paper on which this tutorial is based:

"The unique histidine (His6) was doubly protonated. This residue is solvent-exposed, 10 � away from Asp26, so that its protonation state is not expected to affect Asp26 strongly."

Thus we will doubly protonate HIS6 by changing it's name from HIS to HIP.

Here is the modified pdb file: thioredoxin_chaina_modified.pdb

Next we need to create our prmtop and inpcrd files for the two protonation states of Asp26. We will be using IGB=5 for this calculation so we must remember to set the default PBradii to the mbondi2 set. Here is a leap file for creating the prmtop's:

thio_leap.in

source oldff/leaprc.ff94

set default PBradii mbondi2

thio = loadpdb thioredoxin_chaina_modified.pdb

#Add the disulphide bond
bond thio.32.SG thio.35.SG

saveamberparm thio thio_ash.prmtop thio_ash.inpcrd

#Adjust the charges to match the Asp- charges from the paper
set thio.26.OD2 charge -0.8014
set thio.26.OD1 charge -0.8014
set thio.26.CG  charge  0.5307
set thio.26.HD2 charge  0.0000

saveamberparm thio thio_asp.prmtop thio_asp.inpcrd

quit

$AMBERHOME/exe/tleap -s -f thio_leap.in

You should check the output for any errors. Particularly check to make sure that the ASP version has a charge of -1 electron 'more' than ASH.

Here are the files: thio_ash.prmtop, thio_ash.inpcrd, thio_asp.prmtop, thio_asp.inpcrd

Now we are ready to proceed to the next step where we will minimise and equilibrate our structures and then run the thermodynamic integration.

CLICK HERE TO GO TO SECTION 4