Computational scheme for pH-dependent binding free energy calculation with explicit solvent

Constant pH molecular dynamics simulations: Current status and recent applications

Many important protein functions are carried out through proton-coupled conformational dynamics. Thus, the ability to accurately model protonation states dynamically has wide-ranging implications. Over the past two decades, two main types of constant pH methods (discrete and continuous) have been developed to enable proton-coupled molecular dynamics (MD) simulations. In this short review, we discuss the current status of the development and highlight recent applications that have advanced our understanding of protein structure-function relationships. We conclude the review by outlining the remaining challenges in the method development and projecting important areas for future applications.

An overview of the SAMPL8 host-guest binding challenge

The SAMPL series of challenges aim to focus the community on specific modeling challenges, while testing and hopefully driving progress of computational methods to help guide pharmaceutical drug discovery. In this study, we report on the results of the SAMPL8 host–guest blind challenge for predicting absolute binding affinities. SAMPL8 focused on two host–guest datasets, one involving the cucurbituril CB8 (with a series of common drugs of abuse) and another involving two different Gibb deep-cavity cavitands. The latter dataset involved a previously featured deep cavity cavitand (TEMOA) as well as a new variant (TEETOA), both binding to a series of relatively rigid fragment-like guests. Challenge participants employed a reasonably wide variety of methods, though many of these were based on molecular simulations, and predictive accuracy was mixed. As in some previous SAMPL iterations (SAMPL6 and SAMPL7), we found that one approach to achieve greater accuracy was to apply empirical corrections to the binding free energy predictions, taking advantage of prior data on binding to these hosts. Another approach which performed well was a hybrid MD-based approach with reweighting to a force matched QM potential. In the cavitand challenge, an alchemical method using the AMOEBA-polarizable force field achieved the best success with RMSE less than 1 kcal/mol, while another alchemical approach (ATM/GAFF2-AM1BCC/TIP3P/HREM) had RMSE less than 1.75 kcal/mol. The work discussed here also highlights several important lessons; for example, retrospective studies of reference calculations demonstrate the sensitivity of predicted binding free energies to ethyl group sampling and/or guest starting pose, providing guidance to help improve future studies on these systems.

O-Acetyl Migration within the Sialic Acid Side Chain: A Mechanistic Study Using the Ab Initio Nanoreactor

Many disease-causing viruses target sialic acids on the surface of host cells. Some viruses bind preferentially to sialic acids with O-acetyl modification at the hydroxyl group of C7, C8, or C9 on the glycerol-like side chain. Studies of proteins binding to sialosides containing O-acetylated sialic acids are crucial in understanding the related diseases but experimentally difficult due to the lability of the ester group. We recently showed that O-acetyl migration among hydroxyl groups of C7, C8, and C9 in sialic acids occurs in all directions in a pH-dependent manner. In the current study, we elucidate a full mechanistic pathway for the migration of O-acetyl among C7, C8, and C9. We used an ab initio nanoreactor to explore potential reaction pathways and density functional theory, pK_a calculations, and umbrella sampling to investigate elementary steps of interest. We found that when a base is present, migration is easy in any direction and involves three key steps: deprotonation of the hydroxyl group, cyclization between the two carbons, and the migration of the O-acetyl group. This dynamic equilibrium may play a defensive role against pathogens that evolve to gain entry to the cell by binding selectively to one acetylation state.

Obtaining QM/MM binding free energies in the SAMPL8 drugs of abuse challenge: indirect approaches

Accurately predicting free energy differences is essential in realizing the full potential of rational drug design. Unfortunately, high levels of accuracy often require computationally expensive QM/MM Hamiltonians. Fortuitously, the cost of employing QM/MM approaches in rigorous free energy simulation can be reduced through the use of the so-called “indirect” approach to QM/MM free energies, in which the need for QM/MM simulations is avoided via a QM/MM “correction” at the classical endpoints of interest. Herein, we focus on the computation of QM/MM binding free energies in the context of the SAMPL8 Drugs of Abuse host–guest challenge. Of the 5 QM/MM correction coupled with force-matching submissions, PM6-D3H4/MM ranked submission proved the best overall QM/MM entry, with an RMSE from experimental results of 2.43 kcal/mol (best in ranked submissions), a Pearson’s correlation of 0.78 (second-best in ranked submissions), and a Kendall \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document}τ correlation of 0.52 (best in ranked submissions).

Approach to Study pH-Dependent Protein Association Using Constant-pH Molecular Dynamics: Application to the Dimerization of β-Lactoglobulin

Protein-protein association is often mediated by electrostatic interactions and modulated by pH. However, experimental and computational studies have often overlooked the effect of association on the protonation state of the protein. In this work, we present a methodological approach based on constant-pH molecular dynamics (MD), which aims to provide a detailed description of a pH-dependent protein-protein association, and apply it to the dimerization of β-lactoglobulin (BLG). A selection of analyses is performed using the data generated by constant-pH MD simulations of monomeric and dimeric forms of bovine BLG, in the pH range 3-8. First, we estimate free energies of dimerization using a computationally inexpensive approach based on the Wyman-Tanford linkage theory, calculated in a new way through the use of thermodynamically based splines. The individual free energy contribution of each titratable site is also calculated, allowing for identification of relevant residues. Second, the correlations between the proton occupancies of pairs of sites are calculated (using the Pearson coefficient), and extensive networks of correlated sites are observed at acidic pH values, sometimes involving distant pairs. In general, strongly correlated sites are also slow proton exchangers and contribute significantly to the pH-dependency of the dimerization free energy. Third, we use ionic density as a fingerprint of protein charge distribution and observe electrostatic complementarity between the monomer faces that form the dimer interface, more markedly at the isoionic point (where maximum dimerization occurs) than at other pH values, which might contribute to guide the association. Finally, the pH-dependent dimerization modes are inspected using PCA, among other analyses, and two states are identified: a relaxed state at pH 4-8 (with the typical alignment of the crystallographic structure) and a compact state at pH 3-4 (with a tighter association and rotated alignment). This work shows that an approach based on constant-pH MD simulations can produce rich detailed pictures of pH-dependent protein associations, as illustrated for BLG dimerization.

Efficient Alchemical Intermediate States in Free Energy Calculations Using λ-Enveloping Distribution Sampling

Alchemical free energy calculations generally require intermediate states along a coupling parameter λ to establish sufficient phase space overlap for obtaining converged results. Such intermediate states can also be engineered to lower the energy barriers and, consequently, reduce the required sampling time. The recently introduced λ-enveloping distribution sampling (λ-EDS) scheme combines the properties of the minimum variance pathway and the EDS methods to improve sampling and allow for larger steps along the alchemical pathway compared to conventional approaches. This scheme also eliminates the need for soft-core potentials and retains the behavior of conventional λ-intermediate states as a limiting case. In this study, an automated procedure is developed to select the parameters of λ-EDS for optimal performance. The underlying theory is illustrated based on simulations of simple test systems (bond length changes in harmonic oscillators, mutations of dihedral angles, and charge creation in water), as well as on the calculation of the absolute hydration free energies of 12 small organic molecules.

Computational predictions of metal-macrocycle stability constants require accurate treatments of local solvent and pH effects

Rational design of molecular chelating agents requires a detailed understanding of physicochemical ligand-metal interactions in solvent phase. Computational quantum chemistry methods should be able to provide this, but computational reports have shown poor accuracy when determining absolute binding constants for many chelating molecules. To understand why, we compare and benchmark static- and dynamics-based computational procedures for a range of monovalent and divalent cations binding to a conventional cryptand molecule: 2.2.2-cryptand ([2.2.2]). The benchmarking comparison shows that dynamics simulations using standard OPLS-AA classical potentials can reasonably predict binding constants for monovalent cations, but these procedures fail for divalent cations. We also consider computationally efficient static procedure using Kohn-Sham density functional theory (DFT) and cluster-continuum modeling that accounts for local microsolvation and pH effects. This approach accurately predicts binding energies for monovalent and divalent cations with an average error of 3.2 kcal mol-1 compared to experiment. This static procedure thus should be useful for future molecular screening efforts, and high absolute errors in the literature may be due to inadequate modeling of local solvent and pH effects.

An Alternative to Conventional λ-Intermediate States in Alchemical Free Energy Calculations: λ-Enveloping Distribution Sampling

Alchemical free energy calculations typically rely on intermediate states to bridge between the relevant phase spaces of the two end states. These intermediate states are usually created by mixing the energies or parameters of the end states according to a coupling parameter λ. The choice of the procedure has a strong impact on the efficiency of the calculation, as it affects both the encountered energy barriers and the phase space overlap between the states. The present work builds on the connection between the minimum variance pathway (MVP) and enveloping distribution sampling (EDS). It is shown that both methods can be regarded as special cases of a common scheme referred to as λ-EDS, which can also reproduce the behavior of conventional λ-intermediate states. A particularly attractive feature of λ-EDS is its ability to emulate the use of soft core potentials (SCP) while avoiding the associated computational overhead when applying efficient free energy estimators such as the multistate Bennett's acceptance ratio (MBAR). The method is illustrated for both relative and absolute free energy calculations considering five benchmark systems. The first two systems (charge inversion and cavity creation in a dipolar solvent) demonstrate the use of λ-EDS as an alternative coupling scheme in the context of thermodynamic integration (TI). The three other systems (change of bond length, change of dihedral angles, and cavity creation in water) investigate the efficiency and optimal choice of parameters in the context of free energy perturbation (FEP) and Bennett's acceptance ratio (BAR). It is shown that λ-EDS allows larger steps along the alchemical pathway than conventional intermediate states.

Exploring pH Dependent Host/Guest Binding Affinities

When the electrostatic environment surrounding binding partners changes between unbound and bound states, the net uptake or release of a proton is possible by either binding partner. This process is pH-dependent in that the free energy required to uptake or release the proton varies with pH. This pH-dependence is typically not considered in conventional free energy methods where the use of fixed protonation states is the norm. In the present paper, we apply a simple two-step approach to calculate the pH-dependent binding free energy of a model cucubit[7]uril host/guest system. By use of λ-dynamics with an enhanced sampling protocol, adaptive landscape flattening, pK_a shifts and reference binding free energies upon complexation were determined. This information enables the construction of pH-dependent binding profiles that accurately capture the pK_a shifts and reproduce binding free energies at the different pH conditions that were observed experimentally. Our calculations illustrate a general framework for computing pH-dependent binding free energies but also point to some issues in modeling the molecular charge distributions within this series of molecules with CGenFF. However, by introducing some minor charge modifications to the CGenFF force field, we saw significant improvement in accuracy of the calculated pK_a shifts.

Kinetics of pH-dependent interactions between PD-1 and PD-L1 immune checkpoint proteins from molecular dynamics

Immune checkpoint blockade of signaling pathways such as PD-1/PD-L1 has recently opened up a new avenue for highly efficient immunotherapeutic strategies to treat cancer. Since tumor microenvironments are characterized by lower pH (5.5-7.0), pH-dependent protein-ligand interactions can be exploited as efficient means to regulate drug affinity and specificity for a variety of malignancies. In this article, we investigate the mechanism and kinetics of pH-dependent binding and unbinding processes for the PD-1/PD-L1 checkpoint pair employing classical molecular dynamics simulations. Two representative pH levels corresponding to circumneutral physiological conditions of blood (pH 7.4) and acidic tumor microenvironment (pH 5.5) are considered. Our calculations demonstrate that pH plays a key role in protein-ligand interactions with small pH changes leading to several orders of magnitude increase in binding affinity. By identifying the binding pocket in the PD-1/PD-L1 complex, we show a pivotal role of the His68 protonation state of PD-1in the complex stabilization at low pH. The results on the reaction rate constants are in qualitative agreement with available experimental data. The obtained molecular details are important for further engineering of binding/unbinding kinetics to formulate more efficient immune checkpoint blockade strategies.

Molecular Dynamics Simulations of Membrane Permeability

This Review illustrates the evaluation of permeability of lipid membranes from molecular dynamics (MD) simulation primarily using water and oxygen as examples. Membrane entrance, translocation, and exit of these simple permeants (one hydrophilic and one hydrophobic) can be simulated by conventional MD, and permeabilities can be evaluated directly by Fick's First Law, transition rates, and a global Bayesian analysis of the inhomogeneous solubility-diffusion model. The assorted results, many of which are applicable to simulations of nonbiological membranes, highlight the limitations of the homogeneous solubility diffusion model; support the utility of inhomogeneous solubility diffusion and compartmental models; underscore the need for comparison with experiment for both simple solvent systems (such as water/hexadecane) and well-characterized membranes; and demonstrate the need for microsecond simulations for even simple permeants like water and oxygen. Undulations, subdiffusion, fractional viscosity dependence, periodic boundary conditions, and recent developments in the field are also discussed. Last, while enhanced sampling methods and increasingly sophisticated treatments of diffusion add substantially to the repertoire of simulation-based approaches, they do not address directly the critical need for force fields with polarizability and multipoles, and constant pH methods.

Tracing the Pathways of Waters and Protons in Photosystem II and Cytochrome c Oxidase

Photosystem II (PSII) uses water as the terminal electron donor, producing oxygen in the Mn4CaO5 oxygen evolving complex (OEC), while cytochrome c oxidase (CcO) reduces O2 to water in its heme–Cu binuclear center (BNC). Each protein is oriented in the membrane to add to the proton gradient. The OEC, which releases protons, is located near the P-side (positive, at low-pH) of the membrane. In contrast, the BNC is in the middle of CcO, so the protons needed for O2 reduction must be transferred from the N-side (negative, at high pH). In addition, CcO pumps protons from N- to P-side, coupled to the O2 reduction chemistry, to store additional energy. Thus, proton transfers are directly coupled to the OEC and BNC redox chemistry, as well as needed for CcO proton pumping. The simulations that study the changes in proton affinity of the redox active sites and the surrounding protein at different states of the reaction cycle, as well as the changes in hydration that modulate proton transfer paths, are described.

Overview of the SAMPL6 host-guest binding affinity prediction challenge

Accurately predicting the binding affinities of small organic molecules to biological macromolecules can greatly accelerate drug discovery by reducing the number of compounds that must be synthesized to realize desired potency and selectivity goals. Unfortunately, the process of assessing the accuracy of current computational approaches to affinity prediction against binding data to biological macromolecules is frustrated by several challenges, such as slow conformational dynamics, multiple titratable groups, and the lack of high-quality blinded datasets. Over the last several SAMPL blind challenge exercises, host-guest systems have emerged as a practical and effective way to circumvent these challenges in assessing the predictive performance of current-generation quantitative modeling tools, while still providing systems capable of possessing tight binding affinities. Here, we present an overview of the SAMPL6 host-guest binding affinity prediction challenge, which featured three supramolecular hosts: octa-acid (OA), the closely related tetra-endo-methyl-octa-acid (TEMOA), and cucurbit[8]uril (CB8), along with 21 small organic guest molecules. A total of 119 entries were received from ten participating groups employing a variety of methods that spanned from electronic structure and movable type calculations in implicit solvent to alchemical and potential of mean force strategies using empirical force fields with explicit solvent models. While empirical models tended to obtain better performance than first-principle methods, it was not possible to identify a single approach that consistently provided superior results across all host-guest systems and statistical metrics. Moreover, the accuracy of the methodologies generally displayed a substantial dependence on the system considered, emphasizing the need for host diversity in blind evaluations. Several entries exploited previous experimental measurements of similar host-guest systems in an effort to improve their physical-based predictions via some manner of rudimentary machine learning; while this strategy succeeded in reducing systematic errors, it did not correspond to an improvement in statistical correlation. Comparison to previous rounds of the host-guest binding free energy challenge highlights an overall improvement in the correlation obtained by the affinity predictions for OA and TEMOA systems, but a surprising lack of improvement regarding root mean square error over the past several challenge rounds. The data suggests that further refinement of force field parameters, as well as improved treatment of chemical effects (e.g., buffer salt conditions, protonation states), may be required to further enhance predictive accuracy.

An explicit-solvent hybrid QM and MM approach for predicting pKa of small molecules in SAMPL6 challenge

In this work we have developed a hybrid QM and MM approach to predict pKa of small drug-like molecules in explicit solvent. The gas phase free energy of deprotonation is calculated using the M06-2X density functional theory level with Pople basis sets. The solvation free energy difference of the acid and its conjugate base is calculated at MD level using thermodynamic integration. We applied this method to the 24 drug-like molecules in the SAMPL6 blind pKa prediction challenge. We achieved an overall RMSE of 2.4 pKa units in our prediction. Our results show that further optimization of the protocol needs to be done before this method can be used as an alternative approach to the well established approaches of a full quantum level or empirical pKa prediction methods.

A fast and accurate computational method for the linear-combination-based isotropic periodic sum

An isotropic periodic sum (IPS) is a powerful technique to reasonably calculate intermolecular interactions for wide range of molecular systems under periodic boundary conditions. A linear-combination-based IPS (LIPS) has been developed to attain computational accuracy close to an exact lattice sum, such as the Ewald sum. The algorithm of the original LIPS method has a high computational cost because it needs long-range interaction calculations in real space. This becomes a performance bottleneck for long-time molecular simulations. In this work, the combination of an LIPS and fast Fourier transform (FFT) was developed, and evaluated on homogeneous and heterogeneous molecular systems. This combinational approach of LIPS/FFT attained computational efficiency close to that of a smooth particle mesh Ewald while maintaining the same high accuracy as the original LIPS. We concluded that LIPS/FFT has great potential to extend the capability of IPS techniques for the fast and accurate computation of many types of molecular systems.

Critical test of isotropic periodic sum techniques with group-based cut-off schemes

Truncation is still chosen for many long-range intermolecular interaction calculations to efficiently compute free-boundary systems, macromolecular systems and net-charge molecular systems, for example. Advanced truncation methods have been developed for long-range intermolecular interactions. Every truncation method can be implemented as one of two basic cut-off schemes, namely either an atom-based or a group-based cut-off scheme. The former computes interactions of “atoms” inside the cut-off radius, whereas the latter computes interactions of “molecules” inside the cut-off radius. In this work, the effect of group-based cut-off is investigated for isotropic periodic sum (IPS) techniques, which are promising cut-off treatments to attain advanced accuracy for many types of molecular system. The effect of group-based cut-off is clearly different from that of atom-based cut-off, and severe artefacts are observed in some cases. However, no severe discrepancy from the Ewald sum is observed with the extended IPS techniques.

Development of constant-pH simulation methods in implicit solvent and applications in biomolecular systems

pH is a critical parameter for biological and technological systems directly related with electrical charges. It can give rise to peculiar electrostatic phenomena, which also makes them more challenging. Due to the quantum nature of the process, involving the forming and breaking of chemical bonds, quantum methods should ideally by employed. Nevertheless, due to the very large number of ionizable sites, different macromolecular conformations, salt conditions, and all other charged species, the CPU time cost simply becomes prohibitive for computer simulations, making this a quite complex problem. Simplified methods based on Monte Carlo sampling have been devised and will be reviewed here, highlighting the updated state-of-the-art of this field, advantages, and limitations of different theoretical protocols for biomolecular systems (proteins and nucleic acids). Following a historical perspective, the discussion will be associated with the applications to protein interactions with other proteins, polyelectrolytes, and nanoparticles.

Isotropic periodic sum for multipole interactions and a vector relation for calculation of the Cartesian multipole tensor

Isotropic periodic sum (IPS) is a method to calculate long-range interactions based on the homogeneity of simulation systems. By using the isotropic periodic images of a local region to represent remote structures, long-range interactions become a function of the local conformation. This function is called the IPS potential; it folds long-ranged interactions into a short-ranged potential and can be calculated as efficiently as a cutoff method. It has been demonstrated that the IPS method produces consistent simulation results, including free energies, as the particle mesh Ewald (PME) method. By introducing the multipole homogeneous background approximation, this work derives multipole IPS potentials, abbreviated as IPSMm, with m being the maximum order of multipole interactions. To efficiently calculate the multipole interactions in Cartesian space, we propose a vector relation that calculates a multipole tensor as a dot product of a radial potential vector and a directional vector. Using model systems with charges, dipoles, and/or quadrupoles, with and without polarizability, we demonstrate that multipole interactions of order m can be described accurately with the multipole IPS potential of order 2 or m - 1, whichever is higher. Through simulations with the multipole IPS potentials, we examined energetic, structural, and dynamic properties of the model systems and demonstrated that the multipole IPS potentials produce very similar results as PME with a local region radius (cutoff distance) as small as 6 Å.

Origin of pKa Shifts of Internal Lysine Residues in SNase Studied Via Equal-Molar VMMS Simulations in Explicit Water

Protein internal ionizable groups can exhibit large shifts in pK_a values. Although the environment and interaction changes have been extensively studied both experimentally and computationally, direct calculation of pK_a values of these internal ionizable groups in explicit water is challenging due to energy barriers in solvent interaction and in conformational transition. The virtual mixture of multiple states (VMMS) method is a new approach designed to study chemical state equilibrium. This method constructs a virtual mixture of multiple chemical states in order to sample the conformational space of all states simultaneously and to avoid crossing energy barriers related to state transition. By applying VMMS to 25 variants of staphylococcal nuclease with lysine residues at internal positions, we obtained the pK_a values of these lysine residues and investigated the physics underlining the pK_a shifts. Our calculation results agree reasonably well with experimental measurements, validating the VMMS method for pK_a calculation and providing molecular details of the protonation equilibrium for protein internal ionizable groups. Based on our analyses of protein conformation relaxation, lysine side chain flexibility, water penetration, and the microenvironment, we conclude that the hydrophobicity of the microenvironment around the lysine side chain (which affects water penetration differently for different protonation states) plays an important role in the pK_a shifts.

Absolute binding free energy calculations of CBClip host-guest systems in the SAMPL5 blind challenge

Herein, we report the absolute binding free energy calculations of CBClip complexes in the SAMPL5 blind challenge. Initial conformations of CBClip complexes were obtained using docking and molecular dynamics simulations. Free energy calculations were performed using thermodynamic integration (TI) with soft-core potentials and Bennett's acceptance ratio (BAR) method based on a serial insertion scheme. We compared the results obtained with TI simulations with soft-core potentials and Hamiltonian replica exchange simulations with the serial insertion method combined with the BAR method. The results show that the difference between the two methods can be mainly attributed to the van der Waals free energies, suggesting that either the simulations used for TI or the simulations used for BAR, or both are not fully converged and the two sets of simulations may have sampled difference phase space regions. The penalty scores of force field parameters of the 10 guest molecules provided by CHARMM Generalized Force Field can be an indicator of the accuracy of binding free energy calculations. Among our submissions, the combination of docking and TI performed best, which yielded the root mean square deviation of 2.94 kcal/mol and an average unsigned error of 3.41 kcal/mol for the ten guest molecules. These values were best overall among all participants. However, our submissions had little correlation with experiments.

Calculating distribution coefficients based on multi-scale free energy simulations: an evaluation of MM and QM/MM explicit solvent simulations of water-cyclohexane transfer in the SAMPL5 challenge

One of the central aspects of biomolecular recognition is the hydrophobic effect, which is experimentally evaluated by measuring the distribution coefficients of compounds between polar and apolar phases. We use our predictions of the distribution coefficients between water and cyclohexane from the SAMPL5 challenge to estimate the hydrophobicity of different explicit solvent simulation techniques. Based on molecular dynamics trajectories with the CHARMM General Force Field, we compare pure molecular mechanics (MM) with quantum-mechanical (QM) calculations based on QM/MM schemes that treat the solvent at the MM level. We perform QM/MM with both density functional theory (BLYP) and semi-empirical methods (OM1, OM2, OM3, PM3). The calculations also serve to test the sensitivity of partition coefficients to solute polarizability as well as the interplay of the quantum-mechanical region with the fixed-charge molecular mechanics environment. Our results indicate that QM/MM with both BLYP and OM2 outperforms pure MM. However, this observation is limited to a subset of cases where convergence of the free energy can be achieved.

Continuum Electrostatics Approaches to Calculating pKas and Ems in Proteins

Proteins change their charge state through protonation and redox reactions as well as through binding charged ligands. The free energy of these reactions is dominated by solvation and electrostatic energies and modulated by protein conformational relaxation in response to the ionization state changes. Although computational methods for calculating these interactions can provide very powerful tools for predicting protein charge states, they include several critical approximations of which users should be aware. This chapter discusses the strengths, weaknesses, and approximations of popular computational methods for predicting charge states and understanding the underlying electrostatic interactions. The goal of this chapter is to inform users about applications and potential caveats of these methods as well as outline directions for future theoretical and computational research.