Lysozyme pKa example for PDB2PQR
Hen egg white lysozyme (HEWL) is a very popular system for pKa calculations as it has a number of interesting values for its titratable residues. Early pKa work on this enzyme is presented in Tanford C and Roxby R (Interpretation of protein titration curves. Application to lysozyme. Biochemistry. 11 (11), 2192-8, 1972e) which also contains the pKa values used in this example. More recent pKa calculations and a review of some of the methodology can be found in Nielsen JE and Vriend G (Optimizing the hydrogen-bond network in Poisson-Boltzmann equation-based pKa calculations. Proteins. 43 (4), 403-12, 2001). Finally, the biological relevance of lysozyme is briefly reviewed at Wikipedia.
HEWL has two active site residues GLU 35 and ASP 52 whose titration states determine the catalytic competency of the enzyme:
In particular, the enzyme is only active when ASP 52 is ionized (pKa ≈ 1.2) and GLU 35 is neutral (pKa = 6.3). Additionally, lysozyme has ASP 66 with a pKa of 2.0 and HIS 15 with a pKa of 5.7.
Overview and disclaimer
In what follows, we’ll calculate the intrinsic pKa of ASP 66. This is the pKa calculated without regard for the titration state changes in other lysozyme residues. Calculating actual pKa changes requires analysis of the coupled titration state energetics of the entire system. This is an extremely labor-intensive process which is best left to computer programs such as WHATIF, PDB2PQR, MCCE, and related software.
Preparing the PDB file
Warning: Problems with explicit water! It is very important that you remove all explicit water from the PDB file before proceeding. (Why?)
We’re going to need to generate a protonated form of ASP 66 to perform our pKa calculations. We will do this with the PDB2PQR web server. Unless otherwise directed, PDB2PQR adds hydrogens to residues based on model pKa values. Therefore, we will need to specify the titration state of ASP 66 for our pKa calculation by changing the residue name from ASP to ASH using your favorite text editor. Save the result as 2LZT-ASH66.pdb.
We’re now ready to run PDB2PQR to generate protonated versions of our PDB files. Use the command line version of PDB2PQR or one of the web servers listed on the PDB2PQR homepage to generate protonoated PQR files for for 2LZT-GLU35.pdb and 2LZT-GLH35.pdb. Name your results 2LZT-GLU35.pqr and 2LZT-GLH35.pqr, respectively. Although it is always important to test sensitivity to various force fields, I’d recommend starting with PARSE.
Note: You can use PDB2PQR to assign titration states with PROPKA but don’t do it for the above steps since we need to set the titration states explicitly for our calculations.
Recall that we’re also going to need the isolated residue for our electrostatics calculations of intrinsic pKa. Use your favorite text editor to extract the entire ASP 66 and ASH 66 residue from 2LZT-ASP66.pqr and 2LZT-ASH66.pqr, respectively. Save the results as separate PQR files containing only the residue of interest: ASP66.pqr and ASH66.pqr, respectively.
Finally, we’ll also need to perform electrostatics calculations on HEWL with and uncharged residue 66. Use your favorite text editor to zero out the charges in 2LZT-ASP66.pqr and 2LZT-ASH66.pqr to create 2LZT-noASP66.pqr and 2LZT-noASH66.pqr. This can be done by setting the second-to-last column in the PQR file to zero; e.g.
in the lysozyme PQR file would become
Setting up the total electrostatic energy calculations
We will be using focusing calculations to calculate the electrostatic potential and free energies for the systems of interest.
Warning: In what follows, we are evaluating total electrostatic free energies – e.g., energies which contain charge self-interaction terms. We will cancel these self-interaction terms in subsequent steps when we calculate solvation or transfer free energies. Therefore, it is very important that you use the same grid parameters (grid centers, dimensions, spacings, etc.) for every calculation.
Here is a template input that we will use for each of the solvation energy calculations:
There are a number of aspects to this input file which are worth noting:
- In general, compound.pqr will change for each calculation but ref.pqr will not. Choose one molecule to be ref.pqr (2LZT-ASP66.pqr is a good choice) and use it in every calculation.
- We are using a solute dielectric constant (εp) of 20 (see pdie). This is a common choice for pKa since
it is thought to implicitly represent internal relaxation and rearrangement of the solute </end vigorous waving of hands>.
- We are using a molecular surface (srfm smol) with a reasonably high density of surface discretization points (sdens 40.0]). pKa and other electrostatics results can be very sensitive to surface choice.
You now have enough information to calculate total electrostatic energies for all of the relevant molecules so far: 2LZT-ASP66.pqr, 2LZT-ASH66.pqr, 2LZT-noASP66.pqr, 2LZT-noASH66.pqr, ASP66.pqr, and ASH66.pqr. You should be able to construct APBS input files for each of these systems by modifying the template above. Once these input files are constructed, you can run the PB calculation by$ apbs foo.in | tee foo.out
where foo.in is the input file of interest and the output is saved in foo.out.
Setting up the transfer free energy calculations
Recall the transfer free energies can be evaluated by direct subtraction of total electrostatic energies for the different dielectric environments and components. This is usually the most stable route, assuming identical grids and solute conformations are used for all calculations. You should always check for a lack of convergence in the calculations and can be resolved by decreasing the grid spacing (e.g., increasing the number of grid points).
Putting it all together
At this point, you should have everything you need to calculate the intrinsic pKa of interest. However, if you get stuck, I’ve attached some example files below that might be helpful:
- 2LZT-ASP66.pqr and 2LZT-ASP66.in PQR and input files for ASP 66 in the protein
- ASP66.pqr and ASP66.in PQR and input files for isolated ASP 66 in solution
- 2LZT-noASP66.pqr and 2LZT-noASP66.in PQR and input files for the protein with an uncharged ASP 66
- 2LZT-ASH66.pqr and 2LZT-ASH66.in PQR and input files for ASH 66 in the protein
- ASH66.pqr and ASH66.in PQR and input files for isolated ASH 66 in solution
- 2LZT-noASH66.pqr and 2LZT-noASH66.in PQR and input files for the protein with an uncharged ASH 66
- un-apbs.sh A Bash shell script for running the various APBS calculations
- process.sh A Bash shell script for transforming the output from run-apbs.sh into a pKa shift
Based on this brief introduction, you should be in a good position to go back and try to evaluate intrinsic pKas for HIS 15 and GLU 35. How well do your results for those residues agree with experiment? What’s different with those residues?
So far, we’ve only examined intrinsic pKas in this example and have ignored coupling between titratable groups. Jens Nielsen has developed a very nice software package called pKaTool which allows you to explore coupling between titratable sites and its impact on titration events in a protein system. He has provided a tutorial (PDF) which you can use to explore coupled titration states and use to familiarize yourself with pKaTool.