PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of plosonePLoS OneView this ArticleSubmit to PLoSGet E-mail AlertsContact UsPublic Library of Science (PLoS)
 
PLoS One. 2010; 5(9): e12722.
Published online 2010 September 29. doi:  10.1371/journal.pone.0012722
PMCID: PMC2947494

APBSmem: A Graphical Interface for Electrostatic Calculations at the Membrane

Darren R. Flower, Editor

Abstract

Electrostatic forces are one of the primary determinants of molecular interactions. They help guide the folding of proteins, increase the binding of one protein to another and facilitate protein-DNA and protein-ligand binding. A popular method for computing the electrostatic properties of biological systems is to numerically solve the Poisson-Boltzmann (PB) equation, and there are several easy-to-use software packages available that solve the PB equation for soluble proteins. Here we present a freely available program, called APBSmem, for carrying out these calculations in the presence of a membrane. The Adaptive Poisson-Boltzmann Solver (APBS) is used as a back-end for solving the PB equation, and a Java-based graphical user interface (GUI) coordinates a set of routines that introduce the influence of the membrane, determine its placement relative to the protein, and set the membrane potential. The software Jmol is embedded in the GUI to visualize the protein inserted in the membrane before the calculation and the electrostatic potential after completing the computation. We expect that the ease with which the GUI allows one to carry out these calculations will make this software a useful resource for experimenters and computational researchers alike. Three examples of membrane protein electrostatic calculations are carried out to illustrate how to use APBSmem and to highlight the different quantities of interest that can be calculated.

Introduction

The relationship between the electric field and the charge in a system is determined by Maxwell's equations; however, several factors contribute to making these equations difficult to solve in a heterogeneous, condensed phase. The most popular method for carrying out electrostatic calculations in a biological setting is to solve the Poisson-Boltzmann (PB) equation. Starting from a known protein structure, this method treats the protein and water as distinct dielectric environments, and the charges on the protein give rise to the electric field. Additionally, PB theory implicitly accounts for counter-ions in solution via a non-linear term that depends on the bulk counter-ion concentration and the electrostatic potential. The PB equation for a one-to-one electrolyte solution is:

equation image
(1)

where An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e002.jpg is the reduced electrostatic potential and An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e003.jpg is the electrostatic potential, An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e004.jpg is the Debye-Hückel screening parameter, which accounts for ionic shielding, An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e005.jpg is the dielectric constant for each of the distinct microscopic regimes in the system, and An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e006.jpg is the density of charge within the protein moiety. Since the 1980s, researchers have studied the electrostatic properties of protein and nucleic acid systems by numerically solving the PB equation using finite difference and finite element methods [1][4]. Today, there are several popular software packages available to perform PB calculations such as DelPhi [5], the Adaptive Poisson-Boltzmann Solver (APBS) [6], MIBPB [7], ZAP [8], and the PBEQ module in CHARMM [9]. Unfortunately, there is a dearth of programs that allow researchers to carry out these calculations at or near a membrane. Nonetheless, over the last two decades the number of high-resolution membrane protein structures has dramatically increased. The membrane has several unique electrical properties. For instance, the core of the membrane is extremely hydrophobic giving rise to a desolvation penalty for moving charged molecules into this region. This property is essential to the membrane's ability to control the flow of materials into and out of the cell. Additionally, most cells have a substantial membrane potential that coordinates the action of voltage-dependent membrane proteins such as voltage-gated ion channels. Without including the effects of the membrane dielectric and the transmembrane potential, there is a huge class of molecules whose electrical properties cannot easily be explored.

Groups have carried out simulations using several different levels of theory to include the effects of the membrane such as fully atomistic calculations (for an incomplete list of references see [10][16]), implicit membrane calculations using Generalized Born models (for an incomplete list see [17][29]), and continuum approaches employing the PB equation (for an incomplete list see [30][41]); however, all of these studies require the user to have a high level of computational sophistication as highlighted by the relatively few papers from non-computational laboratories. Thus, it is desirable to have a program that removes many of the technical and bookkeeping aspects from these calculations. Toward this end, very recently, an online web server was created to facilitate PB calculations on membrane proteins using the PBEQ module [42], and here, we present our program, APBSmem, which combines an easy to use interface with APBS to allow non-experts to calculate the electrostatic properties of membrane proteins. APBSmem can be downloaded, easily installed, and run locally on Windows, Mac, and Linux platforms. APBSmem has several pull down templates that allow researchers to carry out specific membrane related calculations, and it has a built-in graphical window that provides quick visual feedback to make sure that the system is set up correctly and to view results.

Methods

User interface

Though the majority of the calculations described here may be performed using APBS input files, keeping track of the files and parameters can become quite difficult. To improve this process we built a Java-based GUI that writes the input files and runs the calculations (Figure 1). The GUI has an embedded Jmol [43] viewer that allows users to visualize the protein-membrane system and the electrostatic potential. Here we explain the necessary parameters and use of the interface in a step-by-step fashion.

Figure 1
A screenshot of the user interface.

Calculation Type

To perform a calculation, the user first selects a type (Protein Solvation, Ion Solvation, or Gating Charge) from the drop down menu. Each type is described in more detail in the case studies below. The user then selects coordinate files, in PQR format, for all of the protein configurations of interest. PQR files contain the atomic positions of all of the atoms in the system in addition to their charge and radii. PQR files can be generated from PDB files with the PDB2PQR tool [44], which allows the user to choose from several commonly used charge and radii parameter sets: PARSE [45], CHARMM27 [46], Swanson [47], AMBER99 [48], along with several other sets and user defined values. This choice is crucial since calculations can be sensitive to the parameter set [49] especially when performing solvation energy calculations. File locations may be entered manually or found and selected from the filesystem with the Browse button. At present, APBSmem does not allow the spatial orientation and placement of the protein to be altered once read in through the GUI, and external software must be employed if a different orientation is desired. For Protein Solvation calculations, the user should provide a coordinate file with only the membrane protein. Ion Solvation calculations require a PQR file with only the protein and a PQR file with only the ion. Two files corresponding to an open and a closed channel are needed in the case of a Gating Charge calculation.

Grid Spacing

Next, the user must specify the desired level of discretization, which is related to the fidelity with which the underlying equations will be solved. It is necessary to apply the appropriate boundary conditions far from the protein to accurately solve electrostatics calculations, and this requires large grid lengths to remain computationally tractable. However, coarse discretization does not capture the correct electrostatic behavior near the protein, where small grid spacing is needed. A technique known as focusing is employed to rectify this problem by solving the equations in a series of steps starting at the largest length scale and focusing into the smallest length scale [5]. When using multiple levels of focusing, the user can set the desired level in the Focus menu to enable fields for additional grid lengths.

With any numeric calculation, the accuracy of the solution is directly related to the degree of discretization. It is important to check the convergence properties of your solution. This is typically done by recalculating with increasing numbers of grid points without changing any of the other parameters. The exact convergence properties depend on the numeric algorithm and the details of the system. In Figure 2, we calculated the convergence for each test case in the Results section. However, all systems behave differently, and users should not assume that discretization schemes that give accurate results here will also give accurate values for other protein-membrane systems.

Figure 2
Convergence properties of test cases I–III.

Dielectric Parameters

In continuum electrostatics, different regions of the system are defined by unique dielectric values. These values are related to the polarizability of each region in response to an applied electric field. The protein dielectric value is assigned to points within the protein's solvent accessible surface. All points outside the protein, but within the membrane region defined by the geometry settings, are assigned a dielectric value corresponding to the membrane. While the core of the membrane often has a very low dielectric value, the head group region may be characterized by a much higher value. If desired, this physical feature can be included in calculations by increasing the default head group thickness and setting the head group dielectric value. All other points in the system are assigned the solvent dielectric value.

Proteins are heterogeneous, and it is not always appropriate to describe the entire molecule with a single uniform dielectric value [50]. Nonetheless, uniformity is a common assumption of PB solvers. Experiments indicate that the protein interior is modeled best by dielectric values between 2 and 20 [51]. With this in mind, we recommend that researchers test the robustness of their results by repeating calculations with several different dielectric constant values within this range.

Boundary Conditions

Several options are offered for Dirichlet boundary conditions when solving the PB equation in APBS. The user may set all boundaries to zero, use a single Debye-Hückel model, multiple Debye-Hückel model, or focusing, in which the boundaries are determined by a previous calculation. When the Gating Charge calculation type is chosen, the boundary condition is set to impose a range of membrane potentials across the membrane as described in Case III. The user provides a membrane potential value in milliVolts, An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e012.jpg, and the interface performs a sweep of calculations from An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e013.jpg to An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e014.jpg to determine the valence of the gating motion. At present, calculations with a membrane potential are only carried out in the linearized limit of Eq. 1. Additionally, application of the boundary conditions ignores differences between the dielectric values of the head group and the membrane core, which has been included in a recent study [52].

Protein Surface Representation

An accurate representation of the protein surface is important in constructing the dielectric and ion-accessibility maps. A probe-based dielectric function is used to construct the protein surface in APBS. The solvent probe radius specifies the size of water spheres for the determination of the solvent space and is typically set to 1.4 Å for water. The surface sphere density determines the resolution at which this boundary is calculated and is typically set to 10 grid points/ÅAn external file that holds a picture, illustration, etc.
Object name is pone.0012722.e015.jpg.

System Geometry

The system geometry parameters determine the shape and location of the membrane. The membrane thickness and vertical position should be adjusted for the protein of interest so that the bilayer interfaces with the protein as expected. It is difficult to determine this placement, and it is an ongoing area of research in the Grabe lab [53]. In practice, this placement is often done ad hoc based on the location of the hydrophobic residues making up the membrane spanning region. A better alternative is to first estimate the orientation of the membrane protein and extent of the membrane spanning region by using the Orientations of Proteins in Membranes database [54]. A second point of concern is that ion channels and transporters often have hydrophilic, water-filled cavities essential for transport. Users must employ the interface exclusion radii to prevent APBSmem from rewriting these cavities as membrane. These radii should be adjusted so that central cavities, if present, are filled with water and the membrane is arranged flush with the outside of the protein. We intend to provide automated cavity detection in future releases. Figure 3 compares correct and incorrect configurations of the membrane geometry with the KcsA potassium channel.

Figure 3
Top view of the KcsA channel (green) and the An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e016.jpg = 2.01 isocontour highlighting the membrane interface (gray).

System Preview

After the user has entered parameters for the dielectric environment and membrane geometry, the Preview button can be used to visualize the system. This Preview action performs a quick “dummy” calculation with coarse grid dimensions to generate the numeric representation of the membrane and display it graphically. If the system and parameters appear to be correct, the user can click Run to perform the calculation with APBS. When the calculations have completed, the total energy is given in kJ/mol, kcal/mol and kAn external file that holds a picture, illustration, etc.
Object name is pone.0012722.e018.jpgT, and the most focused dielectric map of the membrane is displayed in the Jmol panel. The electrostatic potential may also be viewed in the Jmol panel by selecting the Draw Potential option. The user provides an isocontour value of interest and the interface displays the positive (red) and negative (blue) surfaces over the protein. The exterior bulk and cavity (if any) at the interior of the protein are modeled into APBS as coefficient maps (openDX-format). These maps include dielectric maps (diel), ion-accessibility coefficient maps (kappa) and charge distribution maps (charge) for different regions of the protein-membrane complex. All input and output files, including the potential and DX maps are saved for later use and reference.

Membrane potential boundary conditions

In a typical cell, electrogenic transporters create a difference between the electrical potential inside the cell, An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e019.jpg, and outside the cell, An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e020.jpg. A small violation in electroneutrality near the membrane gives rise to this potential difference; however, more than a Debye length from the membrane, electroneutrality is restored. It is possible to model this behavior with the Poisson-Boltzmann equation as outlined by Roux in his seminal work on this topic [55]. Most researchers will want to compute the membrane potential in the absence of protein charges to understand the profile across the protein, and in some cases, they will be interested in computing the interaction energy of the membrane electric field with the charges on the protein. The field due to the protein charges and the membrane potential are only separable when solving the linearized form of the equation. Thus, in order to be self consistent, APBSmem only solves the linearized Poisson-Boltzmann equation when employing membrane potential boundary conditions. In future releases, we will extend this to the full non-linear equation. The presentation in this section and the next largely follows the supporting text found in Grabe et al. [38], but the essence is similar to Roux's [55]. We start by rewriting Eq. 1 in the linearized form:

equation image
(2)

However, this equation does not satisfy the asymptotic boundary condition: An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e022.jpg. This oversight can be corrected by adding a constant term to the equation for all positions in the inner solution space:

equation image
(3)

where An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e024.jpg is An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e025.jpg for all points in the inner solution space and zero otherwise (see Figure 4). Now far from the protein where An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e026.jpg is no longer changing, An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e027.jpg as desired. Eq. 3 can be rewritten as:

equation image
(4)
Figure 4
A cartoon representation of the distinct dielectric environments in each calculation.

Thus, the modified Poisson-Boltzmann equation above takes the form of Eq. 2 with the membrane potential arising from a term that looks like a uniform source charge. The spatial dependence of An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e038.jpg is carried by An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e039.jpg on the right hand side. Since Eq. 4 is linear, it is possible to separate the total reduced electrostatic potential, An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e040.jpg, into contributions from the membrane potential, An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e041.jpg, and contributions from the protein, An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e042.jpg, as An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e043.jpg. Each field is the solution to corresponding equations as shown:

equation image
(5)

Far away from the protein, An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e045.jpg approaches zero. Poisson-Boltzmann solvers typically set zero boundary conditions at the outer boundary to account for this behavior, or they use some asymptotic approximation to the field based on the protein's total charge. In the case of a membrane potential, the behavior of An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e046.jpg far from the membrane protein is required so that far field boundary conditions can be imposed on the system.

An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e047.jpg on the boundary is determined by considering a planar slab of low-dielectric material with symmetric electrolyte solution in the half-spaces above and below the slab. This follows the work of Roux with a slight change in geometry [55]. By symmetry An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e048.jpg, and we assign An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e049.jpg to the center of the membrane. The slab has a length An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e050.jpg, and there are three distinct regions of space: An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e051.jpg (out); An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e052.jpg (membrane); An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e053.jpg (in). The dielectric of water is An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e054.jpg, and the dielectric constant of the membrane is assigned An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e055.jpg. We assume that ions cannot enter the membrane so An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e056.jpg is set to An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e057.jpg in this region, while the inner and outer spaces have the same value of the screening parameter. According to Eq. 5 the An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e058.jpg satisfies the following equations in each region

equation image
(6)

From elementary electrostatics, we know that the potential is continuous at the membrane boundaries but the z-component of the electric field is discontinuous due to the jump in dielectric value:

equation image
(7)

The potential profile can be determined from Eqs. 6 and 7:

equation image
(8)

where An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e062.jpg. When membrane potential calculations are performed, Eq. 8 is used to set An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e063.jpg on the domain boundary. This requires first providing the z-position of the top and bottom of the membrane and the dielectric constants of the membrane and water.

Addition of the membrane

The influence of the membrane must be included in the calculation. Based on the structure file provided, the program calls on APBS to generate dielectric (An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e064.jpg), charge (An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e065.jpg), and ion-accessibility maps (An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e066.jpg) of the molecule as if it were in solution. The protein dielectric value can be set to any value, and the method for delineating the solvent boundary is also configurable. At present, the GUI allows up to 2 levels of focusing, which corresponds to 3 sets of maps produced at this initial stage. Maps are then modified to add the presence of a low-dielectric slab acting as a surrogate membrane. APBS is run with the finite differencing scheme option; therefore, all map points are associated with a regular grid in 3-space. Next, the initial maps are read by a second routine, and the numeric values of points on the grid are modified based on the spatial position and the user-defined placement of the membrane. The program iterates over every grid position and evaluates each position in the following order:

1) Determine if the point is inside the provided protein

If the initial dielectric map value equals An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e067.jpg, the point is located within the protein. Dielectric map values are not changed for these points.

2) Determine if the point is inside the membrane

If the point does not fall within the protein, it falls within the An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e068.jpg-extent of the membrane determined by An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e069.jpg and An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e070.jpg, and it falls outside the cylinder described by the exclusion radii, the value of the dielectric map is set to An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e071.jpg, the ion-accessibility is set to zero, and the charge map is not changed.

3) Determine if the point is in the inner solution space

If the point is below the membrane and the ion-accessibility is not zero, then the neutral charge map is modified for the calculations of An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e072.jpg. The value assigned to the charge map position is determined from Eq. 5 (bottom equation). The effective charge density, An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e073.jpg, follows from the right hand sides of the upper and lower equations:

equation image

The text maps are written in terms of the number density, An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e075.jpg, and using this along with the definition of the Debye length we arrive at the modified value for the charge map

equation image

where An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e077.jpg is the molar concentration of one of the salt species (assumed balanced) and An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e078.jpg is Avogodro's number. The Debye constant above is twice the value that can be found on page 497 of Jackson's Classical Electrodynamics (Second Edition) [56], since we assume that there are mobile cationic and anionic species, not just one mobile species. Simplifying this equation we arrive at:

equation image

where An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e080.jpg is the reduced inner potential and the counter-ion concentration is given in moles per liter. The effective number density is now in inverse Å ngstroms cubed, which is consistent with the APBS solver.

Results

APBSmem was developed in Java and requires Java Runtime Environment 5.0 or later and APBS version 1.2.0 or later which can be downloaded from http://java.sun.com/ and http://www.poissonboltzmann.org/, respectively. The program can be run from the command-line using java -jar apbsmem.jar. Three case studies are presented here to demonstrate potential calculations. All files necessary to perform these calculations are packaged with the APBSmem program.

CASE I: Protein Solvation

The cell membrane is composed of lipid molecules and hosts membrane proteins which account for a third of all proteins in a cell. The hydrophobic core of the membrane provides a dielectric barrier against polar and charged molecules. The transmembrane segments of membrane proteins are therefore largely composed of hydrophobic residues; but charged and polar residues are also sometimes present, so it is natural to ask how these charged residues can be stably accomodated in the membrane. Choe et al. [53] investigated this question using continuum electrostatics with APBS. Here we revisit this problem to demonstrate the applicability of our graphical interface, and we do this by calculating the solvation energy required to insert a charged helix into the membrane. The total energy of a simple An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e081.jpg-helix in bulk water (Figure 5B) is first computed and then subtracted from the total energy of the helix embedded in the membrane (Figure 5A).

Figure 5
States used to compute protein solvation energies.

Using APBSmem to compute the protein solvation energy requires the protein to be read in as PQR file 1. The system of interest for this case study is an An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e086.jpg-helix composed of 27 residues, aligned along the z-axis and centered at the origin. The helix is composed of nonpolar hydrophobic residues with the exception of a charged arginine at the center. The protein solvation energy calculation is performed on a 161An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e087.jpg grid using two levels of focusing from a cube with side length 200 An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e088.jpg to a cube of side length 50 An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e089.jpg. The bathing solution contains 0.1 M symmetric monovalent salt with 2 Å probe radii. The protein is assigned a dielectric value of An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e090.jpg = 5, bulk water is assigned a value of An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e091.jpg = 80, and membrane is assigned a dielectric of An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e092.jpg = 2. The head group is modeled as a region of high dielectric, An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e093.jpg = 80, with a thickness of 8 Å. Calculations are carried out with the linearized PB equation with a solvent probe radius of 1.4 Å, a surface sphere density of 10 gridpoints/ÅAn external file that holds a picture, illustration, etc.
Object name is pone.0012722.e094.jpg and a temperature of 298.15 K. The total membrane thickness is 42 Å running from z = −21 Å to z = +21 Å. The upper and lower exclusion radii are set to 0 Å since there is no pore. Parameters are summarized in Table 1.

Table 1
Parameters for protein solvation CASE I.

For this system, we obtain a protein solvation energy of 28 kcal/mol, and Figure 2A indicates good convergence with grid spacing smaller than 0.781 Å at the finest level. While this energy is large, it is greatly reduced when non-polar energies are considered. Additionally, a large component of this energy is the cost of inserting the charged arginine. If the arginine is replaced with an alanine, the solvation energy drops to 4 kcal/mol. It has been shown that the electrostatic component of the membrane deformation energy can be considerably reduced by allowing the membrane to bend around the charged residue in the core of the membrane [53]. We will incorporate membrane bending and non-polar energy terms in future releases of APBSmem.

CASE II: Ion Solvation

The primary role of ion channels is to facilitate the movement of ions across the dielectric barrier imposed by the lipid bilayer. The hydrated ions in the bulk water are essentially stripped of water molecules (depending on the channel pore size) upon entering a low-dielectric medium [57][59]. The total ion solvation free energy of an ion consists of a Born solvation term, which corresponds to the removal of water molecules away from the ion and an electrostatic term that corresponds to interaction between protein charges and the ion. APBSmem calculates the ion solvation free energy by first computing the total energy of the protein-ion assembly embedded in the membrane and then subtracting the energies of the membrane-embedded protein without the ion and the energy of the ion in bulk water.

Roux and MacKinnon carried out a classic study using this approach to investigate the transfer energy for a single KAn external file that holds a picture, illustration, etc.
Object name is pone.0012722.e103.jpg from bulk water to the central cavity of the potassium channel KcsA [60]. Here we revisit this calculation. KcsA (PDB ID 1BL8) is aligned along the z-axis and centered at the origin. The ion solvation calculation requires: a PQR file with only the KcsA ion channel and a PQR file consisting of a KAn external file that holds a picture, illustration, etc.
Object name is pone.0012722.e104.jpg ion at the coordinate of interest. The ion transfer free energy is calculated using a finite difference method on a 161An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e105.jpg grid with two levels of focusing from a cubic system of side length 300 An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e106.jpg to a cube of side length 60 An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e107.jpg. The bathing solution contains 0.1 M symmetric monovalent salt with 2 Å probe radii. The protein is assigned a dielectric interior of An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e108.jpg = 2, bulk water above and below the membrane, a dielectric of An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e109.jpg = 80, and a low-dielectric slab of dielectric value An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e110.jpg = 2 represents the membrane. The separate dielectric for the head group region (An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e111.jpg = 80) is not used since its thickness is set to zero. The linearized PB equation is solved using focused boundary conditions (one level of focusing) at 298.15 K in the absence of membrane potential. The solvent probe radius is set to 1.4 Å and a surface sphere density of 10 gridpoints/ÅAn external file that holds a picture, illustration, etc.
Object name is pone.0012722.e112.jpg is used. The z-position of the bottom of the membrane and thickness of the membrane slab are set to −12 Å and 24 Å, respectively. Membrane exclusion radii of 24 Å and 16 Å are used for the channel at the top and bottom, respectively (Figure 3C). Parameters are summarized in Table 2.

Table 2
Parameters for ion solvation free energy CASE II.

APBSmem performs nine calculations: three sequential focusing calculations on the protein-ion system embedded in the membrane (Figure 6A), three sequential focusing calculations on just the protein in the membrane (Figure 6B) and three sequential focusing calculations on the KAn external file that holds a picture, illustration, etc.
Object name is pone.0012722.e122.jpg ion in bulk water (Figure 6C). Note that the system in Figure 6C computes the self energy of KAn external file that holds a picture, illustration, etc.
Object name is pone.0012722.e123.jpg in bulk water. APBSmem obtains the ion solvation energy by subtracting the energy values obtained from the fine grid calculation of systems in Figure 6B and 6C from the system in Figure 6A, and a grid spacing of 0.625 Å or smaller gives well converged values (Figure 2B).

Figure 6
States used to compute ion solvation energies.

Using these parameters, the calculated transfer free energy (from bulk water to the center of the cavity) is 7.5 kcal/mol for a single KAn external file that holds a picture, illustration, etc.
Object name is pone.0012722.e125.jpg ion when protein charges are turned off. When two KAn external file that holds a picture, illustration, etc.
Object name is pone.0012722.e126.jpg ions (blue spheres in Figure 6) are present in the selectivity filter, the calculated transfer free energy increases to 16.2 kcal/mol. This is due to electrostatic repulsion between the KAn external file that holds a picture, illustration, etc.
Object name is pone.0012722.e127.jpg ions in the selectivity filter and the incoming KAn external file that holds a picture, illustration, etc.
Object name is pone.0012722.e128.jpg ion. Upon turning protein charges on and in the presence of two KAn external file that holds a picture, illustration, etc.
Object name is pone.0012722.e129.jpg ions in the selectivity filter, the transfer free energy drops to −8.3 kcal/mol. Four pore helices (residues 62–74) along with the two An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e130.jpg ions in the selectivity filter account for an ion transfer free energy of −3.5 kcal/mol. While there are minor differences between some of our calculated values and those of Roux and MacKinnon (see Table 3), we believe that the same conclusions can be drawn from our values.

Table 3
Ion solvation free energy (kcal/mol).

CASE III: Gating Charge

Voltage-gated ion channels are sensitive to changes in membrane potential. The charged residues of the channel experience a force due to the electric field across the membrane-channel complex, and this force drives the channel to open and closed conformations as the membrane potential changes. The voltage dependence of conformational changes can be described by an equivalent “gating charge” or “sensor valence” that is defined as the fraction of the membrane electric field traversed by charges on the protein during the gating process. Thus, a gating charge of 1 indicates that a unit charge has moved through the entire membrane electric field. The gating charge often adopts non-integer values, and the higher the gating charge of a channel, the steeper its voltage dependence. The theory for using continuum electrostatic calculations to determine sensor valence was developed previously [55]. Briefly, the modified PB equation considers the transmembrane potential and calculates the interaction energy of protein charges with the field.

Here we use the murine voltage dependent anion channel 1 (mVDAC1) to illustrate gating charge calculations using APBSmem. The X-ray crystal structure of mVDAC1 shows that it is a 19-stranded An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e141.jpg-barrel with an N-terminal An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e142.jpg-helix thought to be mVDAC1's primary voltage sensor [61]. Both PB and Poisson-Nernst-Planck (PNP) electrostatic calculations on mVDAC1 suggested that the structure represents the open state of the channel [62]. This case study examines the plausibility of a hypothetical gating motion of the channel, ruled out by Choudhary and co-workers [62]. We consider a gating motion in which the N-terminal helix moves out of the channel and into the outer bath, as shown in Figure 7.

Figure 7
Hypothetical gating motion involving movement of N-terminal helix (green) out of the pore and into the outer bath.

The gating charge calculation for this gating motion requires two PQR files - mVDAC1 (PDB ID 3EMN) and a hypothetical closed state structure, to be read in as PQR file 1 and PQR file 2, respectively. We first align mVDAC1 and the hypothetical closed state structure along the z-axis and center them at the origin. The gating charge calculations in this study are carried out on a 161An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e143.jpg grid with two levels of focusing from a cubic system with side length 300 An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e144.jpg to a smaller cubic system of side length 60 An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e145.jpg. The bathing solution contains 0.1 M symmetric monovalent salt with 2 Å probe radii. The influence of the membrane is included as a dielectric slab of value An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e146.jpg = 2. Water is assigned a dielectric value of An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e147.jpg = 80, and the protein dielectric is set to An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e148.jpg = 5. The head group dielectric (An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e149.jpg = 80) is only a placeholder variable since its thickness is zero.

The linearized PB equation (lpbe) is solved using focused boundary conditions with one level of focusing at 298.15 K. The interface varies the membrane potential of the inner bath from −50 mV to +50 mV, keeping the potential of the outer bath constant at 0 mV, as shown in Figure 7. A solvent probe radius of 1.4 Å and a surface sphere density of 10 gridpoints/ÅAn external file that holds a picture, illustration, etc.
Object name is pone.0012722.e150.jpg is used. The z-position of the bottom of the membrane and thickness of the membrane slab are set to −14 Å and 28 Å, respectively. The upper and lower exclusion radii for the membrane are both set to 18.5 Å. All the parameters used for this case study are summarized in Table 4.

Table 4
Parameters for gating charge calculation CASE III.

APBSmem performs PB calculations to determine the membrane potential's contribution to the energy difference between mVDAC1, An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e160.jpg, and the hypothetical closed structure, An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e161.jpg (An external file that holds a picture, illustration, etc.
Object name is pone.0012722.e162.jpg). The energy difference is due to interaction of the protein charges with the membrane electric field. Note that the N-terminal helix has a net charge of +2. The slope of the voltage dependence curve is a measure of voltage-sensor valence which is 1.58 e in this case. This value is very close to that obtained by Choudhary and co-workers [62]. These calculations are useful for determining the voltage sensitivity of a proposed gating mechanism, and within 2.5% of the best estimate they converge to a coarse grid of 1 Å (Figure 2). As long as researchers can provide models of hypothetical transitions, these gating calculations can be used to help evaluate their biophysical correctness.

Discussion

APBSmem is an easy to use software package that carries out electrostatic calculations in the presence of a membrane. We have provided three common cases of interest to researchers in this field. The first calculates the electrostatic penalty of moving charged proteins into the membrane. This has implications for the stability of membrane proteins and for the design of membrane-permeable molecules. The second example examines the electrostatic energy for moving ions through or into ion channels and transporters. Finally, we showed how APBSmem can be used to determine the voltage dependence of a particular molecular movement. As noted earlier in the protein solvation case study, the membrane is modeled as a dielectric slab of variable thickness. Choe et al. [53] discussed the significant effects of membrane bending and its relationship to charged particles. APBSmem will eventually be expanded to identify optimal membrane deformations near the embedded molecule to provide a more complete picture of membrane protein energetics.

APBSmem has been tested on Linux, Mac OS X, and Windows, and both source code and binaries are available for download at http://mgrabe1.bio.pitt.edu/apbsmem.

Acknowledgments

We thank Katherine Ketchum for her help in constructing the hypothetical structure in Figure 7B. We also thank Seungho Choe for his help throughout this project.

Footnotes

Competing Interests: The authors have declared that no competing interests exist.

Funding: This work was supported by National Science Foundation CAREER Grant MCB0845286. M.G. is a Research Fellow of the Alfred P. Sloan Foundation. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

References

1. Warwicker J, Watson H. Calculation of the electric potential in the active site cleft due to α-helix dipoles. Journal of Molecular Biology. 1982;157:671–679. [PubMed]
2. Zauhar R, Morgan R. A new method for computing the macromolecular electric potential. Journal of Molecular Biology. 1985;186:815–820. [PubMed]
3. Klapper I, Hagstrom R, Fine R, Sharp K, Honig B. Focusing of electric fields in the active site of Cu-Zn superoxide dismutase: effects of ionic strength and amino-acid modification. Proteins: Structure, Function, and Genetics. 1986;1:47–59. [PubMed]
4. Warwicker J. Continuum dielectric modelling of the protein-solvent system, and calculation of the long-range electrostatic field of the enzyme phosphoglycerate mutase*. Journal of Theoretical Biology. 1986;121:199–210. [PubMed]
5. Gilson M, Sharp K, Honig B. Calculating electrostatic interactions in biomolecules: method and error assessment. Proceedings of the National Academy of Sciences. 1987;9:327–335.
6. Baker N, Sept D, Joseph S, Holst M, McCammon J. Electrostatics of nanosystems: application to microtubules and the ribosome. Proceedings of the National Academy of Sciences. 2001;98:10037–10041. [PubMed]
7. Zhou Y, Feig M, Wei G. Highly accurate biomolecular electrostatics in continuum dielectric environments. Journal of Computational Chemistry. 2008;29:87–97. [PubMed]
8. Grant J, Pickup B, Nicholls A. A smooth permittivity function for Poisson-Boltzmann solvation methods. Journal of Computational Chemistry. 2001;22:608–640.
9. Brooks B, Brooks C, III, Mackerell A, Jr, Nilsson L, Petrella R, et al. CHARMM: the biomolecular simulation program. Journal of Computational Chemistry. 2009;30:1545–1614. [PMC free article] [PubMed]
10. Tajkhorshid E, Nollert P, Jensen M, Miercke L, O'Connell J, et al. Control of the selectivity of the aquaporin water channel family by global orientational tuning. Science. 2002;296:525–530. [PubMed]
11. Bernèche S, Roux B. Energetics of ion conduction through the K+ channel. Nature. 2001;414:73–77. [PubMed]
12. MacCallum J, Moghaddam M, Chan H, Tieleman D. Hydrophobic association of α-helices, steric dewetting, and enthalpic barriers to protein folding. Proceedings of the National Academy of Sciences. 2007;104:6206–6210. [PubMed]
13. Bond P, Faraldo-Gómez J, Deol S, Sansom M. Membrane protein dynamics and detergent interactions within a crystal: A simulation study of OmpA. Proceedings of the National Academy of Sciences. 2006;103:9518–9523. [PubMed]
14. Petrache H, Zuckerman D, Sachs J, Killian J, Koeppe R, II, et al. Hydrophobic matching mechanism investigated by molecular dynamics simulations. Langmuir. 2002;18:1340–1351.
15. Krepkiy D, Mihailescu M, Freites J, Schow E, Worcester D, et al. Structure and hydration of membranes embedded with voltage-sensing domains. Nature. 2009;462:473–479. [PMC free article] [PubMed]
16. Allen T, Kuyucak S, Chung S. Molecular dynamics study of the KcsA potassium channel. Biophysical Journal. 1999;77:2502–2516. [PubMed]
17. Khelashvili G, Harries D, Weinstein H. Modeling Membrane Deformations and Lipid Demixing upon Protein-Membrane Interaction: The BAR Dimer Adsorption. Biophysical Journal. 2009;97:1626–1635. [PubMed]
18. Khelashvili G, Weinstein H, Harries D. Protein diffusion on charged membranes: a dynamic mean-field model describes time evolution and lipid reorganization. Biophysical Journal. 2008;94:2580–2597. [PMC free article] [PubMed]
19. Murray D, Arbuzova A, Hangyás-Mihályné G, Gambhir A, Ben-Tal N, et al. Electrostatic properties of membranes containing acidic lipids and adsorbed basic peptides: theory and experiment. Biophysical Journal. 1999;77:3176–3188. [PubMed]
20. Murray D, Hermida-Matsumoto L, Buser C, Tsang J, Sigal C, et al. Electrostatics and the membrane association of Src: theory and experiment. Biochemistry. 1998;37:2145–2159. [PubMed]
21. Murray D, McLaughlin S, Honig B. The role of electrostatic interactions in the regulation of the membrane association of G protein beta gamma heterodimers. Journal of Biological Chemistry. 2001;276:45153–45159. [PubMed]
22. Feig M, Tanizaki S, Sayadi M. Implicit Solvent Simulations of Biomolecules in Cellular Environments. Annual Reports in Computational Chemistry. 2008;4:107–121.
23. Feig M, Im W, Brooks C., III Implicit solvation based on generalized Born theory in different dielectric environments. The Journal of Chemical Physics. 2004;120:903–911. [PubMed]
24. Im W, Feig M, Brooks C., III An implicit membrane generalized born theory for the study of structure, stability, and interactions of membrane proteins. Biophysical Journal. 2003;85:2900–2918. [PubMed]
25. Tanizaki S, Feig M. A generalized Born formalism for heterogeneous dielectric environments: Application to the implicit modeling of biological membranes. The Journal of Chemical Physics. 2005;122:124706. [PubMed]
26. Tanizaki S, Feig M. Molecular dynamics simulations of large integral membrane proteins with an implicit membrane model. The Journal of Physical Chemistry B. 2006;110:548–556. [PubMed]
27. Denning E, Woolf T. Double Bilayers and Transmembrane Gradients: A Molecular Dynamics Study of a Highly Charged Peptide. Biophysical Journal. 2008;95:3161–3173. [PubMed]
28. Zhang J, Lazaridis T. Calculating the free energy of association of transmembrane helices. Biophysical Journal. 2006;91:1710–1723. [PubMed]
29. Ulmschneider M, Ulmschneider J, Sansom M, Di Nola A. A generalized Born implicit-membrane representation compared to experimental insertion free energies. Biophysical Journal. 2007;92:2338–2349. [PubMed]
30. Nina M, Bernèche S, Roux B. Anchoring of a monotopic membrane protein: the binding of prostaglandin H 2 synthase-1 to the surface of a phospholipid bilayer. European Biophysics Journal. 2000;29:439–454. [PubMed]
31. Roux B, Berneche S, Im W. Ion channels, permeation, and electrostatics: insight into the function of KcsA. Biochemistry. 2000;39:13295–13306. [PubMed]
32. Philippsen A, Im W, Engel A, Schirmer T, Roux B, et al. Imaging the electrostatic potential of transmembrane channels: atomic probe microscopy of OmpF porin. Biophysical Journal. 2002;82:1667–1676. [PubMed]
33. Faraldo-Gómez J, Roux B, et al. Electrostatics of ion stabilization in a ClC chloride channel homologue from Escherichia coli. Journal of Molecular Biology. 2004;339:981–1000. [PubMed]
34. Chanda B, Asamoah O, Blunck R, Roux B, Bezanilla F. Gating charge displacement in voltage-gated ion channels involves limited transmembrane movement. Nature. 2005;436:852–856. [PubMed]
35. Jogini V, Roux B. Electrostatics of the intracellular vestibule of K+ channels. Journal of Molecular Biology. 2005;354:272–288. [PubMed]
36. Pathak M, Yarov-Yarovoy V, Agarwal G, Roux B, Barth P, et al. Closing in on the resting state of the Shaker K+ channel. Neuron. 2007;56:124–140. [PubMed]
37. Robertson J, Palmer L, Roux B. Long-pore electrostatics in inward-rectifier potassium channels. Journal of General Physiology. 2008;132:613–632. [PMC free article] [PubMed]
38. Grabe M, Lecar H, Jan Y, Jan L. A quantitative assessment of models for voltage-dependent gating of ion channels. Proceedings of the National Academy of Sciences. 2004;101:17640–17645. [PubMed]
39. Jiang Y, Lee A, Chen J, Cadene M, Chait B, et al. The open pore conformation of potassium channels. Nature. 2002;417:523–526. [PubMed]
40. Hilf R, Dutzler R. X-ray structure of a prokaryotic pentameric ligand-gated ion channel. Nature. 2008;452:375–379. [PubMed]
41. Chatelain F, Alagem N, Xu Q, Pancaroglu R, Reuveny E, et al. The pore helix dipole has a minor role in inward rectifier channel function. Neuron. 2005;47:833–843. [PMC free article] [PubMed]
42. Jo S, Vargyas M, Vasko-Szedlar J, Roux B, Im W. PBEQ-Solver for online visualization of electrostatic potential of biomolecules. Nucleic Acids Research. 2008;36:W270–W275. [PMC free article] [PubMed]
43. Jmol: an open-source Java viewer for chemical structures in 3D. Available: http://www.jmol.org/
44. Dolinsky T, Czodrowski P, Li H, Nielsen J, Jensen J, et al. PDB2PQR: expanding and upgrading automated preparation of biomolecular structures for molecular simulations. Nucleic Acids Research. 2007;35:W522–W525. [PMC free article] [PubMed]
45. Sitkoff D, Sharp K, Honig B. Accurate calculation of hydration free energies using macroscopic solvent models. Journal of Physical Chemistry. 1994;98:1978–1988.
46. MacKerell A, Bashford D, Bellott M, Dunbrack R, Evanseck J, et al. All-atom empirical potential for molecular modeling and dynamics studies of proteins. Journal of Physical Chemistry B. 1998;102:3586–3616. [PubMed]
47. Swanson J, Wagoner J, Baker N, McCammon J. Optimizing the Poisson dielectric boundary with explicit solvent forces and energies: Lessons learned with atom-centered dielectric functions. Journal of Chemical Theory and Computation. 2007;3:170–183. [PubMed]
48. Wang J, Cieplak P, Kollman P. How well does a restrained electrostatic potential(RESP) model perform in calculating conformational energies of organic and biological molecules? Journal of Computational Chemistry. 2000;21:1049–1074.
49. Green D, Tidor B. Evaluation of ab initio charge determination methods for use in continuum solvation calculations. The Journal of Physical Chemistry B. 2003;107:10261–10273.
50. Schutz C, Warshel A. What are the dielectric “constants” of proteins and how to validate electrostatic models? Proteins: Structure, Function, and Bioinformatics. 2001;44:400–417. [PubMed]
51. Cohen B, McAnaney T, Park E, Jan Y, Boxer S, et al. Probing protein electrostatics with a synthetic fluorescent amino acid. Science. 2002;296:1700–1703. [PubMed]
52. Silva J, Pan H, Wu D, Nekouzadeh A, Decker K, et al. A multiscale model linking ion-channel molecular dynamics and electrostatics to the cardiac action potential. Proceedings of the National Academy of Sciences. 2009;106:11102–11106. [PubMed]
53. Choe S, Hecht K, Grabe M. A Continuum Method for Determining Membrane Protein Insertion Energies and the Problem of Charged Residues. Journal of General Physiology. 2008;131:563–573. [PMC free article] [PubMed]
54. Lomize M, Lomize A, Pogozheva I, Mosberg H. OPM: orientations of proteins in membranes database. Bioinformatics. 2006;22:623–625. [PubMed]
55. Roux B. Influence of the membrane potential on the free energy of an intrinsic protein. Biophysical Journal. 1997;73:2980–2989. [PubMed]
56. Jackson J. Classical Electrodynamics, Second Edition. 1975.
57. Hille B. Ionic channels of excitable membranes. Sunderland, MA. Sinauer Associates Inc, 426pp. 1984;174:289–300.
58. Parsegian A. Energy of an ion crossing a low dielectric membrane: solutions to four relevant electrostatic problems. Nature. 1969;221:844–846. [PubMed]
59. Beckstein O, Tai K, Sansom M. Not ions alone: barriers to ion permeation in nanopores and channels. Journal of the American Chemical Society. 2004;126:14694–14695. [PubMed]
60. Roux B, MacKinnon R. The cavity and pore helices in the KcsA K+ channel: electrostatic stabilization of monovalent cations. Science. 1999;285:100–102. [PubMed]
61. Ujwal R, Cascio D, Colletier J, Faham S, Zhang J, et al. The crystal structure of mouse VDAC1 at 2.3 Å resolution reveals mechanistic insights into metabolite gating. Proceedings of the National Academy of Sciences. 2008;105:17742–17747. [PubMed]
62. Choudhary OP, Ujwal R, Kowallis W, Coalson R, Abramson J, et al. The Electrostatics of VDAC: Implications for Selectivity and Gating. Journal of Molecular Biology. 2010;396:580–592. [PMC free article] [PubMed]
63. Humphrey W, Dalke A, Schulten K. VMD – Visual Molecular Dynamics. Journal of Molecular Graphics. 1996;14:33–38. [PubMed]

Articles from PLoS ONE are provided here courtesy of Public Library of Science