PBEQ-Solver for online visualization of electrostatic potential of biomolecules.

PBEQ-Solver provides a web-based graphical user interface to read biomolecular structures, solve the Poisson-Boltzmann (PB) equations and interactively visualize the electrostatic potential. PBEQ-Solver calculates (i) electrostatic potential and solvation free energy, (ii) protein-protein (DNA or RNA) electrostatic interaction energy and (iii) pKa of a selected titratable residue. All the calculations can be performed in both aqueous solvent and membrane environments (with a cylindrical pore in the case of membrane). PBEQ-Solver uses the PBEQ module in the biomolecular simulation program CHARMM to solve the finite-difference PB equation of molecules specified by users. Users can interactively inspect the calculated electrostatic potential on the solvent-accessible surface as well as iso-electrostatic potential contours using a novel online visualization tool based on MarvinSpace molecular visualization software, a Java applet integrated within CHARMM-GUI (http://www.charmm-gui.org). To reduce the computational time on the server, and to increase the efficiency in visualization, all the PB calculations are performed with coarse grid spacing (1.5 A before and 1 A after focusing). PBEQ-Solver suggests various physical parameters for PB calculations and users can modify them if necessary. PBEQ-Solver is available at http://www.charmm-gui.org/input/pbeqsolver.

INTRODUCTION

Implicit solvent treatments are approximate methods that attempt to incorporate the average influence of the molecular environment on a system of interest without having to explicitly simulate the molecules constituting this environment (1). Implicit solvent methods have emerged as a popular strategy to approximate bulk solvent or membrane environments and have been applied successfully to protein–protein or protein–ligand binding thermodynamics, scoring of protein conformations in structure prediction, peptide and protein folding/unfolding studies and ion channels (1–4). In particular, Poisson-Boltzmann (PB) continuum electrostatics, in which the solvent is represented as a featureless dielectric material, is the most rigorous and popular method to estimate the electrostatic solvation energy of a solute with an arbitrary shape, and particular successes in applications to complex biological problems are evident (2,5–7). The characterization of the electrostatic potential on the macromolecular surface by solving the PB equation is becoming a routine practice in structural biology (5).

Over the last two decades, considerable efforts have been made to generalize and enhance the computational methodologies and techniques to solve the PB equation and visualize the calculated electrostatic potential (7–13). As a result, various user-friendly programs that provide numerical solutions of the PB equation using finite-difference or finite-element methods with a discretized grid are now available as standalone software such as APBS (7), MEAD (14), Qnifft (15), DelPhi (16) and Zap (17), or as modules in biomolecular modeling and simulation programs such as PB solver (18) in AMBER (19), PBEQ (12,20–22) in CHARMM (23), APBS (7) in AMBER (7), CHARMM (23), TINKER (24) and NAMD (25), and PB solver (26) in Jaguar (http://www.schrodinger.com). PyMOL (27), VMD (28), GRASP (29), PMV (30) and DINO (http://www.dino3d.org) (31) provide various visualization tools for calculated electrostatic potentials. There are also several useful web-based interfaces to setup and perform PB calculations such as PDB2PQR (32), PCE (33) and PDB|Hydro (34).

We have developed PBEQ-Solver (http://www.charmm-gui.org/input/pbeqsolver) to provide a web-based graphical user interface (GUI) to read biomolecular protein data bank (PDB) structures (35) through PDB Reader at the CHARMM-GUI website (http://www.charmm-gui.org), solve the PB equations using PBEQ module (12,20–22) in CHARMM (23) and interactively visualize the electrostatic potential on the solvent-accessible surface as well as iso-electrostatic potential contours using the MarvinSpace molecular visualization software (http://www.chemaxon.com/product/mspace.html), a Java applet integrated within CHARMM-GUI. In addition to the calculations of electrostatic potential and solvation free energy, PBEQ-Solver also computes protein–protein (DNA or RNA) electrostatic interaction energy and of a selected titratable residue. All the calculations can be performed in both aqueous solvent and membrane environments (with a cylindrical pore in the case of membrane).

METHODS AND ILLUSTRATIVE RESULTS

PB theory

The electrostatic solvation energy ΔGelec is the work required to assemble the charges {q} of the solute in the solvent (1), which can be expressed in terms of the reaction field potential , Based on continuum electrostatics, the reaction field potential, can be computed by solving the PB equation twice for the reference electrostatic potential and the electrostatic potential in the solvent environment , where , and are the dielectric constant, the modified Debye–Hückel screening factor and the fixed charge density of the solute, respectively. For a given dielectric constant in the interior of the solute (), is calculated by setting the dielectric constant to (dielectric constant of reference environment) at all points in the solvent region, while is calculated by setting the dielectric constant to (solvent dielectric constant) and (from input concentration) in the solvent region. The influence of biological membranes and cylindrical pore can be incorporated into by setting the dielectric constants to (membrane dielectric constant) and (pore dielectric constant) in the membrane and pore regions, respectively.

In the case of a biomolecular complex containing proteins, DNA and RNA, the electrostatic interaction energy in solvent between a selected chain [A] and the rest [B] in the complex can be computed as where , and are the electrostatic potentials of the complex, [A] and [B] in solvent or membrane environments, respectively. The of a residue in a protein, , is estimated based on calculating, for both the protonated (p) and unprotonated (u) states of the residue, the difference between its electrostatic free energy when it is in the protein environment () and its electrostatic free energy when it is isolated in solution (). The shift () relative to the of the same amino acid isolated in solution () is given by where

PDB Reader

Since PBEQ-Solver uses the PBEQ module in CHARMM, its first step is to read a PDB file into CHARMM, which is generally considered not to be straightforward due to complexity of PDB files. In general, this is also the first difficulty that any simulation program users may face. It is typically even harder to introduce different protonation states of titratable residues, disulfide bonds or other posttranslational modifications such as phosphorylation. For a seamless, efficient procedure, PBEQ-Solver first uses PDB Reader in CHARMM-GUI that provides a flexible web-based GUI to convert a PDB file [downloaded from RCSB (35), http://www.rcsb.org, or uploaded from user's computer] into CHARMM readable files with the following options; (i) partial selection of protein chains as well as model selection in the case of NMR structures, (ii) modification of engineered residues, (iii) terminal group selection, (iv) protonation selection, (v) disulfide bond selection, (vi) phosphorylation selection, (vii) generation of a biologically functional unit and (viii) generation of a crystal packing. The terminal patch residues available in PDB Reader are listed in Supplementary Table S1. If indicated in a PDB file, PDB Reader automatically detects the disulfide bonds and displays them. Users can always add or remove them in the list. PDB Reader also automatically detects some of engineered residues listed in Supplementary Table S2 and converts them to corresponding natural residues. For example, an engineered residue ‘Sep’ in a PDB file represents a phosphorylated Ser residue. In such case, PDB Reader automatically converts Sep to Ser and turns on phosphorylation of the Ser residue (Supplementary Figure S3). However, due to the complexity of PDB files with various heteroatoms and other engineered residues, it is users that need to check if all the residues to be read do exist in a CHARMM topology file (currently, top|all27|prot|na.rtf). If there are undetermined coordinates in selected chains, PDB Reader simply builds them using a predetermined internal coordinate table (‘IC BUILD’ command in CHARMM). The usages of PDB Reader is well illustrated in a video demo, ‘PDB Reader Tutorial’, available at the CHARMM-GUI website (http://www.charmm-gui.org/?doc=demo).

It should be stressed that RCSB PDB structures do not contain orientation information of a membrane(-bound) protein relative to lipid bilayers. Therefore, in order to perform reasonable PB calculations in a membrane environment, users need to validate if the protein structure is properly oriented with respect to membranes. It is assumed that the membrane normal is parallel to the Z-axis and its center by default is located at Z = 0 (Table 1). Users can preorient the structure in their local machine and upload it, or use preoriented protein structures from the OPM database (http://opm.phar.umich.edu) (36) by selecting OPM as PDB download source (Supplementary Figure S1).

Table 1.

Physical parameters in PBEQ-Solver

Variables	Default values	Physical meaning
EpsP	1.0	Dielectric constant for the solute interior
EpsR	1.0	Dielectric constant for the reference environment
EpsW	80.0	Solvent dielectric constant
Conc	0.15	Salt concentration (in M)
Focus	Yes	Focusing option
Dcel|c	1.5	Coarse grid spacing
Dcel|f	1.0	Finer grid spacing (for focusing)
Ledgea	10.0	Minimum distance between solute and grid boundary
Tmemb	35.0	Thickness of membrane (along the Z-axis)
Zmemb	0.0	Center of membrane (along the Z-axis)
EpsM	2.0	Membrane dielectric constant
Htmembb	0.0	Thickness of headgroup region
EpsH	2.0	Membrane headgroup dielectric constant
Rcyln	0.0	Radius of cylindrical pore
Hcyln	35.0	Height of cylindrical pore
EpsC	80.0	Dielectric constant of cylindrical pore
Xcyln	0.0	Position of cylindrical pore in X
Ycyln	0.0	Position of cylindrical pore in Y
Zcyln	0.0	Position of cylindrical pore in Z
ctom	No	Set the dielectric constant of the overlapped region with membrane to EpsM
ckap	Yes	Make cylinder pore accessible to ions

aLEdge*2 is set to the minimum distance for coarse-gird calculations and LEdge/2 for finer grid calculations.

bThe head group region is defined within the membrane thickness (Tmemb). For example, if Tmemb = 35 Å, Zmemb = 0 Å, and Htmemb = 2.5 Å, EpsM is assigned in −15 Å < Z < 15 Å and EpsH in −17.5 < Z < −15.0 and 15.0 < Z < 17.5.

PBEQ-Solver

After PDB Reader, PBEQ-Solver displays the options for three types of PB calculations as well as adjustable physical parameters summarized in Table 1 (see also Supplementary Figure S4). Note that, to reduce the computational time on the server and to increase the efficiency in visualization, all the PB calculations are presently performed with (unadjustable) coarse grid spacing (Dcel|c = 1.5 Å before and Dcel|f = 1 Å after focusing (37)). All the PB calculations uses the molecular surface (with a probe radius of 1.4 Å) to setup the dielectric boundary based on a set of atomic radii for proteins, DNA and RNA that were optimized to closely reproduce the charging free energies calculated by molecular dynamics free energy simulations for standard amino acids (20) and nucleic acids (22) (filename: step2|radii.str, which is downloadable as shown in Supplementary Figure S5).

PBEQ-Solver currently offers three types of PB calculations for (i) electrostatic potential and solvation free energy [based on Equation 1], (ii) protein–protein (DNA or RNA) electrostatic interaction energy [based on Equation 3] and (iii) pKa of a selected titratable residue [based on Equation 4]. As mentioned before, all the calculations can be performed in both aqueous solvent and membrane environments (with a cylindrical pore in the case of membrane). Based on the user's inputs for the physical parameters listed in Table 1, PBEQ-Solver uses the PBEQ module in CHARMM to solve the finite-difference PB equation of molecules specified by users through PDB Reader. Note that a PDB file must contain at least two chains for the electrostatic interaction calculations, such that the interaction energy is calculated between a selected chain and the rest (Supplementary Figure S8). In the case that only one chain is selected in PDB Reader, the interaction energy option is not available.

Tests and illustrations

We have tested PBEQ-Solver with 360 PDB structures whose number of residues range from about 50 to 4000 to examine the computational time to calculate the electrostatic solvation free energy (see Supplementary Table S3 for the full list of PDB IDs). Each test took about 1–14 min on the server depending on the size of the biomolecule. In addition, to validate if the CHARMM input used in PBEQ-Solver is reliable, we have compared the PB solvation energies calculated based on the PBEQ-Solver input with those from other PB programs that are published by Feig et al. (38). As shown in Supplementary Table S4, the PBEQ-Solver input yields similar results to other PB programs when the same parameters are used. This also demonstrates that the visualization of calculated potentials with the input and the coarse grid in PBEQ-Solver is also reliable.

We have illustrated how to use PBEQ-Solver with an example of PDB:1KDX (complex of KIX and phosphorylated KID domains) in Supplementary Material (Figures S1–S11), including all the snapshots with some useful annotations for the three types of PB calculations. Figures 1 and 2 illustrate various molecular graphics that the users can generate using PBEQ-Solver. As shown in Figure 1, the protein complex (PDB:1KDX) of the KIX and phosphorylated KID domains that play an important role in regulation by posttranslational modification (39) can be viewed with its ribbon representation, solvent-accessible surface representation with electrostatic potential and iso-electrostatic potential contours. Note that the unit for the bounds on electrostatic potential visualization is kcal/(mol·e), where e is the unit charge. The iso-contour map helps users to examine how the electrostatic potentials are distributed in a distal place (other than on the surface of the biomolecule). In Figure 2, the protein–DNA complex (PDB:1E3M) of DNA mismatch repair protein MutS and G–T mismatch DNA (40) is shown in its surface electrostatic potentials with and without DNA. It is clearly shown that the DNA binding site has the high positive electrostatic surface for negatively charged DNA.

Figure 1.

Molecular graphics views of the protein complex (PDB:1KDX) of the KIX domain and one of its co-activators, the phosphorylated KID domain, that play an important role in regulation by posttranslational modification (39). PBEQ-Solver provides a tool for online visualization of its (left) ribbon representation as well as (middle) solvent-accessible surface representation with electrostatic potential [+2 kcal/(mol·e) in blue to −2 kcal/(mol·e) in red] and (right) iso-electrostatic potential contours [+1 kcal/(mol·e) in blue and −1 kcal/(mol·e) in red].

Figure 2.

Molecular graphics views of the protein–DNA complex (PDB:1E3M) of DNA mismatch repair protein MutS and G–T mismatch DNA (40) with the surface electrostatic potentials (left) with and (right) without DNA [+2 kcal/(mol·e) in blue to −2 kcal/(mol·e) in red]. To generate the surface potential without DNA, one has to unselect DNA during the PDB reading step.

The web interfaces for PBEQ-Solver and PDB Reader in CHARMM-GUI have been developed following de facto web standards and extensively tested under Firefox (version 2) and Safari (version 2). Although all the functionalities work just fine with different web browsers, we have so far found some minor inconsistency of user interfaces in Internet Explorer (IE) and the beta version of Safari (version 3). In the case of IE versions 6 and 7, some user interface elements look slightly different than those in Firefox and Safari. In the case of beta version of Safari (version 3), the (job) monitoring panel does not work properly.

CONCLUDING DISCUSSION & FUTURE DIRECTIONS

We have described the functionalities of PBEQ-Solver at the CHARMM-GUI website with illustrations. There are several programs available for visualization of electrostatic potential of biomolecules, calculated by solving the PB equations, so that users can use them in their local machines or generate images through web-based visualization applications. However, we believe that PBEQ-Solver provides unique and interactive web-based GUI for online visualization. We recently found that Honig and coworkers have made similar efforts through the Mark-US server, a functional annotation server (http://luna.bioc.columbia.edu/honiglab/mark-us). In particular, such a GUI development will help non-expert users, especially experimentalist, to perform various PB calculations and visualize calculated electrostatic potentials of their own systems. In a similar way, such an effort will be also useful for educational purposes.

The development of PBEQ-Solver is an ongoing project. Its present drawback is to use relatively coarse grid spacing [1.5 Å before and 1.0 Å after focusing (37)] for efficient PB calculations and visualization. However, it is recommended to use 1–1.5 Å before and at least 0.5 Å after focusing in practical applications for solvation energy calculations. One way to overcome this problem is to download all the CHARMM input files and redo the calculations with finer grid spacing in users’ local machines. We also expect to provide the finer grid calculations as the capability of our server increases. The calculation of transmembrane potential using the modified PB equation (21) and its visualization will be incorporated into PBEQ-Solver.

SUPPLEMENTARY DATA

Supplementary Data are available at NAR Online.