A boundary-integral approach for the Poisson-Boltzmann equation with polarizable force fields

A boundary-integral approach for the Poisson-Boltzmann equation with polarizable force fields

Christopher D. Cooper Departmento de Ingeniería Mecánica and Centro Científico Tecnológico de Valparaíso (CCTVal), Universidad Técnica Federico Santa María, Valparaíso, Chile

Implicit-solvent models are widely used to study the electrostatics in dissolved biomolecules, which are parameterized using force fields. Standard force fields treat the charge distribution with point charges, however, other force fields have emerged which offer a more realistic description by considering polarizability. In this work, we present the implementation of the polarizable and multipolar force field AMOEBA, in the boundary integral Poisson-Boltzmann solver PyGBe. Previous work from other researchers coupled AMOEBA with the finite-difference solver APBS, and found difficulties to effectively transfer the multipolar charge description to the mesh. A boundary integral formulation treats the charge distribution analytically, overlooking such limitations. We present verification and validation results of our software, compare it with the implementation on APBS, and assess the efficiency of AMOEBA and classical point-charge force fields in a Poisson-Botlzmann solver. We found that a boundary integral approach performs similarly to a volumetric method on CPU, however, it presents an important speedup when ported to the GPU. Moreover, with a boundary element method, the mesh density to correctly resolve the electrostatic potential is the same for stardard point-charge and multipolar force fields. Finally, we saw that polarizability plays an important role to consider cooperative effects, for example, in binding energy calculations.

Keywords: Poisson-Boltzmann, Implicit solvent, Polarizable force fields, Boundary element method, Electrostatics  

Here, we model the electrostatics of biomolecular systems using a continuum approach, while describing the charge distribution inside the molecule with point multipoles that polarize. In particular, we parameterize the biomolecule with the AMOEBA force field and solve the electrostatic equations with a boundary integral formulation, which integrates the charge distribution analytically. The implementation is validated, shows good behavior as the size of the biomolecule increases, and is tested for binding energy calculations.


Implicit-solvent models dramatically reduce the problem size in biomolecular simulations. All of the solvent degrees of freedom are averaged in a continuum dielectric description, whereas the solute is accounted for as a cavity that contains a charge distribution obtained from a force field. These models are long dated and heavily used by the biophysics/biochemistry community in its various versions,1, 2 beginning from Kirkwood’s closed expressions for a single spherical cavity,3 to Generalized Born models,4 to Poisson and Poisson-Boltzmann solvers.5, 6 Among the latter, there are numerical packages based on finite difference,7, 8, 9 finite element,10 boundary element,11, 12, 13, 14, 15 analytical,16 and semi-analytical17 methods, which are usually coupled to point-charge-based force fields.

The accuracy of molecular simulations is greatly influenced by the quality of force fields. The most widely used force fields, termed classical force fields, describe the charge distribution as a set of point charges, and have proven to be successful in a large number of applications.18 Nevertheless, recently there has been an important development of more elaborate force fields that consider not only the monopole, but also higher order multipole components which may polarize. This multipole polarizable description yields a more realistic charge distribution, and hence, improved accuracy in the resulting simulations.18, 19, 20, 21, 22 This, however, comes at a higher computational cost.23 In order to overcome such time limitations, several researchers have coupled polarizable force fields with implicit solvent models. For example, there are extensions to Kirkwood’s solution24 and Generalized Born models25 to account for polarizable multipoles. Moreover, we can find several efforts towards including polarizability to Poisson-Boltzmann26, 27, 28, 20, 29 and Poisson30, 31, 32 solvers.

Recently, there has been an increasing amount of research developing polarizable force fields.33 Among the most popular ones is AMOEBA (atomic multipole optimized energetics for biomolecular applications),34, 35 which is available in several packages, such as Tinker,36, 37, 38 Force Field X,111https://ffx.biochem.uiowa.edu/ OpenMM,39 and Amber.40 The AMOEBA force field uses point multipoles up to the quadrupole moment to describe the charge distribution, and allows for the dipole moment to have a polarizable component. A notable effort in coupling AMOEBA with a Poisson-Boltzmann model is the work by Schnieders and co-workers,27 where the authors extended the widely used APBS solver10 (using a finite difference method) to account for polarizable multipoles. In that work AMOEBA was successfully integrated into APBS, however, some limitations may be found in transferring the charge distribution to the finite difference mesh, due to the high order multipole components (dipole and quadrupole). The number of mesh nodes to map a point multipole adequately increases with the order of the multipole, to the point that a quadrupole needs 5 evenly-spaced points per dimension. Moreover, a very fine mesh is required for multipoles near the dielectric interface, to avoid placing charges outside the cavity that represents the biomolecule.

In this work, we present a boundary-integral Poisson-Boltzmann solver, compatible with the AMOEBA force field. We extended the boundary element method (BEM) code PyGBe,15 which uses the formulation presented by Yoon and Lehoff,41 to account for polarizable multipoles. In a boundary integral framework, the charge distribution is integrated analytically,32 avoiding the limitations encountered in a finite-difference model. Moreover, the BEM stiffness matrix is not affected by the charge description, making it realtively easy to port into existing boundary-element codes.


The boundary integral formulation in the implicit solvent model

The implicit-solvent model divides the domain in a protein () and a solvent () region, interfaced by the solvent excluded surface (, SES), as sketched in Figure 2. In the protein region, the dielectric constant is low and there is a charge distribution which is parameterized using force fields. The solvent is usually water (), with salt. Continuum electrostatic theory leads to a coupled system to solve for the potential (), where the Poisson equation models the protein region and the linearized Poisson-Boltzmann equation governs in the solvent, with appropriate interface conditions.


Here, the subscripts and refer to the protein and solvent regions, respectively, is the inverse of the Debye length, and is a unit normal vector to the SES pointing out of the protein, into the solvent.

If we apply Green’s second identity on Equation (The boundary integral formulation in the implicit solvent model), we get


where is the potential in region computed at the surface (). and are the double- and single-layer potentials, evaluated in :


with the free space Green’s function of the Laplace or linearized Poisson-Boltzmann (Yukawa) equations:


Evaluating Equation (The boundary integral formulation in the implicit solvent model) at the interface gives


Equation (The boundary integral formulation in the implicit solvent model) is the formulation presented by Yoon and Lenhoff.41

The implicit-solvent model with polarizable point multipoles

To consider a polarizable point-multipole description of the charge distribution inside the protein, the model in Equation (The boundary integral formulation in the implicit solvent model) needs to be extended. More specifically, the multipolar description has to be taken into account in the right-hand side of Equation (The boundary integral formulation in the implicit solvent model), and in the calculation of the electrostatic energy. Furthermore, the dipole moment in AMOEBA has an induced component, which can be solved with a self-consistent approach,42 where the electric field and induced dipole are iteratively computed until convergence.

Computation of the right-hand side

The right hand side of the Poisson equation in Equation (The boundary integral formulation in the implicit solvent model) is the electrostatic potential due to the charge distribution inside the protein, which is a collection of point charges in a classical force field. In general, the potential due to a charge distribution is


Considering a point multipole centered at the origin, up to the quadrupole moment, Equation (6) becomes


making the right-hand side of Equation (The boundary integral formulation in the implicit solvent model):


where is the monopole (or total charge), the dipole vector, and the quadrupole tensor of the point multipole. There are other formulations where the in the quadrupole term is absorbed into .43

Induced dipole

The point multipoles of the AMOEBA force field have an induced dipole component. That is, the force field assigns values of permanent monopole (), dipole (), and quadrupole () moments, and a polarizability (), such that there is an induced dipole ()


where is the electric field at the location of the point multipole. The total dipole moment in Equation (8) is then


Calculation of the electric field

To compute the induced dipole moment in Equation (9), it is easier to treat the electrostatic potential (hence, the electric field) as two separate components: one from the reaction of the solvent (), and from the multipoles ():


We can obtain by substracting out the influence of the multipoles to Equation (The boundary integral formulation in the implicit solvent model), which leaves:


The component of the electric field is then


where the derivatives of the Green’s function are detailed in Equation (APPENDIX).

On the other hand, the electric field due to multipoles at the location of multipole is the gradient of the Coulombic-type electrostatic potential from Equation (7), taken at the evaluation point. This is


The detailed expressions for the derivatives of each term are shown in Equation (APPENDIX).

AMOEBA uses a group-based polarization scheme for permanent multipoles44 where, in this case, multipoles of the same group do not polarize each other. Then, we need to consider a masking rule in Equation (Calculation of the electric field) that zeroes out the contribution of permanent multipole when it is in the same polarization group as multipole . Moreover, the field generated by induced dipoles ( in Equation (10)) is damped with a Thole-like scheme.45, 34

Having both components of the electric field from equations (Calculation of the electric field) and (13), we can compute the total electric field at the location of the point multipoles as


which goes into Equation (9) to obtain the induced dipole moment.

Self-consistent induced field

The dipole moment is an input to the implicit-solvent model, however, the electric field has an influence in the induced dipole component. Then, we need to solve for the induced dipole with a self-consistent scheme, which is summarized below:

  1. Guess (for example, ).

  2. Compute and on with Equation (The boundary integral formulation in the implicit solvent model).

  3. Find and with Equation (13) and Equation (Calculation of the electric field).

  4. Calculate the total electric field with (15).

  5. Get with Equation (9).

  6. Go back to step 2 with the new , until convergence.

We update the induced dipole with a successive-over relaxation (SOR) scheme, using a coefficient of . After reaching convergence of the induced dipole, we can move on to compute the solvation energy.

Calculation of the solvation energy

The solvation free energy () is the energy required to dissolve a molecule. In other words, it is the difference in free energy between the molecule in vacuum state and in the solvent, plus the energy spent polarizing the multipoles:


Assuming a linear dielectric, we can calculate the free energy as


where is the potential and the charge distribution. Then, we can write the solvation energy as


Decomposing the potential in dissolved state into solvent and Coulombic components, as in Equation (11), Equation (18) becomes


In standard non-polarizable force fields, the charge distribution does not have an induced component, hence and , and only the first term in the right-hand side survives. However, in the case of polarizable force-fields, like AMOEBA, the induced dipole is different in the solvated and vacuum states, and the second and third terms on the right hand side of Equation (19) do not cancel out.

In the context of this work, the charge distribution is a collection of point multipoles, and the integral becomes a sum. Then, for the case of point multipoles, the energy in Equation (17) becomes


If and corresponds to the total dipole in dissolved state, Equation (20) gives the solvent contribution to energy (first term in Equation (19)). Then, if and , Equation (20) yields the coulombic energy due to the point multipoles in dissolved state (second term in Equation (19)). Finally, if and , one obtains the coulombic energy in vacuum (second-to-last term in Equation (19)). In Equation (20), becomes when the formulation considers the term of Equation (8) inside .

We already derived the expressions for in Equation (12) and in Equation (13). The second derivative of that multiplies the quadrupole in Equation (20) is


where the required derivatives of the Green’s function are detailed in Equation (APPENDIX).

The Coulomb potential from a collection of point multipoles is shown in Equation (7) and its first derivative in Equation (Calculation of the electric field). To compute the energy in Equation (20), the second derivative of this potential is also required for the quadrupole component, and its terms are explicited in Equation (APPENDIX).

The last term in Equation (17) is known as the polarization energy, and corresponds to the amount of energy required to polarize the multipoles. Böttcher 46 derived an expression for through the charging process of a spherical cavity by a single polarizable dipole, arriving to:


where is the induced dipole component and the total electric field that polarizes the multipole. In that derivation, the dipole was polarized starting from , however, in our case, there are several multipoles that are already polarized in the original (vacuum) state. Then, the total polarization energy is


where the superscript denotes the vacuum state. Comparing Equation (Calculation of the solvation energy) with Equation (20), we realize that the polarization energy cancels out with the induced dipole (as ), and we can replace with to obtain


If we use Equation (24), rather than (20), the polarization energy is already being considered, and there is no need to explicitly calculate it.

The permanent multipoles for the calculation of and and their derivatives are not masked with AMOEBA’s group-based scheme in the energy calculation. However, the induced dipole component is damped by the same Thole-like model as the electric field calculation in Equation (Calculation of the electric field).

Calculation of the binding energy

We can compute the electrostatic component of the energy of binding between two dissolved biomolecules according to the thermodynamic cycle detailed in Figure 3. For biomolecules ‘A’ and ‘B’ complexed in ‘AB’, the binding energy in dissolved state is


where is the solvation energy of molecule ‘A’, ‘B’, or the complex ‘AB’, computed with Equation (17). On the other hand, is the binding energy in vacuum, which is the energetic difference between bound and unbound states, plus the energy required to polarize the multipoles:


Here, is the coulombic energy in vacuum state for each case, computed with Equation (20). Note that in the case of standard point-charge force fields, there is no polarization energy (), and the coulombic energy is the same in vacuum and dissolved states. To compute , we can use Equation (24) to consider the polarization energy implicitly.

From the thermodynamic cycle in Figure 3 we can also compute the relative binding energy, given by the difference in binding energy between the dissolved and vacuum states:


Algorithmic and computational details

The boundary element method.

We use a boundary element method (BEM) to solve the system in Equation (The boundary integral formulation in the implicit solvent model) numerically. In it, we discretize the molecular surface in flat triangular panels, assume a piecewise constant distribution of the potential and its derivative, and use collocation to generate a linear system of equations, which we solve with a generalized minimum residual (GMRES) method.47

Integrals are calculated with Gaussian quadrature, however, due to the nature of the Green’s function, there are singular integrals that are difficult to solve numerically. For this reason, depending on the distance between the collocation node and the integration panel, we distiguish three integration regions: singular, near-singular, and far-away. The integral becomes singular when the collocation point is inside the integration panel, and we use a semi-analytical approach48, 49 that places Gauss nodes on the edges of the element. If the collocation point and the integration panel are close-by, the integrand is nearly singular, and we use high order quadrature rules. For example, we commonly place a threshold away from the collocation point (where is the area of the element containing the collocation point) and use Gauss nodes to compute the integrals within that distance. Finally, beyond this threshold the integrand is smooth enough that we can use a low order approximation with , , or Gauss nodes.

The treecode algorithm

The most time consuming part of the solver is a matrix-vector product (an process) inside the GMRES, done once in every iteration. As we are using collocation and Gaussian quadrature with free-space Green’s functions, the matrix-vector product can be seen as a N-body problem where the sources of mass are the Gauss nodes, and we evaluate the Green’s function at the collocation points. There are several ways to accelerate N-body calculations to , and even , for example, FFT-based methods,50, 51 fast multipole method,52 and treecode.53 In this work, we use the treecode algorithm, which has been revised for the Green’s functions of the Laplace54 and linearized Poisson-Boltzmann equations.55

The treecode clusters sources (in this case, Gauss nodes) and targets (collocation points) in boxes of an octree structure, where the boxes in the lowest level of the tree contain less than a critical number of particles. Then, if a box is far enough from a target, the influence of the sources in that box is approximated using a Taylor expansion of order around the center of the box. If the box is not far enough, the algorithm checks for child boxes and performs the same operation, until a lowest-level box is reached, and the source-target interaction is performed directly. To determine if a box is far enough from a target, we use the multipole-acceptance criterion (MAC), which is:


where is the size of the box and the box-target distance. This reduces the computational complexity from to , allowing us to control accuracy with the order of the Taylor expansion () and the value of in the multipole-acceptance criterion. Further details on the implementation of the treecode in the boundary element method can be found elsewhere.14, 56

The most time consuming parts of the algorithm are the source-target and box-target interactions, which are completely target independent. For this reason, the treecode is highly parallelizable, and maps very well to the architecture of the GPU. PyGBe uses GPU acceleration, via PyCUDA, for the computation of the source-target and box-target interactions, and also for the computation of the right-hand side in Equation (8) and the energy in Equation (24).


We implemented the method described above in a forked Github repository222https://github.com/barbagroup/pygbe333https://github.com/cdcooper84/pygbe of the open-source code PyGBe.15, 57 The results presented in this section were obtained on a workstation with two Intel Xeon E5-2680v3 CPUs and one NVIDIA Tesla K40 GPU. The portions of the code that ran on CPU were serial, whereas the GPU performed parallel computations. Protein structures for AMOEBA were prepared using pdbxyz from the Tinker package, and parameterized with amoebapro0418 or amoebapro13,58 whereas for point-charge force fields, we used pdb2pqr.59 From the molecular structure and atomic radii, we generated meshes with the msms software.60

Verification and validation

Spherical cavity

Kong and Ponder24 derived a closed expression for spherical cavities with a random distribution of multipoles, equivalent to the classic Kirkwood solution3 for point charges, and suggested using self-consistent iterations to account for polarizability. We implemented the latter method to verify our numerical implementation in PyGBe. The test case was a spherical cavity of Å, with a Å-thick Stern layer, dielectric constant , and three polarizable multipoles with unit charge, dipole, quadrupole, and polarizability, immersed in water () with salt (Å). One of the multipoles was placed in the center of the sphere and the other two were off-centered by Å, as sketched in Figure 4. The solvation energy of such system, computed with Kong and Ponder’s expression, is kcal/mol. Figure 5 shows the discretization error of the boundary element solution for different mesh densities, which is, as expected, decaying with the average area of the boundary elements. Both of the surfaces in Figure 4 have the same amount of mesh elements for this simulation, and the axis in Figure 5 corresponds to the total number of them. This result verifies the implementation of PyGBe’s extension to work with polarizable multipoles.

Comparison with Tinker-APBS

The work by Schnieders and co-workers 27 validates the implementation of polarizable multipoles in APBS using 1CRN,61 1ENH,62 1FSV,63 1PGB,64 and 1VII,65 and comparing their implicit-solvent results with explicit solvent calculations. Figure 6 shows the solvation energy for the same proteins, computed with APBS amd PyGBe, for various grids — in particular, we used a spacing of and for APBS and a mesh with , and vertices per square angstrom in PyGBe. APBS is a volumetric solver, whereas PyGBe’s mesh runs only on the molecular surface, hence, mesh sizes are not comparable. To overcome this situation, we performed Richardson extrapolation66, 15 and obtained an approximation of the exact solution from the numerical calculations, as the mesh density tends to infinity. For APBS runs, we used meshes with and to compute the extrapolated values, whereas in the case of PyGBe, we used grids with , , and vertices per square angstrom. Table 1 presents the extrapolated values for each case, and they are also marked with a dotted line in Figure 6. Using the extrapolated values as an exact solution, we computed the error in solvation energy for each mesh, and plotted them against the time to solution in Figure 7. Moreover, from Richardson extrapolation we are also able to compute an observed order of convergence, which is the rate of convergence of simulations towards the extrapolated value, and they are also included in Table 1.

We further compared APBS and PyGBe using the total dipole moment of the protein. In our simulations, we saw that the mesh density had a very weak effect on the total dipole moment, hence, Table 2 only shows results for the coarsest meshes in each case ( in APBS and 2 vertices per square angstrom in PyGBe).

For consistency with Schnieders’ work,27 these runs were performed using the amoebapro04 force field, a protein dielectric region of and a solvent with dielectric constant and mM of salt (Å). The tolerance of the self-consistent solver for the induced dipole moment was set to and the exit criterion of the GMRES in the Poisson-Boltzmann linear solver to .

In the APBS simulations, we used a sharp surface definition (keyword SMOL) rather than the smoothed definition based on fourth order splines used by Schnieders and co-workers27 (keyword SPL4), which explains the differences in the solvation energy between the latter work and Figure 6. The SMOL definition is closer to the surface description from a boundary integral formulation.

For the PyGBe simulations, we used 1 Gauss quadrature point per boundary element far-away from the collocation point, however, for nearly-singular integrals (within Å  of the collocation point) we used finer quadrature rules with 19 Gauss nodes. With respect to the treecode acceleration, we set the multipole acceptance criterion to , and used Taylor expansions up to order . For efficiency, the tree was built making sure that no box of the lowest level had more than elements for CPU and elements for GPU.

Influence of size

From Figure 7, we can see that a mesh density with 4 vertices per square angstroms yields a solution that is around 1% away from a converged value. To study the behavior of PyGBe with respect to the size of the protein, we computed the solvation energy of 1PGB,64 1LYZ,67 1A7M,68 1X1U,69 and 1IGT70 using 4 vertices per square angstrom, and report them in Table 3. These runs were done using the amoebapro1358 force field, and with the same parameters of the GPU runs that led to Figure 6.

Comparison with standard force fields

The following simulations compare AMOEBA (with amoebapro13) with standard point-charge force fields, such as AMBER, CHARMM, and PARSE.

Mesh refinement study

We performed a further mesh convergence study of the solvation energy of protein GB1 (1PGB) using AMOEBA, and the fixed-charge force field AMBER, with the same simulation parameters that led to the results in Figure 6. Figure 8 shows the solvation energy with both force fields decaying approximately with the average area, which is the expected behavior for a piecewise constant boundary element solution.15 In particular, the observed order con convergence for AMOEBA is and for AMBER it is , whereas the Richardson extrapolated values are kcal/mol and kcal/mol, respectively. Also, Figure 8 shows the time to solution for different meshes using AMOEBA and AMBER, with the expected scaling from the treecode, which is the dominant part of the algorithm.

Binding energy calculations

We computed the electrostatic component of the absolute and relative binding energies for the HIV-1 GP120 core complexed with CD4 and a neutralizing human antibody (PDB: 1GC171). CD4 is a glycoprotein present on the surface of immune cells, and induces a binding site for the monoclonal antibody on GP120. Here, we compute the binding energy of GP120 (that has CD4 already attached to it, totalling atoms) with the antigen-binding fragment (Fab) of the antibody ( atoms), which generates a complex with atoms.

Table 4 shows the results for binding energy calculations computed with AMOEBA (with and without polarizability), AMBER, CHARMM, and PARSE, using Equation (25) and Equation (27). The mesh density for these runs was vertices per square angstrom, and we used the same simulation parameters that led to the results in Figure 6, only differing by setting the multipole acceptance criterion to , and the threshold between near and far-away integrals at . For the case of point-charge force fields, we varied the dielectric constant of the protein () from to , to test if it was possible to implicitly capture polarization in those cases.


Verification and validation

The verification result for the sphere in Figure 5 shows the numerical value from PyGBe approaching the analytical solution from Kong and Ponder24 with the average area of the surface mesh elements, which is the expected behavior for a piecewise constant boundary element approximation.15 This indicates that the implementation of polarizable multipoles on PyGBe is correctly solving the mathematical model. With that, we further validate PyGBe by comparing our numerical results with the implementation from Schnieders and co-workers27 in Figure 6 and Table 1. For well converged simulations, the observed order of convergence is expected to match the order of the method: for APBS and for PyGBe,15 however, this is not the case in Table 1, specially for APBS. This is probably due to the difficult transfer of high order multipole charges to the finite-difference mesh. In PyGBe, slight deviations to the expected slope are forseeable due to the irregular nature of the surface mesh generated by msms, as it is impossible to refine it homogeneously. Regardless of these shortcomings, the extrapolated values of solvation energy for APBS and PyGBe are in good agreement (less than 1% off), proving that the main features of the electrostatic potential are correctly represented in both cases.

Figure 7 shows how time scales with the error in the simulation. The expected behavior is faster simulations with APBS at low accuracy that scale worse than PyGBe as the mesh refines,15 yet, both codes present similar slopes. This unexpected trend is an artifact of having a high observed order of convergence with APBS, but the plots are still useful to perform a fair comparison of the volumetric and boundary integral solvers. All CPU runs are single-core, hence, APBS and PyGBe-CPU timings in Figure 7 are comparable, yet, it is unfair to compare those with GPU run times. From these results, we can conclude that PyGBe and AMOEBA have equivalent performance under the same CPU conditions.

We still plotted GPU timings for reference, considering that the extension of PyGBe to use AMOEBA is implemented in both CPU and GPU. From here, we can see an advantage of the boundary integral approach, as it performs well on graphics cards.

Influence of size

Table 3 shows that the size of the biomolecule has a weak effect on the number of self-consistent iterations required for the induced dipoles to converge. Between 1PGB and 1IGT the number of atoms increases by , however, only one extra iteraton was required. This suggests that this implicit-solvent Poisson-Boltzmann approach is efficient to analyze large biomolecular systems. Moreover, the solvent-excluded surface mesh grows slower than a volumetric mesh as the molecule size increases, indicating that a boundary integral approach is more appropriate for a large-scale application.

Comparison with standard force fields

Mesh refinement study

Figure 8 is a further mesh convergence study of 1PGB, comparing AMOEBA with AMBER. Results with both force fields are converging as expected (observed order of convergence close to ), proving that the electrostatic potential is correctly resolved regardless of the charge description. This means that the dipole and quadrupole components in AMOEBA do not affect the mesh density required to resolve the potential field appropriately, and the same mesh sizes can be used in both cases.

The Richardson extrapolated value of solvation energy with AMBER is kcal/mol and with AMOEBA is kcal/mol, differring by around %. There are a number of reasons that could explain this difference, for example, the fact that the parameterization of the protein might not be optimal. Also, in point-charge force fields, slightly higher ad-hoc permittivities are usually used inside the protein72 to account for reorentation of dipoles, and the solvation energy with AMBER decreases with a higher dielectric constant, approaching the value for AMOEBA.

Binding energy calculations

Even though the relative differences in solvation energy may not be very large, they can be critical when studying changes of this quantity, for example, in binding energy calculations. We find evidence of this in Table 4, where the solvation energies of the complex with AMBER and AMOEBA differ by less than , however, the absolute binding energy is off by .

Table 4 shows smaller binding energies with AMOEBA, compared to its point-charge counterparts using . This behaviour is expected. In fact, in point-charge force fields polarization is implicitly considered with higher permittivities, which, as shown in Table 4, results in a lower absolute binding energy. This effect is already considered in a polarizable force field with a relative permittivity inside the protein of .

Regardless of these differences, Table 4 consistently shows that the electrostatic component of the absolute binding energy in dissolved state () is positive, hence, repulsive. Furthermore, the values of binding energy are similar, which indicates that the parameterization with all force fields is adequate. From here, we can conclude that the nonpolar component of the solvation energy is key to correctly predict the expected binding of the antibody to the CD4-induced binding cite of GP120.71

Even though all force fields arrive to similar values of binding energy, it is interesting to analyze each component of the thermodynamic cycle in Figure 3 separately. For the case of AMOEBA, we see a negative relative binding energy (), which means that binding is more likely to happen in dissolved state than vacuum (see Equation (27)). Hence, the role of the solvent is to enhance attraction, and the repulsion is given by the energetic difference in vacuum. On the other hand, all point charge force fields predict the opposite physical phenomena: the attraction of the biomolecules is due to coulombic-type interactions (), and the solvent induces repulsion (). Hence, though we are predicting correct absolute binding energies, we may be misled in what is the role of each term of the thermodynamic cycle.

Also, in Table 4, the relative binding energy of AMOEBA drops in half when polarizability is not considered. This shows that an important part of the binding energy comes from polarization, and it is critical to consider it in applications such as this one, where the two biomolecules polarize each other.


This paper presents an extension of the Poisson-Boltzmann solver PyGBe to use a charge distribution with high order multipoles and polarizability, in particular, through the AMOEBA force field. The software is based on a boundary integral representation of the partial differential equations, where the molecular charge is considered on the right-hand side of the resulting linear system, making it relatively easy to implement in an existing code base. We verified this extension against closed expressions valid for spherical inclusions, and validated it by contrasting with a similar implementation that uses APBS. PyGBe and APBS perform equivalently in serial CPU runs, however, the boundary integral approach presents an important speedup on the GPU. We also compared the performance of AMOEBA and standard point-charge force fields, like AMBER, in a Poisson-Boltzmann model. In that comparison, we realized that the same mesh densities can appropriately resolve the electrostatic potential in either case, and that they converge to similar results. Also, considering polarizability can be extremely important in situations that have cooperative effects, such as binding, where the two molecules are mutually polarized. In that case, we saw that even though an appropriate parameterization of a force field may yield the correct values of the absolute binding energy, there are differences in the mechanisms present in the interaction, where for point-charge force fields the solvent induces repulsion, whereas for AMOEBA it favors binding.

We conclude that this boundary-integral implicit-solvent approach can efficiently compute the electrostatic potential in biomolecular systems with polarizable force fields. The fact that the charge distribution is computed analytically on the right-hand side avoids any difficult point-multipole transfer to the mesh, making a boundary-integral representation ideal for these computations, specially looking at large-scale simulations.


The author is grateful of the helpful interactions with Prof. Robert Krasny (UMich), Prof. Weihua Geng (SMU), and Prof. Michael Schnieders (UIowa). Part of this work was done while visiting the University of Michigan, hosted by Prof. Robert Krasny. The research was financially supported by CONICYT though FONDECYT Iniciación N 11160768, and Basal Project FB 0821.


This appendix shows the detailed expressions of the derivatives needed in different calculations.

The first derivatives of the Green’s function required in Equation (13) to obtain the electric field due to the solvent reaction and Equation (20) for solvation energy are


where the normal moves with the integration variable . On the other hand, the second derivatives used in Equation (21) are


The derivatives of the terms in Equation (Calculation of the electric field) are:


The terms required to compute the second derivative of the potential due to a collection of multipoles are



  • Roux and Simonson 1999 B. Roux and T. Simonson, Biophys. Chem. 78, 1 (1999).
  • Bardhan 2012 J. P. Bardhan, Comput. Sci. Disc. 5 (2012).
  • Kirkwood 1934 J. G. Kirkwood, J. Chem. Phys. 2, 351 (1934).
  • Still et al. 1990 W. Still, A. Tempczyk, R. C. Hawley, and T. F. Hendrickson, J. Am. Chem. Soc. 112, 6127 (1990).
  • Fogolari et al. 2002 F. Fogolari, A. Brigo, and H. Molinari, Journal of Molecular Recognition 15, 377 (2002).
  • Baker 2004 N. A. Baker, Meth. Enzymol. 383, 94 (2004).
  • Gilson et al. 1987 M. K. Gilson, K. A. Sharp, and B. Honig, J. Comput. Chem. 9, 327 (1987).
  • Bashford and Gerwert 1992 D. Bashford and K. Gerwert, J. Mol. Biol. 224, 473 (1992).
  • Geng et al. 2007 W. H. Geng, S. N. Yu, and G. W. Wei, J. Chem. Phys. 127, 114106 (2007).
  • Baker et al. 2001 N. A. Baker, D. Sept, M. J. Holst, and J. A. McCammon, P. Natl. Acad. Sci. USA 98, 10037 (2001).
  • Lu et al. 2006 B. Lu, X. Cheng, J. Huang, and J. A. McCammon, P. Natl. Acad. Sci. USA 103, 19314 (2006).
  • Altman et al. 2009 M. D. Altman, J. P. Bardhan, J. K. White, and B. Tidor, J. Comput. Chem. 30, 132 (2009).
  • Bajaj et al. 2011 C. Bajaj, S. C. Chen, and A. Rand, SIAM J. Sci. Comput. 33, 826 (2011).
  • Geng and Krasny 2013 W. H. Geng and R. Krasny, J. Comput. Phys. 247, 62 (2013).
  • Cooper et al. 2014 C. D. Cooper, J. P. Bardhan, and L. A. Barba, Comput. Phys. Commun. 185, 720 (2014), preprint on arXiv:/1309.4018.
  • Felberg et al. 2017 L. E. Felberg, D. H. Brookes, E.-H. Yap, E. Jurrus, N. A. Baker, and T. Head-Gordon, Journal of computational chemistry 38, 1275 (2017).
  • Yap and Head-Gordon 2010 E.-H. Yap and T. Head-Gordon, Journal of chemical theory and computation 6, 2214 (2010).
  • Ponder and Case 2003 J. W. Ponder and D. A. Case, in Advances in protein chemistry (Elsevier, 2003), vol. 66, pp. 27–85.
  • Cieplak et al. 2009 P. Cieplak, F.-Y. Dupradeau, Y. Duan, and J. Wang, Journal of Physics: Condensed Matter 21, 333102 (2009).
  • Tong et al. 2010 Y. Tong, Y. Mei, Y. L. Li, C. G. Ji, and J. Z. Zhang, Journal of the American Chemical Society 132, 5137 (2010).
  • Cisneros et al. 2013 G. A. Cisneros, M. Karttunen, P. Ren, and C. Sagui, Chemical reviews 114, 779 (2013).
  • Demerdash et al. 2014 O. Demerdash, E.-H. Yap, and T. Head-Gordon, Annual review of physical chemistry 65, 149 (2014).
  • Lipparini et al. 2014 F. Lipparini, L. Lagardère, B. Stamm, E. Cancès, M. Schnieders, P. Ren, Y. Maday, and J.-P. Piquemal, Journal of Chemical Theory and Computation 10, 1638 (2014).
  • Kong and Ponder 1997 Y. Kong and J. W. Ponder, The Journal of chemical physics 107, 481 (1997).
  • Schnieders and Ponder 2007 M. J. Schnieders and J. W. Ponder, Journal of chemical theory and computation 3, 2083 (2007).
  • Maple et al. 2005 J. R. Maple, Y. Cao, W. Damm, T. A. Halgren, G. A. Kaminski, L. Y. Zhang, and R. A. Friesner, Journal of Chemical Theory and Computation 1, 694 (2005).
  • Schnieders et al. 2007 M. J. Schnieders, N. A. Baker, P. Ren, and J. W. Ponder, The Journal of chemical physics 126, 124114 (2007).
  • Tan et al. 2008 Y.-H. Tan, C. Tan, J. Wang, and R. Luo, The Journal of Physical Chemistry B 112, 7675 (2008).
  • Aleksandrov et al. 2018 A. Aleksandrov, F.-Y. Lin, B. Roux, and A. D. MacKerell Jr, Journal of computational chemistry (2018).
  • Li and Gordon 2007 H. Li and M. S. Gordon, The Journal of chemical physics 126, 124112 (2007).
  • Lipparini and Barone 2011 F. Lipparini and V. Barone, Journal of chemical theory and computation 7, 3711 (2011).
  • Lipparini et al. 2015 F. Lipparini, L. Lagardère, C. Raynaud, B. Stamm, E. Cancès, B. Mennucci, M. Schnieders, P. Ren, Y. Maday, and J.-P. Piquemal, Journal of chemical theory and computation 11, 623 (2015).
  • Shi et al. 2015 Y. Shi, P. Ren, M. Schnieders, and J.-P. Piquemal, Polarizable force fields for biomolecular modeling, vol. 28 (John Wiley & Sons, Inc.: Hoboken, NJ, 2015).
  • Ren and Ponder 2003 P. Ren and J. W. Ponder, The Journal of Physical Chemistry B 107, 5933 (2003).
  • Ponder et al. 2010 J. W. Ponder, C. Wu, P. Ren, V. S. Pande, J. D. Chodera, M. J. Schnieders, I. Haque, D. L. Mobley, D. S. Lambrecht, R. A. DiStasio Jr, et al., The journal of physical chemistry B 114, 2549 (2010).
  • Ponder and Richards 1987 J. Ponder and F. Richards, J. Comput. Chem 8, 1016 (1987).
  • Kundrot et al. 1991 C. E. Kundrot, J. W. Ponder, and F. M. Richards, Journal of computational chemistry 12, 402 (1991).
  • Pappu et al. 1998 R. V. Pappu, R. K. Hart, and J. W. Ponder, The Journal of Physical Chemistry B 102, 9725 (1998).
  • Friedrichs et al. 2009 M. S. Friedrichs, P. Eastman, V. Vaidyanathan, M. Houston, S. Legrand, A. L. Beberg, D. L. Ensign, C. M. Bruns, and V. S. Pande, Journal of computational chemistry 30, 864 (2009).
  • Case et al. 2005 D. A. Case, T. E. C. III, T. Darden, H. Gohlke, R. Luo, K. M. Merz Jr., A. Onufriev, C. Simmerling, B. Wang, and R. J. Woods, J. Comput. Chem. 26, 1668 (2005).
  • Yoon and Lenhoff 1990 B. J. Yoon and A. M. Lenhoff, J. Comput. Chem. 11, 1080 (1990).
  • Cramer and Truhlar 1999 C. J. Cramer and D. G. Truhlar, Chemical Reviews 99, 2161 (1999).
  • Stone 2013 A. J. Stone, The Theory of Intermolecular Forces (Oxford University Press, 2013), 2nd ed.
  • Ren and Ponder 2002 P. Ren and J. W. Ponder, Journal of computational chemistry 23, 1497 (2002).
  • Thole 1981 B. T. Thole, Chemical Physics 59, 341 (1981).
  • Böttcher 1973 C. J. F. Böttcher, Theory of electric polarization (Elsevier, 1973).
  • Saad and Schultz 1986 Y. Saad and M. Schultz, SIAM J. Sci. Stat. Comput. 7, 856 (1986).
  • Hess and Smith 1967 J. L. Hess and A. Smith, Progress in Aerospace Sciences 8, 1 (1967).
  • Zhu et al. 2001 Z. Zhu, J. Huang, B. Song, and J. White, in Proceedings of the 2001 IEEE/ACM Int. Conf. on Computer-Aided Design (2001), pp. 592–597.
  • Phillips and White 1997 J. R. Phillips and J. K. White, IEEE Trans. Comput. Aid. D. 16, 1059 (1997).
  • Altman et al. 2006 M. D. Altman, J. P. Bardhan, B. Tidor, and J. K. White, IEEE Trans. Comput. Aid. D. 25, 274 (2006).
  • Greengard and Rokhlin 1987 L. Greengard and V. Rokhlin, J. Comput. Phys. 73, 325 (1987).
  • Barnes and Hut 1986 J. Barnes and P. Hut, Nature 324, 446 (1986).
  • Lindsay and Krasny 2001 K. Lindsay and R. Krasny, J. Comput. Phys. 172, 879 (2001).
  • Li et al. 2009 P. Li, H. Johnston, and R. Krasny, J. Comput. Phys. 228, 3858 (2009).
  • Cooper and Barba 2013 C. D. Cooper and L. A. Barba, Validation of the PyGBe code for Poisson-Boltzmann equation with boundary element methods, Technical Report on figshare, CC-BY license (2013).
  • Cooper et al. 2016 C. D. Cooper, N. C. Clementi, G. Forsyth, and L. A. Barba, The Journal of Open Source Software 1 (2016).
  • Shi et al. 2013 Y. Shi, Z. Xia, J. Zhang, R. Best, C. Wu, J. W. Ponder, and P. Ren, Journal of chemical theory and computation 9, 4046 (2013).
  • Dolinsky et al. 2004 T. J. Dolinsky, J. E. Nielsen, J. A. McCammon, and N. A. Baker, Nucleic Acids Research 32, W665 (2004).
  • Sanner 1996 M. F. Sanner (1996), http://mgltools.scripps.edu/packages/MSMS (last checked Feb. 4, 2010).
  • Teeter 1984 M. Teeter, Proceedings of the National Academy of Sciences 81, 6014 (1984).
  • Clarke et al. 1994 N. D. Clarke, C. R. Kissinger, J. Desjarlais, G. L. Gilliland, and C. O. Pabo, Protein Science 3, 1779 (1994).
  • Dahiyat and Mayo 1997 B. I. Dahiyat and S. L. Mayo, Science 278, 82 (1997).
  • Gallagher et al. 1994 T. Gallagher, P. Alexander, P. Bryan, and G. L. Gilliland, Biochemistry 33, 4721 (1994).
  • McKnight et al. 1997 C. J. McKnight, P. T. Matsudaira, and P. S. Kim, Nature structural biology 4, 180 (1997).
  • Roache 1998 P. J. Roache, Verification and validation in computational science and engineering (Hermosa Albuquerque, 1998).
  • Diamond 1974 R. Diamond, Journal of molecular biology 82, 371 (1974).
  • Hinds et al. 1998 M. G. Hinds, T. Maurer, J.-G. Zhang, N. A. Nicola, and R. S. Norton, Journal of Biological Chemistry 273, 13738 (1998).
  • Ikura et al. 2004 T. Ikura, Y. Urakubo, and N. Ito, Chemical Physics 307, 111 (2004).
  • Harris et al. 1997 L. J. Harris, S. B. Larson, K. W. Hasel, and A. McPherson, Biochemistry 36, 1581 (1997).
  • Kwong et al. 1998 P. D. Kwong, R. Wyatt, J. Robinson, R. W. Sweet, J. Sodroski, and W. A. Hendrickson, Nature 393, 648 (1998).
  • Gilson and Honig 1986 M. K. Gilson and B. H. Honig, Biopolymers 25, 2097 (1986).
Figure 1: Place Figure 1 caption here. In the case of reproduced figures in review articles, you must obtain the publisher’s permission and state a suitable notice here along with a citation.
Figure 2: Graphical representation of the implicit-solvent model.
Figure 3: Thermodynamic cycle for the calculation of binding energy.
Figure 4: Graphical representation of the spherical model used for validation.
Figure 5: Mesh refinement study for the system in Figure 4.
Figure 6: Solvation energy for different mesh sizes. The dotted line represents the Richardson extrapolation from Table 1.
Figure 7: Time to solution versus errors for different mesh sizes. Errors were computed using the values of Table 1 as the reference.
Figure 8: Convergence and timing comparison between AMBER and AMOEBA

Figure 1

Christopher D. Cooper

J. Comput. Chem.

Figure 2

Christopher D. Cooper

J. Comput. Chem.

Figure 3

Christopher D. Cooper

J. Comput. Chem.

Figure 4

Christopher D. Cooper

J. Comput. Chem.

Figure 5

Christopher D. Cooper

J. Comput. Chem.

Figure 6

Christopher D. Cooper

J. Comput. Chem.

Figure 7

Christopher D. Cooper

J. Comput. Chem.

Figure 8

Christopher D. Cooper

J. Comput. Chem.

Protein G [kcal/mol] Convergence
1CRN -271.80 -268.85
1ENH -1356.87 -1363.9
1FSV -902.99 -895.46
1PGB -809.94 -812.06
1VII -574.87 -573.12
Table 1: Richardson extrapolated solvation energy and observed order of convergence.
Total dipole [Debye]
PDB code Vacuum APBS PyGBe
1CRN 62.54 82.02 82.60
1ENH 209.94 266.25 267.02
1FSV 184.20 210.70 207.80
1PGB 101.38 130.55 129.37
1VII 158.72 194.57 194.21
Table 2: Total dipole moment in vacuum and dissolved sate, calculated with APBS and PyGBe.
PDB code G [kcal/mol] Time [s]
1PGB 927 24634 -822.42 31.6 4
1LYZ 1961 42544 -1597.47 64.9 4
1A7M 2809 60582 -2587.44 129.1 4
1X1U 9476 184399 -4918.82 1525 4
1IGT 20176 426588 -12261.13 4239 5
Table 3: Solvation energy, time to solution using GPU, and number of self-consistent iterations for different sized biomolecules.
Force G [kcal/mol] G G G
field Complex HIV+CD4 Antibody [kcal/mol] [kcal/mol] [kcal/mol]
AMBER 1 -8754.31 -5181.56 -4128.96 -508.50 556.21 47.71
2 -4291.84 -2543.03 -2023.35 -254.84 274.54 19.70
4 -2062.65 -1225.64 -972.00 -127.13 135.00 7.87
6 -1322.33 -787.81 -622.69 -85.00 88.17 3.17
CHARMM 1 -8799.66 -5187.31 -4169.20 -519.70 556.85 37.15
2 -4313.50 -2546.07 -2043.32 -259.84 275.89 16.05
4 -2073.57 -1227.20 -981.84 -129.93 135.5 5.57
6 -1329.51 -788.98 -629.17 -86.61 88.64 2.03
PARSE 1 -12151.82 -6926.11 -5839.22 -552.72 613.51 60.78
2 -5955.71 -3395.56 -2862.21 -276.35 302.06 25.71
4 -8260.06 -1632.20 -1376.24 -138.18 148.38 10.2
6 -1834.25 -1047.78 -882.49 -92.12 96.02 3.9
AMOEBA 1 -8511.01 -4846.20 -3547.57 152.25 -117.24 35.01
w/ 1 -9724.01 -5541.78 -4121.54 113.32 -60.69 52.63
Table 4: Electrostatic component of binding energy of 1GC1 computed with the thermodynamic cycle in Figure 3.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description