Biomolecular electrostatics—I want your solvation (model)^*

The nonpolar term is often modeled as proportional to the solute surface area, i.e. as

$$\begin{equation} \Delta W^{(\mathrm{np)}}(\mathbf{X}) = \gamma A(\mathbf{X}), \end{equation} \tag{ 8 }$$

where γ is the surface tension and A(X) is the surface area—a function that is heuristic, in the sense that it does not have an exactly defined form which can be approximated numerically to an arbitrary controlled precision. Roux and Simonson [36] present a simple, intuitive explanation of the widespread use of this functional form, and we will not address the nonpolar term further except to note that the subject remains an active research area [66–72]. Wagoner and Baker [73, 74] used MD simulations to illustrate the challenges with existing models, and Isele-Holder et al [75] have discussed calculating surface areas from MD.

Early studies of solvent influence on molecules demonstrated that macroscopic dielectric theory could be surprisingly useful in modeling the electrostatic component of solvation [26]. In this approach, the electrostatic potential φ(r) obeys a Poisson equation, which in its simplest form can be written as

$$\begin{equation} \nabla \cdot \left(\varepsilon(\mathbf{r})\nabla\varphi(\mathbf{r})\right) = -\rho_{\mathrm{solute}}(\mathbf{r})/\varepsilon_0, \end{equation} \tag{ 9 }$$

where r is the spatial position, ε(r) is a spatially varying dielectric constant, ε₀ is the permittivity of free space and ρ_solute(r) is the solute charge distribution (usually a set of discrete point charges situated at the atom positions). Equation (9) is usually presented as the fundamental continuum law, which is unfortunate, because dielectric theory offers substantially more sophistication.

In particular, for a pure-water solvent (i.e. with no mobile ions), equation (9) is just Gauss's law, which is a conservation law for the electric flux,

$$\begin{equation} \nabla \cdot \mathbf{E}(\mathbf{r}) = (\rho_{\mathrm{solute}}(\mathbf{r}) +\rho_{\mathrm{bound}}(\mathbf{r}))/\varepsilon_0, \end{equation} \tag{ 10 }$$

where ρ_bound(r), the second term on the right-hand side, represents the density of induced charge that arises in a dielectric medium in response to an imposed electric field. It is defined by

$$\begin{equation} \nabla \cdot \mathbf{P}(\mathbf{r}) = -\rho_{\mathrm{bound}}(\mathbf{r}), \end{equation} \tag{ 11 }$$

where the polarization density field P(r) is defined by the medium's constitutive relation. Equation (9) assumes that the constitutive relation is

$$\begin{equation} \mathbf{P}(\mathbf{r}) = \varepsilon_0 \chi(\mathbf{r}) \mathbf{E}(\mathbf{r}), \end{equation} \tag{ 12 }$$

where the electric susceptibility χ(r) is related to the dielectric constant by ε(r) = 1 + χ(r). Often these expressions are simplified by the definition of the electric displacement field D and its relationship to the electric field E(r) = −∇φ(r),

$$\begin{equation} \mathbf{D}(\mathbf{r}) = \varepsilon(\mathbf{r})\mathbf{E}(\mathbf{r}). \end{equation} \tag{ 13 }$$

Of course, as with all other forms of differential equations, the mathematical problem is not well posed until the boundary conditions are specified.

Very commonly, a simple dielectric model is used in which space is separated into distinct solute and solvent regions, and then each is modeled as a uniform dielectric. In other words, the solvent is usually modeled with water's macroscopic dielectric constant of about 80, and the solute permittivity is assigned a value between 1 and 10; the definition of the solute dielectric constant remains controversial (e.g. [60]) and depends on the application. For that matter, in reality, the solvent dielectric constant (relative permittivity) also depends on the ion concentrations and temperature. In the PMF view, the appropriate dielectric constant is usually 1 (giving intramolecular interactions the same exact value they would have in the corresponding MD simulation) or 2. The dielectric constant ε(r) is discontinuous as r crosses the boundary Ω separating these two regions.

From the Maxwell equations, the potential is continuous across the boundary, as is the normal component of the displacement:

$$\begin{equation} \varphi(\mathbf{r}_{\Omega}^{-}) = \varphi(\mathbf{r}_\Omega^{+}), \end{equation} \tag{ 14 }$$

$$\begin{equation} \varepsilon_{\mathrm{protein}} \hat{\bf n}(\mathbf{r}_\Omega) \cdot\nabla \varphi(\mathbf{r}_\Omega^{-}) = \varepsilon_{\mathrm{water}} \hat{\bf n}(\mathbf{r}_\Omega) \cdot\nabla \varphi(\mathbf{r}_\Omega^{+}), \end{equation} \tag{ 15 }$$

and defining

$$\begin{equation} \hat{\varepsilon} = \mathrm{\frac{\varepsilon_{\mathrm{p}}-\varepsilon_{\mathrm{w}}}{\frac{1}{2}(\varepsilon_p+\varepsilon_w)},} \end{equation} \tag{ 16 }$$

we can derive a boundary-integral equation (BIE) that is exactly equivalent to the mixed-dielectric Poisson problem of equation (9):

$$\begin{equation} \sigma(\mathbf{r}) + \hat{\varepsilon} \int_{\Omega} \frac{\partial}{\partial \mathbf{n}(\mathbf{r})}\frac{1}{4\pi \Vert\mathbf{r}-\mathbf{r}'\Vert} \sigma(\mathbf{r}') \,\mathrm{d}^2 \mathbf{r}' = -\hat{\varepsilon} \sum_{k=1}^{Q} \frac{\partial}{\partial \mathbf{n}(\mathbf{r})}\frac{1}{4\pi \Vert\mathbf{r}-\mathbf{r}_k\Vert}, \end{equation} \tag{ 17 }$$

where σ(r) is the infinitesimally thin layer of charge that develops at a dielectric boundary in response to an applied field; consider that the polarization field P(r) changes discontinuously across the boundary. Equation (17) is well known in several areas of science and engineering [54, 76–79]. We refer readers unfamiliar with BIEs to [80, 81].

The situation is orders of magnitude more complicated still if the solvent includes mobile ions—which is to say, every biologically relevant solvent. A simplistic approximation assumes that the ions' effect can be modeled using concentration fields that are proportional to the Boltzmann factor for an ion placed at the given point, i.e.

$$\begin{equation} \nabla \cdot (\varepsilon(\mathbf{r}) \nabla \varphi(\mathbf{r})) = \rho(\mathbf{r}) + \sum_{i=1}^{K} c_i q_i \exp(-q_i \varphi(\mathbf{r})/k_{\mathrm{B}} T), \end{equation} \tag{ 18 }$$

where K denotes the number of ionic species, k_B is the Boltzmann constant, T is the temperature and c_i and q_i represent the bulk concentration and charge of the ith ion species. Equation (18) is a particularly simple form of the famous PB equation. Explicit-atom approaches are particularly challenging for dilute salt solutions; in contrast, continuum models capture these effects naturally. Nevertheless, the PB equation has its own severe limitations, neglecting numerous chemically and biologically important phenomena (see section 6). Furthermore, because the PB model is still a nonlinear PDE, it is far from trivial to analyze, and its linearized form (the LPB equation) has become a popular alternative. Sharp and Honig [82] have addressed many subtle aspects, as do numerous reviews [14] . Many implicit-solvent models rely on a still more drastic simplification by linearizing the right-hand side of equation (18), which amounts to assuming that the electrostatic potential is always very small (so that q_iφ(r)/k_BT ≪ 1) or that the bulk ion concentrations are vanishingly small.

The first numerical calculation of the mixed-dielectric Poisson problem equation (9) for a protein was performed in 1982 [86] using a finite-difference method (FDM). Careful work by Honig and co-workers in their development of DelPhi software [87–90] established many important details for numerical accuracy. Important advances addressed the need to approximate the boundary condition at infinity, ways to represent the complicated dielectric boundaries and methods to represent the protein charge distribution (a set of discrete point charges) using the grid points for the FDM [91, 92]. One particularly important advantage of volume-based methods is that they can be extended easily to solve the PB equation in either its linear or nonlinear forms, with the latter requiring more delicate analysis [82]. A second important capability is that the underlying dielectric model is easily changed from a piecewise-constant form (i.e. ε_protein in the protein and ε_water outside) to one with multiple regions or ion-exclusion regions around the molecule (Stern layer) [93, 94], or with a permittivity that varies smoothly from ε_protein to ε_water, e.g. [25].

The general aim of efforts to improve volume-based methods has been, of course, to allow fast and accurate calculations, which means high accuracy at low grid resolutions [91]. Here, we cover improvements that apply equally well to both static and dynamic calculations (Luo et al have studied extensively issues of dynamics [95–98]). An important recent development was the introduction of a post-processing step to compute the induced surface charge σ(r) from the grid potentials [99]; it was found that this approach offers substantially more accurate reaction potentials at the charge locations, for a given grid resolution.

One of the most valuable theoretical developments of the last few years has been the decomposition of the electrostatic potential into separate Coulombic and reaction terms [33, 100–102]:

$$\begin{equation} \varphi(\mathbf{r}) = \varphi^{\mathrm{Coul}}(\mathbf{r}) + \varphi^{\mathrm{reac}}(\mathbf{r}). \end{equation} \tag{ 19 }$$

The first term contains all the singularities due to the use of point charges, and the function φ^reac is harmonic, i.e. it satisfies the Laplace equation. These properties make the reaction potential suitable for rigorous mathematical analysis [33], and calculations exhibit superior accuracy at lower grid resolutions because the (smooth) reaction potential is easier to approximate accurately than the total potential which has sharp peaks corresponding to the singularities [101, 102]. We note that this is the same reason that smoothed dielectric boundaries [25] allow reduced grid resolutions. Wei and co-workers have combined this decomposition with the matched interface and boundary method to develop a second-order accurate solver [100, 101, 103, 104], and the Luo group has characterized several ways to improve the convergence rate of the iterative method [98, 105], including using a smooth dielectric boundary [98]. A decomposition similar to equation (19), proposed by Chern et al [106] and developed thoroughly by Holst et al [107], splits the total potential into three terms rather than two. The extra term derives from further decomposing the reaction potential φ^reac in equation (19), in order to address an instability that arises for volumetric numerical methods. Holst and co-workers [108] have also applied this decomposition to the Poisson–Nernst–Planck model for diffusion–reaction problems.

Early work on FEM simulations of elliptic PDEs [109–113] established the value of FEM for molecular electrostatics, and later work by Baker and Holst provided the mathematical and algorithmic foundations for the widely used APBS software [114, 115]. It is instructive to note the substantial advantages of the software's rigorous and general mathematical basis, which allow APBS to be a scalable and flexible framework for numerous continuum models, and Holst and co-workers [116] continue to advance finite-element modeling in a variety of ways, including meshing (as have several others recently, including in particular [85]). Very recently, adaptive grid FDMs [117] and novel finite-element approaches have also been introduced [118], and it will be interesting to see how their advantages can be leveraged among the many applications of implicit-solvent models.

The simplest example of a surface method is the BIE of equation (17) for the mixed-dielectric Poisson problem. BIEs have several advantages and disadvantages compared to PDE methods. An important advantage is that the boundary conditions at infinity are handled exactly, with all of the unknowns on the (finite) dielectric boundary itself, rather than the potential throughout all space and towards infinity. A closely related advantage, which we will not discuss further here, is that BIEs can analytically account for other semi-infinite domains, such as a planar membrane, using what are essentially image charges. A second related advantage is that the unknowns in BIE problems, being distributions on a surface, are only 'two-dimensional unknowns', rather than the three-dimensional unknown electrostatic potential. This means that, in principle, BIEs can be more efficient due to a smaller number of unknowns. Finally, the potential due to the protein charge distribution is captured exactly, because one computes its influence on the boundary (or boundaries) directly. This property is one reason why surface approaches have been comparatively popular in electronic structure methods, where the charge distribution (i.e. the electron density) is very complex compared to a set of point charges; furthermore, the smooth potential throughout the solute volume is the main subject of interest.

The main disadvantages of BIE are of three types. Firstly, implementing numerical calculations, using what are known as boundary-element methods (BEM), is more challenging. Secondly, the usual BIE approaches are limited to problems with piece-wise-constant dielectric constants, and also to linear problems. The latter issue is severely restrictive when exploring the space of possible theories, and especially given that standard BIE can be applied to linearized PB problems, but not the fully NLPB equation.

Expanding on the first class of challenges: numerical difficulties include primarily the fact that most approaches lead to a dense matrix equation. That is, from BEM one obtains a linear system of equations Ax = b where every entry of A is non-zero; thus, the memory required to store A grows quadratically with the number of unknowns n, and the time required to compute x using LU factorization increases cubically, i.e. as O(n³). In contrast, because finite-difference, finite-volume and finite-element methods approximate the differential form of the Poisson equation, they produce very sparse matrices.

The dense nature of A is easy to conceptualize if one imagines representing the induced boundary charge σ(r) as a set of discrete point charges on the boundary—every point charge produces an electric field everywhere else. Happily, this same conceptual picture leads to a much faster approach, a class of techniques known as fast BEM solvers, which allow these problems to be solved using time and memory that scale as O(n). First of all, the N² interaction can be computed using fast N-body algorithms such as the fast-multipole method [29, 119–121], tree-codes [122], multiple grid methods [123], P³M [124], precorrected-FFT [125], kernel-independent fast multipole method [126] or FFTSVD [127]. These algorithms allow one to perform accurate matrix–vector multiplication using time and memory that scale as O(n) or O(nlogn) rather than the normal O(n²) time and memory. The second component of the fast solver approach is to use the Krylov-subspace iterative methods [128, 129] such as GMRES [130] to compute an approximate solution $\hat {x}\approx A^{-1} b$ that lies in the Krylov subspace 〈b,Ab,A²b,...〉. It should be noted that iterative methods are widely used in volume-based methods as well, with successive over-relaxation among the most popular [90].

The second major challenge for BEM calculations is that the matrix entries are comparatively difficult to compute accurately, because they involve calculating integrals of singular or near singular functions [131]. Calculation accuracy can depend quite strongly on implementation details that might not seem important [132, 133]. The third disadvantage is that dynamics have been much more difficult to address, because most boundary-integral methods require some kind of mesh (see below), which must be regenerated along with the discretized matrix at each time step. This disadvantage is not fundamental, however: Totrov and Abagyan [134] demonstrated that BEM simulations could be fast enough for Monte Carlo simulations of protein folding, and Wang et al [135] introduced methods for stochastic dynamics.

3.2.1. Boundary-element methods.

3.2.2. Generating a discretized boundary representation (meshing).

3.2.3. Higher-order methods and other approaches.

Clearly, there is no shortage of numerical methods for solving the Poisson equation or its relatives. Popular high-quality software packages already offer high-speed calculations, thorough documentation and developer support. Some of these packages are freely available as open-source software (e.g. APBS [32, 114, 231]), whereas others are available to academic institutions for very modest fees. Almost always, such software is used in 'production environments'; in other words, the user has significant trust in the software and its developers, has extensively tested it in his or her own laboratory and finally a community of users (perhaps led by the developers themselves) has created software infrastructure to interface it with her other simulation code. For example, Baker et al have developed interfaces between their APBS application and CHARMM (Chemistry at Harvard Macromolecular Mechanics), as well as visual molecular dynamics (VMD). New algorithms and software will have minimal impact in the biophysical community unless the developers provide a variety of support mechanisms to make their new approach competitive as a production-ready code. We list a few below.

Firstly, new software must offer a reasonable interface to existing file formats in computational chemistry, e.g. files for the MD software CHARMM and AMBER (Assisted Model Building with Energy Refinement) or the PQR (Point, Charge(Q), Radius) format [232]. Secondly, the software source code must be readable and modifiable by other groups, or at least offer a high-quality application programming interface, because an interface that goes deeper than a 'black box' can lead to dramatic performance improvements. For instance, coupling fast BEM codes to electrostatic optimization theory can speed up calculations for drug analysis and design by several orders of magnitude [181, 233]. Thirdly, the numerical algorithm and software should be extensible to other models; that is, it should provide a means to solve other models. For instance, Orland and co-workers have implemented an efficient solver for their dipolar Poisson model [234–237] using APBS. Fourthly, developers must be willing to offer significant user support. Implicit-solvent models are applied to a diverse and growing range of problems in molecular science and engineering, and no single piece of software can cover them all without some modification, or at the least a moderate amount of assistance.

As mentioned in section 2.2, macroscopic dielectric theory assumes that a microscopic polarization charge density P(r) develops in response to an electric field E(r), with a purely local relationship determined by the susceptibility χ(r):

$$\begin{equation} \mathbf{P}(\mathbf{r})=\chi(\mathbf{r}) \mathbf{E}(\mathbf{r}). \end{equation} \tag{ 20 }$$

The electric displacement field D(r) is then

$$\begin{equation} \mathbf{D}(\mathbf{r}) = \varepsilon_0 \mathbf{E}(\mathbf{r}) + \mathbf{P}(\mathbf{r}), \end{equation} \tag{ 21 }$$

so that

$$\begin{equation} \nabla \cdot \mathbf{D}(\mathbf{r})= -\rho_{\mathrm{solute}}(\mathbf{r}). \end{equation} \tag{ 22 }$$

To motivate the CFA, consider the computational cost of the Poisson model as arising from the need to find a self-consistent set of fields P and E, because the polarization field induces an extra potential. The CFA makes an approximation to the self-consistency requirement by assuming that

$$\begin{equation} \mathbf{P}_{\mathrm{CFA}}(\mathbf{r}) = \chi(\mathbf{r}) \mathbf{E}_{\mathrm{Coul}}(\mathbf{r}), \end{equation} \tag{ 23 }$$

and thus the Coulomb field directly determines every other field; that is, the polarization charge density does not itself create any further polarization charge. This approximation is extensively discussed in an early paper by Kharkats et al  [240], who address its use in both local and nonlocal dielectric theories, and note its earlier use in modeling Marcus electron transfer. A later development of the theory by Schaefer and Froemmel emphasized its interpretation as a means of approximating the electrostatic energy density [241], and these authors introduced a quasi-image-charge approach for fast modeling. Like other works before and after [242], they decomposed the energy calculation into two summations, one over the atomic 'self energies' and the other over the pair-wise interactions. Luo et al [243] also developed a CFA-like approximation for Warshel's solvent model, which uses a lattice of Langevin dipoles rather than a macroscopic dielectric.

Ghosh et al [244] noted the possible advantages for surface integration over volume integration, and applied the divergence theorem to convert the CFA volume integral to a surface integral. They noted that the resulting surface integral can be viewed as an operator approximation to the BIE of equation (17), and commented that the CFA's exactness for a sphere with central charge can be attributed to symmetry. However, this explanation is not the entire story for the exactness, as discussed in section 4.4. The authors developed two sets of empirical correction terms that improved agreement with Poisson calculations, both short-range and long-range. An early and important review of Bashford and Case [22] employed a detailed analysis of the semi-infinite planar-boundary case to illustrate the CFA's weaknesses, and discuss the introduction by Srinivasan et al [245] of a GB-style approximation to include salt effects via LPB theory. Bashford and Case note, however, that 'there is no good reason to expect GB models to overcome ... intrinsic limitations of continuum models'. As implicit-solvent models add complexity to mitigate their inaccuracies, it will be imperative to heed to this consideration and aim at developing rigorous approaches to allow simple and easy transfer of new physics from advanced models to faster approximate ones.

We recently observed [246] an interesting common feature of most fast Poisson approximations: the reaction field is assumed to take a separable form in which the dielectric constants and the geometry are decoupled, i.e. the reaction potential due to a charge at r_j obeys the ansatz

$$\begin{equation} \varphi(\mathbf{r}_i) = g(\varepsilon_p,\varepsilon_w) h(V, \mathbf{r}_i,\mathbf{r}_j) q_j. \end{equation} \tag{ 24 }$$

The classic GB model follows this approximation, as do some fast models based on approximating equation (17) [247]. Failures associated with this assumption have been noted repeatedly [248], but not explained from the separability perspective. For example, Lee et al [249] raised the possibility that parameters could be dependent on the dielectric constants, and the subject has been discussed in a recent review [239]. A more detailed analysis of this assumption may enable simple and physically rigorous improvements to GB models and others.

The main idea of the GB model dates back to more than five decades ago (its early history is characterized in a paper detailing the pair-wise-descreening approximation [250]). One separates the electrostatic solvation free energy of the quadratic expression

$$\begin{equation} \Delta G^{\mathrm{solv,es}} = \frac{1}{2} q^T L q , \end{equation} \tag{ 25 }$$

where q is the vector of atomic charges and L is the reaction-potential operator, into contributions from the diagonal and off-diagonal components of L:

$$\begin{equation} \Delta G^{\mathrm{solv,es}} = \frac{1}{2} \sum_{i=1}^N L_{ii} q_i^2 + \frac{1}{2} \sum_{i}\sum_{j\ne i} L_{ij} q_i q_j, \end{equation} \tag{ 26 }$$

where the first sum accounts for the individual charges' 'self-energies', and the double summation includes the screened interactions between different charges. This decomposition then readily leads to the identification of the first term in terms of effective Born radii, which is defined as the radius of a sphere which would have the same solvation free energy (L_ii). Still et al [251] gave a simple empirical functional form to approximate L_ij in terms of the distance r_ij and the effective radii R_i and R_j. Even though more than 20 years have passed since their paper [251], only a few attempts have been made to modify this form, and none have made enormous advances in accuracy. In fact, a insightful mathematical analysis by Onufriev et al [248] showed that for charges in a spherical molecule, a strikingly similar function to Still's may be obtained. This work provided a much stronger basis for the GB approach than had been known previously.

In 2002, two important publications significantly changed the direction of GB research [27, 249]. Firstly, Onufriev et al [27] published a paper of significant importance for both practical and computational reasons; they used extensive PB calculations to calculate the 'perfect Born radii' and assessed the accuracy of the GB model with these ideal parameters. This process, by removing 'effective Born radius error' due to CFA or other heuristics, focuses distinctly on one central question about GB: the adequacy of the underlying Still equation. As a practical matter, the study was important because it demonstrated that GB is extremely accurate for Poisson problems, giving confidence in the method's wide use. Furthermore, their results indicated that the emphasis in model development could focus on designing methods to more accurately estimate radii (not just solvation energies). From the philosophical point of view, their work provided a new understanding of the limitations of GB models, namely that the GB should not, except in the cases of overparameterization (the 'fitting an elephant' problem), perform better than the results indicated in that paper. This work shows how an accurate numerical computation can provide a substitute for 'exact results', and thus give much deeper insights for the general problem than are available for analytically solvable cases. We also note that the authors have commented that the system Green's function is the basic Coulomb contribution plus the reaction-field component; we will return to this point in greater detail in section 4.4.

In a second paper, Lee et al [249] reported a major breakthrough, apparently serendipitous, in finding systematic corrections to the CFA's weaknesses. They noted that the CFA involves the energy-density volume integral $\int 1/r^4\, \mathrm {d}V$, where r is the distance from the charge location, so the electric field falls as 1/r² and the energy density |E|² falls as 1/r⁴. The 1/r⁴ integrand is of the same form as the integration between a monopole and an induced dipole. This observation, combined with the clear need for more accurate self-energies (beyond those computed via the CFA), suggested that higher-order interactions might significantly improve accuracy. They found that a 1/r⁵ volume integral (and the square root of the result) accomplished precisely this effect, although as the authors noted it is empirical, i.e. it lacks the physical interpretation of the CFA. Furthermore, this correction involves only a single parameter, which can be set so that the computed Born radii are exact for the sphere. The authors presented two versions of this much improved model, called GBMV. One was a slower version that employs a grid-based numerical integration (i.e. it actually computes the volume integral). The other was an approximate method based on analytically evaluated integrals. Although the former naturally proved to be more accurate than the latter, the latter's speed advantage and ready implementation in dynamics argued for its adoption in implicit-solvent MD over the grid-based method.

The work of Lee et al spurred several further analyses. The same group developed an improved version (called GBMV2) that used the fourth root of a 1/r⁷ volume integral and a modified version of the Still equation [252]. Also, Grycuk [253] and Wojciechowski and Lesyng [254] analyzed the sphere case extensively. Grycuk suggested a 1/r⁶ volume integral, and explicitly assumed that ε_protein ≪ ε_water (the work of the Onufriev group addressed the limitations of this common assumption [248]). Another development 'beyond the CFA' presented a generalization of the volume-integral methods [255], using what the authors termed an 'expansion of energy' approach (again, rather than potential). In particular, they suggest computing self-energies of several of these volume integrals,

$$\begin{equation} I_n = \left(\int_S \frac{n\cdot r}{r^n}\,\mathrm{d}S \right)^{1/n-3},\quad n\geqslant 4, \end{equation} \tag{ 27 }$$

with the self-energy defined by

$$\begin{equation} G_i^{\mathrm{self}} = -\left(\frac{1}{\varepsilon_{\mathrm{protein}}}-\frac{1}{\varepsilon_{\mathrm{water}}}\right) \cdot \left(A_0 + \sum_{n=4}^{\infty}A_n I_n\right), \end{equation} \tag{ 28 }$$

and A_n are determined by parameter fitting, and obviously the sum in practice has a finite number of terms (in practice they use 4–7). Note that this energy is still separable in epsilons and geometry. However, the authors do discuss the non-separable aspects of dependence and point out that their approach improves accuracy over a range of dielectric constants.

Bajaj and Zhao recently published a sophisticated implementation of surface GB [224], which combines a number of advanced numerical algorithms. Motivated by the computational cost of approximating Born radii from perfect ones [27, 238] and the need to compute them accurately to have stable simulations, their software comprises four key steps. Firstly, they triangulate the surface (they use a Gaussian definition [256]); secondly, they improve the mesh quality; thirdly, they represent the mesh using a splice surface (i.e. approximate curvature with high accuracy); finally, they use a fast, non-uniform FFT to perform the needed surface integrations. Their detailed analysis of numerical error and convergence should be a model for other work in the area, and allows them to discuss with confidence the application of their methods to the recent GBr6NL model [257]. An important point not discussed in that work, however, is that empirical correction terms such as in the original surface GB [244] are not always suitable for fast, scalable computation.

Whereas GB models have been developed to give at the end a simple, analytical functional form for the charge–charge interactions, relatively few publications have addressed the problem as one of mathematical approximation. The exception, in general, has been studies using spheres or other model geometries as limiting cases, and the fast approximations are not always exact even in simple geometries. This is unfortunate, because the formal problem statement is extremely simple, and as we shall see it provides a novel framework for developing fast continuum-model approximations: the challenge is to find a rapidly computed, and hopefully analytical approximation to the system Green's function, which furthermore should be differentiable with respect to the atom positions.

The readers are, of course, familiar with the notion that in a vacuum, a point charge at r₀ gives rise to the Coulomb potential

$$\begin{equation} \varphi(\mathbf{r}) = \frac{1}{4\pi\varepsilon_0 \Vert\mathbf{r}-\mathbf{r}_0\Vert}; \end{equation} \tag{ 29 }$$

this Green's function is known as the free-space Green's function, i.e. φ(r) = G^FS(r,r₀) satisfies

$$\begin{equation} \nabla^2 \varphi(\mathbf{r}) = -\delta(\mathbf{r}-\mathbf{r}_0)/\varepsilon_0, \end{equation} \tag{ 30 }$$

where δ(·) is the Dirac delta function, and the boundary condition is that the potential decays to zero at infinity. This is just a special case of the system Green's function G^sys(r,r₀), which is the response at r to a point source at r₀ in a specified system with boundary conditions, i.e.

$$\begin{equation} \nabla \cdot (\varepsilon(\mathbf{r}) \nabla_{\mathbf{r}} G^{\mathrm{sys}}(\mathbf{r},\mathbf{r}_0)) = -\delta(\mathbf{r}-\mathbf{r}_0)/\varepsilon_0, \end{equation} \tag{ 31 }$$

where ε(r) is just the relative permittivity throughout space. Crucially, the system Green's function is the sum of the direct Coulomb potential and the reaction potential—this is the same decomposition that has proven valuable in volume-based numerical methods (equation (19)). The function G^sys essentially inverts the PDE—including the dielectric constants, the shape of the boundary, the Debye screening parameter κ for LPB problems—so that given an arbitrary charge distribution ρ(r') in the molecule, one obtains the potential via the convolution

$$\begin{equation} \varphi^{\mathrm{total}}(\mathbf{r}) = \int G^{\mathrm{sys}}(\mathbf{r},\mathbf{r}') \rho(\mathbf{r}')\,\mathrm{d}^3 \mathbf{r}'\,. \end{equation} \tag{ 32 }$$

To obtain the reaction free energy, we make use of the fact that G^sys = G^FS + G^reac to obtain the familiar relation

$$\begin{equation} \Delta G^{\mathrm{solv,es}} = \frac{1}{2} \int \rho(\mathbf{r}')\varphi^{\mathrm{reac}}(\mathbf{r}')\,\mathrm{d}^3 \mathbf{r}', \end{equation} \tag{ 33 }$$

which simplifies for point charges to

$$\begin{equation} \Delta G^{\mathrm{solv,es}} = \frac{1}{2} \sum_{i=1}^N \sum_{j=1}^N q_i G^{\mathrm{reac}}(\mathrm{r}_i,\mathrm{r}_j) q_j. \end{equation} \tag{ 34 }$$

Roux and Simonson made essentially this point in early reviews [36, 63], but apparently its import went unrecognized. A likely explanation for this missed opportunity is that the usual decomposition for the electrostatic solvation free energy is as a pair of discrete summations over the atom charges, i.e.

$$\begin{equation} \Delta G^{\mathrm{solv,es}} = \frac{1}{2}\sum_{i=1}^{N} q_i^2 F(\mathbf{r}_i,\mathbf{r}_i) + \frac{1}{2} \sum_{i=1}^N \sum_{j \ne i} q_i q_j F(\mathbf{r}_i,\mathbf{r}_j), \end{equation} \tag{ 35 }$$

and we use Roux and Simonson's notation F(r_i,r_j) to highlight that these two equivalent expressions, equations (34) and (35), do not suggest the same approaches to modeling and computation.

The first term is appealing because it connects with the Born expression for self-energy of an ion; in this author's view, the self-energy is just a red herring, because it obscures the need to approximate the reaction field potential as a function. It should be noted that much early work, particularly by Schaefer et al [241, 242], others (SEDO etc), addressed the problem of the total energy density. However, the singular nature of the energy density (due to the Coulomb term in the total potential) necessitated complicated re-arrangements, such as integrating over volumes outside the atomic one. Our intention in highlighting this observation is not to dispute that self-energy and charge–charge interaction terms can indeed be treated separately, given the appropriate care, as the ACE model [242] demonstrated and as Sigalov and Onufriev showed more recently [248]. However, recent work has demonstrated convincingly [247, 258] that separating the two is not necessary and, in fact, that treating them together in a field approach presents new opportunities for mathematical analysis.

Using the PDE form of the Poisson problem, Borgis et al [258] developed a functional approach to the electrostatic free energy. Their work begins by developing a local approximation of the polarization field; this approximation appears to be equivalent to the BIBEE/CFA method described below, although more analysis is needed. As noted in work on BIBEE [247], this approximation is generally less accurate for problems dominated by the dipolar component of the total charge distribution. Borgis et al developed an effective re-scaling for these contributions using mathematical analysis and found excellent agreement for many conformations of small molecules. The authors then introduce what they term the variational CFA (VCFA), which approximates the polarization field by a set of 'virtual charges' situated at the charge locations; the variational approach implies that they obtain the best approximation possible for the given basis set, and note that the virtual charge distribution could also include multipoles centered at the charge locations. They show, furthermore, that the VCFA is exact for a dipole inside a sphere, which suggests that this approximate model automatically goes beyond the separability approximation. Their framework offers significant modeling power, as they note that their framework can naturally allow dynamics via Car–Parrinello-type methods, that different definitions of the dielectric boundary (sharp and smoothed) pose no difficulty, that they can improve GB radii estimation and that the variational approximation is readily extended to more sophisticated solvation models. In light of these advantages, Borgis's formalism appears to have a lot of potential (pun intended).

The possible value of fictitious charges for reducing computational complexity is not new, of course, and here we detail a connection between Borgis's framework and the DiSCO multiscale model of Beard and Schlick [259]. The DiSCO model was developed to study extremely large macromolecular systems (chromatin folding processes) for which efficient computational methods are an absolute necessity. Citing earlier 'reduced charge' models, including Stigter's pre-HPC line-of-charge approach for long DNA [260] and the work of Gabdoulline and Wade [261], the authors proposed to assemble a set of discrete LPB 'point charges' situated on a fictitious surface surrounding each solute of interest. Charge values are found by numerically minimizing the difference between a computed PB field and the approximated field. The fictitious point charge approach recalls the early CFA work of Schaefer and Froemmel [241], as well as the basic ideas of potential theory and BIEs. However, Borgis's VCFA approach, or a derivative thereof, may provide a rigorous unifying methodology for applications that span in scale from small molecules (ions or drugs) to large macromolecules. Beard and Schlick performed careful studies to validate the accuracy of their model, and a solid theoretical foundation such as VCFA should allow the design of rigorous error estimators to allow seamless adaptivity with well-characterized error control.

Bardhan et al have developed a complementary operator approximation based on the BIE of equation (17) [246, 247, 262]. This approximation, called boundary-integral-based electrostatics estimation (BIBEE), relies on the spectral properties of the boundary-integral operators. The first BIBEE study showed that the CFA can be interpreted in terms of the extremal eigenvalue and eigenvector of the boundary-integral operator, which are associated with the overall monopole of the solute charge distribution [247]. The same boundary-integral analysis illustrated that the CFA is exact not only for a sphere with central charge, but also for any charge distribution that generates a uniform potential (see a related discussion in [240]). A different approximation is called BIBEE/P, because it relies on BEM preconditioning, which is a way to accelerate the convergence of the iterative method for solving a full BEM problem. The BIBEE/P variant employs the other end of the operator spectrum, and later work proved that the two variants provide upper and lower bounds on the electrostatic solvation free energy [262]. More precisely, BIBEE/CFA provides an upper bound for all shapes, and BIBEE/P is a provable bound only for spheres and prolate spheroids; however, tests on hundreds of proteins (the Feig et al test set [263]) did not provide an example of its failure as an effective bound.

One advantage of the BIBEE algorithm is that the main computational expense is in computing the solute charges' Coulomb potentials and electric fields at the dielectric boundaries. As a result, existing scalable algorithms (e.g. parallel fast multipole methods [264]) are directly applicable without modification [187]. By developing an approximate model that uses only 'algorithmic primitives', it becomes possible to use external libraries such that the same BIBEE code can run on a single or multiple CPUs, on one GPU or multiple GPUs on a CPU, or even on a large cluster containing hundreds of GPUs [187]. More recently, the method's accuracy was substantially improved by studying the operator eigenvalues more carefully, using the analytically solved case of a sphere [246]. The new one-parameter model (called BIBEE/I) provided 4% accuracy on the Feig test set, and was competitive with GBMV for an ensemble of peptide conformations. However, to date the BIBEE model has not yet been applied to LPB problems, and there exists no implementation that might be suitable for dynamics.

The most elementary challenge in parameterization is, of course, to avoid overfitting. This problem should rarely, if ever, arise in implicit-solvent model parameterizations based on other computational results, because one can always generate more reference data. The more difficult task is evaluating what data the model must be able to predict and what it should be able to predict and what it should not be able to predict (that is, that success would have to be fortuitous). Because parameter fitting is inherently an ill-posed inverse problem, parameters must be fit to the most discriminating test cases that are available. For example, when gradients of potentials are of interest (forces [273, 274] or electric fields [259]), it is better (more stable) to parameterize against the gradients rather than the potentials themselves, because differentiation amplifies numerical error. Furthermore, semi-phenomenological approaches make transferability a serious problem; loosely speaking, a transferable parameter is one that can be fit for one application and still give accurate answers for others.

Although solvation free energies have been the most popular source of reference data, recent evidence suggests that more challenging data are needed if models are to reliably predict binding free energies. A study on estimating the electrostatic contributions to binding free energies illustrated that parameter variations can have substantial impact on the change in the electrostatic solvation free energy, also known as the 'desolvation penalty' [275]. Scarsi and Caflisch [276] recommended that parameterization should be conducted using binding free energies directly. Their work is especially noteworthy for its use of the Poisson model as the reference for the parameterization of faster approximate models.

Rankin et al [277] demonstrated the severity of the problem by conducting an exhaustive search in a reduced parameter space with three variables: a single scale factor for atomic radii (using AMBER as the reference radii), the solute dielectric constant, and, for the nonpolar solvation component, one coefficient of proportionality for the surface-area term. Using the experimental solvation free energies of over 200 small molecules and their BEM software for the Poisson model [152], they found that a surprising range of parameter values fit the data equally well; in fact, the nonpolar term could have been dropped entirely without substantially altering their results. A clear example of problems with transferability could be found in their search space, which included the AMBER charges with PARSE radii; predictions were substantially poorer than with other models.

More generally, for a fixed value of the nonpolar coefficient, there exists a compensating effect between reduced radii and increased dielectric constant that leads to a meaningful 'valley' in the parameter space, i.e. instead of a deep local minimum, there is a region of parameter space in which every point gives essentially the same agreement with experimental solvation free energies. However, these points give very different predictions for binding free energies, spanning almost 25 kcal mol⁻¹ in predicted affinity. These empirical results showed, following earlier work, that the electrostatic contribution to binding depends much more strongly on the dielectric constant than does the solvation free energy. The authors offered an incisive demonstration of this property using the analytically solvable system of two spherical molecules binding, so the smaller is completely buried in the bound state. Their results led them to state, 'Parameter sets that yield equivalent results for transfer free energies yield highly disparate results for binding free energies. Hence, parameterization of a solvation model against hydration data, although a necessary condition, is not sufficient to validate its use for calculating binding free energies in solution' (p 959 of [277]).

Similar concerns have been raised for GB models [278]. Michel et al used genetic algorithms to search a high-dimensional parameter space and generate multiple candidates with equivalently good fits. These calculations led to parameter sets with noticeably different parameters, but almost equally good fits for solvation free energies, finding that some 'better' models actually achieved their accuracy at the expense of predicting nonphysical positive electrostatic solvation free energies. As the authors note that GB models are often used in molecular binding applications, they also tested how well their optimized parameters reproduced PMFs of the approach between small molecules; again, equally good parameterizations for solvation free energies were insufficient for accuracy in the binding-type application.

Finally, MD force fields continue to evolve—for instance, adding atomic polarizability [279–282]—so it is important that the community develop standardized, automated methods and software for parameterizing the large number of implicit-solvent models. A still more complex challenge, with important applications in drug design and analysis, is the problem of developing general parameter sets suitable for drug-like compounds, which are many times more diverse than the most commonly simulated biological molecules (proteins, nucleic acids, lipids and carbohydrates) [283]. MD force fields for this enormous space of molecules have emerged relatively recently, including OPLS-AA, AMBER and CHARMM force fields [284–287]. Therefore, large-scale optimization methods such as implicit-solvent parameterization ought to be an active area, particularly in testing new and improved solvation models.

As noted in equation (9), the standard Poisson theory treats water as a linear dielectric. Numerous results from both experiment and MD suggest that this assumption is adequate in many cases, with numerous failure modes of varying subtlety and means for remedy. In particular, for monovalent ions, linear response appears to hold quite well [34, 288, 289, 326]; Rick and Berne [288] also comment on the problems with arbitrary cavity definitions, with findings similar to that discussed by Papazyan and Warshel [306].

Aqvist and Hansson approached their analysis of linear response in implicit-solvent models from the PMF perspective [327]. They follow Zwanzig's approach to obtain the free-energy difference between the uncharged and charged states of the solute. They make an important observation that the assumption of linear response implies that the electrostatic potential in the uncharged solute is zero (or negligible), and go beyond this point to derive thermodynamic invariants for linear-response theory and for the case when the potential in the uncharged solute is not zero [327]. To illustrate their theoretical analysis they study electrostatic solvation free energies of numerous solutes (both polar and ionic) in several solvents. Also, in discussing the implications of the linear-response assumption, they note that hydration-layer waters do exhibit net orientation [288, 326], and the arbitrariness of the cavity definition. They also make the important comment, 'Even if violations of the linear response assumption are found, the deviations might turn out to behave in a systematic way which could still be utilized for simplified free-energy calculations' (p 9512). Challenging this assumption, Cerutti et al [328] found via extensive MD simulations that the net electrostatic potential in an uncharged protein is significant, rather than close to zero; their results seem to hold several important implications for implicit-solvent models, and merit much more study.

Lin and Pettitt [329] have also addressed the assumption of zero potential in the uncharged state. Their work provided a complementary source of simulation data, as they studied amino-acid side-chain analogs using highly accurate continuum-model calculations, replica-exchange methods and an explicit-solvent linear-response approximation. They found noticeable deviations for highly polar molecules (analogs for asparagine and histidine), which occurred near the uncharged states. However, the total energies estimated by linear-response methods were not grossly in error, and furthermore, the relative energies were quite good; as the authors noted, these relative energies are more often the relevant quantities for applications involving protein mutations or molecular design. Matyushov et al [330] have also shown that the potential inside an uncharged cavity is not zero, using atomistic simulations and a model spherical solute called a 'Rossky cavity'. The authors suggested that the standard Maxwell dielectric boundary conditions are in error because the solvent molecules are free to rotate (for instance, in water, to improve hydrogen bonding). A small applied field—in their case, applied externally—does not necessarily re-orient the waters as macroscopic dielectric theory holds them. Their MD simulations indicated deviations from the induced surface charge predicted by standard Poisson theory, and the authors indicate that their results better match a 'Lorentz virtual cavity' model. It is unclear whether one might find modified Poisson boundary conditions for capturing these effects. Kimura et al [331] have used a novel discretecharge approximation to address related shortcomings, but more work is needed. One of the notable arguments raised by Matyushov and co-workers [330] is the case for a wider range of macroscopic observables in the validation of models: they show that the standard Poisson model and their Lorentz-cavity approach predict quite different results for dielectric constants of mixtures.

The standard Poisson model is obviously incapable of predicting atomic-level interactions between solvent and solute. For example, it cannot address how water molecules can 'bridge' two protein side chains by hydrogen bonding to both, thereby stabilizing them at a distance beyond direct contact but before continuum-level theory would be appropriate. Still more complex failures can occur in confined spaces such as enzyme active sites; even at nearly planar interfaces such as lipid bilayers, water structure leads to failures of linear response and local dielectric response [332]. Many results indicate that only a few layers of surrounding solvent are needed to substantially improve model realism [332], which has motivated the development of several semi-implicit solvation models to include these few layers. Gouda et al [42] demonstrated the impact of first-solvation-shell errors on binding free-energy estimation. Reduced screening due to solvent structure has motivated the development of nonlocal dielectric theory (see the following section). However, even these theories assume zero potential in the charged state, and as the work of Cerutti et al [328] demonstrates, solvent structure can indeed create a sizable electrostatic potential inside an uncharged protein. We note that solvent structure is captured in a natural way by more complex statistical–mechanical integral-equation theories [84, 333–335], even though such models are outside the scope of this review.

The PB equation's key limitations are well known, e.g. [336, 337], and easily stated. Firstly, the model assumes that the ions are point charges, which they obviously are not; secondly, that the ions are in equilibrium (thus the Boltzmann weighting for their concentrations) with all correlations between them neglected. An important consequence of assuming ions to be point charges is that predicted ion concentrations can be unphysically large [338]. A classical empirical correction, due to Stern, is to enforce an ion-exclusion layer of the ion radius near charged surfaces [339]. Lastly, the linearized form of the PB equation is only a justifiable approximation for very small electrostatic potentials, a condition which is violated near highly charged species such as DNA or protein active sites.

Implicit-solvent models are often used to avoid the large cost of explicit-solvent calculations, which are dominated by thousands of solvent molecules. Relatively few of these solvent molecules actually interact directly with the solute, i.e. the water and ions in the hydration layer immediately surrounding it. As one of the noted weaknesses of implicit-solvent models has been its lack of molecular realism, it is of interest to explore a 'best of both worlds' approach in which a small number of explicit-solvent molecules are modeled near the solute, with an implicit-solvent model employed outside the explicit region [417, 418]. Another advantage is that these models also eliminate the need for periodic boundary conditions in explicit-solvent calculations. Beglov and Roux [418] derived a solvent boundary potential model for a sphere whose radius varies with time to include all of the explicit-solvent molecules. The strong fundamentals in their approach allowed application to a variety of problems, including the computationally expensive parameterization of implicit-solvent models [266]. Later extensions allowed much more general boundary shapes [419, 420] and application to quantum-mechanical simulations [421]. A more recent approach was presented by Wagoner and Pande [422], who developed a semi-implicit model in which the solvent molecules are smoothly decoupled from the system as they approach a spherical boundary; other methods have used GB electrostatics [423] and BIE methods [424]. Recently, Li et al [425] introduced the 'elastic bag' model which allows a complex boundary to vary dynamically. These models and other studies suggest that in many cases, fewer than ten 'solvent shells' are needed to recover the results of fully explicit-solvent calculations [332, 424].

Semi-explicit models for electrolyte solutions have also been an active area [54, 426–430]. In these models, the water is treated implicitly (as a high-dielectric continuum) but the ions in solution are treated explicitly as charged hard spheres. This approach eliminates important weaknesses of the PB theories, namely that the ions' excluded volume is handled properly and ion–ion correlations are included naturally. However, these advantages are balanced by the need to compute statistical averages over the ions' degrees of freedom, often using Monte Carlo sampling [54, 426], but sometimes with Langevin dynamics [427]. These types of electrolyte models have been studied for decades in the physical chemistry community, where they are known as primitive-model electrolytes [431, 432]. These primitive-model calculations have been very successful in predicting the selectivity of some ion-channel proteins [429, 430].

Perhaps the two most widely known successes for massive parallelism in molecular simulations are MD [30, 433] and scalable FEMs [32]. However, promising results have been demonstrated for many other theories, including statistical–mechanical integral equations [388], DFT [411], BEM simulations [187] and GB-type models [434]. Significant effort is being made to port these algorithms to an emerging computing platform called graphics processing units (GPUs), and in this section we address how GPU computing may influence the evolution of implicit-solvent models.

For years, the data-processing capabilities of high-end graphics cards have grown significantly faster than the capabilities of standard microprocessors (CPU). Most graphics applications require fairly simple computations applied independently to a large number of data elements, for instance applying the same rotation to millions of triangles. This independence means that the overall computation is massively parallel, and so GPUs combine a large number of relatively simple compute cores, each of which applies the same computational kernel to a single element. Many problems in computational biology seem to map naturally to this hardware, and demonstrated applications include bioinformatics, estimating molecular binding affinities [435] and MD in explicit and implicit solvent [436–438]. In some problems, a GPU code can be dozens or even hundreds of times faster than the corresponding CPU code. However, before the exciting potential of these computing platforms can be reached for arbitrary applications, a GPU implementation must resolve several practical complications.

The main price for this dramatic performance improvement is code complexity, especially for scalability to clusters of GPU-enabled computers. Consider how many software packages are designed to run on a single CPU or at most a small cluster of CPUs, and how few applications can make effective use of massively parallel computers (due to the need for complicated analysis of what limits scalability). Single GPU computing is considerably more complex than single CPU computing, so very few algorithms can make effective use of even two GPUs attached to the same CPU (node)—much less exploit hundreds of nodes which each have one or more GPUs. Computational scientists specializing in parallel computing should therefore be integral team members for any serious GPU implementation.

Even when implementing an algorithm for a single GPU on a single node, there are significant complications that argue for collaboration with GPU experts and computer scientists. First, modern GPUs are fast enough that data movement is usually a bottleneck, i.e. it is hard to keep the compute cores busy all the time because data cannot be transferred fast enough to them. Algorithm experts can often improve performance without sacrificing calculation accuracy by minimizing data movement. A second complication arises because most compute cores are substantially faster at single-precision floating point arithmetic than they are at double precision. In pursuit of performance, many GPU implementations of existing software have moved from double to single precisions; the user should be careful to check that the applications of interest have been tested thoroughly for subtle issues that can arise (see, e.g., the thorough analysis by Elber and co-workers [439]).

To be sure, at some point, we may arrive at new mathematical models due to novel insights into algorithms that run efficiently on GPUs. In the short run, however, we must resist the temptation to introduce new, modified implicit-solvent models simply because they can be optimized to run extremely quickly on a single GPU. Firstly, unless the algorithm is implemented within an existing code, e.g. [440], coupling codes is a time-consuming and error-prone task. Secondly, GPU hardware continues to evolve rapidly, which risks an early obsolescence for proposed models. Thirdly, non-scalable GPU algorithms will not be suitable for the increasingly large biological systems of current study. Thus, a new model that runs on one GPU but does not scale to multiple GPUs (many workstations hold two or even four), or clusters of GPU-enabled computers, is of limited value. For these reasons, it seems prudent to focus GPU development on solidly understood and accepted models (and algorithms).

Why is defining a new mathematical model of solvation so much easier than solving it numerically? This state of affairs is as lamentable as it is obvious—yet the disparity is not in any sense 'fundamental'. In this final section, I would like to point out an emerging paradigm of automated code generation that may ultimately make testing new implicit-solvent models as simple as writing down the governing equations.

Admittedly, at first glance it appears that each new mathematical model would require a new implementation—usually a slow and complicated business, even if one is only modifying an existing code. Compounding the difficulty, the new implementation must be fast enough to run on problems of current scientific interest, and as always, parallelism significantly increases complexity. Thus, the computational chemist/biophysicist developing more realistic models seems to face an unavoidable reality that the vast majority of her efforts will address mundane but complicated problems of software engineering and numerical analysis. Similar problems are faced by scientists across many disciplines and length scales that range from the molecular to the global, including cardiac electrophysiology and geological dynamics (to name just two examples).

Computational methods built around high-level mathematical interfaces offer a powerful way to address these challenges. Just as high-level programming interfaces streamline the design of sophisticated new algorithms, high-level mathematical interfaces streamline the design of sophisticated new models. Many scientists in the computational biophysics and biochemistry communities are familiar with an excellent example without realizing it: the APBS software [32] is built on top of the extremely general FETk software toolkit developed by M Holst [441]. The FEniCS [442–446] and Sundance [447] projects also address the problem very generally, by automating the generation of finite-element simulation codes for continuum-theory models. In FEniCS, the user needs only to specify the so-called weak form of the PDE, along with basic aspects of the numerical method, e.g. how the finite-element mesh is to be used as input, and what basis functions should be used. These are specified at a high level of mathematical abstraction in either C + + or Python. The 'model compiler' then automatically generates a parallel FEM code for computing the model [448, 449]. Numerous applications show that these solvers already scale to several hundreds of nodes and improved scalability is a priority in ongoing work (J Hake, personal communication); such a parallelism should be more than adequate for many basic modeling studies in the near future. A demonstration calculation of turbulence in computational fluid dynamics illustrates the capabilities of high-level modeling flexibility [450], and recent work extends these tools from traditional FEM to newer discontinuous Galerkin FEM [451].

The main mechanism for the success of FEniCS and related projects such as Sundance is that the developers were willing to circumscribe a class of models (specifically, those based on weak forms) and then developed automatic methods to solve that whole class. Importantly, they did so while retaining a programming interface that allows users to solve closely related problems. The point is not that the present version of any of these libraries will solve all of our implicit-solvent models, of course. Instead, consider these ideas as important steps to let the scientist explore an important part of 'model space' efficiently. Sundance and FEniCS address models based on standard PDEs, but one could imagine similar tools for variational models, which have become increasingly popular recently [72, 375], classical DFT [378–380, 403, 404, 408, 412, 415, 416], or statistical–mechanical integral equations [84, 382, 384, 386]. With such a software framework, the scientist can stay focused on the application science itself, rather than spending time implementing complex numerical algorithms on sophisticated hardware. Many chemists are already familiar with the MADNESS software [452], which uses wavelets for electronic structure methods; it has been constructed to be very general, and could be a suitable starting point.