Applications of MMPBSA to Membrane Proteins I: Efficient Numerical Solutions of Periodic Poisson-Boltzmann Equation

Abstract

Continuum solvent models have been widely used in biomolecular modeling applications. Recently much attention has been given to inclusion of implicit membranes into existing continuum Poisson-Boltzmann solvent models to extend their applications to membrane systems. Inclusion of an implicit membrane complicates numerical solutions of the underlining Poisson-Boltzmann equation due to the dielectric inhomogeneity on the boundary surfaces of a computation grid. This can be alleviated by the use of the periodic boundary condition, a common practice in electrostatic computations in particle simulations. The conjugate gradient and successive over-relaxation methods are relatively straightforward to be adapted to periodic calculations, but their convergence rates are quite low, limiting their applications to free energy simulations that require a large number of conformations to be processed. To accelerate convergence, the Incomplete Cholesky preconditioning and the geometric multigrid methods have been extended to incorporate periodicity for biomolecular applications. Impressive convergence behaviors were found as in the previous applications of these numerical methods to tested biomolecules and MMPBSA calculations.

Figure 1. Discretized Laplacian Operator for 3×4×5 Grid

Left: Illustration of the band structure for the finite-difference discretization of the Laplacian operator under periodic boundary conditions. Right: Illustration of the band structure for the finite-difference discretization of the Laplacian operator under fixed potential boundary conditions.

Figure 2. Differences between Electrostatic Energies by Different Boundary Conditions vs Fill Ratios

Illustration of the grid-dimension dependence in deviations of electrostatic energies between the periodic/conductor boundary condition and the free boundary condition in vacuum (upper) and in water (lower) for three model systems: low dielectric spherical cavities with a radius of Angstrom imbedded with a dipole, quadrapole, or octapole, respectively. The dipole consists of unit positive and negative charges located at positive and negative .5 Å from the cavity center along the x axis. The quadrapole consists of four unit positive and negative charges located at the vertices of a unit square in the x–y plane and centered at the cavity center. The octapole consists of eight unit positive and negative charges located at the vertices of a unit cube centered at the cavity center. The program reports the sum of the Coulomb and reaction field energies as a single electrostatic energy.

Figure 3. Differences between Electrostatic Energies by Different Solvers vs Grid Volumes

Log-Log plots of differences in electrostatic energies by PICCG (top), PMG (middle), and PSOR (bottom) with respect to PCG versus computation grid volumes for the protein test set.

Figure 4. Optimization of Scaling Coefficients for PICCG

Contour plot of the number of iterations required for convergence of the PICCG method on a small model peptide system as a function of the scaling coefficients for periodic and non-periodic bands.

Figure 5. Solver Iteration Scaling: Log Iteration Required vs Log Grid Volume

Log-log plots of iteration steps required to reach convergence versus total number of grid points (grid volume). In the case of the PMG method, “effective” iteration steps are used. This is computed as the sum of iteration steps at each grid level divided by the scaling factor of that grid level relative to the finest grid level, e.g. factors of 1, 8, 64, and 512, respectively. Top Left: PICCG, Top Right: PCG, Bottom Left: 4 Level V-Cycle PMG with SOR/Gauss-Siedel Relaxation, Bottom Right: PSOR.

Figure 6. Solver Time Scaling: Log Solver Time Required vs Log Grid Volume

Log-log plots of the solver computation time versus total number of grid points (grid volume). The solver computation time excludes other time such as boundary potential setup, energy calculation, and molecular surface generation.

Figure 7. Total PB Time Scaling: Log Computation Time vs Log Grid Volume

Log-log plots of the total computation time required versus total number of grid points (grid volume). The total computation time includes all time needed for the PB calculations to finish normally, such as boundary potential setup, energy calculation, and molecular surface generation.

Figure 8. Fractional Solver Time Scaling: Log Relative Solver Time vs Log Grid Volume

Log-log plots of the fractional solver computation time versus total number of grid points (grid volume). The fractional solver computation time is with respect to the total computation time.

Figure 9. Charge Neutrality Achieved: Log Unsigned Relative Net Charge vs Log Absolute Net Solute Charge

Log-log plot of the unsigned relative net charge of the combined solute and mobile ion charge distributions versus the net solute charge for various nucleic acid systems under periodic boundary conditions. The relative net charge is computed with respect to the absolute solute net charge.

Figure 10. Dependence of Solver Performance upon Ionic Strength: Log Solver Time vs Log Ionic Strength

Log-log plot of the solver computation time versus the ionic strength. Note that the solver time for PSOR degrades rapidly, following a seemingly exponential increase in time as the ionic strength approaches zero.

Figure 11. Effect of Implicit Solvent Modeling on Binding Affinity Calculation for P2Y12R

Correlation plots of computed and measured relative binding free energies (ΔΔG) for P2Y12R complexes.