Ammm:PB (APBS)

What is Poisson-Boltzmann Equation(PBE)?

Electrostatics plays a fundamental role in virtually all processes involving biomolecules in solution. The Poisson–Boltzmann equation is a second order nonlinear partial differential equation. It is one of the most fundamental approaches to treat electrostatic effects in solution. Poisson-Boltzmann equation (PBE).

In the Poisson–Boltzmann approach all solute atoms are often considered explicitly as particles with low dielectric constant (typical of organic molecules) with point partial charges at atomic positions. We will assume that the dielectric constant of the solute (typically a protein or a nucleic acid) is in the range 2–4. The solvent surrounding the solute is taken into account through its high dielectric constant(~80).

Pro & Cons:

Advantages

• Compromise between explicit and GB methods
• Reasonably fast and “accurate”
• Linear scaling
• Applicable to very large systems

Disadvantages

• Fails badly with highly-charged systems and/or high salt concentrations
• Neglects molecular details of solvent and solvation

PBE Derivation

The concentration of salt in physiological system is on the order of 150 mM, which corresponds approximately to 350 water molecules for each cation-anion pair. For this reason, investigations of salt effects in biological systems using detailed atomic models and molecular dynamic simulations become rapidly prohibitive. and meanfield treatments based on continuum electrostatics are advantageous.

Numerical Solution of the Poisson-Boltzmann Equation

Finite Difference(FD) discretization

Some of the most popular discretization techniques employ Cartesian meshes to subdivide the domain in which the PB equation is to be solved. Of these, the finite difference (FD) method has been at the forefront of PBE solvers.

Adaptive finite element (FE) discretization

Unlike finite difference methods, adaptive finite element (FE) discretizations offer the ability to place computational effort in specific regions of the problem domain. Finite element meshes (see below) are composed of simplices (e.g., triangles or tetrahedra) which are joined at edges and vertices. The solution is constructed from piecewise polynomial basis functions which are associated with mesh vertices and typically re non-zero only over a small set of neighboring simplices. Solution accuracy can be increased in specific areas by locally increasing the number of vertices through simplex refinement.

5-7 is the conformity loop. A globally conforming simplex mesh is defined as a collection of simplices which meet only at vertices and faces; for example, removing the dotted bisection in the third figure from the left in Figure 2 produces a non-conforming mesh. Non-conforming simplex meshes create several theoretical as well as practical implementation difficulties.

Finite difference and uniform mesh methods	Fast solvers Low memory overhead Cartesian mesh	Non-adaptive Poor solution resolution Previous parallel methods complicated and inefficient
Finite element methods	Highly adaptive Relatively fast solvers	Previous solver and adaptive methods inadequate Previous parallel methods complicated and inefficient

Parallel Multilevel Adaptive Finite Element Methods

Algorithm:

• The first step in the Bank-Holst parallel finite element algorithm is the solution of the equation over the entire problem domain using a coarse resolution basis set by each processor.
• This solution is then used, in conjunction with an error estimator, to partition the problem domain into P subdomains which are assigned to the P processors of a parallel computer. The algorithm is load balanced by equidistributing the error estimates across the subdomains.
• Each processor then solves the partial differential equation over the global mesh, but confines its adaptive refinement to the local subdomain and a small surrounding overlap region. This procedure results in a very accurate representation of the solution over the local subdomain.
• After all processors have completed their local adaptive solution of the equation, a global solution is constructed by the piecewise assembly of the solutions in each subdomain.

Parallel Focusing

• Given the problem data and P processors of a parallel machine:
  • Each processor i = 1, …, P:
    • Subdivides the global domain into P subdomains, each of which is assigned a processor
    • Assigns boundary conditions to a fine discretization of its subdomain using the coarse global solution
    • Solves the equation on its subdomain
• A master processor collects observable data from other processors and controls I/O

Adaptive Poisson-Boltzmann Solver (APBS)

Introduction

APBS is a software package for the numerical solution of the Poisson-Boltzmann equation (PBE) written by Nathan Baker et al. It was designed to efficiently evaluate electrostatic properties for such simulations for a wide range of length scales to enable the investigation of molecules with tens to millions of atoms.

http://apbs.sourceforge.net

What can APBS do?

• Qualitative and quantitative potential analysis
• Biomolecular energetics
  • Calculating the solvation energy
  • Calculating binding energies
  • Calculating pKa
• Biomolecular dynamics

Get Structures and Input Files Ready for Electrostatic Calculations

PQR File

PQR Format:

Field Atm_num Atm_name Res_name Chain_ID Res_numb X Y Z Charge Radius

ATOM      1  N   ASY     1      46.331  15.935  -4.837 -0.470 1.850
ATOM      2  HN  ASY     1      47.159  16.151  -4.326  0.310 1.000
ATOM      3  CA  ASY     1      45.149  15.619  -4.067  0.070 2.275
ATOM      4  HA  ASY     1      44.863  14.628  -4.402  0.090 1.320
ATOM      5  CB  ASY     1      45.552  15.559  -2.613 -0.280 2.175
ATOM      6  HB1 ASY     1      46.234  16.403  -2.388  0.090 1.320
ATOM      7  HB2 ASY     1      44.664  15.645  -1.958  0.090 1.320
ATOM      8  CG  ASY     1      46.250  14.231  -2.494  0.620 2.000
ATOM      9  OD1 ASY     1      45.579  13.233  -2.252 -0.760 1.700
ATOM     10  OD2 ASY     1      47.447  14.211  -2.761 -0.760 1.700
ATOM     11  C   ASY     1      43.945  16.516  -4.272  0.510 2.000
ATOM     12  O   ASY     1      43.470  16.441  -5.403 -0.510 1.700
ATOM     13  CAY ASY     1      47.542  16.266  -6.896 -0.270 2.175
ATOM     14  HY1 ASY     1      48.344  16.474  -6.156  0.090 1.320
ATOM     15  HY2 ASY     1      47.382  17.165  -7.529  0.090 1.320
ATOM     16  HY3 ASY     1      47.844  15.411  -7.539  0.090 1.320
ATOM     17  CY  ASY     1      46.279  15.929  -6.181  0.510 2.000
ATOM     18  OY  ASY     1      45.252  15.661  -6.802 -0.510 1.700
ATOM     19  N   GLU     2      43.440  17.391  -3.372 -0.470 1.850
ATOM     20  HN  GLU     2      43.901  17.558  -2.505  0.310 1.000
ATOM     21  CA  GLU     2      42.201  18.170  -3.552  0.070 2.275
ATOM     22  HA  GLU     2      42.016  18.756  -2.661  0.090 1.320
ATOM     23  CB  GLU     2      42.430  19.133  -4.770 -0.180 2.175
ATOM     24  HB1 GLU     2      42.615  18.509  -5.675  0.090 1.320
ATOM     25  HB2 GLU     2      41.499  19.705  -4.979  0.090 1.320
ATOM     26  CG  GLU     2      43.608  20.160  -4.678 -0.280 2.175
ATOM     27  HG1 GLU     2      43.773  20.622  -5.669  0.090 1.320
...

$ ${PQR2PDB_PREFIX}/pdb2pqr [options] --ff=<forcefield> <path> <output-path>

<forcefield>  : The forcefield to use - currently amber, charmm, parse, and tyl06 are supported.

<path>  : The path to the PDB file or an ID to obtain from the PDB archive.

<output-path> : The desired output name of the PQR file to be generate.

APBS Input

read 
    mol pqr mol1.pqr
    mol pqr mol2.pqr
    mol pqr complex.pqr
end

#CALCULATE POTENTIAL FOR COMPONENT 1

elec name mol1
    mg-para  # Automatic parallel focusing multigrid calculation
    ofrac 0.1 # Specify the amount of overlap to include between the individual processors meshes 
    pdime 2 2 2 # specifies the 8-processor array dimensions
    dime  97  97  97 # Specifies the number of grid points per processor 
    fglen 112 91 116 # Specifies the fine mesh domain lengths in a focusing calculation 
    cglen 156 121 162 # Specifies the coarse mesh domain lengths in a focusing calculation 
    fgcent mol 3 # Specify the center of the fine grid 
    cgcent mol 3 # Specify the center of the coarse grid 
    mol 1 # Center the grid on molecule with integer ID id
    npbe # non-linear PBE
    bcfl sdh # boundary conditions: zero, sdh, mdh, focus
    ion charge 1 conc 0.050 radius 2.0
    ion charge -1 conc 0.050 radius 2.0
    pdie 2.0 # the dielectric constant of the biomolecule
    sdie 78.54 # the dielectric constant of solvent
    srfm mol # Specify the model used to construct the dielectric and ion-accessibility coefficients
    chgm spl0 # Specify the method by which the biomolecular point charges
    srad 1.4 # Specify the radius of the solvent molecules
    swin 0.3 # Specify the size of the support (i.e., the rate of change) for spline-based surface definitions
    sdens 10.0
    temp 298.15
    calcenergy total
    calcforce no
end

#CALCULATE POTENTIAL FOR COMPONENT 2

elec name mol2
    mg-para
    ofrac 0.1
...skipping...
    calcforce no
end

#CALCULATE POTENTIAL FOR COMPLEX

elec name complex
    mg-para
    ofrac 0.1
    pdime 2 2 2
...skipping...
    write pot dx pot  # write out OpenDX-format maps of the potential to 8 files pot#.dx
end

#COMBINE TO GIVE BINDING ENERGY

print energy complex - mol1 - mol2 end

Useful Keywords of ELEC

$ ${MPI_PREFIX}/bin/mpirun -n 8 ${APBS_PREFIX}/apbs inputfile &

$ mergedx 65 65 65 pot0.dx pot1.dx pot2.dx pot3.dx pot4.dx pot5.dx pot6.dx pot7.dx

APBS Output

  # Comments
  object 1 class gridpositions counts nx ny nz
origin xmin ymin zmin
delta hx 0.0 0.0
  delta 0.0 hy 0.0 
  delta 0.0 0.0 hz
object 2 class gridconnections counts nx ny nz
object 3 class array type double rank 0 items n data follows
u(0,0,0) u(0,0,1) u(0,0,2)
  ...
 u(0,0,nz-3) u(0,0,nz-2) u(0,0,nz-1)
 u(0,1,0) u(0,1,1) u(0,1,2)
 ...
 u(0,1,nz-3) u(0,1,nz-2) u(0,1,nz-1)
 ...
 u(0,ny-1,nz-3) u(0,ny-1,nz-2) u(0,ny-1,nz-1)
 u(1,0,0) u(1,0,1) u(1,0,2)
 ...
 attribute "dep" string "positions"
object "regular positions regular connections" class field
component "positions" value 1
component "connections" value 2
component "data" value 3

  # Data from APBS 1.0.0
  # 
  # mergedx
  # 
  object 1 class gridpositions counts 65 65 65
  origin -5.348200e+01 -4.796450e+01 -4.125750e+01
  delta 1.750000e+00 0.000000e+00 0.000000e+00
  delta 0.000000e+00 1.421875e+00 0.000000e+00
  delta 0.000000e+00 0.000000e+00 1.812500e+00
  object 2 class gridconnections counts 65 65 65
  object 3 class array type double rank 0 items 274625 data follows
  -1.509577e-02 -1.669933e-02 -1.848610e-02 
  -2.047263e-02 -2.266871e-02 -2.511585e-02 
  -2.779038e-02 -3.077484e-02 -3.403626e-02 
  -3.763402e-02 -4.158453e-02 -4.591748e-02 
  -5.066451e-02 -5.585899e-02 -6.153566e-02 
  ...skipping...
  -1.619205e-02 -1.463064e-02 -1.322073e-02 
  -1.192318e-02 -1.075077e-02 -9.683742e-03 
  -8.717906e-03 -7.845471e-03 
  attribute "dep" string "positions"
  object "regular positions regular connections" class field
  component "positions" value 1
  component "connections" value 2
  component "data" value 3

Visualization in VMD

Surface potentials

• Load structure, surface mesh, or solvent surface density grid into a molecule in VMD. (e.g. 'mol new structure.pqr' or load the PQR file into a new molecule, using the graphical interface)
• Display the structure with your favorite representation. (e.g. set the drawing method to MSMS, Surf, or if you've loaded an APBS solvent accessibility map, you can use the Isosurface representation instead)
• Add the potential grid (or other volumetric data you wish to color the structure with) into the SAME molecule. (e.g. 'mol addfile pot.dx')
• Once both the structure and potential are loaded into the same molecule, the next step is to set the structure coloring method to "Volume". The volume coloring method requires that you select one of the loaded volumetric datasets (from within the same molecule). If you only have on volumetric dataset loaded, it will be selected for you automatically.
• At this point, you will see the coloring of the surface change, but the the coloring range may not be set to the appriate range for the data, so we must adjust the "color scale data range" to be able to see the potential values we're interested in more clearly. The color scale data range settings are located in the "trajectory" tab of the graphical representations window. A good initial range for APBS potential data is -10 to 10. If you prefer more intense colors, you can set the range lower, to something like -4 to 4. Potential values outside of the minima and maxima defined in the color scale data range will be clamped to solid red and solid blue (in the case of the RWB color scale).
• The last item is to change the active vmd color scale, to something appropriate for viewing potential data, such as the "RWB" color scale which displays low or (in this case negative) values in red, and high or positive values in blue. The color scale is set in the Color menu, in the "color scale tab". In new versions of VMD, the default color scale is already set to "RWB", so users of new versions may not feel the need to change the color scale.
  

Isocontour visualization

One of the most popular visualization methods is the isocontour.

• First, choose Graphics → Representations... from the VMD Main window.
• Hit the Create Rep button and change Drawing Method to "Isosurface".
• Change Draw from "Points" to "Solid Surface" and Material to "Transparent".
• Note that the current isovalue is 0; change this to 1 (or 5 or 10 or whatever) depending on your personal preferences.
• To continue the lonstanding tradition of electrostatic potential coloring, choose "ColorID 0" for the Coloring Method.
• For the negative isocontour, hit "Create Rep" and select the newly created representation. Change the isocontour value to -1 (or 5 or 10...) and the ColorID to 1.
  

 

Reference