The Poisson-Boltzmann equation is widely used for modeling the electrostatics of biomolecules, but the calculation results are sensitive to the choice of the boundary between the low solute dielectric and the high solvent dielectric. The default choice for the dielectric boundary has been the molecular surface, but the use of the van der Waals surface has also been advocated. Here we review recent studies in which the two choices are tested against experimental results and explicit-solvent calculations. The assignment of the solvent high dielectric constant to interstitial voids in the solute is often used as a criticism against the van der Waals surface. However, this assignment may not be as unrealistic as previously thought, since hydrogen exchange and other NMR experiments have firmly established that all interior parts of proteins are transiently accessible to the solvent. 

Many molecular events are associated with small or large conformational changes occurring in the corresponding proteins. Modeling such changes is a challenge and requires significant amount of computing time. From point of view of electrostatics, these changes can be viewed as a reorganization of local charges and dipoles in response to the changes of the electrostatic field, if the cause is insertion or deletion of a charged amino acid. Here we report a large scale investigation of modeling the changes of the folding energy due to single mutations involving charged group. This allows the changes of the folding energy to be considered mostly electrostatics in origin and to be calculated with DelPhi assigning residue-specific value of the internal dielectric constant of protein. The predicted energy changes are benchmarked against experimentally measured changes of the folding energy on a set of 257 single mutations. The best fit between experimental values and predicted changes is used to find out the effective value of the internal dielectric constant for each type of amino acid. The predicted folding free energy changes with the optimal, amino acid specific, dielectric constants are within RMSD=0.86 kcal/mol from experimentally measured changes.

This is a brief review of the computational modeling of protein-ligand interactions using a recently developed fully polarizable continuum model (FPCM) and rational drug design. Computational modeling has become a powerful tool in understanding detailed protein-ligand interactions at molecular level and in rational drug design. To study the binding of a protein with multiple molecular species of a ligand, one must accurately determine both the relative free energies of all of the molecular species in solution and the corresponding microscopic binding free energies for all of the molecular species binding with the protein. In this paper, we aim to provide a brief overview of the recent development in computational modeling of the solvent effects on the detailed protein-ligand interactions involving multiple molecular species of a ligand related to rational drug design. In particular, we first briefly discuss the main challenges in computational modeling of the detailed protein-ligand interactions involving the multiple molecular species and then focus on the FPCM model and its applications. The FPCM method allows accurate determination of the solvent effects in the first-principles quantum mechanism (QM) calculations on molecules in solution. The combined use of the FPCM-based QM calculations and other computational modeling and simulations enables us to accurately account for a protein binding with multiple molecular species of a ligand in solution. Based on the computational modeling of the detailed protein-ligand interactions, possible new drugs may be designed rationally as either small-molecule ligands of the protein or engineered proteins that bind/metabolize the ligand. The computational drug design has successfully led to discovery and development of promising drugs.

The definition of a molecular surface which is physically sound and computationally efficient is a very interesting and long standing problem in the implicit solvent continuum modeling of biomolecular systems as well as in the molecular graphics field. In this work, two molecular surfaces are evaluated with respect to their suitability for electrostatic computation as alternatives to the widely used Connolly-Richards surface: the blobby surface, an implicit Gaussian atom centered surface, and the skin surface. As figures of merit, we considered surface differentiability and surface area continuity with respect to atom positions, and the agreement with explicit solvent simulations. Geometric analysis seems to privilege the skin to the blobby surface, and points to an unexpected relationship between the non connectedness of the surface, caused by interstices in the solute volume, and the surface area dependence on atomic centers. In order to assess the ability to reproduce explicit solvent results, specific software tools have been developed to enable the use of the skin surface in Poisson-Boltzmann calculations with the DelPhi solver. The results indicate that the skin and Connolly surfaces have a comparable performance from this last point of view.

Salt influences protein stability through electrostatic mechanisms as well as through nonpolar Hofmeister effects. In the present work, a continuum solvation based model is developed to explore the impact of salt on protein stability. This model relies on a traditional Poisson-Boltzmann (PB) term to describe the polar or electrostatic effects of salt, and a surface area dependent term containing a salt concentration dependent microscopic surface tension function to capture the non-polar Hofmeister effects. The model is first validated against a series of cold-shock protein variants whose salt-dependent protein fold stability profiles have been previously determined experimentally. The approach is then applied to HIV-1 protease in order to explain an experimentally observed enhancement in stability and activity at high (1M) NaCl concentration. The inclusion of the salt-dependent non-polar term brings the model into quantitative agreement with experiment, and provides the basis for further studies into the impact of ionic strength on protein structure, function, and evolution. 

This paper summarizes the mathematical and numerical theories and computational elements of the adaptive fast multipole Poisson-Boltzmann (AFMPB) solver. We introduce and discuss the following components in order: the Poisson-Boltzmann model, boundary integral equation reformulation, surface mesh generation, the node-patch discretization approach, Krylov iterative methods, the new version of fast multipole methods (FMMs), and a dynamic prioritization technique for scheduling parallel operations. For each component, we also remark on feasible approaches for further improvements in efficiency, accuracy and applicability of the AFMPB solver to large-scale long-time molecular dynamics simulations. The potential of the solver is demonstrated with preliminary numerical results.

We recently developed the Image-Charge Solvation Model (ICSM), which is an explicit/implicit hybrid model to accurately account for long-range electrostatic forces in molecular dynamics simulations [Linet al., J. Chem. Phys., 131, 154103, 2009]. The ICSM has a productive spherical volume within the simulation cell for which key physical properties of bulk water are reproduced, such as density, radial distribution function, diffusion constants and dielectric properties. Although the reaction field (RF) is essential, it typically accounts for less than 2% of the total electrostatic force on a water molecule. This observation motivates investigating further the role of the RF within the ICSM. In this report we focus on distributions of forces and torques on water molecules as a function of distance from the origin and make extensive tests over a range of model parameters where Coulomb forces are decomposed into direct interactions from waters modeled explicitly and the RF. Molecular torques due to the RF typically account for 20% of the total torque, revealing why the RF plays an important role in the dielectric properties of simulated water. Moreover, it becomes clear that the buffer layer in the ICSM is essential to mitigate artifacts caused by the discontinuous change in dielectric constants at the explicit/implicit interface.

We present a solver for the Poisson-Boltzmann equation and demonstrate its applicability for biomolecular electrostatics computation. The solver uses a level set framework to represent sharp, complex interfaces in a simple and robust manner. It also uses non-graded, adaptive octree grids which, in comparison to uniform grids, drastically decrease memory usage and runtime without sacrificing accuracy. The basic solver was introduced in earlier works [16, 27], and here is extended to address biomolecular systems. First, a novel approach of calculating the solvent excluded and the solvent accessible surfaces is explained; this allows to accurately represent the location of the molecule's surface. Next, a hybrid finite difference/finite volume approach is presented for discretizing the nonlinear Poisson-Boltzmann equation and enforcing the jump boundary conditions at the interface. Since the interface is implicitly represented by a level set function, imposing the jump boundary conditions is straightforward and efficient.

A nonlocal continuum electrostatic model, defined as integro-differential equations, can significantly improve the classic Poisson dielectric model, but is too costly to be applied to large protein simulations. To sharply reduce the model's complexity, a modified nonlocal continuum electrostatic model is presented in this paper for a protein immersed in water solvent, and then transformed equivalently as a system of partial differential equations. By using this new differential equation system, analytical solutions are derived for three different nonlocal ionic Born models, where a monoatomic ion is treated as a dielectric continuum ball with point charge either in the center or uniformly distributed on the surface of the ball. These solutions are analytically verified to satisfy the original integro-differential equations, thereby, validating the new differential equation system.

Stochastic walk-on-spheres (WOS) algorithms for solving the linearized Poisson-Boltzmann equation (LPBE) provide several attractive features not available in traditional deterministic solvers: Gaussian error bars can be computed easily, the algorithm is readily parallelized and requires minimal memory and multiple solvent environments can be accounted for by reweighting trajectories. However, previouslyreported computational times of these Monte Carlo methods were not competitive with existing deterministic numerical methods. The present paper demonstrates a series of numerical optimizations that collectively make the computational time of these Monte Carlo LPBE solvers competitive with deterministic methods. The optimization techniques used are to ensure that each atom’s contribution to the variance of the electrostatic solvation free energy is the same, to optimize the bias-generating parameters in the algorithm and to use an epsilon-approximate rather than exact nearest-neighbor search when determining the size of the next step in the Brownian motion when outside the molecule. 

Experimental data shows that the binding of human prolactin (hPRL) to human prolactin receptor (hPRLr-ECD) is strongly pH-dependent, while the binding of the same receptor to human growth hormone (hGH) is pH-independent. Here we carry in silico analysis of the molecular effects causing such a difference and reveal the role of individual amino acids. It is shown that the computational modeling correctly predicts experimentally determined pKa's of histidine residues in an unbound state in the majority of the cases and the pH-dependence of the binding free energy. Structural analysis carried in conjunction with calculated pH-dependence of the binding revealed that the main reason for pH-dependence of the binding of hPRL-hPRLr-ECD is a number of salt-bridges across the interface of the complex, while no salt-bridges are formed in the hGH-hPRlr-ECD. Specifically, most of the salt-bridges involve histidine residues and this is the reason for the pH-dependence across a physiological range of pH. The analysis not only revealed the molecular mechanism of the pH-dependence of the hPRL-hPRLr-ECD, but also provided critical insight into the underlying physic-chemical mechanism.

Recently, it was demonstrated that implicit solvent models were capable of generating stable B-form DNA structures. Specifically, generalized Born (GB) implicit solvent models have improved regarding the solvation of conformational sampling of DNA [1, 2]. Here, we examine the performance of the GBSW and GBMV models in CHARMM for characterizing base flipping free energy profiles of undamaged and damaged DNA bases. Umbrella sampling of the base flipping process was performed for the bases cytosine, uracil and xanthine. The umbrella sampling simulations were carried-out with both explicit (TIP3P) and implicit (GB) solvent in order to establish the impact of the solvent model on base flipping. Overall, base flipping potential of mean force (PMF) profiles generated with GB solvent resulted in a greater free energy difference of flipping than profiles generated with TIP3P. One of the significant differences between implicit and explicit solvent models is the approximation of solute-solvent interactions in implicit solvent models. We calculated electrostatic interaction energies between explicit water molecules and the base targeted for flipping. These interaction energies were calculated over the base flipping reaction coordinate to illustrate the stabilizing effect of the explicit water molecules on the flipped-out state. It is known that nucleic base pair hydrogen bonds also influenced the free energy of flipping since these favorable interactions must be broken in order for a base to flip-out of the helix. The Watson-Crick base pair hydrogen bond fractions were calculated over the umbrella sampling simulation windows in order to determine the effect of base pair interactions on the base flipping free energy. It is shown that interaction energies between the flipping base and explicit water molecules are responsible for the lower base flipping free energy difference in the explicit solvent PMF profiles.

The structural and functional description of protein-protein complexes and their comprehension is a key concept, not only to increase the scientific knowledge in basic terms but also for the application to the biomedical and pharmaceutical industry. The binding association between proteins is nowadays attribute to a few key residues at the interface – the hot-spots. The complex between the RNase inhibitor (RI) and RNaseA protein provides an excellent system to study the role of the functional epitope as it is essential in various molecular recognition processes and constitute one of the tightest complexes known. An energetic pattern of the interface is accomplished by computational alanine scanning mutagenesis and a dynamical characterization is accomplished by a detailed study of the molecular dynamical simulations. A special emphasis is given to the role of solvation across the interface and the shielding of warm- and hot-spots from water.

A new VMD plugin that interfaces with DelPhi to provide ensemble-averaged electrostatic calculations using the Poisson-Boltzmann equation is presented. The general theory and context of this approach are discussed, and examples of the plugin interface and calculations are presented. This new tool is applied to systems of current biological interest, obtaining the ensemble-averaged electrostatic properties of the two major influenza virus glycoproteins, hemagglutinin and neuraminidase, from explicitly solvated all-atom molecular dynamics trajectories. The differences between the ensemble-averaged electrostatics and those obtained from a single structure are examined in detail for these examples, revealing how the plugin can be a powerful tool in facilitating the modeling of electrostatic interactions in biological systems.

Here we report a web server, the DelPhi web server, which utilizes DelPhi program to calculate electrostatic energies and the corresponding electrostatic potential and ionic distributions, and dielectric map. The server provides extra services to fix structural defects, as missing atoms in the structural file and allows for generation of missing hydrogen atoms. The hydrogen placement and the corresponding DelPhi calculations can be done with user selected force field parameters being either Charmm22, Amber98 or OPLS. Upon completion of the calculations, the user is given option to download fixed and protonated structural file, together with the parameter and Delphi output files for further analysis. Utilizing Jmol viewer, the user can see the corresponding structural file, to manipulate it and to change the presentation. In addition, if the potential map is requested to be calculated, the potential can be mapped onto the molecule surface. The DelPhi web server is available from http://compbio.clemson.edu/delphi webserver.

Proton transport is one of the most important and interesting phenomena in living cells. The present work proposes a multiscale/multiphysics model for the understanding of the molecular mechanism of proton transport in transmembrane proteins. We describe proton dynamics quantum mechanically via a density functional approach while implicitly model other solvent ions as a dielectric continuum to reduce the number of degrees of freedom. The densities of all other ions in the solvent are assumed to obey the Boltzmann distribution. The impact of protein molecular structure and its charge polarization on the proton transport is considered explicitly at the atomic level. We formulate a total free energy functional to put proton kinetic and potential energies as well as electrostatic energy of all ions on an equal footing. The variational principle is employed to derive nonlinear governing equations for the proton transport system. Generalized Poisson-Boltzmann equation and Kohn-Sham equation are obtained from the variational framework. Theoretical formulations for the proton density and proton conductance are constructed based on fundamental principles. The molecular surface of the channel protein is utilized to split the discrete protein domain and the continuum solvent domain, and facilitate the multiscale discrete/continuum/quantum descriptions. A number of mathematical algorithms, including the Dirichlet to Neumann mapping, matched interface and boundary method, Gummel iteration, and Krylov space techniques are utilized to implement the proposed model in a computationally efficient manner. The Gramicidin A (GA) channel is used to demonstrate the performance of the proposed proton transport model and validate the efficiency of proposed mathematical algorithms. The electrostatic characteristics of the GA channel is analyzed with a wide range of model parameters. The proton conductances are studied over a number of applied voltages and reference concentrations. A comparison with experimental data verifies the present model predictions and validates the proposed model.