A boundaryintegral approach for the PoissonBoltzmann equation with polarizable force fields
Abstract
Implicitsolvent 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 PoissonBoltzmann solver PyGBe. Previous work from other researchers coupled AMOEBA with the finitedifference 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 pointcharge force fields in a PoissonBotlzmann 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 pointcharge and multipolar force fields. Finally, we saw that polarizability plays an important role to consider cooperative effects, for example, in binding energy calculations.
Keywords: PoissonBoltzmann, Implicit solvent, Polarizable force fields, Boundary element method, Electrostatics
Introduction
Implicitsolvent 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 PoissonBoltzmann 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 semianalytical^{17} methods, which are usually coupled to pointchargebased 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 solution^{24} and Generalized Born models^{25} to account for polarizable multipoles. Moreover, we can find several efforts towards including polarizability to PoissonBoltzmann^{26, 27, 28, 20, 29} and Poisson^{30, 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,^{1}^{1}1https://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 PoissonBoltzmann model is the work by Schnieders and coworkers,^{27} where the authors extended the widely used APBS solver^{10} (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 evenlyspaced 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 boundaryintegral PoissonBoltzmann 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 finitedifference model. Moreover, the BEM stiffness matrix is not affected by the charge description, making it realtively easy to port into existing boundaryelement codes.
Methodology
The boundary integral formulation in the implicit solvent model
The implicitsolvent 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 PoissonBoltzmann equation governs in the solvent, with appropriate interface conditions.
(1) 
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
(2) 
where is the potential in region computed at the surface (). and are the double and singlelayer potentials, evaluated in :
(3) 
with the free space Green’s function of the Laplace or linearized PoissonBoltzmann (Yukawa) equations:
(4) 
Evaluating Equation (The boundary integral formulation in the implicit solvent model) at the interface gives
(5) 
Equation (The boundary integral formulation in the implicit solvent model) is the formulation presented by Yoon and Lenhoff.^{41}
The implicitsolvent model with polarizable point multipoles
To consider a polarizable pointmultipole 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 righthand 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 selfconsistent approach,^{42} where the electric field and induced dipole are iteratively computed until convergence.
Computation of the righthand 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
(6) 
Considering a point multipole centered at the origin, up to the quadrupole moment, Equation (6) becomes
(7) 
making the righthand side of Equation (The boundary integral formulation in the implicit solvent model):
(8) 
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 ()
(9) 
where is the electric field at the location of the point multipole. The total dipole moment in Equation (8) is then
(10) 
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 ():
(11) 
We can obtain by substracting out the influence of the multipoles to Equation (The boundary integral formulation in the implicit solvent model), which leaves:
(12) 
The component of the electric field is then
(13) 
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 Coulombictype electrostatic potential from Equation (7), taken at the evaluation point. This is
(14) 
The detailed expressions for the derivatives of each term are shown in Equation (APPENDIX).
AMOEBA uses a groupbased polarization scheme for permanent multipoles^{44} 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 Tholelike 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
(15) 
which goes into Equation (9) to obtain the induced dipole moment.
Selfconsistent induced field
The dipole moment is an input to the implicitsolvent model, however, the electric field has an influence in the induced dipole component. Then, we need to solve for the induced dipole with a selfconsistent scheme, which is summarized below:

Guess (for example, ).

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

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

Calculate the total electric field with (15).

Get with Equation (9).

Go back to step 2 with the new , until convergence.
We update the induced dipole with a successiveover 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:
(16) 
Assuming a linear dielectric, we can calculate the free energy as
(17) 
where is the potential and the charge distribution. Then, we can write the solvation energy as
(18) 
Decomposing the potential in dissolved state into solvent and Coulombic components, as in Equation (11), Equation (18) becomes
(19) 
In standard nonpolarizable force fields, the charge distribution does not have an induced component, hence and , and only the first term in the righthand side survives. However, in the case of polarizable forcefields, 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
(20) 
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 (secondtolast 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
(21) 
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:
(22) 
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
(23) 
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
(24) 
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 groupbased scheme in the energy calculation. However, the induced dipole component is damped by the same Tholelike 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
(25) 
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:
(26) 
Here, is the coulombic energy in vacuum state for each case, computed with Equation (20). Note that in the case of standard pointcharge 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:
(27) 
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, nearsingular, and faraway. The integral becomes singular when the collocation point is inside the integration panel, and we use a semianalytical approach^{48, 49} that places Gauss nodes on the edges of the element. If the collocation point and the integration panel are closeby, 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 matrixvector product (an process) inside the GMRES, done once in every iteration. As we are using collocation and Gaussian quadrature with freespace Green’s functions, the matrixvector product can be seen as a Nbody 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 Nbody calculations to , and even , for example, FFTbased 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 Laplace^{54} and linearized PoissonBoltzmann 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 lowestlevel box is reached, and the sourcetarget interaction is performed directly. To determine if a box is far enough from a target, we use the multipoleacceptance criterion (MAC), which is:
(28) 
where is the size of the box and the boxtarget 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 multipoleacceptance 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 sourcetarget and boxtarget 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 sourcetarget and boxtarget interactions, and also for the computation of the righthand side in Equation (8) and the energy in Equation (24).
Results
We implemented the method described above in a forked Github repository^{2}^{2}2https://github.com/barbagroup/pygbe^{3}^{3}3https://github.com/cdcooper84/pygbe of the opensource code PyGBe.^{15, 57} The results presented in this section were obtained on a workstation with two Intel Xeon E52680v3 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 amoebapro04^{18} or amoebapro13,^{58} whereas for pointcharge 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 Ponder^{24} derived a closed expression for spherical cavities with a random distribution of multipoles, equivalent to the classic Kirkwood solution^{3} for point charges, and suggested using selfconsistent 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 offcentered 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 TinkerAPBS
The work by Schnieders and coworkers ^{27} validates the implementation of polarizable multipoles in APBS using 1CRN,^{61} 1ENH,^{62} 1FSV,^{63} 1PGB,^{64} and 1VII,^{65} and comparing their implicitsolvent 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 extrapolation^{66, 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 selfconsistent solver for the induced dipole moment was set to and the exit criterion of the GMRES in the PoissonBoltzmann 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 coworkers^{27} (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 faraway from the collocation point, however, for nearlysingular 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 1IGT^{70} using 4 vertices per square angstrom, and report them in Table 3. These runs were done using the amoebapro13^{58} 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 pointcharge 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 fixedcharge 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 HIV1 GP120 core complexed with CD4 and a neutralizing human antibody (PDB: 1GC1^{71}). 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 antigenbinding 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 faraway integrals at . For the case of pointcharge 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.
Discussion
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 Ponder^{24} 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 coworkers^{27} 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 finitedifference 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 singlecore, hence, APBS and PyGBeCPU 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 selfconsistent 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 implicitsolvent PoissonBoltzmann approach is efficient to analyze large biomolecular systems. Moreover, the solventexcluded surface mesh grows slower than a volumetric mesh as the molecule size increases, indicating that a boundary integral approach is more appropriate for a largescale 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 pointcharge force fields, slightly higher adhoc permittivities are usually used inside the protein^{72} 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 pointcharge counterparts using . This behaviour is expected. In fact, in pointcharge 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 CD4induced 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 coulombictype 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.
Conclusions
This paper presents an extension of the PoissonBoltzmann 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 righthand 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 pointcharge force fields, like AMBER, in a PoissonBoltzmann 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 pointcharge force fields the solvent induces repulsion, whereas for AMOEBA it favors binding.
We conclude that this boundaryintegral implicitsolvent 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 righthand side avoids any difficult pointmultipole transfer to the mesh, making a boundaryintegral representation ideal for these computations, specially looking at largescale simulations.
Acknowledgments
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.
Appendix
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
(29) 
where the normal moves with the integration variable . On the other hand, the second derivatives used in Equation (21) are
(30) 
The derivatives of the terms in Equation (Calculation of the electric field) are:
(31) 
The terms required to compute the second derivative of the potential due to a collection of multipoles are
(32) 
References
 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. HeadGordon, Journal of computational chemistry 38, 1275 (2017).
 Yap and HeadGordon 2010 E.H. Yap and T. HeadGordon, 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. HeadGordon, Annual review of physical chemistry 65, 149 (2014).
 Lipparini et al. 2014 F. Lipparini, L. LagardeÌre, B. Stamm, E. CanceÌ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 ComputerAided 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 PoissonBoltzmann equation with boundary element methods, Technical Report on figshare, CCBY 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
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  
APBS  PyGBe  APBS  PyGBe  
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 
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 
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 
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 