Potentials of Mean Force

Poisson–Boltzmann profile for an ion channel

This example is hosted off-site here

The polar solvation potential of mean force for a helix in a dielectric slab membrane


This example will examine the differing polar solvation free energies of an alpha helix as it translates through a low-dielectric slab, a model membrane-like environment. The low dielectric slab is intended to crudely represent the nonpolar membrane environment. Note that the peptide used in this example is largely nonpolar (glycine, alanine, and leucine) with a centrally-located arginine. The starting configuration of this example has the peptide nearly completely buried in the membrane environment:

The helix (arginine in blue, nonpolar residues in white, glycine in green) in the membrane environment (gray planes)

The final configuration of this example has the peptide significantly exposed to solvent:

The helix (arginine in blue, nonpolar residues in white, glycine in green) in the membrane environment (gray planes)

This example will use APBS to solve the Poisson equation for these different configurations to determine the polar solvation energy. The resulting polar solvation energy profile (as a function of helix position) is called a “potential of mean force” for the solvation of this helix through this low dielectric slab membrane mimic.

All of the software and input files needed for this example are available here.


Software requirements

In addition to APBS, you will use the draw_membrane2 program written by Michael Grabe. The source code for this program can be compiled very easily with the following command:

gcc draw_membrane2.c -lm -o draw_membrane2


The basic steps

We will illustrate the basic steps of the polar solvation calculation with the peptide in its membrane-immersed started configuration. The PQR file for this starting configuration is called Membrane-helix-0.pqr.


Generating the peptide dielectric map

The system will have three dielectric regions:

  • The high dielectric solvent exterior (usually 80 for water)
  • The low dielectric interior of the protein (10 for this example)
  • The low dielectric interior of the membrane (2 for this example)

and two regions of ion-accessibility:

  • The ion-accessible solvent exterior
  • The ion-inaccessible interior of the protein and the membrane

We will model these regions in APBS using OpenDX-format coefficient maps read into APBS through a READ statement.

The first step in such modeling is to construct the ion-accessibility and dielectric maps for the isolated peptide in the absence of the membrane environment. This can be accomplished through APBS mg-dummy calculations. These calculations simply construct the coefficient input maps and write them out in the desired format without performing any numerical solutions of the Poisson-Boltzmann equation.

In the latter steps of this example, we’ll be calculating solvation energies using sequential focusing to solve the Poisson-Boltzmann equation. We’ll need three sets of maps for the focusing procedure:

  • A large map of lengths 300 × 300 × 300 Å3.
  • A medium map of lengths 200 × 200 × 200 Å3.
  • A small map of lengths 100 × 100 × 100 Å3.

All of these maps will have 97 × 97 × 97 grid points. This coarse grid resolution is necessary for a relatively quick example; however, it is not recommended for calculations where quantitative accuracy is desired. The APBS mg-dummy input file for these calculations is called Apbs_dummy.in. Run APBS with this input file to generate the following dielectric maps:

  • diel{x,y,z}_{S,M,L}.dx The nine dielectric maps for the x-, y-, and z- components of the Poisson operator for the small (S), medium (M), and large (L) problem domains
  • kappa_{S,M,L}.dx The three ion accessibility coefficient maps for the small (S), medium (M), and large (L) problem domains
  • charge_{S,M,L}.dx The three peptide charge distribution maps for the small (S), medium (M), and large (L) problem domains


Adding the membrane environment

The next step is to incorporate the membrane environment through the draw_membrane2 program, described above. This program is invoked separately for the small, medium and large systems as follows:

./draw_membrane2 dielx_S.dx zmem Lmem pdie 0.0 I R_top R_bottom
./draw_membrane2 dielx_M.dx zmem Lmem pdie 0.0 I R_top R_bottom
./draw_membrane2 dielx_L.dx zmem Lmem pdie 0.0 I R_top R_bottom

and assumes a common naming scheme for all dielectric, kappa, and charge maps. It also assumes that the membrane is oriented such that its plane is in the xy-plane and the bilayer normal is aligned along the z-axis. Finally, it sets the solvent dielectric to 80.0 and the membrane dielectric to 2.0. The parameters for this program are

  • zmem The location of the bottom of the membrane slab
  • Lmem The thickness of the membrane slab
  • pdie Internal dielectric of the protein
  • I Value of the symmetric ion concentration in molar
  • R_top Membrane exclusion radius in case your protein is a pore
  • R_bottom Membrane exclusion radius in case your protein is a pore

For the current example, we will invoke the draw_membrane2 program as

./draw_membrane2 dielx_S.dx -20 40 10.0 0.0 0.1 0.0 0.0
./draw_membrane2 dielx_M.dx -20 40 10.0 0.0 0.1  0.0 0.0
./draw_membrane2 dielx_L.dx -20 40 10.0 0.0 0.1 0.0 0.0

This will produce a number of *m.dx files (e.g., dielx_Lm.dx) that are the modified coefficient maps incorporating the membrane environment.


Performing the electrostatic calculation

Finally, we will use these input coefficient maps in an APBS calculation to solve the linearized Poisson-Boltzmann equation using sequential focusing. The APBS input file, Apbs-solv.in, performs six calculations:

  • Three sequential focusing calculations for the problem with the peptide in solution (“solvated”)
  • Three sequential focusing calculations for the problem with the peptide in the membrane dielectric (“reference”)

The difference between the finest “solvated” calculation and the finest “reference” calculation is an estimate of the polar solvation energy associated with transferring the peptide from solution into the membrane.


Visualizing the result

The result of these electrostatic potential calculations can be visualized with VMD as described elsewhere in this documentation. Perhaps one of the most interesting VMD features to use with this example is the FieldLines representation of the electrostatic potential which demonstrates how the dielectric interface between the membrane and bulk solution significantly perturbs the electrostatic field:


An arginine-containing alpha helix in a low dielectric (membrane) slab. The gray surfaces represent the surface of the membrane. The blue and red surfaces are +1 and -1 kT/e electrostatic isocontours, respectively. The lines show the direction of the electric field in the system.

This image was generated using VMD 1.8.6. The membrane was represented as an isocontour of the dielx_Sm.dx map with a value of 72. The field lines were colored by potential values from pot_S.dx and displayed with a GradientMag value of 1.61, a minimum length of 7.24, and a maximum length of 50.00. Finally, the isocontours were displayed from the pot_S.dx file with values of +1 kT/e (blue) and -1 kT/e (red).


Putting it all together into a PMF

The steps above should be repeated for every structure in the potential of mean force:

  • membrane-helix-0.pqr
  • membrane-helix-4.pqr
  • membrane-helix-8.pqr
  • membrane-helix-12.pqr
  • membrane-helix-16.pqr

and doing so demonstrates an interesting trend in solvation energies:

Displacement (A) Polar solvation energy (kJ/mol) Relative change in energy (kJ/mol)
0 127.1 68.31
4 127.2 68.39
8 119.9 61.16
12 85.74 26.96
16 58.78 0.000

What causes the shape of this curve?


Automating the PMF calculation

Repeating the above calculations by hand can be very tedious; the bash script, Run_membrane-helix.sh, was designed to simplify the process. Note that, in addition to the software described above, this script requires:

  • A template for the APBS coefficient map construction named apbs_dummy-TEMPLATE.in.
  • A template for the APBS coefficient map construction named apbs_solv-TEMPLATE.in.
  • All of the PQR files from above named as membrane-helix-0.pqr, membrane-helix-4.pqr, …

You will also need to modify the bash script to specify the correct locations of your APBS and draw_membrane2 executables.


Caveats and comments

As mentioned above, this calculation was performed using a grid with 97 × 97 × 97 points for the sake of efficiency. However, this is almost certainly too coarse for quantitative work. For real applications, one should increase the number of grid points until the desired observable (e.g., solvation energy) stops changing with increasing grid points.

Additionally, this calculation in no way represents a real membrane. Instead, it is a model of the low dielectric environment commonly assumed for membrane systems. Of course, real membranes deform, include some water permeation, form defects of large and small scale, and have very inhomogeneous dielectric properties. Additionally, real peptides can also deform and rearrange sidechains to maximize solvent exposure of polar or charged groups.