Edited by: Zoran Nikoloski, Max-Planck Institute of Molecular Plant Physiology, Germany
Reviewed by: Zhaoyong Yang, Chinese Academy of Medical Sciences, China; Daehee Lee, Korea Research Institute of Bioscience and Biotechnology, South Korea
Specialty section: This article was submitted to Synthetic Biology, a section of the journal Frontiers in Bioengineering and Biotechnology
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
Knowledge about the electrostatic potential on the surface of biomolecules or biomembranes under physiological conditions is an important step in the attempt to characterize the physico-chemical properties of these molecules and, in particular, also their interactions with each other. Additionally, knowledge about solution electrostatics may also guide the design of molecules with specified properties. However, explicit water models come at a high computational cost, rendering them unsuitable for large design studies or for docking purposes. Implicit models with the water phase treated as a continuum require the numerical solution of the Poisson–Boltzmann equation (PBE). Here, we present a new flexible program for the numerical solution of the PBE, allowing for different geometries, and the explicit and implicit inclusion of membranes. It involves a discretization of space and the computation of the molecular surface. The PBE is solved using finite differences, the resulting set of equations is solved using a Gauss–Seidel method. It is shown for the example of the sucrose transporter ScrY that the implicit inclusion of a surrounding membrane has a strong effect also on the electrostatics within the pore region and, thus, needs to be carefully considered, e.g., in design studies on membrane proteins.
Electrostatic interactions govern the physical–chemical interactions in and between biomolecules (Perutz,
Coulombic forces are modulated by the environment (Warshel et al.,
The electrostatic contribution to free energy differences between two states of a biomolecular system – e.g., two proteins bound to each other vs. the two proteins separated in space – is in many cases difficult to directly access; however, it may be computed from ensembles of microscopic structures. The simulation-based free energy perturbation (FEP) approach or the thermodynamic integration (TI) method is frequently used to analyze free energy difference of biomolecular states. The computational cost of such methods is, however, at variance with the need for high throughput, e.g., in protein design or protein–ligand docking and in the pKa analysis of titratable sites as well. The latter problem requires the accurate estimation of the free energy differences between protonated and de-protonated states of all titratable groups in proteins, either on crystal structures or on trajectories obtained from molecular dynamics simulations to better grasp the protein flexibility (Narzi et al.,
Due to the computational efficiency required to tackle the above problems, the solvation free energy is usually addressed in implicit models by solving the Poisson–Boltzmann equation (PBE) for the different states [see, e.g., Ullmann and Bombarda (
A number of program packages for the numerical solution of the Poisson–Boltzmann equation have been developed in the past, namely APBS (Baker et al.,
The analysis of the electrostatic potential Φ of a biomolecule in an implicitly treated solvent environment requires the numerical solution of the Poisson–Boltzmann equation (PBE). For the most simple case of a solute molecule in a homogeneous medium, three different domains (Holst, Inside the molecule (green) (i.e., within the van der Waals radii of the atoms of the solute) Φ can be computed using the Poisson equation and Green’s identity resulting in
Outside the molecule (light blue), the charge density of the ions in the solvent is assumed to follow a Boltzmann distribution that yields
In the ion exclusion layer (dark blue) around the molecule no mobile ions are present. The Poisson equation reads accordingly
Combining these conditions into one equation, results in the non-linear PBE
The framework provided by equation (
The above partial differential equation cannot be solved analytically for objects shaped more complex than, e.g., a single sphere or cylinder. Therefore, the PBE has to be solved numerically.
The first step in the numerical solution of the PBE is to map all physical quantities (charges at atom centers, dielectric values, etc.) onto a three-dimensional uniform grid. Such a grid allows to replace differential operators by grid value differences. This approach is facilitated by linearizing the PBE [LPBE (Holst,
Discretization of this equation yields for every grid point
Application of equation (
Several alternatives were suggested for the treatment of the grid boundaries: setting the electrostatic potential at the boundaries to 0, application of distance-dependent quasi-Coulombic boundary conditions, or periodic extensions of the system box in one or more directions.
The set of linear equations (
The grid-based electrostatic potential is used to compute the solvation free energy of the system. This is achieved by summing the product of the potential and the charge at each grid point, followed by subtraction of the corresponding energy as obtained for the potential using a uniform dielectric inside and outside of the solute. This approach eliminates the self-energy terms that are physically not meaningful. The disadvantage of this approach is the required double computation of the electrostatic potential.
A different approach can be used if the molecular surface is known (see below). The reaction field effects due to a dielectric boundary are replaced by the computation of the induced charges at this boundary (Rocchia et al.,
Some simplifications of equation (
In each iteration step of the SOR first this stencil is used, followed by analysis of the correction for dielectric discontinuities and charges if necessary. This approach is also a first step for a parallelization as the uniform stencil in equation (
In order to define which grid points lie inside and outside of the solute, the assignment of dielectric constants to the midpoints and, in particular for the computation of the free energy in a more sophisticated way (see above), the solute surface needs to be defined. The molecular surface is determined by its implicit description from the atom positions and radii of the solute. This surface is defined as the contact surface between the van der Waals surface and the surface of a spherical probe representing the solvent.
The first step is to determine a van der Waals map by mapping the atoms onto the grid. Every midpoint is categorized as inside or outside of the solute. This is done by examining all midpoints inside a box surrounding each atom and comparing their distance to the atom center with the radius of the atom. The same approach is used to determine which grid points are in solution, which is important for knowing which grid points’ stencil has to be modified by salt. Additionally, the grid points are classified as internal points if all surrounding midpoints are inside, external points if all surrounding midpoints are outside, or boundary points if some midpoints are in solution and some are not.
Taking into account the probe radius, some of these points have to be reclassified as sketched in the red marked region in Figure
The actual surface points are finally constructed by projecting the boundary grid points onto the molecular surface either directly by moving them along the line connecting the grid point and the nearest atom center or by projecting them first onto the closest point of the solvent accessible surface and then back on the molecular surface in a similar way (Rocchia et al.,
Figure
The approach described in the Section “
The above-described input files and additional parameters like grid size, percentage of filling, probe radius, etc. are listed in one parameter file that serves as input. Table
in(tpr, |
Gromacs input file | Can also be given in the command line |
in(pdb, |
Protein positions | Alternative to tpr |
in(siz, |
Protein sizes | Necessary with pdb, optional with tpr |
in(crg, |
Protein charges | Necessary with pdb, optional with tpr |
in(sph, |
Positions in pdb format | Optional, to compute the potential |
Along a path or a specific postilions | ||
Filling | Percentage of box that is filled with protein | Default: 80% |
Spacing | Spacing (in Å) between two grid points | Default: 1 Å |
gridS | Grid size, i.e., number of grid points in each direction (odd) | Default: is computed |
Two of these three parameters | ||
Can be chosen | ||
rmsc | Convergence criterion | Default: 0.0001 |
maxc | Convergence criterion | Default: 0.0001 |
maxit | Maximal iteration | Default: 500 |
epsIn | Dielectric constant inside the molecule | Default: 2.0 |
epsOut | Dielectric constant in the solution | Default: 80.0 |
salt | Concentration in moles/liter | Default: 0.0 |
temp | Temperature | Default: 273.15 |
bc = {1,2} | Boundary conditions | 1: 0-boundary, 2: quasi-Coulombic |
pb = |
Periodic boundary conditions in |
Default: false |
nomem | No membrane input | Default: no information |
mem = zmin, zmax, eps | Minimal and maximal spread of membrane | Default: no information |
In |
||
tprmem = res, atm, eps | Residue and atoms to be considered membrane in tpr | Default: no information |
Its dielectric constant |
In the following, some of these parameters are explained in more detail.
Membranes are modeled as a slab in the If the exact values are known beforehand, the mem = zmin,zmax,eps option in the parameter file can be used. These values can also be entered in the interactive mode after the If a membrane is included in the “.tpr” file, with the tprmem = res,atm,eps option in the parameter file (or in the interactive mode) the associated residue names and atoms (e.g., POPC as residuum and P as atom) can be used to specify the membrane. The chosen atoms are separated into upper and lower parts of the membrane and the according center of mass is used as zmin and zmax.
There are three possibilities to specify boundary conditions. These have to be fixed, as the grid points at the grid’s boundary do not have enough neighbors to be computed directly. The first and easiest way is to set the boundary grid points to 0. Another possibility is to approximate the potential at the boundary using quasi-Coulombic dipole conditions:
Alternatively, periodic boundary conditions can be used in each grid dimension separately. Then for the missing neighbors, the corresponding points at the opposite side of the grid are used. This provides the opportunity to model infinite boxes.
There are three possible ways to define the convergence of the iterative procedure. The first is to set a threshold on the root mean squared change (rmsc), defined as
The implementation was tested using molecules of different sizes, ranging from small proteins with 49 residues to proteins as large as 1,070 residues (1,260–16,260 atoms, respectively). The grid size was chosen such that the filling rate was ≈80% using a grid size of 1 Å. Different results regarding the correlation between the number of atoms and grid size, number of iterations, and time consumption for different parts of the program are shown in Figures
The results show that the number of iterations and, therefore, the time to solve the linear system strongly depends on the number of grid points. The surface computation and energy calculation, however, depend on the size of the molecule.
The electrostatic potential computed using GroPBS may be mapped onto the surface of biomolecules using, e.g., PyMOL (Schrödinger,
In order to evaluate the influence of the low dielectric of a lipid membrane on the electrostatic potential of an embedded membrane protein, the potential was compared for the sucrose-specific porin ScrY (Forst et al.,
The obtained electrostatic potential through the ScrY pore (see Figure
A new program for the numerical solution of the Poisson–Boltzmann equation around biological macromolecules is presented (GroPBS). Apart from soluble proteins, GroPBS may as well be used to analyze the electrostatic potential of integral membrane proteins. The low-dielectric membrane environment may be modeled implicitly or explicitly. Additionally, GroPBS is shown to efficiently perform such computations both on pdb files and using the Gromacs input file format. This significantly simplifies the fast calculation of, e.g., the solvation free energies of biomolecules for different force fields or on ensembles of structures obtained from molecular dynamics simulations.
On the example of the sucrose-specific porin ScrY, we show that the inclusion of a membrane may have a substantial influence also on the potential inside the protein, and thus should not be neglected in PB calculations of membrane proteins.
In a subsequent step, GroPBS will be parallelized for multi-core architectures, and in particular, GPUs to enable, e.g., the fast analysis of changes in pKa values on the fly during biomolecular simulations. Parallelization may be easily achieved by the so-called
The program will be made available free of charge on the following website:
Research was designed by FB and RB, performed by FB and LS, and the manuscript was written by all.
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
We acknowledge support by the Deutsche Forschungsgemeinschaft within the Research Training Group grant No. 1962/1, Dynamic Interactions at Biological Membranes – From Single Molecules to Tissue (FB, LS, and RB). This work was also supported by the Emerging Field Initiative Synthetic Biology of the Friedrich-Alexander University of Erlangen-Nürnberg. Liping Sun was supported by the China Scholarship Council (CSC).