Coupling finite and boundary element methods to solve the Poisson-Boltzmann equation for electrostatics in molecular solvation.

The Poisson-Boltzmann equation is widely used to model electrostatics in molecular systems. Available software packages solve it using finite difference, finite element, and boundary element methods, where the latter is attractive due to the accurate representation of the molecular surface and partial charges, and exact enforcement of the boundary conditions at infinity. However, the boundary element method is limited to linear equations and piecewise constant variations of the material properties. In this work, we present a scheme that couples finite and boundary elements for the linearised Poisson-Boltzmann equation, where the finite element method is applied in a confined solute region and the boundary element method in the external solvent region. As a proof-of-concept exercise, we use the simplest methods available: Johnson-Nédélec coupling with mass matrix and diagonal preconditioning, implemented using the Bempp-cl and FEniCSx libraries via their Python interfaces. We showcase our implementation by computing the polar component of the solvation free energy of a set of molecules using a constant and a Gaussian-varying permittivity. As validation, we compare against well-established finite difference solvers for an extensive binding energy data set, and with the finite difference code APBS (to 0.5%) for Gaussian permittivities. We also show scaling results from protein G B1 (955 atoms) up to immunoglobulin G (20,148 atoms). For small problems, the coupled method was efficient, outperforming a purely boundary integral approach. For Gaussian-varying permittivities, which are beyond the applicability of boundary elements alone, we were able to run medium to large-sized problems on a single workstation. The development of better preconditioning techniques and the use of distributed memory parallelism for larger systems remains an area for future work. We hope this work will serve as inspiration for future developments that consider space-varying field parameters, and mixed linear-nonlinear schemes for molecular electrostatics with implicit solvent models.

Monovalent cation binding to model systems and the macrocyclic depsipeptide, emodepside.

This study focuses on the systematic exploration of the emodepside conformations bound to monovalent K+ ion using quantum mechanical density functional theory (DFT) calculations at the M06-2X/6-31+G(d,p) level of theory. Nine conformers of emodepside and their complexes with K+ ion were characterized as stationary points on the potential energy surface. The conformational isomers were examined for their 3D structures, bonding, energetics, and interactions with the cation. A cavitand-like structure (CC) is identified to be the energetically most stable arrangement. To arrive at a better understanding of the K+ ion binding, calculations were initially performed on complexes formed by the K+ and Na+ ions with model ligands (methyl ester and N,N-dimethyl acetamide). Both the natural bond orbital (NBO) method and the block-localized wavefunction (BLW) energy decomposition approach was employed to assess the bonding and energetic contributions stabilizing the ion-bound model complexes. Finally, the solvent effect was evaluated through complete geometry optimizations and energy minimizations for the model ion-ligand complexes and the emodepside-K+ bound complexes using an implicit solvent model mimicking water and DMSO.

Prediction of Strong Solvatochromism in a Molecular Photocatalyst.

Based on quantum chemical calculations, we predict strong solvatochromism in a light-driven molecular photocatalyst for hydrogen generation, that is we show that the electronic and optical properties of the photocatalyst strongly depend on the solvent it is dissolved in. Our calculations in particular indicate a solvent-dependent relocation of the highest occupied molecular orbital (HOMO). Ground-state density functional theory and linear response time-dependent density functional theory calculations were applied in order to investigate the influence of implicit solvents on the structural, electronic and optical properties of a molecular photocatalyst. Only at high dielectric constants of the solvent, is the HOMO located at the metal center of the photosensitizer, whereas at low dielectric constants the HOMO is centered at the metal atom of the catalytically active complex. We elucidate the electronic origins of this strong solvatochromic effect and sketch the consequences of these insights for the use of photocatalysts in different environments.

Comparing the influence of explicit and implicit solvation models on site-specific thermodynamic stability of proteins.

Understanding the molecular basis for protein stability requires a thermodynamic analysis of protein folding. Thermodynamic analysis is often performed by sampling many atomistic conformations using molecular simulations that employ either explicit or implicit water models. However, it remains unclear to what extent thermodynamic results from different solvation models are reliable at the molecular level. In this study, we quantify the influence of both solvation models on folding stability at the individual backbone and side chain resolutions. We assess the residue-specific folding free energy components of a ß-sheet protein and a helical protein using trajectories resulting from TIP3P explicit and generalized Born/surface area implicit solvent simulations of model proteins. We found that the thermodynamic discrepancy due to the implicit solvent mostly originates from charged side chains, followed by the under-stabilized hydrophobic ones. In contrast, the contributions of backbone residue in both proteins were comparable for explicit and implicit water models. Our study lays out the foundation for detailed thermodynamic assessment of solvation models in the context of protein simulation.

Short solvent model for ion correlations and hydrophobic association.

Coulomb interactions play a major role in determining the thermodynamics, structure, and dynamics of condensed-phase systems, but often present significant challenges. Computer simulations usually use periodic boundary conditions to minimize corrections from finite cell boundaries but the long range of the Coulomb interactions generates significant contributions from distant periodic images of the simulation cell, usually calculated by Ewald sum techniques. This can add significant overhead to computer simulations and hampers the development of intuitive local pictures and simple analytic theory. In this paper, we present a general framework based on local molecular field theory to accurately determine the contributions from long-ranged Coulomb interactions to the potential of mean force between ionic or apolar hydrophobic solutes in dilute aqueous solutions described by standard classical point charge water models. The simplest approximation leads to a short solvent (SS) model, with truncated solvent-solvent and solute-solvent Coulomb interactions and long-ranged but screened Coulomb interactions only between charged solutes. The SS model accurately describes the interplay between strong short-ranged solute core interactions, local hydrogen-bond configurations, and long-ranged dielectric screening of distant charges, competing effects that are difficult to capture in standard implicit solvent models.

The Oxygen Evolution Reaction at MoS2 Edge Sites: The Role of a Solvent Environment in DFT-Based Molecular Simulations.

Density functional theory (DFT) calculations are employed to study the oxygen evolution reaction (OER) on the edges of stripes of monolayer molybdenum disulfide. Experimentally, this material has been shown to evolve oxygen, albeit with low efficiency. Previous DFT studies have traced this low catalytic performance to the unfavourable adsorption energies of some reaction intermediates on the MoS2 edge sites. In this work, we study the effects of the aqueous liquid surrounding the active sites. A computational approach is used, where the solvent is modeled as a continuous medium providing a dielectric embedding of the catalyst and the reaction intermediates. A description at this level of theory can have a profound impact on the studied reactions: the calculated overpotential for the OER is lowered from 1.15 eV to 0.77 eV. It is shown that such variations in the reaction energetics are linked to the polar nature of the adsorbed intermediates, which leads to changes in the calculated electronic charge density when surrounded by water. These results underline the necessity to computationally account for solvation effects, especially in aqueous environments and when highly polar intermediates are present.

Identifying Systematic Force Field Errors Using a 3D-RISM Element Counting Correction.

Hydration free energies of small molecules are commonly used as benchmarks for solvation models. However, errors in predicting hydration free energies are partially due to the force fields used and not just the solvation model. To address this, we have used the 3D reference interaction site model (3D-RISM) of molecular solvation and existing benchmark explicit solvent calculations with a simple element count correction (ECC) to identify problems with the non-bond parameters in the general AMBER force field (GAFF). 3D-RISM was used to calculate hydration free energies of all 642 molecules in the FreeSolv database, and a partial molar volume correction (PMVC), ECC, and their combination (PMVECC) were applied to the results. The PMVECC produced a mean unsigned error of 1.01±0.04kcal/mol and root mean squared error of 1.44±0.07kcal/mol, better than the benchmark explicit solvent calculations from FreeSolv, and required less than 15 s of computing time per molecule on a single CPU core. Importantly, parameters for PMVECC showed systematic errors for molecules containing Cl, Br, I, and P. Applying ECC to the explicit solvent hydration free energies found the same systematic errors. The results strongly suggest that some small adjustments to the Lennard-Jones parameters for GAFF will lead to improved hydration free energy calculations for all solvent models.

Accelerating the 3D reference interaction site model theory of molecular solvation with treecode summation and cut-offs.

The 3D reference interaction site model (3D-RISM) of molecular solvation is a powerful tool for computing the equilibrium thermodynamics and density distributions of solvents, such as water and co-ions, around solute molecules. However, 3D-RISM solutions can be expensive to calculate, especially for proteins and other large molecules where calculating the potential energy between solute and solvent requires more than half the computation time. To address this problem, we have developed and implemented treecode summation for long-range interactions and analytically corrected cut-offs for short-range interactions to accelerate the potential energy and long-range asymptotics calculations in non-periodic 3D-RISM in the AmberTools molecular modeling suite. For the largest single protein considered in this work, tubulin, the total computation time was reduced by a factor of 4. In addition, parallel calculations with these new methods scale almost linearly and the iterative solver remains the largest impediment to parallel scaling. To demonstrate the utility of our approach for large systems, we used 3D-RISM to calculate the solvation thermodynamics and density distribution of 7-ring microtubule, consisting of 910 tubulin dimers, over 1.2 million atoms.

Towards optimal boundary integral formulations of the Poisson-Boltzmann equation for molecular electrostatics.

The Poisson-Boltzmann equation offers an efficient way to study electrostatics in molecular settings. Its numerical solution with the boundary element method is widely used, as the complicated molecular surface is accurately represented by the mesh, and the point charges are accounted for explicitly. In fact, there are several well-known boundary integral formulations available in the literature. This work presents a generalized expression of the boundary integral representation of the implicit solvent model, giving rise to new forms to compute the electrostatic potential. Moreover, it proposes a strategy to build efficient preconditioners for any of the resulting systems, improving the convergence of the linear solver. We perform systematic benchmarking of a set of formulations and preconditioners, focusing on the time to solution, matrix conditioning, and eigenvalue spectrum. We see that the eigenvalue clustering is a good indicator of the matrix conditioning, and show that they can be easily manipulated by scaling the preconditioner. Our results suggest that the optimal choice is problem-size dependent, where a simpler direct formulation is the fastest for small molecules, but more involved second-kind equations are better for larger problems. We also present a fast Calderón preconditioner for first-kind formulations, which shows promising behavior for future analysis. This work sets the basis towards choosing the most convenient boundary integral formulation of the Poisson-Boltzmann equation for a given problem.

Variational implicit-solvent predictions of the dry-wet transition pathways for ligand-receptor binding and unbinding kinetics.

Ligand-receptor binding and unbinding are fundamental biomolecular processes and particularly essential to drug efficacy. Environmental water fluctuations, however, impact the corresponding thermodynamics and kinetics and thereby challenge theoretical descriptions. Here, we devise a holistic, implicit-solvent, multimethod approach to predict the (un)binding kinetics for a generic ligand-pocket model. We use the variational implicit-solvent model (VISM) to calculate the solute-solvent interfacial structures and the corresponding free energies, and combine the VISM with the string method to obtain the minimum energy paths and transition states between the various metastable ("dry" and "wet") hydration states. The resulting dry-wet transition rates are then used in a spatially dependent multistate continuous-time Markov chain Brownian dynamics simulation and the related Fokker-Planck equation calculations of the ligand stochastic motion, providing the mean first-passage times for binding and unbinding. We find the hydration transitions to significantly slow down the binding process, in semiquantitative agreement with existing explicit-water simulations, but significantly accelerate the unbinding process. Moreover, our methods allow the characterization of nonequilibrium hydration states of pocket and ligand during the ligand movement, for which we find substantial memory and hysteresis effects for binding vs. unbinding. Our study thus provides a significant step forward toward efficient, physics-based interpretation and predictions of the complex kinetics in realistic ligand-receptor systems.

Modulation of free energy landscapes as a strategy for the design of antimicrobial peptides.

Computational design of antimicrobial peptides (AMPs) is a promising area of research for developing novel agents against drug-resistant bacteria. AMPs are present naturally in many organisms, from bacteria to humans, a time-tested mechanism that makes them attractive as effective antibiotics. Depending on the environment, AMPs can exhibit α-helical or ß-sheet conformations, a mix of both, or lack secondary structure; they can be linear or cyclic. Prediction of their structures is challenging but critical for rational design. Promising AMP leads can be developed using essentially two approaches: traditional modeling of the physicochemical mechanisms that determine peptide behavior in aqueous and membrane environments and knowledge-based, e.g., machine learning (ML) techniques, that exploit ever-growing AMP databases. Here, we explore the conformational landscapes of two recently ML-designed AMPs, characterize the dependence of these landscapes on the medium conditions, and identify features in peptide and membrane landscapes that mediate protein-membrane association. For both peptides, we observe greater conformational diversity in an aqueous solvent than in a less polar solvent, and one peptide is seen to alter its conformation more dramatically than the other upon the change of solvent. Our results support the view that structural rearrangement in response to environmental changes is central to the mechanism of membrane-structure disruption by linear peptides. We expect that the design of AMPs by ML will benefit from the incorporation of peptide conformational substates as quantified here with molecular simulations.

Efficient mesh refinement for the Poisson-Boltzmann equation with boundary elements.

The Poisson-Boltzmann equation is a widely used model to study electrostatics in molecular solvation. Its numerical solution using a boundary integral formulation requires a mesh on the molecular surface only, yielding accurate representations of the solute, which is usually a complicated geometry. Here, we utilize adjoint-based analyses to form two goal-oriented error estimates that allow us to determine the contribution of each discretization element (panel) to the numerical error in the solvation free energy. This information is useful to identify high-error panels to then refine them adaptively to find optimal surface meshes. We present results for spheres and real molecular geometries, and see that elements with large error tend to be in regions where there is a high electrostatic potential. We also find that even though both estimates predict different total errors, they have similar performance as part of an adaptive mesh refinement scheme. Our test cases suggest that the adaptive mesh refinement scheme is very effective, as we are able to reduce the error one order of magnitude by increasing the mesh size less than 20% and come out to be more efficient than uniform refinement when computing error estimations. This result sets the basis toward efficient automatic mesh refinement schemes that produce optimal meshes for solvation energy calculations.

ARGOS: An adaptive refinement goal-oriented solver for the linearized Poisson-Boltzmann equation.

An adaptive finite element solver for the numerical calculation of the electrostatic coupling between molecules in a solvent environment is developed and tested. At the heart of the solver is a goal-oriented a posteriori error estimate for the electrostatic coupling, derived and implemented in the present work, that gives rise to an orders of magnitude improved precision and a shorter computational time as compared to standard finite difference solvers. The accuracy of the new solver ARGOS is evaluated by numerical experiments on a series of problems with analytically known solutions. In addition, the solver is used to calculate electrostatic couplings between two chromophores, linked to polyproline helices of different lengths and between the spike protein of SARS-CoV-2 and the ACE2 receptor. All the calculations are repeated by using the well-known finite difference solvers MEAD and APBS, revealing the advantages of the present finite element solver.

Rationalizing the Mechanism of Peroxyformate Decomposition: Computational Insights To Understand Solvent Influence.

The heterolytic decomposition of tert-butyl peroxyformate to tert-butanol and carbon dioxide, catalyzed by pyridine, is a long-known example of a reaction whose kinetics are strongly affected by solvent polarity. From DFT and ab initio methods together with the SMD implicit solvation model, an extension on the formerly accepted mechanism is proposed. This novel proposal involves the formation of a carbonic acid ester intermediate and its further decomposition, through an unreported pyridine-mediated stepwise route. Computed barriers for this mechanism at DLPNO/CCSD(T)-def2-TZVP are in excellent agreement with experimental kinetic data across different solvents. Furthermore, the strong relationships between activation energies, geometric parameters in the transition state and the characteristics of the different solvents are also analyzed in depth.

An Effective MM/GBSA Protocol for Absolute Binding Free Energy Calculations: A Case Study on SARS-CoV-2 Spike Protein and the Human ACE2 Receptor.

The binding free energy calculation of protein-ligand complexes is necessary for research into virus-host interactions and the relevant applications in drug discovery. However, many current computational methods of such calculations are either inefficient or inaccurate in practice. Utilizing implicit solvent models in the molecular mechanics generalized Born surface area (MM/GBSA) framework allows for efficient calculations without significant loss of accuracy. Here, GBNSR6, a new flavor of the generalized Born model, is employed in the MM/GBSA framework for measuring the binding affinity between SARS-CoV-2 spike protein and the human ACE2 receptor. A computational protocol is developed based on the widely studied Ras-Raf complex, which has similar binding free energy to SARS-CoV-2/ACE2. Two options for representing the dielectric boundary of the complexes are evaluated: one based on the standard Bondi radii and the other based on a newly developed set of atomic radii (OPT1), optimized specifically for protein-ligand binding. Predictions based on the two radii sets provide upper and lower bounds on the experimental references: -14.7(ΔGbindBondi)<-10.6(ΔGbindExp.)<-4.1(ΔGbindOPT1) kcal/mol. The consensus estimates of the two bounds show quantitative agreement with the experiment values. This work also presents a novel truncation method and computational strategies for efficient entropy calculations with normal mode analysis. Interestingly, it is observed that a significant decrease in the number of snapshots does not affect the accuracy of entropy calculation, while it does lower computation time appreciably. The proposed MM/GBSA protocol can be used to study the binding mechanism of new variants of SARS-CoV-2, as well as other relevant structures.

Prediction of protein-protein complexes using replica exchange with repulsive scaling.

The realistic prediction of protein-protein complex structures is import to ultimately model the interaction of all proteins in a cell and for the design of new protein-protein interactions. In principle, molecular dynamics (MD) simulations allow one to follow the association process under realistic conditions including full partner flexibility and surrounding solvent. However, due to the many local binding energy minima at the surface of protein partners, MD simulations are frequently trapped for long times in transient association states. We have designed a replica-exchange based scheme employing different levels of a repulsive biasing between partners in each replica simulation. The bias acts only on intermolecular interactions based on an increase in effective pairwise van der Waals radii (repulsive scaling (RS)-REMD) without affecting interactions within each protein or with the solvent. For a set of five protein test cases (out of six) the RS-REMD technique allowed the sampling of near-native complex structures even when starting from the opposide site with respect to the native binding site for one partner. Using the same start structures and same computational demand regular MD simulations sampled near native complex structures only for one case. The method showed also improved results for the refinement of docked structures in the vicinity of the native binding geometry compared to regular MD refinement.

A Boundary-Integral Approach for the Poisson-Boltzmann Equation with Polarizable Force Fields.

Implicit-solvent models are widely used to study the electrostatics in dissolved biomolecules, which are parameterized using force fields. Standard force fields treat the charge distribution with point charges; however, other force fields have emerged which offer a more realistic description by considering polarizability. In this work, we present the implementation of the polarizable and multipolar force field atomic multipole optimized energetics for biomolecular applications (AMOEBA), in the boundary integral Poisson-Boltzmann solver PyGBe. Previous work from other researchers coupled AMOEBA with the finite-difference solver APBS, and found difficulties to effectively transfer the multipolar charge description to the mesh. A boundary integral formulation treats the charge distribution analytically, overlooking such limitations. This becomes particularly important in simulations that need high accuracy, for example, when the quantity of interest is the difference between solvation energies obtained from separate calculations, like happens for binding energy. We present verification and validation results of our software, compare it with the implementation on APBS, and assess the efficiency of AMOEBA and classical point-charge force fields in a Poisson-Boltzmann solver. We found that a boundary integral approach performs similarly to a volumetric method on CPU. Also, we present a GPU implementation of our solver. Moreover, with a boundary element method, the mesh density to correctly resolve the electrostatic potential is the same for standard point-charge and multipolar force fields. Finally, we saw that for binding energy calculations, a boundary integral approach presents more consistent results than a finite difference approximation for multipolar force fields. © 2019 Wiley Periodicals, Inc.

Unravelling the covalent binding of zampanolide and taccalonolide AJ to a minimalist representation of a human microtubule.

Many natural products target mammalian tubulin but only a few can form a covalent bond and hence irreversibly affect microtubule function. Among them, zampanolide (ZMP) and taccalonolide AJ (TAJ) stand out, not only because they are very potent antitumor agents but also because the adducts they form with ß-tubulin have been structurally characterized in atomic detail. By applying model building techniques, molecular orbital calculations, molecular dynamics simulations and hybrid QM/MM methods, we have gained insight into the 1,2- and 1,4-addition reactions of His229 and Asp226 to ZMP and TAJ, respectively, in the taxane-binding site of ß-tubulin. The experimentally inaccessible precovalent complexes strongly suggest a water-mediated proton shuttle mechanism for ZMP adduct formation and a direct nucleophilic attack by the carboxylate of Asp226 on C22 of the C22R,C23R epoxide in TAJ. The M-loop, which is crucially important for interprotofilament interactions, is structured into a short helix in both types of complexes, mostly as a consequence of the fixation of the phenol ring of Tyr283 and the guanidinium of Arg284. As a side benefit, we obtained evidence supporting the existence of a commonly neglected intramolecular disulfide bond between Cys241 and Cys356 in ß-tubulin that contributes to protein compactness and is absent in the ßIII isotype associated with resistance to taxanes and other drugs.

Breaking the polar-nonpolar division in solvation free energy prediction.

Implicit solvent models divide solvation free energies into polar and nonpolar additive contributions, whereas polar and nonpolar interactions are inseparable and nonadditive. We present a feature functional theory (FFT) framework to break this ad hoc division. The essential ideas of FFT are as follows: (i) representability assumption: there exists a microscopic feature vector that can uniquely characterize and distinguish one molecule from another; (ii) feature-function relationship assumption: the macroscopic features, including solvation free energy, of a molecule is a functional of microscopic feature vectors; and (iii) similarity assumption: molecules with similar microscopic features have similar macroscopic properties, such as solvation free energies. Based on these assumptions, solvation free energy prediction is carried out in the following protocol. First, we construct a molecular microscopic feature vector that is efficient in characterizing the solvation process using quantum mechanics and Poisson-Boltzmann theory. Microscopic feature vectors are combined with macroscopic features, that is, physical observable, to form extended feature vectors. Additionally, we partition a solvation dataset into queries according to molecular compositions. Moreover, for each target molecule, we adopt a machine learning algorithm for its nearest neighbor search, based on the selected microscopic feature vectors. Finally, from the extended feature vectors of obtained nearest neighbors, we construct a functional of solvation free energy, which is employed to predict the solvation free energy of the target molecule. The proposed FFT model has been extensively validated via a large dataset of 668 molecules. The leave-one-out test gives an optimal root-mean-square error (RMSE) of 1.05 kcal/mol. FFT predictions of SAMPL0, SAMPL1, SAMPL2, SAMPL3, and SAMPL4 challenge sets deliver the RMSEs of 0.61, 1.86, 1.64, 0.86, and 1.14 kcal/mol, respectively. Using a test set of 94 molecules and its associated training set, the present approach was carefully compared with a classic solvation model based on weighted solvent accessible surface area. © 2017 Wiley Periodicals, Inc.

Combining the polarizable Drude force field with a continuum electrostatic Poisson-Boltzmann implicit solvation model.

In this work, we have combined the polarizable force field based on the classical Drude oscillator with a continuum Poisson-Boltzmann/solvent-accessible surface area (PB/SASA) model. In practice, the positions of the Drude particles experiencing the solvent reaction field arising from the fixed charges and induced polarization of the solute must be optimized in a self-consistent manner. Here, we parameterized the model to reproduce experimental solvation free energies of a set of small molecules. The model reproduces well-experimental solvation free energies of 70 molecules, yielding a root mean square difference of 0.8 kcal/mol versus 2.5 kcal/mol for the CHARMM36 additive force field. The polarization work associated with the solute transfer from the gas-phase to the polar solvent, a term neglected in the framework of additive force fields, was found to make a large contribution to the total solvation free energy, comparable to the polar solute-solvent solvation contribution. The Drude PB/SASA also reproduces well the electronic polarization from the explicit solvent simulations of a small protein, BPTI. Model validation was based on comparisons with the experimental relative binding free energies of 371 single alanine mutations. With the Drude PB/SASA model the root mean square deviation between the predicted and experimental relative binding free energies is 3.35 kcal/mol, lower than 5.11 kcal/mol computed with the CHARMM36 additive force field. Overall, the results indicate that the main limitation of the Drude PB/SASA model is the inability of the SASA term to accurately capture non-polar solvation effects. © 2018 Wiley Periodicals, Inc.