Problems of robustness in Poisson-Boltzmann binding free energies

Some Challenges of Diffused Interfaces in Implicit-Solvent Models

The standard Poisson-Boltzmann (PB) model for molecular electrostatics assumes a sharp variation of the permittivity and salt concentration along the solute-solvent interface. The discontinuous field parameters are not only difficult numerically, but also are not a realistic physical picture, as it forces the dielectric constant and ionic strength of bulk in the near-solute region. An alternative to alleviate some of these issues is to represent the molecular surface as a diffuse interface, however, this also presents challenges. In this work we analyzed the impact of the shape of the interfacial variation of the field parameters in solvation and binding energy. However we used a hyperbolic tangent functiontanh k p x $$ \left(\tanh \left({k}_px\right)\right) $$ to couple the internal and external regions, our analysis is valid for other definitions. Our methodology, restricted to the linear PB, was based on a coupled finite element (FEM) and boundary element (BEM) scheme that allowed us to have a special treatment of the permittivity and ionic strength in a bounded FEM region near the interface, while maintaining BEM elsewhere. Our results suggest that the shape of the function (represented byk p $$ {k}_p $$ ) has a large impact on solvation and binding energy. We saw that high values ofk p $$ {k}_p $$ induce a high gradient on the interface, to the limit of recovering the sharp jump whenk p → ∞ $$ {k}_p\to \infty $$ , presenting a numerical challenge where careful meshing is key. Using the FreeSolv database to compare with molecular dynamics, our calculations indicate that an optimal value ofk p $$ {k}_p $$ for solvation energies was around 3. However, more challenging binding free energy tests make this conclusion more difficult, as binding showed to be very sensitive to small variations ofk p $$ {k}_p $$ . In that case, optimal values ofk p $$ {k}_p $$ ranged from 2 to 20.

Optimal Dielectric Boundary for Binding Free Energy Estimates in the Implicit Solvent

Accuracy of binding free energy calculations utilizing implicit solvent models is critically affected by parameters of the underlying dielectric boundary, specifically, the atomic and water probe radii. Here, a multidimensional optimization pipeline is used to find optimal atomic radii, specifically for binding calculations in the implicit solvent. To reduce overfitting, the optimization target includes separate, weighted contributions from both binding and hydration free energies. The resulting five-parameter radii set, OPT_BIND5D, is evaluated against experiment for binding free energies of 20 host-guest (H-G) systems, unrelated to the types of structures used in the training. The resulting accuracy for this H-G test set (root mean square error of 2.03 kcal/mol, mean signed error of -0.13 kcal/mol, mean absolute error of 1.68 kcal/mol, and Pearson's correlation of r = 0.79 with the experimental values) is on par with what can be expected from the fixed charge explicit solvent models. Best agreement with the experiment is achieved when the implicit salt concentration is set equal or close to the experimental conditions.

Computational modeling of protein conformational changes - Application to the opening SARS-CoV-2 spike

We present a new approach to compute and analyze the dynamical electro-geometric properties of proteins undergoing conformational changes. The molecular trajectory is obtained from Markov state models, and the electrostatic potential is calculated using the continuum Poisson-Boltzmann equation. The numerical electric potential is constructed using a parallel sharp numerical solver implemented on adaptive Octree grids. We introduce novel a posteriori error estimates to quantify the solution's accuracy on the molecular surface. To illustrate the approach, we consider the opening of the SARS-CoV-2 spike protein using the recent molecular trajectory simulated through the Folding@home initiative. We analyze our results, focusing on the characteristics of the receptor-binding domain and its vicinity. This work lays the foundation for a new class of hybrid computational approaches, producing high-fidelity dynamical computational measurements serving as a basis for protein bio-mechanism investigations.

Computational Design of PDZ-Peptide Binding

This chapter describes two computational methods for PDZ-peptide binding: high-throughput computational protein design (CPD) and a medium-throughput approach combining molecular dynamics for conformational sampling with a Poisson-Boltzmann (PB) Linear Interaction Energy for scoring. A new CPD method is outlined, which uses adaptive Monte Carlo simulations to efficiently sample peptide variants that tightly bind a PDZ domain, and provides at the same time precise estimates of their relative binding free energies. A detailed protocol is described based on the Proteus CPD software. The medium-throughput approach can be performed with standard MD and PB software, such as NAMD and Charmm. For 40 complexes between Tiam1 and peptide ligands, it gave high a2ccuracy, with mean errors of around 0.5 kcal/mol for relative binding free energies and no large errors. It requires a moderate amount of parameter fitting before it can be applied, and its transferability to other protein families is still untested.

Multidimensional Global Optimization and Robustness Analysis in the Context of Protein-Ligand Binding

Accuracy of protein-ligand binding free energy calculations utilizing implicit solvent models is critically affected by parameters of the underlying dielectric boundary, specifically, the atomic and water probe radii. Here, a global multidimensional optimization pipeline is developed to find optimal atomic radii specifically for protein-ligand binding calculations in implicit solvent. The computational pipeline has these three key components: (1) a massively parallel implementation of a deterministic global optimization algorithm (VTDIRECT95), (2) an accurate yet reasonably fast generalized Born implicit solvent model (GBNSR6), and (3) a novel robustness metric that helps distinguish between nearly degenerate local minima via a postprocessing step of the optimization. A graph-based "kT-connectivity" approach to explore and visualize the multidimensional energy landscape is proposed: local minima that can be reached from the global minimum without exceeding a given energy threshold (kT) are considered to be connected. As an illustration of the capabilities of the optimization pipeline, we apply it to find a global optimum in the space of just five radii: four atomic (O, H, N, and C) radii and water probe radius. The optimized radii, ρW = 1.37 Å, ρC = 1.40 Å, ρH = 1.55 Å, ρN = 2.35 Å, and ρO = 1.28 Å, lead to a closer agreement of electrostatic binding free energies with the explicit solvent reference than two commonly used sets of radii previously optimized for small molecules. At the same time, the ability of the optimizer to find the global optimum reveals fundamental limits of the common two-dielectric implicit solvation model: the computed electrostatic binding free energies are still almost 4 kcal/mol away from the explicit solvent reference. The proposed computational approach opens the possibility to further improve the accuracy of practical computational protocols for binding free energy calculations.

Accuracy Comparison of Generalized Born Models in the Calculation of Electrostatic Binding Free Energies

The need for accurate yet efficient representation of the aqueous environment in biomolecular modeling has led to the development of a variety of generalized Born (GB) implicit solvent models. While many studies have focused on the accuracy of available GB models in predicting solvation free energies, a systematic assessment of the quality of these models in binding free energy calculations, crucial for rational drug design, has not been undertaken. Here, we evaluate the accuracies of eight common GB flavors (GB-HCT, GB-OBC, GB-neck2, GBNSR6, GBSW, GBMV1, GBMV2, and GBMV3), available in major molecular dynamics packages, in predicting the electrostatic binding free energies ( ΔΔ G_el) for a diverse set of 60 biomolecular complexes belonging to four main classes: protein-protein, protein-drug, RNA-peptide, and small complexes. The GB flavors are examined in terms of their ability to reproduce the results from the Poisson-Boltzmann (PB) model, commonly used as accuracy reference in this context. We show that the agreement with the PB of ΔΔ G_el estimates varies widely between different GB models and also across different types of biomolecular complexes, with R² correlations ranging from 0.3772 to 0.9986. A surface-based "R6" GB model recently implemented in AMBER shows the closest overall agreement with reference PB ( R² = 0.9949, RMSD = 8.75 kcal/mol). The RNA-peptide and protein-drug complex sets appear to be most challenging for all but one model, as indicated by the large deviations from the PB in ΔΔ G_el. Small neutral complexes present the least challenge for most of the GB models tested. The quantitative demonstration of the strengths and weaknesses of the GB models across the diverse complex types provided here can be used as a guide for practical computations and future development efforts.

Grid-Based Surface Generalized Born Model for Calculation of Electrostatic Binding Free Energies

Fast and accurate calculation of solvation free energies is central to many applications, such as rational drug design. In this study, we present a grid-based molecular surface implementation of "R6" flavor of the generalized Born (GB) implicit solvent model, named GBNSR6. The speed, accuracy relative to numerical Poisson-Boltzmann treatment, and sensitivity to grid surface parameters are tested on a set of 15 small protein-ligand complexes and a set of biomolecules in the range of 268 to 25099 atoms. Our results demonstrate that the proposed model provides a relatively successful compromise between the speed and accuracy of computing polar components of the solvation free energies (ΔGpol) and binding free energies (ΔΔGpol). The model tolerates a relatively coarse grid size h = 0.5 Å, where the grid artifact error in computing ΔΔGpol remains in the range of kBT ∼ 0.6 kcal/mol. The estimated ΔΔGpols are well correlated (r2 = 0.97) with the numerical Poisson-Boltzmann reference, while showing virtually no systematic bias and RMSE = 1.43 kcal/mol. The grid-based GBNSR6 model is available in Amber (AmberTools) package of molecular simulation programs.

A Simple PB/LIE Free Energy Function Accurately Predicts the Peptide Binding Specificity of the Tiam1 PDZ Domain

PDZ domains generally bind short amino acid sequences at the C-terminus of target proteins, and short peptides can be used as inhibitors or model ligands. Here, we used experimental binding assays and molecular dynamics simulations to characterize 51 complexes involving the Tiam1 PDZ domain and to test the performance of a semi-empirical free energy function. The free energy function combined a Poisson-Boltzmann (PB) continuum electrostatic term, a van der Waals interaction energy, and a surface area term. Each term was empirically weighted, giving a Linear Interaction Energy or “PB/LIE” free energy. The model yielded a mean unsigned deviation of 0.43 kcal/mol and a Pearson correlation of 0.64 between experimental and computed free energies, which was superior to a Null model that assumes all complexes have the same affinity. Analyses of the models support several experimental observations that indicate the orientation of the α₂ helix is a critical determinant for peptide specificity. The models were also used to predict binding free energies for nine new variants, corresponding to point mutants of the Syndecan1 and Caspr4 peptides. The predictions did not reveal improved binding; however, they suggest that an unnatural amino acid could be used to increase protease resistance and peptide lifetimes in vivo. The overall performance of the model should allow its use in the design of new PDZ ligands in the future.

Electrostatic component of binding energy: Interpreting predictions from poisson-boltzmann equation and modeling protocols

Macromolecular interactions are essential for understanding numerous biological processes and are typically characterized by the binding free energy. Important component of the binding free energy is the electrostatics, which is frequently modeled via the solutions of the Poisson-Boltzmann Equations (PBE). However, numerous works have shown that the electrostatic component (ΔΔGelec ) of binding free energy is very sensitive to the parameters used and modeling protocol. This prompted some researchers to question the robustness of PBE in predicting ΔΔGelec . We argue that the sensitivity of the absolute ΔΔGelec calculated with PBE using different input parameters and definitions does not indicate PBE deficiency, rather this is what should be expected. We show how the apparent sensitivity should be interpreted in terms of the underlying changes in several numerous and physical parameters. We demonstrate that PBE approach is robust within each considered force field (CHARMM-27, AMBER-94, and OPLS-AA) once the corresponding structures are energy minimized. This observation holds despite of using two different molecular surface definitions, pointing again that PBE delivers consistent results within particular force field. The fact that PBE delivered ΔΔGelec values may differ if calculated with different modeling protocols is not a deficiency of PBE, but natural results of the differences of the force field parameters and potential functions for energy minimization. In addition, while the absolute ΔΔGelec values calculated with different force field differ, their ordering remains practically the same allowing for consistent ranking despite of the force field used. © 2016 Wiley Periodicals, Inc.

Reconciling the understanding of 'hydrophobicity' with physics-based models of proteins

The idea that a 'hydrophobic energy' drives protein folding, aggregation, and binding by favoring the sequestration of bulky residues from water into the protein interior is widespread. The solvation free energies (ΔGsolv) of small nonpolar solutes increase with surface area (A), and the free energies of creating macroscopic cavities in water increase linearly with A. These observations seem to imply that there is a hydrophobic component (ΔGhyd) of ΔGsolv that increases linearly with A, and this assumption is widely used in implicit solvent models. However, some explicit-solvent molecular dynamics studies appear to contradict these ideas. For example, one definition (ΔG(LJ)) of ΔGhyd is that it is the free energy of turning on the Lennard-Jones (LJ) interactions between the solute and solvent. However, ΔG(LJ) decreases with A for alanine and glycine peptides. Here we argue that these apparent contradictions can be reconciled by defining ΔGhyd to be a near hard core insertion energy (ΔGrep), as in the partitioning proposed by Weeks, Chandler, and Andersen. However, recent results have shown that ΔGrep is not a simple function of geometric properties of the molecule, such as A and the molecular volume, and that the free energy of turning on the attractive part of the LJ potential cannot be computed from first-order perturbation theory for proteins. The theories that have been developed from these assumptions to predict ΔGhyd are therefore inadequate for proteins.

Examining the assumptions underlying continuum-solvent models

Continuum-solvent models (CSMs) have successfully predicted many quantities, including the solvation-free energies (ΔG) of small molecules, but they have not consistently succeeded at reproducing experimental binding free energies (ΔΔG), especially for protein-protein complexes. Several CSMs break ΔG into the free energy (ΔGvdw) of inserting an uncharged molecule into solution and the free energy (ΔGel) gained from charging. Some further divide ΔGvdw into the free energy (ΔGrep) of inserting a nearly hard cavity into solution and the free energy (ΔGatt) gained from turning on dispersive interactions between the solute and solvent. We show that for 9 protein-protein complexes neither ΔGrep nor ΔGvdw was linear in the solvent-accessible area A, as assumed in many CSMs, and the corresponding components of ΔΔG were not linear in changes in A. We show that linear response theory (LRT) yielded good estimates of ΔGatt and ΔΔGatt, but estimates of ΔΔGatt obtained from either the initial or final configurations of the solvent were not consistent with those from LRT. The LRT estimates of ΔGel differed by more than 100 kcal/mol from the explicit solvent model's (ESM's) predictions, and its estimates of the corresponding component (ΔΔGel) of ΔΔG differed by more than 10 kcal/mol. Finally, the Poisson-Boltzmann equation produced estimates of ΔGel that were correlated with those from the ESM, but its estimates of ΔΔGel were much less so. These findings may help explain why many CSMs have not been consistently successful at predicting ΔΔG for many complexes, including protein-protein complexes.

Protein-Ligand Electrostatic Binding Free Energies from Explicit and Implicit Solvation

Accurate yet efficient computational models of solvent environment are central for most calculations that rely on atomistic modeling, such as prediction of protein-ligand binding affinities. In this study, we evaluate the accuracy of a recently developed generalized Born implicit solvent model, GBNSR6 (Aguilar et al. J. Chem. Theory Comput. 2010, 6, 3613-3639), in estimating the electrostatic solvation free energies (ΔG(pol)) and binding free energies (ΔΔG(pol)) for small protein-ligand complexes. We also compare estimates based on three different explicit solvent models (TIP3P, TIP4PEw, and OPC). The two main findings are as follows. First, the deviation (RMSD = 7.04 kcal/mol) of GBNSR6 binding affinities from commonly used TIP3P reference values is comparable to the deviations between explicit models themselves, e.g. TIP4PEw vs TIP3P (RMSD = 5.30 kcal/mol). A simple uniform adjustment of the atomic radii by a single scaling factor reduces the RMS deviation of GBNSR6 from TIP3P to within the above "error margin" - differences between ΔΔG(pol) estimated by different common explicit solvent models. The simple radii scaling virtually eliminates the systematic deviation (ΔΔG(pol)) between GBNSR6 and two out of the three explicit water models and significantly reduces the deviation from the third explicit model. Second, the differences between electrostatic binding energy estimates from different explicit models is disturbingly large; for example, the deviation between TIP4PEw and TIP3P estimates of ΔΔG(pol) values can be up to ∼50% or ∼9 kcal/mol, which is significantly larger than the "chemical accuracy" goal of ∼1 kcal/mol. The absolute ΔG(pol) calculated with different explicit models could differ by tens of kcal/mol. These discrepancies point to unacceptably high sensitivity of binding affinity estimates to the choice of common explicit water models. The absence of a clear "gold standard" among these models strengthens the case for the use of accurate implicit solvation models for binding energetics, which may be orders of magnitude faster.