Numerical Diculties Computing Electrostatic Potentials Near Interfaces

Numerical Difficulties Computing Electrostatic. Potentials Near Interfaces with the. Poisson-Boltzmann Equation. Robert C. Harris,† Alexander H. Bos...
0 downloads 0 Views 3MB Size
Subscriber access provided by University of Newcastle, Australia

Article

Numerical Diculties Computing Electrostatic Potentials Near Interfaces with the Poisson-Boltzmann Equation Marcia Oliveira Fenley, Alexander H Boschitsch, and Robert C Harris J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00487 • Publication Date (Web): 22 Jun 2017 Downloaded from http://pubs.acs.org on June 25, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 21

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Numerical Difficulties Computing Electrostatic Potentials Near Interfaces with the Poisson-Boltzmann Equation Robert C. Harris,† Alexander H. Boschitsch,‡ and Marcia O. Fenley∗,¶ †Department of Pharmaceutical Sciences, University of Maryland School of Pharmacy, Baltimore, Maryland 21201 ‡Continuum-Dynamics Inc., 34 Lexington Avenue, Ewing, NJ 08618 ¶Institute of Molecular Biophysics, Florida State University, Tallahassee, FL, 32306 E-mail: [email protected] Phone: (850)644-7961. Fax: (850)644-7244 Abstract Many researchers compute surface maps of the electrostatic potential (ϕ) with the Poisson-Boltzmann (PB) equation to relate the structural information obtained from X-ray and NMR experiments to biomolecular functions. Here we demonstrate that the usual method of obtaining these surface maps of ϕ, by interpolating from neighboring grid points on the solution grid generated by a PB solver, generates large errors because of the large discontinuity in the dielectric constant (and thus in the normal derivative of ϕ) at the surface. The Cartesian Poisson-Boltzmann solver contains several features that reduce the numerical noise in surface maps of ϕ: First, CPB introduces additional mesh points at the Cartesian grid/surface intersections where the PB equation is solved. This procedure ensures that the solution for interior mesh points only references nodes on the interior or on the surfaces; similarly for exterior points. Second, for added points

1

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

on the surface, a second order least squares reconstruction (LSR) is implemented that analytically incorporates the discontinuities at the surface. LSR is used both during the solution phase to compute ϕ at the surface and during post-processing to obtain ϕ, induced charges, and ionic pressures. Third, it uses an adaptive grid where the finest grid cells are located near the molecular surface.

1

Introduction

Many researchers use surface maps of the electrostatic potential (ϕ) predicted by the PoissonBoltzmann (PB) equation to rationalize biomolecular functions from structural data. 1 Such maps have been used for many tasks, including locating binding sites 2–4 and classifying molecules of unknown functions. 5–7 Visualization codes, such as PyMol 8 and VMD, 9 read in solutions for ϕ on Cartesian meshes generated by PB solvers and generate surface maps by interpolating to triangulated surfaces, usually with trilinear interpolation. However, we demonstrate in the Results that maps generated with this protocol can contain significant amounts of numerical noise. Because the interior dielectric constant (ǫin ) is much smaller than the exterior dielectric constant (ǫout ) in biomolecular systems, the PB equation predicts a large discontinuity in the normal derivative of ϕ at the surface of the molecule. As a result, if one or more of the neighboring points from which an estimate of ϕ will be obtained lie within the molecule (which is almost certain to happen when interpolating to a surface), the resulting prediction can be very different than the correct solution. In particular, if the molecule is predominantly one charge (such as DNA, which is negatively charged), then the estimates of ϕ would be much larger in magnitude than the correct value. Perhaps this observation explains some of the unphysically large magnitudes of ϕ reported on surface maps in the literature. 10–13 One potential method of reducing interpolation errors is to solve for ϕ on the normal solvent-excluded (SE) molecular surface but display the resulting potential values on the solvent-accessible (SA) surface. Such a procedure should reduce the problem with interpolation because most of the points on the SA surface are not at the discontinuity 2

ACS Paragon Plus Environment

Page 2 of 21

Page 3 of 21

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

in the derivative of ϕ at the SE surface, except where the two surfaces nearly coincide (mainly at intersections and cusps). However, although this method worked well for the protein we examined, it did not work well for DNA. In short, this fix is not sufficient to guarantee accurate surface maps of ϕ. The Cartesian PB (CPB) solver 14,15 contains a number of features that reduce the errors discussed above. First, it adds additional mesh points where the solution grid intersects the molecular surface and uses a least-squares reconstruction (LSR, defined in Methods), 16 eliminating the need to interpolate from a solution grid to obtain surface maps of ϕ. Second, it does not use a fixed grid spacing, allowing it to focus more grid points in the vicinity of the molecular surface, where ϕ varies most quickly and is hardest to converge (Figure 1). Third, it builds in the discontinuity in the dielectric constant at the interface, reducing the numerical noise in its estimates of ϕ near the molecular surface.

Figure 1: a) The solution grid used by the Cartesian Poisson-Boltzmann solver around an RNA helix (RCSB Protein Data Bank 17 id 36m7). 18

3

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

2

Methods

2.1

Features of the Cartesian Poisson-Boltzmann solver

The CPB uses an octree-based partition of the domain, called here an adaptive Cartesian grid(Figure 1), where the term ‘adaptive’ refers to the use of a finer mesh spacing near the surface, where the solution gradients and gradient changes are largest. By using a coarser mesh away from the surface it reduces the number of nodes on the solution grid and thus the computational time of the calculation. Additionally, over the interior of the molecule CPB solves for the reaction field potential, rather than ϕ, so that singularities at the interior point charges are eliminated and local solution gradients are small. Mesh refinement can also be requested at select locations specified by the user to zoom in on a region of interest or, in the case of interacting molecules, to enhance mesh resolution at the binding sites. CPB presently supports a variety of surface definitions including the van der Waals, SE and SA surfaces and several Gaussian representations (Figure 2). These surfaces are defined analytically so that an accurate representation is maintained as grid spacing is reduced. 14–16 As with regular lattices, the underlying mesh structure does not conform to the molecular surface. This poses numerical challenges at the molecular surface, where the dielectric discontinuity induces a corresponding jump the in electric field along the normal direction. At points where the surface itself is smooth, the solutions on the interior and exterior sides remain smooth and analytic (i.e., can be expanded into multidimensional Taylor series) and various numerical methods are available to connect the two domains in an accurate and consistent manner. However, at surface discontinuities (e.g., cusps or sharp edges), the electric field becomes singular and the local solution is not analytic. It is possible in some instances to represent the local singular solution analytically and thus effectively separate the full solution into analytical and regular parts, the latter being represented numerically. Such an approach was adopted, for example, in Case 8 of Xia et al. 20 However, for general geometries, particularly near re-entrant corners of van der Waals surfaces (arising at the exposed

4

ACS Paragon Plus Environment

Page 4 of 21

Page 5 of 21

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

does not conform to the molecular surface; some edges of the ACG intersect the surface. The first step in the improved interpolation strategy is to identify these intersection points and add them to the collection of mesh nodes where the PB equation solution is to be developed. The nodes on the solution grid can then be divided into a set of interior nodes (Ni ) comprised of all nodes strictly inside the surface, a set of exterior nodes (Ne ) located outside the molecular surface, and the remaining nodes (Ns ) on the surface corresponding to edge/surface intersections. For each interior node, the Laplacian operator governing the reaction field potential is approximated using neighboring points that are either interior or surface nodes; specifically, no exterior nodes contribute. Similarly for the exterior points the Laplacian is developed using only exterior or surface nodes. No special measures are needed for developing the finite difference approximations at these points, 14 since the dielectric constant at all points involved in the Laplacian stencil is the same. For the nodes on the surface, a LSR is used to solve for ϕ. Since ϕ is analytic in both the interior and exterior domains one can develop separate Taylor series expansions of the interior and exterior solutions, centered on the surface node. The coefficients for each expansion correspond to ϕ, its gradient, and its second order derivatives evaluated at the surface node. Moreover, the coefficients for the interior expansion can be directly related to those for the exterior expansions using the continuity and jump conditions across the surface and the local surface curvature. For a second-order expansion in each domain (note that a second-order representation is required to consistently represent the Laplacian operator using a finite difference method), the total number of coefficients is 10 (the potential, 3 gradients, and 6 second-order derivatives). However, the governing PB equation is used to eliminate one of the second-order derivatives, thereby reducing the number of coefficients to 9 while ensuring that the reconstruction explicitly satisfies the PB equation. Having specified the Taylor series expansion in the interior and exterior domains, one can now determine the 9 unknown coefficients using a least-squares fit to the surrounding points. At least 9 such neighboring points are needed to effect a unique fit; however, generally

6

ACS Paragon Plus Environment

Page 6 of 21

Page 7 of 21

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

more points are included for redundancy and immunity against potentially degenerate or ill-conditioned node arrangements. The complete algebraic details of these relations and the numerical evaluation of the LSR have been given elsewhere. 16 Numerical tests of the LSR for smooth surfaces, such as spheres, have demonstrated consistent convergence to the analytical solutions–i.e., the solutions at the grid nodes (including surface intersection points) converge to exact results as the grid spacing is reduced. These tests have also shown the importance of accounting for surface curvature in the LSR. Additional tests involving surface discontinuities, such as a sharp edge or cusp, have confirmed that the scheme breaks down locally, but the correct behavior is recovered 2-3 grid points away from the discontinuity. As described above, this breakdown is expected since the solution is nonanalytic and a Taylor series expansion is not valid. However, the solution convergence away from the discontinuity is reassuring and can be explained on the basis that the radius of convergence of the Taylor series expansion at such locations is approximately equal to the distance to the nearest discontinuity. If the nodes used to estimate the Laplacian are within this radius, convergence can be expected. Post-processing of the solution to obtain the surface potential and gradients proceeds by first triangulating the intersection points. This is done on a cell-by-cell basis where the edge-surface intersection points for each intersected cell are identified and triangulated. The potential at each intersection node is already known, and the surface gradients are estimated by LSR. For off-surface locations, conventional trilinear interpolation can be used when the point lies inside a cell that is not intersected by the surface. For evaluation points located inside an intersected cell, the potential can be calculated by extrapolating from the nearest surface point using the second order Taylor series expansion already developed for that point during its LSR.

7

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

2.2

Page 8 of 21

Calculation details

All PB calculations used ǫin = 1, ǫout = 80, a temperature of 298 K, a 1:1 salt concentration of 0.1 M, an SE surface generated with a probe radius of 1.4 ˚ A, and no Stern layer. The structure of the nucleosome core particle (1KX5 23 (Figure 3)) was taken from the RCSB protein databank, 17 all atoms but the DNA and protein were removed, and charges and radii were added with the PDB2PQR program. 24,25 The AMBER 94 force field was used A unless otherwise stated. 26 The spherical model problem was a sphere with a radius of 5 ˚ and a central charge of 1 e, where e is the charge of a proton, centered at the origin, (0,0,0). The surface maps of ϕ from CPB 14,15 were generated from a solution mesh with a finest mesh spacing of 0.5 ˚ A and additional points added where the mesh intersected the molecular surface. 16 The equivalent maps from the adaptive Poisson-Boltzmann Solver (APBS) 27 were generated by computing the potentials at the same points in space by trilinear interpolation A grid from a ϕ grid generated by APBS with autofocusing with a fine grid with a 0.5 ˚ A on a side. The potential maps were created in spacing and a coarse grid 180 x 185 x 110 ˚ Tecplot.

3 3.1

Results and discussion Spherical model problem

The difficulty in computing surface maps of ϕ can be illustrated by considering a central point charge in a low-dielectric sphere embedded in a high-dielectric medium without any salt. In such a system ϕ is given by the Poisson rather than Poisson-Boltzmann equation, and the analytical solution to this problem is well known:

ϕ (r) =

   q/ǫ

out r

if r ≥ R

  −q (ǫout − ǫin ) /ǫout ǫin R + q/ǫin r, if r ≤ R 8

ACS Paragon Plus Environment

(1)

Page 9 of 21

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

where q is the central charge, r is the distance from the center of the sphere, and R is the radius of the sphere. A plot of ϕ given by Equation 1 as a function of r for the spherical model problem defined in the Methods is shown in Figure 4. Clearly, the estimate of ϕ obtained with linear interpolation using one point inside the sphere and one point outside the sphere will only converge slowly towards the correct potential. Figure 4 shows the estimate of ϕ that would be obtained with linear interpolation between two points chosen so as to precisely straddle the surface as a function of the distance between the points.

3.2

Electrostatic surface potential maps of nucleic acids

Surface maps of ϕ are used widely in biophysical applications. 1 Researchers use them to rationalize protein function from structural data and to locate potential binding sites. 2–4 However, as we showed in the previous section, computing ϕ near, much less at, the molecular surface can be problematic because of the large discontinuity in the dielectric constant. In particular, simple trilinear interpolation is prone to overestimating the magnitude of ϕ. The resulting overestimation is clearly visible in of Figure 5, where surface maps of ϕ generated by CPB for a nucleosomal DNA (pdbid: 1KX5) are compared to maps generated by computing the same potentials in APBS and using trilinear interpolation to obtain ϕ at the vertices on the SE surface generated in CPB. The ϕ obtained with trilinear interpolation are much larger in general. Indeed, at many points trilinear interpolation predicts unrealistic ϕ, whereas CPB’s estimates of ϕ are never larger than ∼-7 kT/e, which is much more realistic than trilinear interpolation’s largest predictions of ∼-350 kT/e. Figure 6 shows a comparison ϕ for these maps. As can be seen from this figure, CPB produces estimates of ϕ that are at least physically realistic (|ϕ| < 10). In contrast, trilinear interpolation produces estimates of |ϕ| > 100, and the correlation between these two predictions is R2 = 0.01. Because in the nonlinear PB equation the ion concentration is proportional to exp (−β ∗ q ∗ ϕ), where 10

ACS Paragon Plus Environment

Page 10 of 21

8 7 6 φ (kT/e)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

5 4 3 2 1

4.4

4.6

4.8

5

5.2

5.4

5.6

r (Å)

(a) 4.5 4 3.5 φ (kT/e)

Page 11 of 21

3 2.5 2 analytical value

1.5 1

0

0.2

0.4

0.6

0.8

1

δ (Å)

(b) Figure 4: a) The electrostatic potential (ϕ) as a function of distance (r) from the central 1 e charge in the 5 ˚ A radius sphere. b) The estimate of ϕ obtained by linear interpolation between a point located a distance (δ/2) inside the spherical surface and one located a distance (δ/2) outside the spherical surface as a function of the distance (δ) between the points for the spherical model problem. The dashed line represents the analytical value given by Equation 1.

11

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

β = 1/kT where k is Boltzmann’s constant, T is the temperature, and q is the charge of the ion, the values of ϕ produced by trilinear interpolation are clearly unphysical, compared to CPB’s predictions, which are much more reasonable. One strategy that is sometimes employed to produce better surface maps of ϕ is to display ϕ on an SA surface despite computing ϕ with an SE surface. As can be seen in Figure 5, doing this produces maps of ϕ that contain values of ϕ that are smaller in magnitude and vary more smoothly than those displayed on the SE surface. Additionally, as can be seen from Figure 6, this strategy does improve the agreement between the predictions of trilinear interpolation and CPB, as R2 = 0.27 between these two sets of predictions. However, trilinear interpolation still predicts many unrealistic values of ϕ, while CPB does not. Even when this strategy is pursued the predictions of CPB appear to be superior to those given by trilinear interpolation.

3.3

Electrostatic surface potential maps of proteins

Figures 7 and 8 are the corresponding figures to Figures 5 and 6 for the nucleosomal protein. These figures show better agreement between the results given by CPB and those given by trilinear interpolation, probably because of the lower charge density in the protein than in the DNA. The correlation coefficients for the plots in Figure 8 are R2 = 0.41 and R2 = 0.93.

4

Conclusions

Although surface maps of ϕ are widely used in research, 1 the current protocols for generating such maps produce maps with a large amount of numerical noise. Performing trilinear interpolations from neighboring grid points will yield results significantly different from those sought if one or more of the points lie in the molecular interior. Moving the points away from the SE surface by using an SA surface reduces the numerical noise, but the resulting maps can still be noisy. In particular, if the molecule is predominantly negatively or positively

12

ACS Paragon Plus Environment

Page 12 of 21

Page 13 of 21

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(a)

(b)

(c)

(d)

Figure 5: Surface maps of the electrostatic potential (ϕ) around a nucleosomal DNA (pdbid: 1KX5). Figures a. and b. show CPB’s predictions of ϕ displayed on solvent-excluded and solvent-accessible surfaces, respectively. Figures c. and d. are the same figures with trilinear interpolation’s predictions of ϕ.

13

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(a)

Page 14 of 21

(b)

Figure 6: (a) The electrostatic potential (ϕ) at the vertices of the solvent-excluded (SE) surface map of ϕ of the nucleosomal DNA computed with trilinear interpolation plotted against the same values from CPB. (b) The corresponding figure for ϕ plotted on the solventaccessible (SA) surface. In all cases ϕ was computed with an SE molecular surface.

14

ACS Paragon Plus Environment

Page 15 of 21

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(a)

(b)

(c)

(d)

Figure 7: Surface maps of the electrostatic potential (ϕ) around a nucleosomal protein (pdbid: 1KX5). Figures a. and b. show CPB’s predictions of ϕ displayed on solventexcluded and solvent-accessible surfaces, respectively. Figures c. and d. are the same figures with trilinear interpolation’s predictions of ϕ.

15

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 16 of 21

Page 17 of 21

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

charged, as DNA is, then the resulting maps will systematically overestimate the magnitude of ϕ. This observation may help explain some of the unphysically large magnitudes of ϕ on molecular surfaces quoted in the literature. 10–13 CPB has a variety of features that allow it to generate surface maps of ϕ and estimates of ϕ near the molecular surface that are superior to other PB solvers. Because it adds additional points to the solution mesh where the mesh intersects the molecular surface and a LSR, it can generate surface maps of ϕ without needing to perform any of the interpolations that produce numeric noise in the maps produced by other solvers. Because it uses an adaptive grid it can use a finer mesh spacing in the vicinity of the surface than other PB solvers before encountering memory limitations. Also, building the dielectric discontinuity into the finite-difference solution generates higher accuracy than many other PB solvers.

Acknowledgement This research was supported by the NIH grant 5R44GM073391-03.

References (1) Honig, B.; Nicholls, A. Classical electrostatics in biology and chemistry. Science 1995, 268, 1144–1149. (2) Getzoff, E. D.; Tainer, J. A.; Weiner, P. K.; Kollman, P. A.; Richardson, J. S.; Richardson, D. C. Electrostatic recognition between superoxide and copper, zinc superoxide dismutase. Nature 1983, 306, 287–290. (3) Jones, S.; Shanahan, H. P.; Berman, H. M.; Thornton, J. M. Using electrostatic potentials to predict DNA-binding sites on DNA-binding proteins. Nucl. Acids Res. 2003, 31, 7189–7198. (4) Tsuchiya, Y.; Kinoshita, K.; Nakamura, H. Structure-based prediction of DNA-binding

17

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

sites on proteins Using the empirical preference of electrostatic potential and the shape of molecular surfaces. Proteins: Struct., Funct., Bioinf. 2004, 55, 885–894. (5) Shanahan, H. P.; Garcia, M. A.; Jones, S.; Thornton, J. M. Identifying DNA-binding proteins using structural motifs and the electrostatic potential. Nucl. Acids Res. 2004, 32, 4732–4741. (6) Bradford, J. R.; Westhead, D. R. Improved prediction of protein–protein binding sites using a support vector machines approach. Bioinformatics 2005, 21, 1487–1494. (7) Lee, D.; Redfern, O.; Orengo, C. Predicting protein function from sequence and structure. Nat. Rev. Mol. Cell Biol. 2007, 8, 995–1005. (8) Schr¨odinger, LLC, The PyMOL Molecular Graphics System, Version 1.8 ; 2015. (9) Humphrey, W.; Dalke, A.; Schulten, K. VMD – Visual Molecular Dynamics. J. Mol. Graph. 1996, 14, 33–38. (10) Chin, K.; Sharp, K. A.; Honig, B.; Pyle, A. M. Calculating the electrostatic properties of RNA provides new insights into molecular interactions and function. Nat. Struct. Mol. Biol. 1999, 6, 1055–1061. (11) Draper, D. E.; Grilley, D.; Soto, A. M. Ions and RNA folding. Annu. Rev. Biophys. Biomol. Struct. 2005, 34, 221–243. (12) Stahley, M. R.; Adams, P. L.; Wang, J.; Strobel, S. A. Structural metals in the group I intron: A ribozyme with a multiple metal ion core. J. Mol. Biol. 2007, 372, 89–102. (13) Rychlik, M. P.; Chon, H.; Cerritelli, S. M.; Klimek, P.; Crouch, R. J.; Nowotny, M. Crystal structures of RNase H2 in complex with nucleic acid reveal the mechanism of RNA-DNA junction recognition and cleavage. Mol. Cell 2010, 40, 658–670. (14) Boschitsch, A. H.; Fenley, M. O. A fast and robust Poisson–Boltzmann solver based on adaptive Cartesian grids. J. Chem. Theory Comput. 2011, 7, 1524–1540. 18

ACS Paragon Plus Environment

Page 18 of 21

Page 19 of 21

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(15) Fenley, M. O.; Harris, R. C.; Mackoy, T.; Boschitsch, A. H. Features of CPB: A Poisson– Boltzmann solver that uses an adaptive Cartesian grid. J. Comput. Chem. 2015, 36, 235–243. (16) Boschitsch, A. H.; Fenley, M. O. In Computational Electrostatics for Biological Applications; Rocchia, W., Spagnuolo, M., Eds.; Springer: London, 2015; pp 73–110. (17) Berman, H. M.; Westbrook, J.; Feng, Z.; Gilliland, G.; Bhat, T. N.; Weissig, H.; Shindyalov, I. N.; Bourne, P. E. The protein data bank. Nucleic Acids Res. 2000, 28, 235–242. (18) Kiliszek, A.; Kierzek, R.; Krzyzosiak, W. J.; Rypniewski, W. Structural insights into CUG repeats containing the ‘stretched U–U wobble’: Implications for myotonic dystrophy. Nucl. Acids Res. 2009, 37, 4149–4156. (19) Kieft, J. S.; Tinoco, I. Solution structure of a metal-binding site in the major groove of RNA complexed with cobalt (III) hexammine. Structure 1997, 5, 713–721. (20) Xia, K.; Zhan, M.; Wei, G.-W. MIB method for elliptic equations with multi-material interfaces. J. Comput. Phys. 2011, 230, 4588–4615. (21) Van Bladel, J. G. In Electromagnetic fields; Dudley, D. G., Ed.; IEEE Press Series on Electromagnetic Wave Theory; John Wiley & Sons: Hoboken, NJ, 2007; Vol. 19. (22) Grant, J. A.; Pickup, B. T.; Nicholls, A. A smooth permittivity function for Poisson– Boltzmann solvation methods. J. Comput. Chem. 2001, 22, 608–640. (23) Davey, C. A.; Sargent, D. F.; Luger, K.; Maeder, A. W.; Richmond, T. J. Solvent mediated interactions in the structure of the nucleosome core particle at 1.9 ˚ A resolution. J. Mol. Biol. 2002, 319, 1097–1113.

19

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(24) Dolinsky, T. J.; Nielsen, J. E.; McCammon, J. A.; Baker, N. A. PDB2PQR: An automated pipeline for the setup of Poisson-Boltzmann electrostatics calculations. Nucleic Acids Res. 2004, 32, W665–W667. (25) Dolinsky, T. J.; Czodrowski, P.; Li, H.; Nielsen, J. E.; Jensen, J. H.; Klebe, G.; Baker, N. A. PDB2PQR: Expanding and upgrading automated preparation of biomolecular structures for molecular simulations. Nucleic Acids Res. 2007, 35, W522–W525. (26) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. A second generation force field for the simulation of proteins, nucleic acids, and organic molecules. J. Am. Chem. Soc. 1995, 117, 5179–5197. (27) Baker, N. A.; Sept, D.; Joseph, S.; Holst, M. J.; McCammon, J. A. Electrostatics of nanosystems: Application to microtubules and the ribosome. Proc. Natl. Acad. Sci. USA 2001, 98, 10037–10041.

20

ACS Paragon Plus Environment

Page 20 of 21

Page 21 of 21

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Graphical TOC Entry

21

ACS Paragon Plus Environment