pKa Calculations

This example was prepared in collaboration with David Sept for a biomolecular modeling class.


Although pKa calculations may not seem like the most routine applications for demonstrating continuum electrostatics concepts, they have important scientific and teaching value. From a scientific standpoint, pKa values are important determinants of biomolecular (particularly enzymatic) function and can be used to assess functional activity and identify active sites. From a pedagogical perspective, pKa calculations require all of the important continuum electrostatics concepts and therefore serve as a “self-contained” introduction to solvation and binding energies.

This tutorial covers Poisson-Boltzmann methods for determining biomolecule pKa values. Other methods for pKa and titration state determination are given in the PDB2PQR examples. Comparison with PDB2PQR results can provide hours of additional entertainment.

This is a very brief introduction to the concepts behind biomolecular pKas and titration states. More information can be obtained from most biochemistry or biophysics textbooks as well as some of the original articles on pKa evaluation:

Recall that the acid dissociation constant Ka describes the dissocation of an acid into its components

\[ HA \overset{K_{\textrm{a}}}\rightleftharpoons H^+ + A^- \]

where the equilibrium constant is defined related to the activities of the species:

\[ K_{\textrm{a}} = \frac{a_{\textrm{H}} + a_{\textrm{A-}}}{a_{\textrm{HA}}} \]

which, under conditions of “ideality” (never realized in a biological system…), can be replaced with concentrations to give

\[ K_{\textrm{a}}\approx\frac{C_{\textrm{H}} + C_{\textrm{A-}}}{C_{\textrm{HA}}} \]

You should also recall that chemical equilibrium constants can be related to free energies by

\[ -RTlnK_{\textrm{a}} = \Delta_{\textrm{a}} = G_{\textrm{HA}} - G_{\textrm{H}} + - G_{\textrm{A-}} \]

However, chemists found it easier to use quantities in base-10, therefore this energetic quantity is often referred to as a pKa

\[ \mathrm p K_a = -\log_{10}{K_a} = -\frac{K_a}{\ln(10)} = \frac{\Delta_a G}{RT \ln(10)} \approx \frac{\Delta_a G}{2.303RT} \]


Amino acid model pKa values

In many calculations, pKa values are assigned based on model values for amino acid side chains to mimic the reaction for an “isolated” amino acid in solution. Some of these model pKa values are given in the following table:

Amino acid pKa
Arginine 13.0
Aspartic acid 4.0
Cysteine 8.7
C-terminus 3.8
Glutamic acid 4.4
Histidine 6.3
Lysine 10.4
N-terminus 8.0
Tyrosine 9.6

This data is from Nielsen JE, Vriend G. Optimizing the hydrogen-bond network in Poisson-Boltzmann equation-based pKa calculations. Proteins. 43 (4), 403-12, 2001. As we’ll see in the next section, these model values provide the basis for calculating pKa values in proteins.


pKa values in proteins

The role of the model pKa values introduced in the previous sections is to move all the chemical (bond making and breaking) complexity of protonation into the model values. In particular, pKa values in proteins are calculated as pertubations of the model compounds according to the free energy cycles shown below. The pKa for the amino acid in the context of the protein is given by the free energy cycle:

\[ HA(aq) \overset{^{\Delta_a G_{\mathrm{HAmodel}}}}{\rightarrow} H^+(aq) + A^-(protein) \]

\[ HA(protein) \overset{^{\Delta_a G_{\mathrm{HA}}}}{\rightarrow} H^+(aq) + A^-(aq) \]

We are interested in determining the unknown $\Delta_{a}G_{\mathrm{HA}}$ from the known model $\Delta_{a}G_{\mathrm{HA,model}}$, and the unknown $\Delta_{\mathrm{xfer}}G_{\mathrm{HA}}$ and $\Delta_{\mathrm{xfer}}G_{\mathrm{A}^-}$ according to:

\[ \Delta_{\textrm{a}}G_{\textrm{HA,model}} + \Delta_{\textrm{xfer}}G_{\textrm{A\overset{-}}} - \Delta_{\textrm{xfer}}G_{\textrm{HA}} - \Delta_{\textrm{a}}G_{\textrm{HA}} = 0 \]


\[ \Delta_{\textrm{a}}G_{\textrm{HA}} = \Delta_{\textrm{a}}G_{\textrm{HA,model}} - \Delta_{\textrm{A\overset{-}}} - \Delta_{\textrm{xfer}}G_{\textrm{HA}} = 0 \]

In general, the quantities $\Delta_{\mathrm{xfer}}G_{\mathrm{HA}}$ and $\Delta_{\mathrm{xfer}}G_{\mathrm{A}^-}$ are obtained from a computational approach. Nearly any free energy calculation method could be used to obtain these energies according to a scheme where the (de)solvation energies of the charged and uncharged amino acids are calculated according to:


Continuum electrostatics methods for pKa calculations in proteins

Although nearly any free energy method could be used to evaluate the energies of transferring the protonated and unprotonated amino acids from solution into the protein environment, continuum electrostatics offer a (usually) satisfying compromise between accuracy and computational efficiency.

The transfer free energies to be calculated, $\Delta_{\mathrm{xfer}}G_{\mathrm{HA}}$ and $\Delta_{\mathrm{xfer}}G_{\mathrm{A}^-}$, can be determined from Poisson-Boltzmann (PB) energies. In particular, these energies can be calculated as effective “binding energy calculations” similar to those covered in the Binding energies section:

\[ \Delta_{\textrm{xfer}}G_{\textrm{X}} = G_{\textrm{protein with charged X}} \]

\[ -G_{\textrm{protein with uncharged X}} - G_{\textrm{charged X in solution}} \]


  • $G_{\mathrm{protein~with~charged~}X}$ is the electrostatic energy of the protein with group X bound and all charges on X set to their normal values
  • $G_{\mathrm{protein~with~uncharged~}X}$ is the electrostatic energy of the protein with group X bound and all charges on X set to zero
  • $G_{\mathrm{charged~}X\mathrm{~in~solution}}$ is the electrostatic energy of group X in solution with all charges set to their normal values.
  • (I think this is failed latex code?)

Note that, as with binding energies, $\Delta_{\mathrm{xfer}}G_{X}$ can be evaluated two ways:

  • Directly in a PB equation by computing the difference of total electrostatic energies (including self energies) from PB calculations for each of the states. In order for this to work, all conformations/grid positions/charge states must be the same in each PB calculation.
  • Indirectly through PB calculations of solvation free energies and Coulomb’s law calculations (in a dielectric εp) of electrostatic interaction energies. For a sufficiently fine grid, this method is much more stable.

Both of these methods can be combined, via free energy cycles, to give the desired . However, when care is taken to use the same grids and conformations for all calculations, the direct method using total electrostatic energies is usually the most efficient.

Note that none of the methods discussed above have explicitly allowed for changes in titration state of other groups in the protein during protonation/deprotonation of the acid group of interest. Additionally, none of these methods explicitly provide for conformational changes in the protein coupled to protonation/deprotonation. As such, we are not true pKa values with this method. Instead, we are calculating so-called “intrinsic pKa” values.



We currently provide a lysozyme intrinsic pKa example to illustrate these principles.