Poisson-Boltzmann model for protein-surface electrostatic interactions and grid-convergence study using the PyGBe code

Poisson-Boltzmann model for protein-surface electrostatic interactions and grid-convergence study using the PyGBe code

Christopher D. Cooper christopher.cooper@usm.cl Lorena A. Barba labarba@gwu.edu Department of Mechanical Engineering, Boston University, Boston, MA. Department of Mechanical Engineering, Universidad Técnica Federico Santa María, Valparaíso, Chile. Department of Mechanical & Aerospace Engineering, The George Washington University, Washington, D.C.

Interactions between surfaces and proteins occur in many vital processes and are crucial in biotechnology: the ability to control specific interactions is essential in fields like biomaterials, biomedical implants and biosensors. In the latter case, biosensor sensitivity hinges on ligand proteins adsorbing on bioactive surfaces with a favorable orientation, exposing reaction sites to target molecules. Protein adsorption, being a free-energy-driven process, is difficult to study experimentally. This paper develops and evaluates a computational model to study electrostatic interactions of proteins and charged nanosurfaces, via the Poisson-Boltzmann equation. We extended the implicit-solvent model used in the open-source code PyGBe to include surfaces of imposed charge or potential. This code solves the boundary integral formulation of the Poisson-Boltzmann equation, discretized with surface elements. PyGBe has at its core a treecode-accelerated Krylov iterative solver, resulting in scaling, with further acceleration on hardware via multi-threaded execution on gpus. It computes solvation and surface free energies, providing a framework for studying the effect of electrostatics on adsorption. We then derived an analytical solution for a spherical charged surface interacting with a spherical molecule, then completed a grid-convergence study to build evidence on the correctness of our approach. The study showed the error decaying with the average area of the boundary elements, i.e., the method is , which is consistent with our previous verification studies using PyGBe. We also studied grid-convergence using a real molecular geometry (protein G B1 D4), in this case using Richardson extrapolation (in the absence of an analytical solution) and confirmed the scaling in this case. PyGBe is open-source under an MIT license and is hosted under version control at https://github.com/barbagroup/pygbe. In addition, we prepared “reproducibility packages” to supplement this paper, consisting of running and post-processing scripts in Python to allow replication of the grid-convergence studies, all the way to generating the final plots, with a single command.

biomolecular electrostatics, protein surface interaction, implicit solvent, Poisson-Boltzmann, boundary element method, treecode, Python, CUDA
journal: Computer Physics Communications

1 Introduction

Proteins interacting with solid surfaces fundamentally appear in many biological processes. Adsorption serves a key function in natural activities, like blood coagulation, and in biotechnologies like tissue engineering, biomedical implants and biosensors. A full understanding of protein-surface interactions has remained elusive Gray2004 (); RabeVerdesSeegel2011 (), but adsorption mechanisms are governed by surface energy and often the dominant effect is electrostatics. As a free-energy-driven process, protein-surface interaction is difficult to study experimentally MijajlovicETal2013 (), and thus simulations offer a good alternative. Full atomistic molecular dynamics simulations demand large amounts of computing effort, however, so we often must resort to other methods.

Protein electrostatics can be studied via modeling approaches using the Poisson-Boltzmann equation and implicit-solvent representations. These models are popular for computing solvation energies in protein systems RouxSimonson1999 (); Bardhan2012 (), but few studies have included the effect of surfaces. Lenhoff and co-workers studied surface-protein interactions using continuum models discretized with boundary-element YoonLenhoff1992 (); RothLenhoff1993 (); AsthagiriLenhoff1997 () and finite-difference methods YaoLenhoff2004 (); YaoLenhoff2005 (), in the context of ion-exchange chromatography. They realized that van der Waals effects can be neglected for realistic molecular geometries RothNealLenhoff1996 () and that the model is adequate as long as conformational changes in the protein are slight YaoLenhoff2004 (); YaoLenhoff2005 ().

The aim of this work is to develop and assess a computational model to simulate proteins near engineered surfaces of fixed charge, using implicit-solvent electrostatics. We have added the capability of modeling a protein near a charged surface to our code PyGBe, an open-source code111https://github.com/barbagroup/pygbe that solves the Poisson-Boltzmann equations via an integral formulation, using a fast multipole algorithm and gpu hardware acceleration. Previously, we verified and validated PyGBe in its use to obtain solvation and binding energies, by comparing with analytical solutions of the equations and with results obtained using the well-known apbs software CooperBarba-share154331 (); CooperBardhanBarba2013 (). In the present work, we derived an analytical solution for a spherical molecule interacting with a spherical surface of prescribed charge, and used it to verify the code in its new application and study numerical convergence. Using the newly extended code, we also studied the interaction between protein G B1 D4 and a solid surface of imposed charge, and conduct a grid-convergence study using this more realistic surface geometry.

We intend our new modeling tool to be useful in studying the behavior of proteins as they adsorb on surfaces that have been functionalized with self-assembled monolayers (sam), which are modeled within an implicit-solvent framework as surfaces with prescribed charge. One application is biosensing, where the target molecules are captured on the sensor via ligand molecules (for which antibodies are a common choice). Favorable orientations of ligand molecules lead to greatly enhanced sensitivity of biosensors TajimaTakaiIshihara2011 (); TrillingBeekwilderZuilhof2013 (), because binding sites need to be physically accessible to the targets. Studies of protein orientation near charged surfaces might look at how orientation can be influenced by engineering decisions regarding surface preparation, to aid the design of better biosensors. We explore this application in a companion publication that obtains probability of orientations for an antibody near a surface, as function of changing conditions on charge and ionic strength CooperBarba2015b (). In this paper, we present the details of a new analytical solution for spherical charged surfaces and molecules, grid-convergence studies for the interaction free energy in this case, and grid-convergence studies for protein G B1 D4 alone and interacting with a charged surface. The detailed analysis of the model is complemented with a diligent effort for reproducibility and we deposit both input and results data in accessible and permanent archival storage, in addition to the open-source code.

2 Implicit-solvent model for proteins near charged surfaces

The implicit-solvent model uses continuum electrostatics to describe the mean-field potential in a molecular system. A typical system consists of a protein in a solvent, defining two regions: inside and outside the protein, with an interface marked by the solvent-excluded surface (ses). The ses, beyond which a water molecule cannot penetrate into the protein, can be generated by rolling a (virtual) spherical probe of the size of a water molecule around the protein (see Figure 1). Inside the protein, the domain has low permittivity () and there are point charges located at the positions of the atoms. The solvent region, representing water with salt, has a permittivity of . A system of partial differential equations models this situation, with a Poisson equation governing inside the protein and a linearized Poisson-Boltzmann equation governing in the solvent region. Appropriate interface conditions on the ses express the continuity of the potential and electric displacement, completing the mathematical formulation.

Figure 1: Sketch of the process for generating a solvent-excluded surface (ses): a protein molecule contains a set of atoms that define a radius upon applying a force field and a probe the size of a water molecule is rolled to define the ses. is the protein region and the solvent region.

This model has been widely applied to investigate interactions between molecules, such as in protein-ligand binding. We are interested here in an extension of the model to consider interactions between proteins and surfaces with an imposed potential or charge. This new setup is sketched in Figure 2, and is described mathematically by the following equations:


Here, is the potential corresponding to the region with permittivity , and and are the set potential or charge on the nanosurface.

Figure 2: Sketch of a molecule interacting with a surface: is the protein, the solvent region, is the ses and a nanosurface with imposed charge or potential.
Boundary integral formulation

We express the system of partial-differential equations in (2) by the corresponding integral equations along the interface and the nanosurface, and . Many authors have used the boundary-integral representation of the implicit-solvent model to compute solvation energies of proteins YoonLenhoff1990 (); Juffer1991a (); LuETal2006 (); BajajETal2011 (); AltmanBardhanWhiteTidor09 (); GengKrasny2013 (); CooperBardhanBarba2013 (), but apart from work led by Lenhoff YoonLenhoff1992 (), we know of no studies that account for interacting surfaces in the system.

Consider the setting in Figure 2 with prescribed potential at . The application of Green’s second identity on the first two equations of (2) yields:


where is the potential in region evaluated at the surface . and are defined as


corresponding to the double- and single-layer potentials of and evaluated in the region . The functions and are the free-space Green’s functions of the Poisson (Laplace kernel) and linearized Poisson-Boltzmann (Yukawa kernel) equations, respectively:


We then take the limits , , on Equation (2), and apply the boundary conditions: , and to get the following system of boundary equations:


Rearranging terms, we write Equation (2) in matrix form, as follows:


If the surface has prescribed charge, corresponding to a Neumann boundary condition, , the equivalent derivation yields


The formulation detailed in this section differs from the work by Lenhoff and co-workers YoonLenhoff1992 (); RothLenhoff1993 () because they consider an infinite charged surface, modeled using a modified Green’s function to account for the half-space domain. Lenhoff’s approach has the advantage that the charged surface does not require a mesh, but cannot be applied to any situation where the surface has non-planar geometry. An infinite surface may not be a good model if the size of a device is comparable to the protein’s, like it happens with nano-structures. In that case, one needs to be able to represent the detailed geometry of a device via a surface mesh, and satisfy boundary conditions there, as detailed above.

The system in (2) can be extended to account for Stern layers and solvent-filled cavities by adding more surfaces or interfaces. In our code, we deal with multiple surfaces in the manner presented by Altman and co-workers AltmanBardhanWhiteTidor09 (), as described in our previous paper CooperBardhanBarba2013 ().

3 Methods

3.1 Discretization

To numerically solve the system in (2), we discretize the boundaries into flat triangular panels and assume that and are constant within those panels. The discretized form of the integral operators is as follows:


where is the number of discretization elements on , and and are the constant values of and on panel (we are somewhat abusing the nomenclature here by reusing the symbol , which previously referred to the complete surface). By collocating on the center of each panel, we get a linear system of equations that look just like (2) or (2), but the coefficient matrix is formed by sub-matrices of size rather than integral operators. Each element of a sub-matrix is an integral over one panel , with located at the center of the collocation panel , as follows:


The terms on the right-hand side and the unknown vectors in the discretized form of Equation (2) are sub-vectors of size . In this case, each element is the evaluation on the collocation panel , written as


where is located at the center of panel .

In our numerical solution, integrals are calculated in three possible ways, depending on how close the panel is to the collocation point. When the collocation point is inside the element being integrated, we use a semi-analytical technique, with Gauss points placed along the edges of the element (HessSmith1967, , p. 49, ff.)ZhuHuangSongWhite2001 (). If the integrated element is closer than from the collocation point —where for the area of panel — we use a fine Gauss quadrature rule, with 19 or more points per element. Beyond a distance of , elements have only 1, 3, 4 or 7 Gauss points, depending on the case.

3.2 Treecode-accelerated boundary element method

Most modern implementations of the boundary element method (bem) use Krylov methods to solve the linear system, usually a general minimal residual method (gmres), which is agnostic to the structure of the matrix. In practice, Krylov solvers for bem require operations to obtain the unknown vector, where is the number of iterations to get a desired residual, and is much smaller than . The scaling is given by a matrix-vector product (with a dense matrix) done in every iteration; this is the most time-consuming part of the algorithm, and makes bem prohibitive for more than a few thousand discretization elements.

But when we inspect the approximation of the integrals in (3.1) with Gauss quadrature rules, we see that the matrix-vector product has the form of an -body problem, similar to gravitational potential calculations in planetary systems. In this case, the Gauss quadrature points act analogously to planets (sources of mass) and the collocation points are analogous to the locations where the gravitational potential is computed (targets points). There are several ways to accelerate this kind of computations, for example fast-multipole methods GreengardRokhlin1987 (), treecodes BarnesHut1986 (), and fast-Fourier-transform methods PhillipsWhite1997 (). In our numerical solution (developed as the open-source code PyGBe), we accelerate the -body calculation with a treecode BarnesHut1986 (); LiJohnstonKrasny2009 (), making this part of the algorithm scale as rather than .

The treecode algorithm groups the sources and targets in a tree-structured set of boxes and approximates interactions between far-away boxes using a series expansion—a Taylor series, in our case. This allows for controllable accuracy that depends on the number of terms used in the expansion and the multipole-acceptance criterion that defines the threshold where the distance between source and target is far enough to approximate the interactions with expansions. Details of our implementation of the treecode in PyGBe can be found in our previous work CooperBarba-share154331 ().

3.3 Energy calculation

Figure 2 shows a system with three types of free energy: Coulombic energy from the point charges, surface energy due to and solvation energy. The Coulombic energy arises simply from the Coulomb interactions of all point charges. This section describes how we compute the other two components of free energy in the boundary-element framework.

Solvation free energy

When a protein is in a solvated state, surrounded by water molecules that have become polarized, its free energy differs from its state in vacuo by an amount known as the solvation energy. Its free energy again differs in the presence of other structures in the solvent, e.g., other proteins or charged surfaces. In this work we use the term solvation energy to more broadly mean the change in free energy of the protein from its state in a vacuum, to its state in the solvent with any other components or structures. In single-molecule settings, this definition of solvation energy coincides with the energy required to solvate the molecule.

To calculate the solvation energy, the total minus the Coulomb potential is applied inside the protein, i.e.,


where is the charge distribution, consisting of point charges (which transforms the integral into a sum). The total minus Coulomb potential includes the reaction potential—representing the response of the solvent by polarization and rearrangement of free ions—and any effects from the immersed surface. We can also interpret it as the potential generated by the boundary of the molecular region . Taking the first expression of Equation (2) and subtracting out the Coulombic effect yields


Equation (11) requires evaluating for each point-charge location . We obtain this by discretizing Equation (13) and using the solution of the linear system in Equation (2) or Equation (2) as inputs.

Surface free energy

Chan and co-workers ChanMitchell1983 (); CarnieChan1993 () derived the free energy for a surface with a set charge or potential. They describe the free energy on a surface as


where and are the prescribed potential and surface charge, respectively. The potential is given by for the first expression and the surface charge by for the second one. This is valid because we are using a linearized Poisson-Boltzmann model.

Using constant values of and per panel, the discretized version of Equation (3.3) takes the form


where is the area of panel , and . To obtain the surface free energy, we can plug in the solution of the system in Equation (2) or (2) to Equation (3.3).

Interaction free energy

When there are two or more bodies in the solvent, they will interact electrostatically. In order to compute the energy of interaction, we need to take the difference between the total energy of the interacting system and the total energy of each isolated component, where the total free energy is given by


The interaction free energy is


where is the number of components in the system and is calculated over the isolated component .

4 Analytical solution

It is possible to derive a closed-form expression for the free energy of interaction between a spherical molecule with a centered charge and a spherical surface with imposed potential or charge, like the one sketched in Figure 3. There are such analytical expressions for interacting charged surfaces CarnieChanGunning1994 (), and interacting spherical molecules with multiple point charges inside LotanHead-Gordon2006 (), but not for a situation where surfaces and molecules interact. Having such an analytical solution is of great utility in the development of a computational model for protein-surface interaction, because it will allow for proper code verification.

4.1 Expansion in Legendre polynomials

The system of partial differential equations from Equation (2) models the electrostatic potential field in the setting of Figure 3. Following Carnie and co-workers CarnieChanGunning1994 (), the axial symmetry lets us formulate the solution of Equation (2) as an expansion in Legendre polynomials:


being the -degree Legendre polynomial and the modified spherical Bessel function of the second kind.

Figure 3: Sketch of system solved with Legendre polynomials expansions.

We make use of the following addition formula MarceljaMitchellNinhamSculley1977 (),


to reformulate the expression for in Equation (4.1) as


Here, is the modified spherical Bessel function of the first kind; is defined by


where is the center-to-center distance; and is given by the following expression, with (in this context only) representing the gamma function:


Legendre polynomials are orthogonal to each other, and is independent of . Thus, taking the inner product of the expressions in Equations (4.1) and (4.1) with , where or , yields


for the first expression of Equation (4.1), and


for Equation (4.1).

Applying the interface conditions for on Equation (23) and the first expression of Equation (4.1), produces


where is the radius of surface .

Constant potential on .

The application of the boundary condition on , , where is independent on , gives


Combining Equations (4.1) and (26) yields the following system of equations for the coefficients and




Constant surface charge on .

In this case, the application of the boundary condition on , , where is independent on , gives


Combining Equations (4.1) and (26) produces a system of equations for the coefficients and




4.2 Energy calculation

Solvation free energy of the molecule

According to Equation (11), the solvation free energy of a molecule with a centered charge is given by


and using Equation (4.1), the reaction potential from Equation (13) is:


Applying the boundary conditions at on Equation (23), we can rewrite in terms of the already computed and :


Because the charge is located at , only the terms of Equation (33) will survive, and the potential at this location is:


The result from Equation (4.2) in Equation (32) yields the solvation free energy.

For the isolated molecule, makes , which nullifies the sum in Equation (4.2) and for , from the system in Equation (4.1), is

Surface free energy with set potential

We can expand from Equation (3.3) in Legendre polynomials as


Applying Equation (4.2) in Equation (3.3) gives


If the surface is isolated, makes , and the free energy in this case is


where is taken from the system in (4.1) considering , which results in

Surface free energy with set charge

We can expand from Equation (3.3) in Legendre polynomials as


Applying Equation (4.2) into Equation (3.3) gives


For the isolated surface, and , and the free energy is


where is calculated from the system in (4.1) considering , which results in


5 Results of the grid-convergence study

To obtain the following results, we extended the PyGBe code to consider surfaces with prescribed charge or potential. For all runs, we used a workstation with Intel Xeon X5650 cpus and one nvidia Tesla C2075 gpu card (2011 Fermi). We used the free msms software SannerOlsonSpehner1995 () to generate meshes, and pdb2pqr Dolinsky04 () with an amber force field to determine the charges and van der Waals radii.

5.1 Verification against analytical solution

Using the analytical solution detailed in Section 4, we carried out a grid-convergence study of PyGBe extended to treat interacting surfaces with biomolecules. The setup consists of a spherical molecule with a Å radius and a centered charge of , interacting with a spherical surface of Å radius and an imposed potential of . The center-to-center distance between the spheres is Å, and they are dissolved in water with salt at 145mM, which gives a Debye length of 8 (), and permittivity . The permittivity inside the spherical protein is . Figure 4 shows a sketch of this system.

Figure 4: Sketch of system used in the convergence study of Figure 5.

Figure 5 presents the results of the grid-convergence analysis, where the error is the relative difference in interaction free energy between the analytical result from Section 4 and the numerical solution computed with PyGBe. The observed order of convergence of the three finest meshes was 1.007. Table 1 presents the numerical parameters used in this case. Recall from section 3 that we calculate the boundary-element integrals differently for close-by and far-away elements, and use a semi-analytical method for the element that contains the collocation point. The fine Gauss quadrature rule is used for elements closer than from the collocation point, where . For the treecode, is the maximum number of boundary elements per box, is the Taylor expansion truncation parameter and is the multipole-acceptance criterion. The final numerical parameter is the exit tolerance of the gmres solver.

# Gauss points: Treecode: gmres:
in-element close-by far-away tol.
9 per side 37 3 300 15 0.5
Table 1: Numerical parameters used in the code-verification runs with the analytical solution.
Figure 5: Grid-convergence study for the interaction free energy between a spherical molecule with a centered charge and a sphere with potential . Data sets, figure files plus running/plotting scripts are available under cc-by CooperBarba2015-share1348841 ().

As seen in Figure 5, the error decays with the average area of the boundary elements (), which is the expected behavior and consistent with our previous verification exercises CooperBarba-share154331 (). This proves that the extension of PyGBe to treat charged surfaces is solving the mathematical model correctly.

Obtaining the interaction free energy involves three separate calculations: one with both bodies (molecule and interacting surface with set potential) and one for each isolated body. Figure 6 shows the total time to solution for each mesh, including all three cases that need to be computed. The most time-consuming part of the algorithm is the matrix-vector product within the Krylov solver, which scales as thanks to the treecode acceleration. However, the total time-to-solution scales slightly worse than because the condition number of the system increases with , and we need more iterations of the Krylov solver to converge to the desired tolerance, as shown in Figure 7.

Figure 6: Total runtime for obtaining interaction free energy (requiring three separate runs: one for each surface in isolation and one for both together), for the two-sphere system of 5.1.
Figure 7: Number of iterations to converge for the two-sphere system of 5.1. Iteration count increases with problem size due to the conditioning properties of the boundary integral formulation in our model.

5.2 Protein G B1 D4

We computed the electrostatic field of protein G B1 D4 interacting with a 100Å100Å10Å block with surface charge density C/m. The protein was centered with respect to a 100Å100Å face, a distance 2Å above it. Since we did not consider any Stern layers or solvent-filled cavities, these tests contain only two surfaces: the protein’s ses and the charged surface. We also computed the electrostatic field generated by protein G B1 D4 and the surface by themselves.

The angle between the dipole moment of the protein and the vector normal to the surface was . The dipole-moment vector placed at the center of mass of the protein generates an axis, and we used the line of shortest distance between the outermost atom and this axis as a reference vector . The rotation angle is the angle between the normal vector to a 100Å10Å side face of the block and when , and is equal to in these tests. Figure 8 is a sketch of this arrangement.

Figure 8: Orientation of a protein near a charged surface. is the dipole moment vector, the vector between and the atom that is the furthest, and and are normal to their corresponding surfaces. is the angle between and , and the angle between and when and are aligned.

In this test, the solvent has no salt, i.e., , and its relative permittivity was 80. The region inside the protein had a relative permittivity of 4. We computed the solvation and surface energy using meshes with 1, 2, 4 and 8 elements per square Angstrom with the parameters detailed in Table 2. Using Richardson extrapolation and the result of the three finest meshes we calculated an approximate exact solution, shown in Table 3, which we consider as a reference to calculate estimated errors. The errors plotted in Figures 9 and 10 are the relative differences between the energy obtained with Richardson extrapolation and the results computed with each mesh. They decay as in both the solvation and surface energies for the finest three meshes, indicating that the calculations are in the asymptotic region and the geometry is well resolved in these cases. The observed order of convergence is 0.95 for the solvation energy and 1.12 for the surface energy in the isolated cases, and 0.96 for the solvation energy and 0.94 for the surface energy when the protein and surface were interacting. Using the extrapolated values from Table 3, we obtain an interaction free energy of [kcal/mol]. For details on the Richardson-extrapolation method for performing grid-convergence analysis, see our previous work CooperBardhanBarba2013 ().

# Gauss points: Treecode: gmres:
in-element close-by far-away tol.
9 per side 19 7 500 15 0.5
Table 2: Numerical parameters used in the convergence runs with protein G B1 D4.
Energy [kcal/mol]
Solvation Surface
Table 3: Extrapolated values of energy for protein G B1 D4.
Figure 9: Grid-convergence study of the solvation energy for an isolated protein G B1 D4 mutant, and the surface energy of an isolated surface with charge density of 0.05C/m.
Figure 10: Grid-convergence study of the solvation and surface energy for protein G B1 D4 mutant, interacting with a surface with a charge density of 0.05C/m. Data sets, figure files and running/plotting scripts available under cc-by CooperBarba2015-share1348803 ()


5.3 Reproducibility and data management

To facilitate the replication of our work, we consistently release code and data associated with every publication. In that context, PyGBe was released under an MIT open-source license with our previous publication CooperBardhanBarba2013 (), and continues to be available via its version-control repository. Supplementing this paper, we prepared “reproducibility packages” containing running and post-processing scripts in Python to generate Figures 5 and 10. The packages invoke PyGBe with the parameters and meshes reported here, and then produce the plots, all with a single command. The reproducibility packages are hosted on figshare, and are referenced in the respective captions.

6 Discussion

In order to study the interaction of proteins and charged surfaces, we extended PyGBe to account for surfaces with prescribed charge or potential. Unfortunately, there was no analytical solution available in the literature to compare and verify PyGBe’s extension. Section 4 derives a closed-form expression for a spherical molecule with a centered charge interacting with a spherical nanosurface with imposed charge or potential. We used this new analytical solution to conduct a grid-convergence study of the interaction energy (Figure 5). The error decays with the area, which is the expected behavior for a boundary element method with constant elements CooperBardhanBarba2013 (); CooperBarba-share154331 (). Discretization error is very small for a spherical geometry. To make sure that the errors due to integration, the treecode approximation and the gmres solver were even smaller, we chose all the numerical parameters for high accuracy. This allows us to observe the convergence with respect to the discretization only and extract the order. With more realistic molecular geometries, however, discretization errors will be larger and the requirements can be relaxed in the other numerical parameters of PyGBe, resulting in lower runtimes.

The results in Figures 9 and 10 show the applicability of this approach in more realistic situations. The setting in Figure 8 can model a nanostructure coated with a self-assembled monolayer (sam) interacting with a protein that will adsorb on that surface. Our application of interest in developing this model is the field of biosensors, where it can assist design through studies of electrostatic adsorption affecting protein orientation near the biosensor. With antibody-based biosensors, orientation determines the accessibility of reaction sites and is critical for sensitivity. In a companion publication CooperBarba2015b (), we present the first studies of protein orientation near charged surfaces using our modeling framework.

Figures 9 and 10 show the expected 1/N convergence with a simple protein. Just like in the case of the sphere, the simple geometry of the charged surface forced us to use very fine parameters in order to extract the order of convergence. If we needed to run this computation many times—for example, to sample different protein orientations—we might relax these parameters to obtain shorter computation times.

The extrapolated values in Table 3 are useful to find more relaxed parameters that still give acceptable results. For example, using a mesh density of 2 elements per square Angstrom on the charged surface and 4 elements per square Angstrom on the protein, and the parameters detailed in Table 4, we get the results in Table 5. These results are less than 2% away from those in Table 3, and each run takes less than one minute. Moreover, using the energy values from Table 5, the interaction free energy is [kcal/mol], which is very close to the extrapolated case ( [kcal/mol]).

# Gauss points: Treecode: gmres:
in-element close-by far-away tol.
9 per side 19 1 300 4 0.5
Table 4: Numerical parameters for relaxed runs with protein G B1 D4.
Energy [kcal/mol]
Solvation Surface
Table 5: Values of energy for protein G B1 D4 using the parameters in Table 4, and a mesh density of 4 elements per square angstrom in the protein and 2 elements per square angstrom on the charged surface

7 Conclusion

In this work, we used an implicit-solvent model to study protein-surface interaction. We present for the first time and apply an extension of our open-source PyGBe code to account for the presence of surfaces with imposed potential or charge. The new feature of the code was verified against an analytical solution, which we derived for that purpose.

To demonstrate the power of this approach in a more realistic setting, we performed tests of protein G B1 D4 near a brick-shaped surface with an imposed charge. The error in energy scaling with the area of boundary elements demonstrates that this extension of PyGBe is capable of resolving the mathematical model correctly. This test was motivated by the biosensing application, where a ligand molecule is adsorbed on a sam-coated nanoparticle, which can be represented by the brick-shaped surface.

The addition of a surface with imposed charge or potential in the implicit-solvent model falls naturally in a boundary integral approach. In this case, the region enclosed by the surface is not part of the domain, then, this surface only adds one equation to the linear system, rather than two, which is the case with the molecular solvent-excluded surface.

We conclude that this implicit-solvent model can offer a valuable approach in protein-surface interaction studies. This tool can be useful for orientation studies of ligand molecules in biosensors, either to find optimal adsorption conditions of salt concentration and surface charge, or to guide the design of better ligand molecules.


This work was supported by ONR via grant #N00014-11-1-0356 of the Applied Computational Analysis Program. LAB also acknowledges support from NSF CAREER award OCI-1149784 and from NVIDIA, Inc. via the CUDA Fellows Program. We are grateful for many helpful conversations with members of the Materials and Sensors Branch of the Naval Research Laboratory, especially Dr. Jeff M. Byers and Dr. Marc Raphael.


  • (1) J. J. Gray, The interaction of proteins with solid surfaces, Curr. Opin. Struct. Biol. 14 (2004) 110–115.
  • (2) M. Rabe, D. Verdes, S. Seeger, Understanding protein adsorption phenomena at solid surfaces, Adv. Colloid Interface Sci. 162 (2011) 87–106.
  • (3) M. Mijajlovic, M. J. Penna, M. J. Biggs, Free energy of adsorption for a peptide at a liquid/solid interface via nonequilibrium molecular dynamics, Langmuir 29 (9) (2013) 2919–2926, pMID: 23394469. arXiv:http://dx.doi.org/10.1021/la3047966, doi:10.1021/la3047966.
    URL http://dx.doi.org/10.1021/la3047966
  • (4) B. Roux, T. Simonson, Implicit solvent models, Biophys. Chem. 78 (1999) 1–20.
  • (5) J. P. Bardhan, Biomolecular electrostatics — I want your solvation (model), Comput. Sci. Disc. 5 (013001).
  • (6) B. J. Yoon, A. Lenhoff, Computation of the electrostatic interaction energy between a protein and a charged surface, J. Phys. Chem. 96 (1992) 3130–3134.
  • (7) C. M. Roth, A. M. Lenhoff, Electrostatic and van der Waals contributions to protein adsorption: computation of equilibrium constants, Langmuir 9 (1993) 962–972.
  • (8) D. Asthagiri, A. M. Lenhoff, Influence of structural details in modeling electrostatically driven protein adsorption, Langmuir 13 (1997) 6761–6768.
  • (9) Y. Yao, A. M. Lenhoff, Electrostatic contributions to protein retention in ion-exchange chromatography. 1. cytochrome c variants, Anal. Chem. 76 (2004) 6743–6752.
  • (10) Y. Yao, A. M. Lenhoff, Electrostatic contributions to protein retention in ion-exchange chromatography. 2. proteins with vaious degrees of structural differences., Anal. Chem. 77 (2005) 2157–2165.
  • (11) C. M. Roth, B. L. Neal, A.