Accurate solution of multi-region continuum biomolecule electrostatic problems using the linearized Poisson-Boltzmann equation with curved boundary elements

Abstract

We present a boundary-element method (BEM) implementation for accurately solving problems in biomolecular electrostatics using the linearized Poisson-Boltzmann equation. Motivating this implementation is the desire to create a solver capable of precisely describing the geometries and topologies prevalent in continuum models of biological molecules. This implementation is enabled by the synthesis of four technologies developed or implemented specifically for this work. First, molecular and accessible surfaces used to describe dielectric and ion-exclusion boundaries were discretized with curved boundary elements that faithfully reproduce molecular geometries. Second, we avoided explicitly forming the dense BEM matrices and instead solved the linear systems with a preconditioned iterative method (GMRES), using a matrix compression algorithm (FFTSVD) to accelerate matrix-vector multiplication. Third, robust numerical integration methods were employed to accurately evaluate singular and near-singular integrals over the curved boundary elements. Fourth, we present a general boundary-integral approach capable of modeling an arbitrary number of embedded homogeneous dielectric regions with differing dielectric constants, possible salt treatment, and point charges. A comparison of the presented BEM implementation and standard finite-difference techniques demonstrates that for certain classes of electrostatic calculations, such as determining absolute electrostatic solvation and rigid-binding free energies, the improved convergence properties of the BEM approach can have a significant impact on computed energetics. We also demonstrate that the improved accuracy offered by the curved-element BEM is important when more sophisticated techniques, such as nonrigid-binding models, are used to compute the relative electrostatic effects of molecular modifications. In addition, we show that electrostatic calculations requiring multiple solves using the same molecular geometry, such as charge optimization or component analysis, can be computed to high accuracy using the presented BEM approach, in compute times comparable to traditional finite-difference methods.

Figure 1

A one-surface problem in molecular electrostatics. The molecular interior (Region I), containing point charges q_i, is surrounded by a salt solution with high dielectric constant and inverse Debye length k (Region II).

Figure 2

A two-surface problem in molecular electrostatics. The molecular interior (Region I), containing point charges q_i, is surrounded by an ion-exclusion layer with solvent dielectric and no salt (Region II), which in turn is surrounded by solvent with a salt treatment (Region III).

Figure 3

A three-surface problem in molecular electrostatics. This geometry is analogous to the two-surface problem (Figure 2) except that a solvent-filled cavity has been added within the molecular interior (Region IV). Note that in contrast to previous examples, the regions and surfaces have been labeled in reverse order.

Figure 4

Tree representation of a general surface problem. The example molecular geometry shown in (A) might correspond to an encounter complex between two associating proteins (Regions III_a and III_b) with point charges q_i, surrounded by a single ion-exclusion layer (Region II), which in turn is surrounded by solvent with salt (Region I). The binding partners contain several solvent filled cavities (Regions IV_a–c), and one cavity is large enough to contain a small ion-exclusion layer (Region V). The tree representation for this example multi-surface geometry is shown in (B).

Figure 5

An overview of the FFTSVD matrix compression algorithm. FFTSVD uses a multi-level octree spatial decomposition to separate element–evaluation point interactions into near- and far-field components at multiple length scales. When two cubes at the finest length scale are nearby, interactions are computed through direct integration. However, when two interacting cubes are well separated, dominant sources are projected onto a cubic grid and translated to a grid surrounding the recipient cube. The FFT is used to accelerate this translation operation. Finally, the grid potentials can be interpolated back onto the dominant responses of the element centroids. This Figure has been adapted from reference .

Figure 6

The two types of curved elements used to discretize accessible and molecular surfaces. A generalized spherical triangle (GST) (A), is a three-sided region on the surface of a sphere bounded by three circular arcs. These arcs are not necessarily geodesic arcs. Torus patches on molecular surfaces are discretized using toroidal elements (B), which are isomorphic to a rectangle.

Figure 7

A rendering of a curved element discretization for the molecular surface of the barnase–barstar protein complex. Red regions indicate convex spherical patches, green regions are re-entrant spherical patches, and blue regions are re-entrant toroidal patches. Black lines indicate the boundaries between elements. The graphic depicts an approximation to the discretized geometry used for calculation. Every GST and torus element has been approximated by a very large number of flat triangles for the purpose of visualization only, and the true surface normal in conjunction with Phong shading have been used to render the image. Surface patches with higher curved element densities were regions where finer discretization was necessary to exactly represent the molecular geometry.

Figure 8

Convergence plots for the solvation free energy for a protein-sized sphere with a near-surface charge and ion-exclusion layer. Results are compared between the curved BEM and FDM implementations in terms of absolute electrostatic solvation free energy per unit compute time (A), relative error from the analytical solution per unit compute time (B), relative error per unit computer memory (C), and relative error per unit inverse discretization length scale (where the length scale is average mesh curved edge length for BEM and distance between grid points for FDM) (D). In addition the effect of several common improvements to the FDM on the relative error is shown (E). For (C), dashed lines have been added at 1 GB, 2 GB, and 4 GB of computer memory for convenience.

Figure 9

Computed solvation free energies, using curved BEM and FDM, for an HIV-1 protease substrate peptide (A, B) and the barnase–barstar complex (C, D). The absolute electrostatic solvation free energy is plotted as a function of compute time (A, C) or computer memory usage (B, D). In (B, D) dashed lines have been added at 1 GB, 2 GB, and 4 GB of computer memory for convenience. The slight decrease in memory usage with increasing refinement observed in B is an artifact of the spatial partitioning scheme used in the FFTSVD algorithm.

Figure 10

Comparison of preconditioning strategies when solving for the electrostatic solvation free energy of an HIV-1 protease substrate peptide discretized with 18,657 and 7,089 elements on the dielectric and ion-exclusion surfaces, respectively (A) or the barnase–barstar complex discretized with 61,493 and 28,800 elements (B). In both cases, the block-diagonal preconditioner significantly reduced the number of GMRES iterations required to solve the linear system of BEM equations to a relative residual of 10⁻6.

Figure 11

Comparison between curved BEM and FDM for computing the electrostatic component of the rigid binding free energy between the wild-type barnase–barstar complex (A), and three mutant complexes, E73Q in barnase (B), D39A in barstar (C), and T42A in barstar (D). The binding free energy obtained is plotted as a function of the compute time required. In (A), several FDM and BEM results are labeled with their inverse discretization length scale in ångströms⁻1 (inverse grid spacing or inverse average curved element edge length, respectively).

Figure 12

Comparison between curved BEM and FDM for computing relative rigid electro-static binding energies between mutant and wild-type barnase–barstar complexes. Results are shown for the mutations E73Q in barnase (A), D39A in barstar (B), and T42A in barstar (C). The relative binding free energy is plotted as a function of the compute time for the mutant complex rigid binding free energy.

Figure 13

Comparison between curved BEM and FDM for computing relative non-rigid electrostatic binding energies between mutant and wild-type barnase–barstar complexes. Results are shown for the mutations E73Q in barnase (A,B), D39A in barstar (C,D), and T42A in barstar (E,F). The relative binding free energy is plotted as a function of compute time (A,C,E) or computer memory usage (B,D,F) when calculating the mutant complex non-rigid binding free energy. In (B,D,F), dashed lines have been added at 1 GB, 2 GB, and 4 GB of computer memory for convenience.