(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 and T. Dwight McGee, 2013

---

Analyzing the Results:

This section describes how to analyze the results obtained from the CpHMD simulations.

Using cphstats

The various protonation states that are sampled throughout the course of the CpHMD simulations are written to a cpout file. The program cphstats can be used to parse the cpout file and for calculating the predicted pK_a values. Use the -h flag to get a full listing of the available command-line flags for cphstats.

cphstats -i cpin [cpout1 cpout2 cpout3 cpoutN] [-o output] [--population population_ouput]

For example, the command

cphstats -i 4LYT.equil.cpin 0/4LYT.md1.cpout -o pH0_calcpka.dat --population pH0_populations.dat

will create the following two files for the simulations at pH 0: pH0_calcpka.dat and pH0_populations.dat

You can use a small shell command to analyze all results at once. For instance:

for ph in 0 1 2 3 4 5 6 7; do
   cphstats -i 4LYT.equil.cpin ${ph}/4LYT.md1.cpout -o pH${ph}_calcpka.dat --population pH${ph}_populations.dat
done

will create files with the titration statistics for all simulations. We recommend you generate these files yourself, but a tarball containing these files are available here: pH_statistics.tar.bz2

 

Offset: is the difference between the predicted pK_a and the system pH
Pred: is the predicted pK_a
Frac Prot: is the fraction of time the residue spends protonated
Transitions: gives the number of accpeted protonations state transitions
Average total molecular protonation: is the sum of the fractional protonations

cphstats also printed another file with the populations of every state for every titratable residue, pH0_populations.dat. This file prints the fraction of snapshots that the system spent in each state for each residue. For example, His 15 was protonated 100% of the time, whereas Asp 119 was deprotonated 1.42% of the time, syn- protonated ~48% of the time on each oxygen (this is good since they are rotationally degenerate!), and anti-protonated ~0.7% of the time.

NOTE: The parentheses (*), indicates the number protons present in that state.

The cpout file can be divided into chunks of time using calcpka. This feature is enabled when a dump_interval, is specified (using the -t flag on the command-line). The information is printed to the file pKa_evolution.dat if the -ao flag is not given.

Calculating the pK_as

Here we will look at a couple ways of computing pK_a values using the output from calcpka. The predicted pK_as printed in the calcpka output assumes ideal, Hendersen-Hasselbalch titration behavior obeying the equation:

In proteins, titration behavior for each titratable residue is not always this straight-forward, and two or more residues may impact each others' titration curves. Therefore, the Hill equation, shown below, is typically used to fit titration curves and calculate pK_as.

The trick now becomes to plot our titration curve—typically plotted as the deprotonated fraction (which is directly related to [A^-]) versus pH. Knowing that f_d = \frac{[A^-]} {[HA] + [A^-]}, where f_d is the fraction deprotonated, we can rearrange Equation 2 to obtain:

Now the trick is to extract the deprotonated fraction from the calcpka output shown above (which is just 1 − fraction protonated) for each residue as a function of the pH. The data is shown below, where the first column is the pH and the subsequent columns are the deprotonated fraction for each residue.

In this case, I use gnuplot to fit all of this data to Equation 3 using gnuplot's built-in fitting functionality. This function needs to be fit using a non-linear, least-squares fitting algorithm (e.g., the Marquardt-Levenberg algorithm). As a result, you should choose a reasonable starting guess for the fit (the predicted pK_a printed out by calcpka for the pH with the lowest offset value is a good starting point). The gnuplot script shown below will fit all of the titration curves, plot them, and print the calculated pK_a and Hill coefficients. I wrote a gnuplot tutorial (available here) if you wish to better understand what each line in the script is doing.

The calculated pK_as here show very good qualitative agreement with experimental measurements for HEWL pK_as. Glu 35, whose pK_a has been measured near 6.1, is correctly predicted to have a large shift away from the model compound pK_a of 4.4 for an isolated Glutamate. In each case, this method correctly predicts the direction of the pK_a shift from the model compound, and even gives quantitatively reasonable results for most residues as well. Keep in mind, though, that we ran only 1 ns simulations at each pH, which is most likely insufficient for extensive conformational sampling. It remains, however, illustrative in demonstrating how to carry out CpHMD simulations and calculate resulting pK_as.

Comparing the Root Mean Squared Deviations of the CpHMD Simulations

Next, we calculate the rmsd of each the CpHMD simulations and histogram the data using cpptraj.

The input shown above is using the trajectory for pH 0 and can be modified to work with any of the other pH values. The coordinates from 4LYT.equil.rst7 were used as the reference structure and the rmsd of only the backbone atoms were considered. The hist function in cpptraj was used to histogram the data. The bin width is set to .1, and the wildcards (*) tell cpptraj to use the smallest and largest values of the data set to set the smallest and highest bin, respectively. The last wild-card tells cpptraj to use the number of bins necessary to go from the data set minimum to the data set maximum using a step size of 0.1 Angstroms. See the AmberTools manual for more detailed instructions on using cpptraj.

The data can be generated by executing the following command:

cpptraj 4LYT.parm7 cpptraj.rmsd.in

A tarball containing all of the rmsd and histogram files can be downloaded here: rmsd_data.tar.bz2

Because our simulations were relatively short and we only have 1000 data points, the histograms appear pretty noisy. However, we can qualitatively say that the high pH and low pH simulations were more flexible and explored conformations farther from the crystal structure. This makes sense, since HEWL is most active near pH 4, and 4LYT was solved at pH 4.8. [2]

Visualizing the trajectories in VMD

If you do not already know how to use VMD with Amber files, please see Basic Tutorial 2 for instructions

The trajectories of the CpHMD simulations can be viewed using VMD. The program VMD_pH.py (along with getatms.py) get_creates a tcl script to load into VMD with your topology and trajectory file. This will script will allow you to see the different protonation states sampled by the protein during the CpHMD simulations. VMD_pH.py is inefficient for large trajectories, and was written quickly with limited knowledge of VMD's scripting environment.

VMD_pH.py should be used as:

VMD_pH.py cpin cpout1 cpout2 cpout3 cpoutN...

where the cpin file is just parsed for which residues are titrating (so it should typically be the cpin file written by cpinutil.py). The cpout files must be listed in the same order as the trajectories are loaded into VMD (the full records in the cpout file are matched to each frame of the trajectory). For large numbers of snapshots, this script becomes very inefficient.

Perfoming the command below will create a file named pHVMDscript.tcl that can be loaded into VMD.

VMD_pH.py 4LYT.equil.cpin 0/4LYT.md1.cpout

This particular script identifies the titratable protons that are "gone" for each residue in each frame, then moves those hydrogen atoms 99 Angstroms away in each direction. Therefore, we need to use the Dynamic Bonds representation in VMD to avoid drawing all of the bonds to the H-atoms we have moved away from the protein.

Loading the files in VMD

1. Open VMD
2. In the console titled VMD Main, click File, followed by New Molecule
3. Load the topology file, 4LYT.parm7, (AMBER 7 Parm) format, followed by the trajectory file, 0/4LYT.md1.nc, (NetCDF AMBER,MMTK) format
4. In the vmd terminal type source pHVMDscript.tcl
5. In the console titled VMD Main click Graphics, followed by Representations
6. Change the drawing method in Graphical Representations console to New Cartoon
7. Click the tab Create Rep and under Selected Atoms type resid 7 15 18 35 48 52 66 87 101 119
8. Change the drawing method to Dynamic Bonds
9. Click the tab Create Rep and change the drawing method to VDW and reduce the Sphere Scale to 0.4

---

Click here to go back to the Introduction

---

Click here to go back to Section 3

---

References

[2] A. C. M. Young, J. C. Dewan, C. Nave, and R. F. Tilton, "Comparison of radiation-induced decay and structure refinement from X-ray data collected from lysozyme crystals at low and ambient temperatures", J. Appl. Cryst., 1993, 26, pp. 309-319

---

(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 and T. Dwight McGee, 2013