Constant pH Molecular Dynamics Example - Calculating pKas for titratable side chains in HEWL

(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 Jason Swails & T. Dwight McGee Jr 2013

Section 1: Building the Input structure and input files

This section describes the steps required to prepare a PDB file for use with CpHMD simulations as well as the additional input files required to run simulations at constant pH. We will be using the 4LYT HEWL structure from the PDB.

Download the original version of the 4LYT PDB file (you can download it here: 4LYT.pdb, or from the PDB itself).

Modify the residue names of the residues you wish to titrate, making the following substitutions (the reference pK_as shown are for the capped, neutral amino acid in solution):

Amino acid Titratable Residue Name Reference pK_a

Aspartate AS4 4.0

Glutamate GL4 4.4

Histidine HIP 6.6

Cysteine CYS 8.5

Lysine LYS 10.4

Tyrosine TYR 9.6

In this example, we are only titrating the acidic residues, so the only changes we are making in the PDB file are the ASP->AS4, GLU->GL4, and HIS->HIP substitutions. Furthermore, HEWL contains 4 disulfide bonds (8 Cysteine residues forming 4 Cystines linked via disulfides). Therefore, we must change all CYS residues into CYX so that LEaP chooses the correct residue template for Cysteines linked via disulfide bonds. Finally, 4LYT contains two conformations, of which we will only keep the first structure. The resulting PDB file should look like this with all waters and extra data removed: 4LYT_Fixed.pdb

Prepare the topology file (prmtop) and input coordinate (inpcrd) file using tleap. A specific leaprc file, leaprc.constph, has been prepared to set up all variables necessary for your simulations. The underlying force field is ff10 (equivalent to ff99SB for proteins), and should not be changed without verifying that the model compounds titrate to the correct pK_a. The leaprc.constph file also sets the appropriate PBRadii set (mbondi2).

The following leap script will generate the necessary topology and coordinate files:

tleap.in

source leaprc.constph # Load the PDBs 4lyt = loadPDB 4LYT_Fixed.pdb # Make the disulfide bonds bond 4lyt.6.SG 4lyt.127.SG bond 4lyt.30.SG 4lyt.115.SG bond 4lyt.64.SG 4lyt.80.SG bond 4lyt.76.SG 4lyt.94.SG # Save topology files saveAmberParm 4lyt 4LYT.parm7 4LYT.rst7 # Quit quit

The command

tleap -f tleap.in

will generate the following files: 4LYT.parm7 and 4LYT.rst7. (Ignore the warning that the net charge is not neutral---this is not necessary for implicit solvent calculations).

Now we will generate the constant pH input file (cpin file) using the cpinutil.py program. This step is important, as it will tell sander which residues it should titrate during the course of the simulation. You can use the --help flag to obtain a list of available options. (Run cpinutil.py --help in your terminal now to see).

Keep in mind here that we only want to titrate glutamate, aspartate, and histidine residues. The flags that are of general interest are listed below:

Flag Meaning

-p Input topology file to search for titrating residues

-resnames Only titratable residues in the given list of names will be titrated.

-notresnames All titratable residues except those in the given list of names will be titrated

-resnums All titratable residues in the list of numbers will be titrated.

-notresnums All titratable residues except those in the list of numbers will be titrated.

-states Defines the initial protonation states of all residues. See --describe below for info on what the state numbers mean. You must supply trescnt values here.

-l,--list List all available titratable residues. Unless restrictive flags (above) are used, all residues matching these names will be marked for titration.

--describe List the charge vectors and reference energies for each state of the described residue(s). This can help you decide what starting states to choose.

By default, all protonation states start at state 0 (AS4, GL4 deprotonated, all others are protonated). We will generate our cpin file with the default starting states, but only for residues AS4, GL4, and HIP.

cpinutil.py -resnames HIP GL4 AS4 -p 4LYT.parm7 -o 4LYT.cpin

generates the following cpin file: 4LYT.cpin. Looking at the cpin file, we can see that 10 residues will be titrated (TRESCNT=10), whose names and numbers are given in the RESNAME section above.

At this point, we have all the files we need to run a constant pH simulation. It is always a good idea to visualize your system at this point, making sure that everything appears correct (especially the disulfide bonds!).

Click here to go to section 2

Click here to go back to the Introduction