Exploring a charge-central strategy in the solution of Poisson's equation for biomolecular applications

Solvent-scaling as an alternative to coarse-graining in adaptive-resolution simulations: The adaptive solvent-scaling (AdSoS) scheme

A new approach termed Adaptive Solvent-Scaling (AdSoS) is introduced for performing simulations of a solute embedded in a fine-grained (FG) solvent region itself surrounded by a coarse-grained (CG) solvent region, with a continuous FG ↔ CG switching of the solvent resolution across a buffer layer. Instead of relying on a distinct CG solvent model, the AdSoS scheme is based on CG models defined by a dimensional scaling of the FG solvent by a factor s, accompanied by an s-dependent modulation of the atomic masses and interaction parameters. The latter changes are designed to achieve an isomorphism between the dynamics of the FG and CG models, and to preserve the dispersive and dielectric solvation properties of the solvent with respect to a solute at FG resolution. This scaling approach offers a number of advantages compared to traditional coarse-graining: (i) the CG parameters are immediately related to those of the FG model (no need to parameterize a distinct CG model); (ii) nearly ideal mixing is expected for CG variants with similar s-values (ideal mixing holding in the limit of identical s-values); (iii) the solvent relaxation timescales should be preserved (no dynamical acceleration typical for coarse-graining); (iv) the graining level N_G (number of FG molecules represented by one CG molecule) can be chosen arbitrarily (in particular, N_G = s³ is not necessarily an integer); and (v) in an adaptive-resolution scheme, this level can be varied continuously as a function of the position (without requiring a bundling mechanism), and this variation occurs at a constant number of particles per molecule (no occurrence of fractional degrees of freedom in the buffer layer). By construction, the AdSoS scheme minimizes the thermodynamic mismatch between the different regions of the adaptive-resolution system, leading to a nearly homogeneous scaled solvent density s³ρ. Residual density artifacts in and at the surface of the boundary layer can easily be corrected by means of a grid-based biasing potential constructed in a preliminary pure-solvent simulation. This article introduces the AdSoS scheme and provides an initial application to pure atomic liquids (no solute) with Lennard-Jones plus Coulomb interactions in a slab geometry.

Recent Developments in Free Energy Calculations for Drug Discovery

The grand challenge in structure-based drug design is achieving accurate prediction of binding free energies. Molecular dynamics (MD) simulations enable modeling of conformational changes critical to the binding process, leading to calculation of thermodynamic quantities involved in estimation of binding affinities. With recent advancements in computing capability and predictive accuracy, MD based virtual screening has progressed from the domain of theoretical attempts to real application in drug development. Approaches including the Molecular Mechanics Poisson Boltzmann Surface Area (MM-PBSA), Linear Interaction Energy (LIE), and alchemical methods have been broadly applied to model molecular recognition for drug discovery and lead optimization. Here we review the varied methodology of these approaches, developments enhancing simulation efficiency and reliability, remaining challenges hindering predictive performance, and applications to problems in the fields of medicine and biochemistry.

Improved Poisson-Boltzmann Methods for High-Performance Computing

Implicit solvent models based on the Poisson-Boltzmann equation (PBE) have been widely used to study electrostatic interactions in biophysical processes. These models often treat the solvent and solute regions as high and low dielectric continua, leading to a large jump in dielectrics across the molecular surface which is difficult to handle. Higher order interface schemes are often needed to seek higher accuracy for PBE applications. However, these methods are usually very liberal in the use of grid points nearby the molecular surface, making them difficult to use on high-performance computing platforms. Alternatively, the harmonic average (HA) method has been used to approximate dielectric interface conditions near the molecular surface with surprisingly good convergence and is well suited for high-performance computing. By adopting a 7-point stencil, the HA method is advantageous in generating simple 7-banded coefficient matrices, which greatly facilitate linear system solution with dense data parallelism, on high-performance computing platforms such as a graphics processing unit (GPU). However, the HA method is limited due to its lower accuracy. Therefore, it would be of great interest for high-performance applications to develop more accurate methods while retaining the simplicity and effectiveness of the 7-point stencil discretization scheme. In this study, we have developed two new algorithms based on the spirit of the HA method by introducing more physical interface relations and imposing the discretized Poisson's equation to the second order, respectively. Our testing shows that, for typical biomolecules, the new methods significantly improve the numerical accuracy to that comparable to the second-order solvers and with ∼65% overall efficiency gain on widely available high-performance GPU platforms.

An efficient second-order poisson-boltzmann method

Immersed interface method (IIM) is a promising high-accuracy numerical scheme for the Poisson-Boltzmann model that has been widely used to study electrostatic interactions in biomolecules. However, the IIM suffers from instability and slow convergence for typical applications. In this study, we introduced both analytical interface and surface regulation into IIM to address these issues. The analytical interface setup leads to better accuracy and its convergence closely follows a quadratic manner as predicted by theory. The surface regulation further speeds up the convergence for nontrivial biomolecules. In addition, uncertainties of the numerical energies for tested systems are also reduced by about half. More interestingly, the analytical setup significantly improves the linear solver efficiency and stability by generating more precise and better-conditioned linear systems. Finally, we implemented the bottleneck linear system solver on GPUs to further improve the efficiency of the method, so it can be widely used for practical biomolecular applications. © 2019 Wiley Periodicals, Inc.

Robustness and Efficiency of Poisson-Boltzmann Modeling on Graphics Processing Units

Poisson-Boltzmann equation (PBE) based continuum electrostatics models have been widely used in modeling electrostatic interactions in biochemical processes, particularly in estimating protein-ligand binding affinities. Fast convergence of PBE solvers is crucial in binding affinity computations as numerous snapshots need to be processed. Efforts have been reported to develop PBE solvers on graphics processing units (GPUs) for efficient modeling of biomolecules, though only relatively simple successive over-relaxation and conjugate gradient methods were implemented. However, neither convergence nor scaling properties of the two methods are optimal for large biomolecules. On the other hand, geometric multigrid (MG) has been shown to be an optimal solver on CPUs, though no MG have been reported for biomolecular applications on GPUs. This is not a surprise as it is a more complex method and depends on simpler but limited iterative methods such as Gauss-Seidel in its core relaxation procedure. The robustness and efficiency of MG on GPUs are also unclear. Here we present an implementation and a thorough analysis of MG on GPUs. Our analysis shows that robustness is a more pronounced issue than efficiency for both MG and other tested solvers when the single precision is used for complex biomolecules. We further show how to balance robustness and efficiency utilizing MG's overall efficiency and conjugate gradient's robustness, pointing to a hybrid GPU solver with a good balance of efficiency and accuracy. The new PBE solver will significantly improve the computational throughput for a range of biomolecular applications on the GPU platforms.

Exploring a multi-scale method for molecular simulation in continuum solvent model: Explicit simulation of continuum solvent as an incompressible fluid

We explored a multi-scale algorithm for the Poisson-Boltzmann continuum solvent model for more robust simulations of biomolecules. In this method, the continuum solvent/solute interface is explicitly simulated with a numerical fluid dynamics procedure, which is tightly coupled to the solute molecular dynamics simulation. There are multiple benefits to adopt such a strategy as presented below. At this stage of the development, only nonelectrostatic interactions, i.e., van der Waals and hydrophobic interactions, are included in the algorithm to assess the quality of the solvent-solute interface generated by the new method. Nevertheless, numerical challenges exist in accurately interpolating the highly nonlinear van der Waals term when solving the finite-difference fluid dynamics equations. We were able to bypass the challenge rigorously by merging the van der Waals potential and pressure together when solving the fluid dynamics equations and by considering its contribution in the free-boundary condition analytically. The multi-scale simulation method was first validated by reproducing the solute-solvent interface of a single atom with analytical solution. Next, we performed the relaxation simulation of a restrained symmetrical monomer and observed a symmetrical solvent interface at equilibrium with detailed surface features resembling those found on the solvent excluded surface. Four typical small molecular complexes were then tested, both volume and force balancing analyses showing that these simple complexes can reach equilibrium within the simulation time window. Finally, we studied the quality of the multi-scale solute-solvent interfaces for the four tested dimer complexes and found that they agree well with the boundaries as sampled in the explicit water simulations.

Recent Developments and Applications of the MMPBSA Method

The Molecular Mechanics Poisson-Boltzmann Surface Area (MMPBSA) approach has been widely applied as an efficient and reliable free energy simulation method to model molecular recognition, such as for protein-ligand binding interactions. In this review, we focus on recent developments and applications of the MMPBSA method. The methodology review covers solvation terms, the entropy term, extensions to membrane proteins and high-speed screening, and new automation toolkits. Recent applications in various important biomedical and chemical fields are also reviewed. We conclude with a few future directions aimed at making MMPBSA a more robust and efficient method.

A Continuum Poisson-Boltzmann Model for Membrane Channel Proteins

Membrane proteins constitute a large portion of the human proteome and perform a variety of important functions as membrane receptors, transport proteins, enzymes, signaling proteins, and more. Computational studies of membrane proteins are usually much more complicated than those of globular proteins. Here, we propose a new continuum model for Poisson-Boltzmann calculations of membrane channel proteins. Major improvements over the existing continuum slab model are as follows: (1) The location and thickness of the slab model are fine-tuned based on explicit-solvent MD simulations. (2) The highly different accessibilities in the membrane and water regions are addressed with a two-step, two-probe grid-labeling procedure. (3) The water pores/channels are automatically identified. The new continuum membrane model is optimized (by adjusting the membrane probe, as well as the slab thickness and center) to best reproduce the distributions of buried water molecules in the membrane region as sampled in explicit water simulations. Our optimization also shows that the widely adopted water probe of 1.4 Å for globular proteins is a very reasonable default value for membrane protein simulations. It gives the best compromise in reproducing the explicit water distributions in membrane channel proteins, at least in the water accessible pore/channel regions. Finally, we validate the new membrane model by carrying out binding affinity calculations for a potassium channel, and we observe good agreement with the experimental results.

Acceleration of Linear Finite-Difference Poisson-Boltzmann Methods on Graphics Processing Units

Electrostatic interactions play crucial roles in biophysical processes such as protein folding and molecular recognition. Poisson-Boltzmann equation (PBE)-based models have emerged as widely used in modeling these important processes. Though great efforts have been put into developing efficient PBE numerical models, challenges still remain due to the high dimensionality of typical biomolecular systems. In this study, we implemented and analyzed commonly used linear PBE solvers for the ever-improving graphics processing units (GPU) for biomolecular simulations, including both standard and preconditioned conjugate gradient (CG) solvers with several alternative preconditioners. Our implementation utilizes the standard Nvidia CUDA libraries cuSPARSE, cuBLAS, and CUSP. Extensive tests show that good numerical accuracy can be achieved given that the single precision is often used for numerical applications on GPU platforms. The optimal GPU performance was observed with the Jacobi-preconditioned CG solver, with a significant speedup over standard CG solver on CPU in our diversified test cases. Our analysis further shows that different matrix storage formats also considerably affect the efficiency of different linear PBE solvers on GPU, with the diagonal format best suited for our standard finite-difference linear systems. Further efficiency may be possible with matrix-free operations and integrated grid stencil setup specifically tailored for the banded matrices in PBE-specific linear systems.

Numerical interpretation of molecular surface field in dielectric modeling of solvation

Continuum solvent models, particularly those based on the Poisson-Boltzmann equation (PBE), are widely used in the studies of biomolecular structures and functions. Existing PBE developments have been mainly focused on how to obtain more accurate and/or more efficient numerical potentials and energies. However to adopt the PBE models for molecular dynamics simulations, a difficulty is how to interpret dielectric boundary forces accurately and efficiently for robust dynamics simulations. This study documents the implementation and analysis of a range of standard fitting schemes, including both one-sided and two-sided methods with both first-order and second-order Taylor expansions, to calculate molecular surface electric fields to facilitate the numerical calculation of dielectric boundary forces. These efforts prompted us to develop an efficient approximated one-dimensional method, which is to fit the surface field one dimension at a time, for biomolecular applications without much compromise in accuracy. We also developed a surface-to-atom force partition scheme given a level set representation of analytical molecular surfaces to facilitate their applications to molecular simulations. Testing of these fitting methods in the dielectric boundary force calculations shows that the second-order methods, including the one-dimensional method, consistently perform among the best in the molecular test cases. Finally, the timing analysis shows the approximated one-dimensional method is far more efficient than standard second-order methods in the PBE force calculations. © 2017 Wiley Periodicals, Inc.

Modeling Membrane Protein-Ligand Binding Interactions: The Human Purinergic Platelet Receptor

Membrane proteins, due to their roles as cell receptors and signaling mediators, make prime candidates for drug targets. The computational analysis of protein-ligand binding affinities has been widely employed as a tool in rational drug design efforts. Although efficient implicit solvent-based methods for modeling globular protein-ligand binding have been around for many years, the extension of such methods to membrane protein-ligand binding is still in its infancy. In this study, we extended the widely used Amber/MMPBSA method to model membrane protein-ligand systems, and we used it to analyze protein-ligand binding for the human purinergic platelet receptor (P2Y12R), a prominent drug target in the inhibition of platelet aggregation for the prevention of myocardial infarction and stroke. The binding affinities, computed by the Amber/MMPBSA method using standard parameters, correlate well with experiment. A detailed investigation of these parameters was conducted to assess their impact on the accuracy of the method. These analyses show the importance of properly treating the nonpolar solvation interactions and the electrostatic polarization in the binding of nucleotide agonists and non-nucleotide antagonists to P2Y12R. On the basis of the crystal structures and the experimental conditions in the binding assay, we further hypothesized that the nucleotide agonists lose their bound magnesium ion upon binding to P2Y12R, and our computational study supports this hypothesis. Ultimately, this work illustrates the value of computational analysis in the interpretation of experimental binding reactions.

Calculating protein-ligand binding affinities with MMPBSA: Method and error analysis

Molecular Mechanics Poisson-Boltzmann Surface Area (MMPBSA) methods have become widely adopted in estimating protein-ligand binding affinities due to their efficiency and high correlation with experiment. Here different computational alternatives were investigated to assess their impact to the agreement of MMPBSA calculations with experiment. Seven receptor families with both high-quality crystal structures and binding affinities were selected. First the performance of nonpolar solvation models was studied and it was found that the modern approach that separately models hydrophobic and dispersion interactions dramatically reduces RMSD's of computed relative binding affinities. The numerical setup of the Poisson-Boltzmann methods was analyzed next. The data shows that the impact of grid spacing to the quality of MMPBSA calculations is small: the numerical error at the grid spacing of 0.5 Å is already small enough to be negligible. The impact of different atomic radius sets and different molecular surface definitions was further analyzed and weak influences were found on the agreement with experiment. The influence of solute dielectric constant was also analyzed: a higher dielectric constant generally improves the overall agreement with experiment, especially for highly charged binding pockets. The data also showed that the converged simulations caused slight reduction in the agreement with experiment. Finally the direction of estimating absolute binding free energies was briefly explored. Upon correction of the binding-induced rearrangement free energy and the binding entropy lost, the errors in absolute binding affinities were also reduced dramatically when the modern nonpolar solvent model was used, although further developments were apparently necessary to further improve the MMPBSA methods. © 2016 Wiley Periodicals, Inc.

Charge Central Interpretation of the Full Nonlinear PB Equation: Implications for Accurate and Scalable Modeling of Solvation Interactions

Continuum solvation modeling based upon the Poisson-Boltzmann equation (PBE) is widely used in structural and functional analysis of biomolecules. In this work, we propose a charge-central interpretation of the full nonlinear PBE electrostatic interactions. The validity of the charge-central view or simply charge view, as formulated as a vacuum Poisson equation with effective charges, was first demonstrated by reproducing both electrostatic potentials and energies from the original solvated full nonlinear PBE. There are at least two benefits when the charge-central framework is applied. First the convergence analyses show that the use of polarization charges allows a much faster converging numerical procedure for electrostatic energy and forces calculation for the full nonlinear PBE. Second, the formulation of the solvated electrostatic interactions as effective charges in vacuum allows scalable algorithms to be deployed for large biomolecular systems. Here, we exploited the charge-view interpretation and developed a particle-particle particle-mesh (P3M) strategy for the full nonlinear PBE systems. We also studied the accuracy and convergence of solvation forces with the charge-view and the P3M methods. It is interesting to note that the convergence of both the charge-view and the P3M methods is more rapid than the original full nonlinear PBE method. Given the developments and validations documented here, we are working to adapt the P3M treatment of the full nonlinear PBE model to molecular dynamics simulations.

Net charge changes in the calculation of relative ligand-binding free energies via classical atomistic molecular dynamics simulation

The calculation of binding free energies of charged species to a target molecule is a frequently encountered problem in molecular dynamics studies of (bio-)chemical thermodynamics. Many important endogenous receptor-binding molecules, enzyme substrates, or drug molecules have a nonzero net charge. Absolute binding free energies, as well as binding free energies relative to another molecule with a different net charge will be affected by artifacts due to the used effective electrostatic interaction function and associated parameters (e.g., size of the computational box). In the present study, charging contributions to binding free energies of small oligoatomic ions to a series of model host cavities functionalized with different chemical groups are calculated with classical atomistic molecular dynamics simulation. Electrostatic interactions are treated using a lattice-summation scheme or a cutoff-truncation scheme with Barker-Watts reaction-field correction, and the simulations are conducted in boxes of different edge lengths. It is illustrated that the charging free energies of the guest molecules in water and in the host strongly depend on the applied methodology and that neglect of correction terms for the artifacts introduced by the finite size of the simulated system and the use of an effective electrostatic interaction function considerably impairs the thermodynamic interpretation of guest-host interactions. Application of correction terms for the various artifacts yields consistent results for the charging contribution to binding free energies and is thus a prerequisite for the valid interpretation or prediction of experimental data via molecular dynamics simulation. Analysis and correction of electrostatic artifacts according to the scheme proposed in the present study should therefore be considered an integral part of careful free-energy calculation studies if changes in the net charge are involved.

Exploring accurate Poisson-Boltzmann methods for biomolecular simulations

Accurate and efficient treatment of electrostatics is a crucial step in computational analyses of biomolecular structures and dynamics. In this study, we have explored a second-order finite-difference numerical method to solve the widely used Poisson-Boltzmann equation for electrostatic analyses of realistic bio-molecules. The so-called immersed interface method was first validated and found to be consistent with the classical weighted harmonic averaging method for a diversified set of test biomolecules. The numerical accuracy and convergence behaviors of the new method were next analyzed in its computation of numerical reaction field grid potentials, energies, and atomic solvation forces. Overall similar convergence behaviors were observed as those by the classical method. Interestingly, the new method was found to deliver more accurate and better-converged grid potentials than the classical method on or nearby the molecular surface, though the numerical advantage of the new method is reduced when grid potentials are extrapolated to the molecular surface. Our exploratory study indicates the need for further improving interpolation/extrapolation schemes in addition to the developments of higher-order numerical methods that have attracted most attention in the field.