PoissonBoltzmann model for proteinsurface electrostatic interactions and gridconvergence study using the PyGBe code
Abstract
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 freeenergydriven 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 PoissonBoltzmann equation. We extended the implicitsolvent model used in the opensource code PyGBe to include surfaces of imposed charge or potential. This code solves the boundary integral formulation of the PoissonBoltzmann equation, discretized with surface elements. PyGBe has at its core a treecodeaccelerated Krylov iterative solver, resulting in scaling, with further acceleration on hardware via multithreaded 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 gridconvergence 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 gridconvergence 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 opensource 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 postprocessing scripts in Python to allow replication of the gridconvergence studies, all the way to generating the final plots, with a single command.
keywords:
biomolecular electrostatics, protein surface interaction, implicit solvent, PoissonBoltzmann, boundary element method, treecode, Python, CUDA1 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 proteinsurface interactions has remained elusive Gray2004 (); RabeVerdesSeegel2011 (), but adsorption mechanisms are governed by surface energy and often the dominant effect is electrostatics. As a freeenergydriven process, proteinsurface 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 PoissonBoltzmann equation and implicitsolvent 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 coworkers studied surfaceprotein interactions using continuum models discretized with boundaryelement YoonLenhoff1992 (); RothLenhoff1993 (); AsthagiriLenhoff1997 () and finitedifference methods YaoLenhoff2004 (); YaoLenhoff2005 (), in the context of ionexchange 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 implicitsolvent electrostatics. We have added the capability of modeling a protein near a charged surface to our code PyGBe, an opensource code^{1}^{1}1https://github.com/barbagroup/pygbe that solves the PoissonBoltzmann 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 wellknown apbs software CooperBarbashare154331 (); 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 gridconvergence 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 selfassembled monolayers (sam), which are modeled within an implicitsolvent 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, gridconvergence studies for the interaction free energy in this case, and gridconvergence 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 opensource code.
2 Implicitsolvent model for proteins near charged surfaces
The implicitsolvent model uses continuum electrostatics to describe the meanfield 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 solventexcluded 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 PoissonBoltzmann 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.
This model has been widely applied to investigate interactions between molecules, such as in proteinligand 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:
(1) 
Here, is the potential corresponding to the region with permittivity , and and are the set potential or charge on the nanosurface.
Boundary integral formulation
We express the system of partialdifferential equations in (2) by the corresponding integral equations along the interface and the nanosurface, and . Many authors have used the boundaryintegral representation of the implicitsolvent 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:
(2) 
where is the potential in region evaluated at the surface . and are defined as
(3) 
corresponding to the double and singlelayer potentials of and evaluated in the region . The functions and are the freespace Green’s functions of the Poisson (Laplace kernel) and linearized PoissonBoltzmann (Yukawa kernel) equations, respectively:
(4) 
We then take the limits , , on Equation (2), and apply the boundary conditions: , and to get the following system of boundary equations:
(5) 
Rearranging terms, we write Equation (2) in matrix form, as follows:
(6) 
If the surface has prescribed charge, corresponding to a Neumann boundary condition, , the equivalent derivation yields
(7) 
The formulation detailed in this section differs from the work by Lenhoff and coworkers YoonLenhoff1992 (); RothLenhoff1993 () because they consider an infinite charged surface, modeled using a modified Green’s function to account for the halfspace 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 nonplanar 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 nanostructures. 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 solventfilled cavities by adding more surfaces or interfaces. In our code, we deal with multiple surfaces in the manner presented by Altman and coworkers 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:
(8) 
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 submatrices of size rather than integral operators. Each element of a submatrix is an integral over one panel , with located at the center of the collocation panel , as follows:
(9) 
The terms on the righthand side and the unknown vectors in the discretized form of Equation (2) are subvectors of size . In this case, each element is the evaluation on the collocation panel , written as
(10) 
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 semianalytical 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 Treecodeaccelerated 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 matrixvector product (with a dense matrix) done in every iteration; this is the most timeconsuming 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 matrixvector 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 fastmultipole methods GreengardRokhlin1987 (), treecodes BarnesHut1986 (), and fastFouriertransform methods PhillipsWhite1997 (). In our numerical solution (developed as the opensource 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 treestructured set of boxes and approximates interactions between faraway 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 multipoleacceptance 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 CooperBarbashare154331 ().
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 boundaryelement 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 singlemolecule 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.,
(11)  
(12) 
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
(13) 
Surface free energy
Chan and coworkers ChanMitchell1983 (); CarnieChan1993 () derived the free energy for a surface with a set charge or potential. They describe the free energy on a surface as
(14) 
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 PoissonBoltzmann model.
Using constant values of and per panel, the discretized version of Equation (3.3) takes the form
(15) 
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
(16) 
The interaction free energy is
(17) 
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 closedform 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 LotanHeadGordon2006 (), 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 proteinsurface 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 coworkers CarnieChanGunning1994 (), the axial symmetry lets us formulate the solution of Equation (2) as an expansion in Legendre polynomials:
(18) 
being the degree Legendre polynomial and the modified spherical Bessel function of the second kind.
We make use of the following addition formula MarceljaMitchellNinhamSculley1977 (),
(19) 
to reformulate the expression for in Equation (4.1) as
(20) 
Here, is the modified spherical Bessel function of the first kind; is defined by
(21) 
where is the centertocenter distance; and is given by the following expression, with (in this context only) representing the gamma function:
(22) 
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
(23) 
for the first expression of Equation (4.1), and
(24) 
for Equation (4.1).
Applying the interface conditions for on Equation (23) and the first expression of Equation (4.1), produces
(25) 
where is the radius of surface .
Constant potential on .
The application of the boundary condition on , , where is independent on , gives
(26) 
Combining Equations (4.1) and (26) yields the following system of equations for the coefficients and
(27) 
where
(28) 
Constant surface charge on .
In this case, the application of the boundary condition on , , where is independent on , gives
(29) 
where
(31) 
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
(32) 
Applying the boundary conditions at on Equation (23), we can rewrite in terms of the already computed and :
(34) 
Because the charge is located at , only the terms of Equation (33) will survive, and the potential at this location is:
(35) 
For the isolated molecule, makes , which nullifies the sum in Equation (4.2) and for , from the system in Equation (4.1), is
(36) 
Surface free energy with set potential
We can expand from Equation (3.3) in Legendre polynomials as
(37) 
(38) 
If the surface is isolated, makes , and the free energy in this case is
(39) 
where is taken from the system in (4.1) considering , which results in
(40) 
Surface free energy with set charge
We can expand from Equation (3.3) in Legendre polynomials as
(41) 
(42) 
For the isolated surface, and , and the free energy is
(43) 
where is calculated from the system in (4.1) considering , which results in
(44) 
5 Results of the gridconvergence 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 gridconvergence 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 centertocenter 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 5 presents the results of the gridconvergence 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 boundaryelement integrals differently for closeby and faraway elements, and use a semianalytical 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 multipoleacceptance criterion. The final numerical parameter is the exit tolerance of the gmres solver.
# Gauss points:  Treecode:  gmres:  

inelement  closeby  faraway  tol.  
9 per side  37  3  300  15  0.5 
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 CooperBarbashare154331 (). 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 timeconsuming part of the algorithm is the matrixvector product within the Krylov solver, which scales as thanks to the treecode acceleration. However, the total timetosolution 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.
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 solventfilled 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 dipolemoment 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.
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 Richardsonextrapolation method for performing gridconvergence analysis, see our previous work CooperBardhanBarba2013 ().
# Gauss points:  Treecode:  gmres:  

inelement  closeby  faraway  tol.  
9 per side  19  7  500  15  0.5 
Energy [kcal/mol]  

Solvation  Surface  
Isolated  
Interacting 
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 opensource license with our previous publication CooperBardhanBarba2013 (), and continues to be available via its versioncontrol repository. Supplementing this paper, we prepared “reproducibility packages” containing running and postprocessing 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 closedform 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 gridconvergence 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 (); CooperBarbashare154331 (). 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 selfassembled 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 antibodybased 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:  

inelement  closeby  faraway  tol.  
9 per side  19  1  300  4  0.5 
Energy [kcal/mol]  

Solvation  Surface  
Isolated  
Interacting 
7 Conclusion
In this work, we used an implicitsolvent model to study proteinsurface interaction. We present for the first time and apply an extension of our opensource 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 brickshaped 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 samcoated nanoparticle, which can be represented by the brickshaped surface.
The addition of a surface with imposed charge or potential in the implicitsolvent 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 solventexcluded surface.
We conclude that this implicitsolvent model can offer a valuable approach in proteinsurface 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.
Acknowledgments
This work was supported by ONR via grant #N000141110356 of the Applied Computational Analysis Program. LAB also acknowledges support from NSF CAREER award OCI1149784 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.
References
 (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 ionexchange chromatography. 1. cytochrome c variants, Anal. Chem. 76 (2004) 6743–6752.
 (10) Y. Yao, A. M. Lenhoff, Electrostatic contributions to protein retention in ionexchange chromatography. 2. proteins with vaious degrees of structural differences., Anal. Chem. 77 (2005) 2157–2165.
 (11) C. M. Roth, B. L. Neal, A.