Copyright © 2002, 2003, 2004, 2005, 2006, 2007 Washington University in St. Louis
This file is part of APBS.
This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
Linking APBS statically or dynamically with other modules is making a combined work based on APBS. Thus, the terms and conditions of the GNU General Public License cover the whole combination.
SPECIAL GPL EXCEPTION: In addition, as a special exception, the copyright holders of APBS give you permission to combine the APBS program with free software programs and libraries that are released under the GNU LGPL or with code included in releases of ISIM, Ion Simulator Interface, PMV, PyMOL SMOL, VMD, and Vision. Such combined software may be linked with APBS and redistributed together in original or modified form as mere aggregation without requirement that the entire work be under the scope of the GNU General Public License. This special exception permission is also extended to any software listed in the SPECIAL GPL EXCEPTION clauses by the PMG, FEtk, MC, or MALOC libraries.
Note that people who make modified versions of APBS are not obligated to grant this special exception for their modified versions; it is their choice whether to do so. The GNU General Public License gives permission to release a modified version without this exception; this exception also makes it possible to release a modified version which carries forward this exception.
Table of Contents
List of Figures
List of Tables
List of Examples
List of Equations
This tutorial is designed as a "how to" document to get users acquainted with computational electrostatics calculations using APBS. In order to perform the provided examples, you will need the latest version of APBS (http://apbs.sf.net). Other requirements are listed in the individual sections.
![]() | Important |
---|---|
Note that many of the examples in this tutorial can also be run through the web-based Gemstone service without the need for local installations of the APBS software. More information is available in Chapter 10, Can I run APBS via the web? (Gemstone) . |
![]() | Warning |
---|---|
This document is still under construction and should be completed before the next APBS release. Many of the topics covered in the incomplete sections are demonstrated in the APBS |
In order to perform electrostatics calculations on your biomolecular structure of interest, you need to provide atomic charge and radius information to APBS. Charges are used to form the biomolecular charge distribution for the Poisson-Boltzmann (PB) equation while the radii are used to construct the dielectric and ionic accessibility functions. This charge and radius information can be provided to APBS in a few different formats which are described in more detail below.
The PQR format provides a very simple way to include parameter information by by replacing the occupancy and temperature columns of a PDB-format structure file with charge ("Q") and radius ("R") information. Unfortunately, the simplicity of this format also limits its extensibility: it can be very difficult to add new atom types and parameters in a generic format without the use of external software such as PDB2PQR. The XML parameter format described below is much easier to modify.
The XML format provides a somewhat more complicated format for including parameter information with increased flexibility in formatting, extension, and other modifications. As in the PQR format, atom coordinates are supplemented with charge and radius information. Please see the APBS user guide for complete format specifications.
![]() | Note |
---|---|
This tutorial was prepared using PDB2PQR 1.2.1 |
The PDB2PQR web service and software will convert most PDB files into PQR format with some caveats, in particular the inability to "repair" large regions of missing coordinates. Although PDB2PQR can fix some missing heavy atoms in sidechains, it does not currently have the (nontrivial) capability to model in large regions of missing backbone and sidechain coordinates.
PDB2PQR will also perform hydrogen bond optimization, sidechain rotamer search, limited titration state assignment, ligand parameterization, and APBS input file preparation. Please see http://pdb2pqr.sf.net/ for more details.
As mentioned above, PDB2PQR is discussed in more detailed on the PDB2PQR homepage . Therefore, we will review the minimal steps required to produce a PQR file from a PDB file here. To begin, choose a server from http://pdb2pqr.sf.net/#availability.
Either enter the 4-character PDB ID or accession number (e.g., 1FAS, 1MAH, 1LYS, etc.) or upload your own PDB file.
![]() | Note |
---|---|
Note that, if you choose to enter a 4-character PDB file, PDB2PQR will process all recognizable chains of PDB file as it was deposited in the PDB (e.g., not the biological unit, any related transformations, etc.). |
For most applications, the choice is easy: PARSE. This forcefield has been optimized for implicit solvent calculation and is probably the best choice for visualization of protein electrostatics and many common types of energetic calculations for proteins. However, AMBER and CHARMM may be more appropriate if you are attempting to compare directly to simulations performed with those force fields, require nucleic acid support, are simulating ligands parameterized with those force fields, etc.
It is also possible to upload a user-defined forcefield (e.g., to define radii and charges for ligands or unusual residues). Please see the PDB2PQR documentation for more information.
This is largely irrelevant to electrostatics calculations but may be important for visualization. When in doubt, choose the "Internal naming scheme" which attempts to conform to IUPAC standards.
These options fall into two categories: how to build missing atoms (including hydrogens) onto the structure and additional output configuration. Please see http://pdb2pqr.sf.net/ for more details.
Table of Contents
There are a number of excellent molecular graphics programs which can be used to visualize electrostatic properties in the context of the biomolecule. Since there are far too many to review in any depth, we will focus on the 3 packages with which we are most familiar.
The VMD molecular graphics software package provides support for both the execution of APBS and the visualization of the resulting electrostatic potentials. Documentation on the APBS interface has been provided by the VMD developers at http://www.ks.uiuc.edu/Research/vmd/plugins/apbsrun/. We will supplement this with a basic demonstration of how go from a PDB entry to a plot of structure & potential in VMD using APBS.
![]() | Note |
---|---|
This tutorial was written using VMD 1.8.5. |
We'll perform this example with arc repressor (PDB ID 1MYK), a DNA binding protein. The first step in the visualization process is generating a PQR file for use with VMD and APBS. Please follow the directions in Section 2.3, “Generating PQR files from PDB files (PDB2PQR)” to generate a PQR file for this structure.
Load the PQR file you just created into VMD (1] Depending on the force field[2] you used in PDB2PQR, you may notice that VMD is displaying some strange bonds. When a PQR file is loaded into VMD, bond distances are inferred from the PQR radii. These radii were chosen for continuum electrostatics calculations, not visualization. Don't worry if some of the bonds appear a bit strange (hide the hydrogens if it really bothers you).
→ in the VMD Main window). Adjust the molecule to the desired view. [You're now ready for electrostatics calculations.
Choose
→ → in the VMD Main window. A new APBS Tool window should appear.Choose
→ in the APBS Tool window and set the path to your local copy of the APBS executable. Hit to close this window.Select the "0" calculation from the "Individual PB calculations (ELEC):" window. Hit the button. Adjust the APBS settings here as you desire. The defaults are usually fine, although it may be useful to edit the ions to adjust the ionic strength used in the calculation. Hit when you're done.
You're now ready to run the calculation. Hit
in the APBS Tool window. Check the VMD Console window for information about the job while it's running. When the job is finished, a "APBSRun: Load APBS Maps" window should appear. Select " " and hit .
One of the most popular visualization methods is the isocontour.
First, choose
→ from the VMD Main window.Hit the Drawing Method to "Isosurface".
button and changeChange Draw from "Points" to "Solid Surface" and Material to "Transparent".
Note that the current isovalue is 0; change this to 1 (or 5 or 10 or whatever) depending on your personal preferences.
To continue the lonstanding tradition of electrostatic potential coloring, choose "ColorID 0" for the Coloring Method.
For the negative isocontour, hit "
" and select the newly created representation. Change the isocontour value to -1 (or 5 or 10...) and the ColorID to 1.At this point, you probably have an image that looks something like this (note that I changed the surface to points to make the figure a bit cleaner):
Another popular method of electrostatic potential visualization maps the electrostatic potential to the biomolecular surface. Before proceeding, you may want to delete the two isocontours you just created using the
in the Graphical Representations window.Go to the Graphical Representations window (
→ from the VMD Main window) and create a new representation of the molecule with the button.Go to the "Draw style" tab of the Graphical Representations window and change Drawing Method to "Surf"[3] and Coloring Method to "Volume".
Go to the "Trajectory tab" of the Graphical Representations window and change the "Color Scale Data Range:" to -10 to 10 (or whatever your favorite values are).
Based on your version of VMD and your personal preferences, you might want to change the color scale for this image. Go to Method menu).
→ in the VMD Main window and select the "Color Scale" tab from the new "Color Controls" window. The traditional coloring scheme for electrostatics is "RWB" (in theYour molecule now probably looks like this:
The PyMOL molecular graphics software package provides support for both the execution of APBS and the visualization of the resulting electrostatic potentials. We will provide a basic demonstration of how go from a PDB entry to a plot of structure & potential in PyMOL using APBS.
![]() | Note |
---|---|
This tutorial was written using PyMOLX11Hybrid 0.99 on a Mac. |
We'll perform this example with fasciculin-2 (PDB ID 1FAS), a snake neurotoxin which binds the negatively-charged acetylcholinesterase. Please generate the PQR file using the steps outlined in Section 2.3, “Generating PQR files from PDB files (PDB2PQR)”.
Load the PQR file you created into PyMOL (
→ ) and choose your favorite graphical representation.Go to the
→ to open the APBS calculation plugin.Under the "Main" tab of the PyMOL APBS Tools window, select and either browse to (via the button) or input the path to your PQR file. This step is necessary to ensure you use the radii and charges assigned by PDB2PQR.
Under the "APBS Location" tab of the PyMOL APBS Tools window, either browse to (via the button) or input the path to your local APBS binary. It is not necessary to provide a path to the APBS psize.py binary for most biomolecules.
Under the "Temporary File Locations" tab of the PyMOL APBS Tools window, customize the locations of the various temporary files created during the run. This can be useful if you want to save the generated files for later use.
Under the "Configuration" tab of the PyMOL APBS Tools window, hit the to set the grid spacings. The default values are usually sufficient for all but the most highly charged biomolecules.
Under the "Configuration" tab of the PyMOL APBS Tools window, customize the remaining parameters; the defaults are usually OK[4].
Under the "Configuration" tab of the PyMOL APBS Tools window, hit the button to start the APBS calculation. Depending on the speed of your computer, this could take a few minutes. The button will become unselected when the calcaultion is finished.
Before proceeding with the remaining steps, you must load the electrostatic potential data into PyMOL. Under the "Visualization" tab of the PyMOL APBS Tools window, hit the button.
PyMOL makes this step very easy: adjust the positive and negative "Contour" fields to the desired values (usually +/-1, +/-5, or +/-10 kT/e[5]) and hit the Positive Isosurface and Negative Isosurface and buttons.
At this point, you probably have a figure that looks something like:
If you haven't already, hide the isocontours by hitting Positive Isosurface and Negative Isosurface and buttons.
The surface potential is also straightforward to visualize. Set the "Low" and "High" values to the desired values (usually +/-1, +/-5, or +/-10 kT/e) at which the surface colors are clamped at red (-) or blue (+). Check the " " and " " buttons to plot the potential on the solvent-accessible (probe-inflated or Lee-Richards) surface[6]. Hit the "Molecular Surface" button to load the surface potential.
![]() | Warning |
---|---|
Under construction -- surfaces are currently broken with some PQR structure/force field combinations. |
[1] At this point you may notice (to many users' chagrin) that PDB2PQR removed the chain IDs from the structure. Sorry about this; we're trying to fix it in future releases of APBS and PDB2PQR.
[2] Visualization seems to work best with PARSE
[3] FYI, this appears to generate a molecular or Connolly type of surface rather than the solvent-accessible or Lee-Richards surface I personally prefer for this type of visualization...
[4] Although I usually prefer 0.150 concentrations for the +1 and -1 ion species so electrostatic properties aren't overly exaggerated.
[5] Note that PyMOL has a slight units error on the contour values.
[6] In my opinion, the solvent-accessible surface tends to reveal more global features of the surface potential. Tighter surfaces (e.g., van der Waals and molecular or Connolly surfaces) tend to simply map atomic surface onto the biomolecular surface.
Table of Contents
Solvation energies are usually decomposed into a free energy cycle as shown in Figure 4.1, “Full solvation energy cycle”. Note that such solvation energies often performed on fixed conformations; as such, they are more correctly called "potentials of mean force". More details on using APBS for the polar and nonpolar portions of such a cycle are given in the following sections.
Figure 4.1. Full solvation energy cycle. This cycle incorporates several processes into the solvation energy (step 1). Step 2 indicates charging of the solute in solution (e.g., inhomogeneous dielectric, ions present). Step 3 indicates the introduction of attractive solute-solvent dispersive interaction interactions (e.g., an integral of Weeks-Chandler-Andersen interactions over the solvent-accessible volume). Step 4 indicates the introduction of repulsive solute-solvent interaction (e.g., cavity formation). Steps 5 and 6 are basically null steps although they could be used to offset unwanted energies added in Steps 3 and 4 above. Finally, Step 6 represents the charging of the solute in a vacuum or homogeneous dielectric environment in the absence of mobile ions.
The full free energy cycle in Figure 4.1, “Full solvation energy cycle” is usually decomposed into polar and nonpolar parts. The polar portion is usually represented by the charging energies in Steps 2 and 6:
Energies returned from APBS electrostatics calculations are charging free energies. Therefore, to calculate the polar contribution to the solvation free energy, we simply need to setup two calculations corresponding to Steps 2 and 6 in Figure 4.1, “Full solvation energy cycle”. However, the electrostatic charging free energies returned by APBS include self-interaction terms. These are the energies of a charge distribution interacting with itself. Such self-interaction energies are typically very large and extremely sensitive to the problem discretization (grid spacing, location, etc.). Therefore, it is very important that the two calculations are performed with identical grid spacings, lengths, and centers, in order to ensure appropriate matching (or "cancellation") of self-energy terms.
One of the canonical examples for polar solvation is the Born ion: a nonpolarizable sphere with a single charge at its center surrounded by an aqueous medium. In the absence of external ions, the polar solvation energy for this sytem is given by Equation 4.1, “Born ion polar solvation energy”.
Equation 4.1. Born ion polar solvation energy. q is the ion charge, a is the ion radius, and the two epsilons denote the internal and external (solution) dielectric constants. This model assumes zero ionic strength.
We can setup a PQR file for the Born ion for
use with APBS as shown in Example 4.1, “Born ion PQR” and downloaded
here.
Example 4.1. Born ion PQR
REMARK This is an ion with a 3 A radius and a +1 e charge ATOM 1 I ION 1 0.000 0.000 0.000 1.00 3.00
We're interested in performing two APBS
calculations for the charging free energies in
homogeneous and heterogeneous dielectric
coefficients. We'll assume the internal
dielectric coefficient is 1 (e.g., a vacuum)
and the external dielectric coefficient is
78.54 (e.g., water).
for these settings, the polar solvation energy
expression has the form
where R is the ion radius in Angstroms. One possible APBS input file for this calculation is shown in Example 4.2, “Example born input file” and can be downloaded here.
Example 4.2. Example born input file
# READ IN MOLECULES read mol pqr born.pqr end elec name solv # Electrostatics calculation on the solvated state mg-manual # Specify the mode for APBS to run dime 97 97 97 # The grid dimensions nlev 4 # Multigrid level parameter grid 0.33 0.33 0.33 # Grid spacing gcent mol 1 # Center the grid on molecule 1 mol 1 # Perform the calculation on molecule 1 lpbe # Solve the linearized Poisson-Boltzmann equation bcfl mdh # Use all multipole moments when calculating the potential pdie 1.0 # Solute dielectric sdie 78.54 # Solvent dielectric chgm spl2 # Spline-based discretization of the delta functions srfm mol # Molecular surface definition srad 1.4 # Solvent probe radius (for molecular surface) swin 0.3 # Solvent surface spline window (not used here) sdens 10.0 # Sphere density for accessibility object temp 298.15 # Temperature calcenergy total # Calculate energies calcforce no # Do not calculate forces end elec name ref # Calculate potential for reference (vacuum) state mg-manual dime 97 97 97 nlev 4 grid 0.33 0.33 0.33 gcent mol 1 mol 1 lpbe bcfl mdh pdie 1.0 sdie 1.0 chgm spl2 srfm mol srad 1.4 swin 0.3 sdens 10.0 temp 298.15 calcenergy total calcforce no end # Calculate solvation energy print energy solv - ref end quit
Running this example with a recent version of
APBS should give an answer of -229.59 kJ/mol
which is in good agreement with the -230.62
kJ/mol predicted by the analytic formula above.
Note that the Born example above can be easily
generalized to other polar solvation energy
calculations. For example, ions could be added
to the solv
ELEC
section, dielectric
constants could be modified, surface
definitions could be changed (in both
ELEC
sections!), or more
complicated molecules could be examined. Many
of the examples included with APBS (e.g.,
solv
and
ionize
) also demonstrate
solvation energy calculations.
![]() | Note |
---|---|
Note that, as molecules get larger, it is important to examine the sensitivity of the calculated solvation energies with respect to grid spacings and dimensions. |
Referring back to the full free energy cycle in Figure 4.1, “Full solvation energy cycle”, the nonpolar solvation free energy is usually represented by the energy changes in Steps 3 through 5:
where Step 4 represents the energy of creating a cavity in solution and Steps 3-5 is the energy associated with dispersive interactions between the solute and solvent. There are many possible choices for modeling this nonpolar solvation process. APBS implements a relatively general model described by Wagoner and Baker (PNAS 2006) and references therein. The implementation and invocation of this model is described in more detail in the APBS user guide.
Our basic model for the cavity creation term (Step 4) is motivated by scaled particle theory and has the form
where p is the solvent pressure (press
keyword in APBS), Δ V is the solvent-accessible volume of the solute, γ is the solvent surface tension (gamma
keyword in APBS), and Δ A is the solvent-accessible surface area of the solute.
Our basic model for the dispersion terms (Steps 3 and 5) follow a Weeks-Chandler-Anderson framework as proposed by Levy, Zhang, and Gallicchio (J Comput Chem, 2002):
where ρ is the bulk solvent density (bconc
keyword in APBS), Ω is the problem domain, u(att)(y) is the attractive dispersion interaction between the solute and the solvent at point y with dispersive Lennard-Jones parameters specified in APBS parameter files, and θ(y) describes the solvent accessibility of point y.
![]() | Note |
---|---|
The ability to independently adjust |
Unlike the Born ion, there are no good simple examples for demonstrating these types of nonpolar calculations. APBS includes several examples of calculations using the apolar model above. Interested readers should examine the alkanes
apolar example provided with APBS.
Table of Contents
In general, implicit solvent models are used to calculation the contribution of solvation to binding free energies. Additional binding free energy contributions (molecular mechanics energies, entropic changes, etc.) must be calculated separately and are not discussed in this tutorial. One exception is the inclusion of intermolecular Coulombic interations; we'll discuss how these can be calculated in APBS below.
Our framework for calculating solvation contributions to binding free energies is shown in Figure 5.1, “ Binding free energy cycle ”.
Figure 5.1. Binding free energy cycle illustrating binding in terms of transfer free energies from a homogeneous dielectric environment (where interactions are described by Coulomb's law) to an inhomogeneous dielectric environment with differing internal (green) and external (cyan) dielectric constants. The binding (dissociation) free energy is depicted in Step 3.
The binding free energy is given by Equation 5.1, “ Binding free energy formula
”.
The following sections provide more detail on
calculating individual terms of this equation.
If we're just interested in calculating the solvation contributions to binding (steps 4 and 2 in Figure 5.1, “ Binding free energy cycle ”), then we simply need to follow the instructions from Chapter 4, How do I calculate a solvation energy? for the complex and isolated components. The solvation energy contribution to the binding is then given by Equation 5.2, “ Solvation contributions to the binding free energy ”.
Equation 5.2. Solvation contributions to the binding free energy for a two-component complex of molecules "mol 1" and "mol 2". Refer to Figure 5.1, “ Binding free energy cycle ” for numbering.
This solvation energy contribution can be decomposed
into polar and nonpolar parts as outlined in Chapter 4, How do I calculate a solvation
energy?.
To complete the binding free energy cycle (see Figure 5.1, “ Binding free energy cycle ”), we need to add intermolecular Coulombic contributions to the solvation energy change upon binding to get the total electrostatic/solvent contribution to the binding free energy. In particular, we're interested in the change in Coulombic electrostatic energy upon binding, as given by Equation 5.3, “ Coulombic contributions to the binding free energy ”.
Equation 5.3. Coulombic contributions to the binding free energy for a two-component complex of molecules "mol 1" and "mol 2". See Equation 5.1, “ Binding free energy formula ” for numbering.
Each of the ΔGcoul
quantities in Equation 5.3, “ Coulombic contributions to the binding
free energy ” is the
sum of pairwise Coulombic interactions between
all atoms in the molecule (or complex)
for a particular uniform
dielectric. In order to
combine these Coulombic binding energies with
the solvation energies described above,
we need to make sure consistent dielectric
constants are used. In particular, Coulombic
interactions should be calculated using the
same uniform dielectric constant as the
reference state of the solvation energy above.
For example, if solvation energies are
calculated for transfering a protein from a
medium with uniform dielectric of
εin
to a medium with internal dielectric
εin
and external dielectric
εout,
then Coulombic energies should
be calculated using a dielectric of
εin.
The APBS accessory program
tools/manip/coulomb
was
created to help with the calculation of these
individual per-molecule Coulombic energies.
Given a PQR file as input, the
tools/manip/coulomb
program calculates Coulombic energies for a
vacuum dielectric (e.g., a uniform dielectric
of 1). If the reference dielectric is
εin, then all
energies returned by
tools/manip/coulomb
need to be divided by
εin.
With the appropriate dielectric constant of εin used for the Coulombic binding energy, a total electrostatic/solvation free energy can be calculated from Equation 5.4, “ Binding free energy ”. Equation 5.1, “ Binding free energy formula ”, Equation 5.2, “ Solvation contributions to the binding free energy ”, and Equation 5.3, “ Coulombic contributions to the binding free energy ” as
APBS provides force calculations for both polar and nonpolar solvation following the same procedures used in Chapter 4, How do I calculate a solvation
energy?.
In general, forces can be obtained by modifying input files used for solvation energy calculations to include calcforce total
(for total forces on the entire solute) or calcforce comps
to obtain detailed force information for each atom.
It is important to note that, like solvation energy calculations, "self-interaction" terms must be removed (cf. the Example 4.2, “Example born input file” example).
Table of Contents
![]() | Important |
---|---|
This tutorial was prepared with Dave Sept as part of a lab for a biomolecular simulation class. |
![]() | Note |
---|---|
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. |
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[7].
Recall that the acid dissociation constant Ka describes the dissocation of an acid into its components:
by way of the activities
Under “conditions of ideality,”[8] activities can be replaced with concentrations to give
You should also recall that chemical equilibrium constants can be related to free energies by
However, chemists found it easier to use base-10 logarithms than natural logarithms for measuring pH, therefore the pKa is defined as
In many calculations, pKa values are assigned based on model values for amino acid side chains to mimic
for an isolated amino acid in solution. Some model pKa values are given in Table 7.1, “Model amino acid pKa values for common titratable groups; data from Nielsen et al (see footnote) ”.
Table 7.1. Model amino acid pKa values for common titratable groups; data from Nielsen et al (see footnote)
Amino acid | Model 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 |
As we'll see in the next section, these model values provide the basis for calculating 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:
We are interested in determining the unknown ΔaGHA from the known model ΔaGHA,model and the unknown ΔxferGHA and ΔxferGA- according to:
In general, the quantities ΔxferGHA and ΔxferGA- are obainted 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:
Note that these energies have all assumed the same background state
of the protein for the calculations; in other words, the same states
are used for other titratable groups in the protein for charged and
uncharged amino acid states. We'll discuss the implications of this
assumption a bit later in this tutorial.
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, ΔxferGHA and ΔxferGA-, 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 APBS examples and tutorial:
where
Gprotein with charged X is the electrostatic energy of the protein with group X bound and all charges on X set to their normal values.
Gprotein with uncharged X is the electrostatic energy of the protein with group X bound and all charges on X set to zero.
Gcharged X in solution is the electrostatic energy of group X in solution with all charges set to their normal values.
Note that, as with binding energies, Equation 7.2, “Transfer free energy” can be evaluated two ways:
Ways to evaluate transfer free energies
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. This method is much more stable.
Both of these methods can be combined, via free energy cycles, to give the desired ΔxferGX. 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 stable.
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/deprotonoation. As such, we cannot claim to be computing true pKa values with this method. Instead, we are calculating so-called intrinsic pKa values.
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, Roxby R. Interpretation of protein titration curves. Application to lysozyme. Biochemistry. 11 (11), 2192-8, 1972 which also contains the pKa values used in this lab. More recent pKa calculations and a review of some of the methodology can be found in Nielsen JE, Vriend G. Optimizing the hydrogen-bond network in Poisson-Boltzmann equation-based pk(a) calculations. Proteins-Structure Function and Genetics. 43 (4), 403-12, 2001. Finally, the biological relevance of lysozyme is briefly reviewed at the PDB Molecule of the Month site.
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 6.6 and HIS 15 with a pKa of 5.7.
Download the PDB entry 2LZT
from the PDB; save it as 2LZT-GLU35.pdb
. If
you have time, you should also visit the PDBSum
analysis page as well for additional information about the structure.
![]() | 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 GLU 35 to perform our pKa calculations. We will do this with the PDB2PQR web service. PDB2PQR “cleans up” PDB files for electrostatics calculations by performing a number of operations described on the PDB2PQR homepage.
Unless otherwise directed, PDB2PQR adds hydrogens to residues based
on model pKa values. Therefore, we will need to specify the
titration state of GLU 35 for our pKa calculation by changing the
residue name from GLU
to GLH
using your favorite text
editor. Save the result as
2LZT-GLH35.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 explicit 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 GLU 35 or GLH 35 residue from
2LZT-GLU35.pqr
and
2LZT-GLH35.pqr
. Save the results as
GLU35.pqr
and GLH35.pqr
,
respectively.
Finally, we'll also need to perform electrostatics calculations on
HEWL with the uncharged sidechains. Use your favorite text editor
to zero out the charges in 2LZT-GLU35.pqr
and
2LZT-GLH35.pqr
to create
2LZT-noGLU35.pqr
and
2LZT-noGLH35.pqr
. This can be done by setting
the second-to-last column in the PQR file to zero; e.g.
ATOM 534 N GLU 35 6.456 16.408 22.487 -0.5163 1.8240 ATOM 535 CA GLU 35 6.354 15.205 23.349 0.0397 1.9080 ATOM 536 C GLU 35 6.889 13.941 22.695 0.5366 1.9080 ATOM 537 O GLU 35 7.705 13.192 23.277 -0.5819 1.6612 ATOM 538 CB GLU 35 4.906 14.973 23.796 0.0560 1.9080 ATOM 539 CG GLU 35 4.476 15.932 24.948 0.0136 1.9080 ATOM 540 CD GLU 35 5.171 15.736 26.255 0.8054 1.9080 ATOM 541 OE1 GLU 35 5.844 14.786 26.570 -0.8188 1.6612 ATOM 542 OE2 GLU 35 5.038 16.682 27.056 -0.8188 1.6612 ATOM 543 H GLU 35 5.671 16.799 22.056 0.2936 0.6000 ATOM 544 HA GLU 35 6.856 15.382 24.215 0.1105 1.3870 ATOM 545 HB2 GLU 35 4.318 15.158 23.042 -0.0173 1.4870 ATOM 546 HB3 GLU 35 4.832 14.066 24.142 -0.0173 1.4870 ATOM 547 HG2 GLU 35 4.654 16.857 24.629 -0.0425 1.4870 ATOM 548 HG3 GLU 35 3.500 15.796 25.084 -0.0425 1.4870
might become
ATOM 534 N GLU 35 6.456 16.408 22.487 0.0000 1.8240 ATOM 535 CA GLU 35 6.354 15.205 23.349 0.0000 1.9080 ATOM 536 C GLU 35 6.889 13.941 22.695 0.0000 1.9080 ATOM 537 O GLU 35 7.705 13.192 23.277 0.0000 1.6612 ATOM 538 CB GLU 35 4.906 14.973 23.796 0.0000 1.9080 ATOM 539 CG GLU 35 4.476 15.932 24.948 0.0000 1.9080 ATOM 540 CD GLU 35 5.171 15.736 26.255 0.0000 1.9080 ATOM 541 OE1 GLU 35 5.844 14.786 26.570 0.0000 1.6612 ATOM 542 OE2 GLU 35 5.038 16.682 27.056 0.0000 1.6612 ATOM 543 H GLU 35 5.671 16.799 22.056 0.0000 0.6000 ATOM 544 HA GLU 35 6.856 15.382 24.215 0.0000 1.3870 ATOM 545 HB2 GLU 35 4.318 15.158 23.042 0.0000 1.4870 ATOM 546 HB3 GLU 35 4.832 14.066 24.142 0.0000 1.4870 ATOM 547 HG2 GLU 35 4.654 16.857 24.629 0.0000 1.4870 ATOM 548 HG3 GLU 35 3.500 15.796 25.084 0.0000 1.4870
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 (download it here):
read mol pqr compound.pqr # This is the compound for which we will calculate solvation energies mol pqr ref.pqr # This is a compound used as a reference for grid centering end elec name inhom mg-auto # Focusing calculations dime 129 129 129 # This is a good grid spacing for this system cglen 52.0 66.0 79.0 # These are reasonable coarse grid settings for this system (PDB2PQR-recommended) fglen 51.0 59.0 67.0 # These are reasonable fine grid settings for this system (PDB2PQR-recommended) cgcent mol 2 # Center the grid on the reference molecule fgcent mol 2 # Center the grid on the reference molecule mol 1 lpbe bcfl sdh pdie 20.00 sdie 78.54 srfm smol sdens 20.0 chgm spl2 srad 1.40 swin 0.30 temp 298.15 calcenergy total calcforce no end # Print the final energy print energy inhom end quit
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-GLU35.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[9]
it is thought to implicitly represent
internal relaxation and rearrangement of the
solute[10].
You now have enough information to calculate total electrostatic
energies for all of the relevant molecules so far:
2LZT-GLU35.pqr
,
2LZT-GLH35.pqr
,
2LZT-noGLU35.pqr
,
2LZT-noGLH35.pqr
,
GLU35.pqr
, and GLH35.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
| teefoo.out
where foo.in
is the input file of interest and
the output is saved in foo.out
.
Recall from our earlier discussion that the transfer free energies can be evaluated in two ways. It is up to you to decide which method to use; both should give comparable results for a sufficiently fine grid. As mentioned above, the direct subtraction of total electrostatic energies is usually the most stable route, assuming identical grids and solute conformations are used for all calculations.
Any differences in results between these two transfer free energy methods are typically due to a lack of convergence in the calculations and can be resolved by decreasing the grid spacing (e.g., increasing the number of grid points).
If you choose to decompose the transfer free energies into solvation and Coulombic energy changes, then you will need to perform two additional electrostatics calculations beyond the APBS calculations outlined in the previous section:
First, you will need to re-run all of your APBS calculations with a homogeneous dielectric εs = εp to calculate solvation energies. This can be accomplished by appropriate modification of the input files from the previous section.
Second, you will need to calculate Coulombic interaction energies. The APBS program coulomb will evaluate the Coulomb's law electrostatic energy of a set of charges in a vacuum. The program's usage is
$
coulomb foo.pqr
where foo.pqr
is the particular PQR file you are
interested in. Note that this program can
also be run with no arguments to obtain an
expanded list of options. Don't forget
to scale by εp!
At this point, you should have everything you need to calculate the intrinsic pKa of interest. However, if you get stuck, I've included some example files here that might be helpful:
GLU35.pqr
, GLH35.pqr
, 2LZT-GLH35.pqr
, 2LZT-GLU35.pqr
, 2LZT-noGLH35.pqr
, 2LZT-noGLU35.pqr
PQR files described above
2LZT-noGLH35.in
, 2LZT-GLH35.in
, GLH35.in
APBS input files needed to calculate the transfer free energy for GLH directly using total electrostatic energies.
2LZT-noGLH35-vac.in
, 2LZT-GLH35-vac.in
, GLH35-vac.in
Additional APBS input files needed to calculate the transfer free energy for GLH using polar solvation energies.
2LZT-noGLU35.in
, 2LZT-GLU35.in
, GLU35.in
APBS input files needed to calculate the transfer free energy for GLU directly using total electrostatic energies.
2LZT-noGLU35-vac.in
, 2LZT-GLU35-vac.in
, GLU35-vac.in
Additional APBS input files needed to calculate the transfer free energy for GLU using polar solvation energies.
run-apbs.sh
Bash
script to run all the APBS input
files listed above. Assumes the
apbs
binary is
in your
PATH
.
run-coul.sh
Bash
script to run supporting Coulombic
energy calculations for solvation
energy-based evaluation of the
transfer free energies. Assumes the
APBS coulomb
binary (usually found in
apbs/tools/manip
)
is in your
PATH
.
process.sh
Bash
script to process APBS output and
calculate pKa shifts directly
using total electrostatic energies.
Assumes results were obtained by
running
run-apbs.sh
,
linked above.
[7] The following references provide a good introduction to biomolecular pKa calculations:
Bashford D, Karplus M. pKa's of ionizable groups in proteins: Atomic detail from a continuum electrostatic model. Biochemistry. 29 (44), 10219-25, 1990. http://www.ncbi.nlm.nih.gov/entrez/query.fcgi?cmd=Retrieve&db=PubMed&dopt=Citation&list_uids=2271649
Antosiewicz J, McCammon JA, Gilson MK. The determinants of pKas in proteins. Biochemistry. 35 (24), 7819-33, 1996. http://pubs.acs.org/cgi-bin/pagelookup?bichaw/35/7819
Nielsen JE, Vriend G. Optimizing the hydrogen-bond network in Poisson-Boltzmann equation-based pKa calculations. Proteins-Structure Function and Genetics. 43 (4), 403-12, 2001. http://www3.interscience.wiley.com/cgi-bin/abstract/78505732/ABSTRACT
[8] Such conditions might happen at extremely low ion and coion concentrations in the absence of other interacting species, etc. In other words, for real biological systems, these conditions are rarely encountered.
[9] Start vigorous waving of hands...
[10] Stop vigorous waving of hands...
APBS requires approximately 200 B memory per grid point.
Memory usage can be predicted before performing the
calculations using the tools/manip/psize.py
Python script provided with APBS.
If it appears your calculation is going to require more memory than is currently available on your system, you have a few options:
APBS calculations can be performed in
parallel across multiple processors
(hopefully, sharing distributed
memory!). This functionality is
provided by the
mg-para
keyword
which is described in more detail in
the APBS user manual and below.
APBS calculations can be broken into a
series of smaller, asynchronous runs
which (individually) require less
memory. This functionality is provided
by the mg-para
async
keyword
which is described in more detail in
the APBS user manual and below.
Submit your calculations through Gemstone to external computational resources.
The actin dimer example provided with the APBS
source code
(examples/actin-dimer/complex.pqr
)
is a fairly large system and a good example for
parallel focusing calculations. This example
will use the actin dimer complex PQR file
complex.pqr
.
We're going to use an 8-processor parallel calculation to write out the electrostatic potential map for this complex. Each processor will solve a portion of the overall problem using the parallel focusing method (see Baker et al, PNAS, 2001) on a 973 mesh with 20% overlap between meshes for neighboring processors. An example input file for this calculation might look like (download here):
read mol pqr complex.pqr end elec name complex mg-para ofrac 0.1 pdime 2 2 2 dime 97 97 97 fglen 150 115 160 cglen 156 121 162 cgcent mol 1 fgcent mol 1 mol 1 npbe bcfl sdh ion 1 0.150 2.0 ion -1 0.150 2.0 pdie 2.0 sdie 78.54 srfm mol chgm spl0 srad 1.4 swin 0.3 sdens 10.0 temp 298.15 calcenergy total calcforce no write pot dx pot end quit
where the pdime 2 2 2
specifies the 8-processor array dimensions, the
ofrac 0.1
specifies the 20%
overlap between processor calculations, and the
dime 97 97 97
specifies the
size of each processor's
calculation. The write pot dx
pot
instructs APBS to write
out OpenDX-format
maps of the potential to 8 files
pot#.dx
, where
#
is the number of the
particular processor.
Running this input file with an MPI-compiled version
of APBS runs 8 parallel focusing calculations,
with each calculation generating fine-scale
solutions on a different region of the
(fglen
) problem domain.
Note that 8 separate OpenDX files are written
by the 8 processors used to perform the
calculation. Writing separate OpenDX files
allows us to avoid communication in the
parallel run and keeps individual file sizes
(relatively) small. Additionally, if a user is
interested in a specific portion of the problem
domain, only a few files are needed to get
local potential information.
However, most users are interested in global
potentials. For some programs (OpenDX,
DataTank),
the individual potential files can simply be
read into the program separately and the
program will reconstruct the global map.
Most other programs will require the user to
reassemble the global map first; APBS provides
the mergedx
program for this purpose.
mergedx
is a simple
program that allows users to combine several
OpenDX files from a parallel focusing
calculation into a single map. This map can be
downsampled from the original resolution to
provide coarser datasets for fast
visualization, etc.
For example, the command
$
mergedx 65 65 65 pot0.dx pot1.dx pot2.dx pot3.dx pot4.dx pot5.dx pot6.dx pot7.dx
will generate a file
gridmerged.dx
which has downsampled the much larger dataset
contained in the 8 OpenDX files into a
653 file which
would be suitable for rough visualization. An
example of the mergedx
output (see Section 3.1.4.1, “Isocontour visualization”
for methods) is shown in
Note that downsampling isn't necessary -- and
often isn't desirable for high quality
visualization or quantitative analysis.
The steps described in the previous section can also be performed for systems or binaries which are not equipped for parallel calculations via MPI. In particular, you can add
async
n
to the ELEC
mg-para
APBS input file to
make the single-processor calculation
masquerade as processor
n
of a parallel
calculation.
Scalar maps from asynchronous APBS calculations
can be combined using the
mergedx
program as
described above. Currently, energies and
forces from asynchronous APBS calculations need
to merged manually (e.g., summed) from the
individual asynchronous calculation output.
This can be accomplished by simple shell
scripts.
As a specific example, we can modify the input file
above to include a
async 0
command in the ELEC
statement and thus cause APBS to perform the
operations of the first processor in the
parallel focusing calculation. The modified
input file (available for download here)
should look like:
read mol pqr complex.pqr end elec name complex mg-para ofrac 0.1 pdime 2 2 2 async 0 dime 97 97 97 fglen 150 115 160 cglen 156 121 162 cgcent mol 1 fgcent mol 1 mol 1 npbe bcfl sdh ion 1 0.150 2.0 ion -1 0.150 2.0 pdie 2.0 sdie 78.54 srfm mol chgm spl0 srad 1.4 swin 0.3 sdens 10.0 temp 298.15 calcenergy total calcforce no write pot dx pot end quit
This should create an OpenDX-format potential
map called pot.dx
,
corresponding to the output from processor 0 in
a parallel focusing calculation. Performing
additional APBS calculations with
async 1
,
async 2
,
...,
async 7
will generate the corresponding OpenDX maps for
the remaining processors of the parallel
focusing calculations. These can then be
reassembled with mergedx
as described above.
Robert Konecny (McCammon group) has developed the iAPBS package which provides a C/C++/FORTRAN interface to APBS for use with AMBER, NAMD, and CHARMM. More information is available from the iAPBS homepage.
APBS also links against developmental versions of the TINKER software package. Public versions of TINKER with APBS support should be available in Fall 2007.
Table of Contents
Yes, you can run APBS via the web, thanks to the Gemstone service. This section demonstrates how Gemstone can be used to perform most APBS calculations commonly done through the command-line interface. Better yet, it allows you to use NBCR's resources instead of your own!
You need to download and install the Gemstone extension for the Firefox web browser from http://gemstone.mozdev.org before starting this tutorial.
We'll start by preparing a structure with PDB2PQR. First, you need to download a PDB file from the Protein Data Bank. For the purposes of this example, we'll use 1FAS in PDB format (download directly here).
Open the Gemstone interface through the Firefox menu
→ . Select the PDB2PQR application through the Gemstone menu → on the righthand side of the Gemstone screen. Select the PDB file you just downloaded and any other options for PDB2PQR of interest. Your screen should look something like:
Go ahead and "Run" the Gemstone
calculation, check the
stdOut
and
stdErr
for any
warning messages, and download the
resulting PQR file.
We're now ready to perform an APBS calculation of the 1FAS electrostatic potential through Gemstone. Open the
→ menu on the righthand side of the Gemstone screen. We'll step through each of the input tabs in the middle pane in the following sections.Move to the Calculation tab in the middle pane.
Most of the "Mathematics" settings can be left as-is; however, feel free to play with the discretization and surface definitions to see how it affects your results.
Under "Molecules", load the 1FAS PQR file you just prepared in the previous section.
Don't select anything under "Energy and Force" for this calculation.
At this point, your screen should look something like this:
Next, move to the Grid tab in the middle pane.
Under "General", enter 97 for X, Y, and Z "Grid points".
Under "Coarse Grid", choose 100 for all X, Y, and Z "Lengths". Center the mesh on the molecule.
Under "Fine Grid", choose 80 for all X, Y, and Z "Lengths". Center the mesh on the molecule.
Leave the "Parallel Execution" alone for this calculation.
At this point, your screen should look something like this:
Next, move to the Physics tab in the middle pane.
Choose any values you like here; the defaults should be fine for most visualization calculations.
At this point, your screen should look something like this:
Next, move to the File I/O tab in the middle pane.
Unless you have other properties you wan to visualize, you can leave these settings as-is.
At this point, your screen should look something like this:
Finally, return to the Calculation tab and select "Run APBS". After a while, the results will appear under the Output tab in the middle pane. When the run is complete, you should see something like this:
In particular, you should have the
option to save the potential file
calculated in this APBS run.
In addition to the relatively simple applications presented here, there are a number of additional examples distributed with the APBS software package (http://apbs.sf.net) which include ionization energies, protein-protein interactions, parallel calculations, etc.
Please visit the APBS-users mailing list. After you've looked for the answer to your question in the archive, please post it to the mailing list.