Protocols utilizing constant pH molecular dynamics to compute pH-dependent binding free energies

Resolving coupled pH titrations using alchemical free energy calculations

In a protein, nearby titratable sites can be coupled: the (de)protonation of one may affect the other. The degree of this interaction depends on several factors and can influence the measured p K a . Here, we derive a formalism based on double free energy differences ( Δ Δ G ) for quantifying the individual site p K a values of coupled residues. As Δ Δ G values can be obtained by means of alchemical free energy calculations, the presented approach allows for a convenient estimation of coupled residue p K a s in practice. We demonstrate that our approach and a previously proposed microscopic p K a formalism, can be combined with alchemical free energy calculations to resolve pH-dependent protein p K a values. Toy models and both, regular and constant-pH molecular dynamics simulations, alongside experimental data, are used to validate this approach. Our results highlight the insights gleaned when coupling and microstate probabilities are analyzed and suggest extensions to more complex enzymatic contexts. Furthermore, we find that naïvely computed p K a values that ignore coupling, can be significantly improved when coupling is accounted for, in some cases reducing the error by half. In short, alchemical free energy methods can resolve the p K a values of both uncoupled and coupled residues.

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.

TLR7 gain-of-function genetic variation causes human lupus

Although circumstantial evidence supports enhanced Toll-like receptor 7 (TLR7) signalling as a mechanism of human systemic autoimmune disease1-7, evidence of lupus-causing TLR7 gene variants is lacking. Here we describe human systemic lupus erythematosus caused by a TLR7 gain-of-function variant. TLR7 is a sensor of viral RNA8,9 and binds to guanosine10-12. We identified a de novo, previously undescribed missense TLR7Y264H variant in a child with severe lupus and additional variants in other patients with lupus. The TLR7Y264H variant selectively increased sensing of guanosine and 2',3'-cGMP10-12, and was sufficient to cause lupus when introduced into mice. We show that enhanced TLR7 signalling drives aberrant survival of B cell receptor (BCR)-activated B cells, and in a cell-intrinsic manner, accumulation of CD11c+ age-associated B cells and germinal centre B cells. Follicular and extrafollicular helper T cells were also increased but these phenotypes were cell-extrinsic. Deficiency of MyD88 (an adaptor protein downstream of TLR7) rescued autoimmunity, aberrant B cell survival, and all cellular and serological phenotypes. Despite prominent spontaneous germinal-centre formation in Tlr7Y264H mice, autoimmunity was not ameliorated by germinal-centre deficiency, suggesting an extrafollicular origin of pathogenic B cells. We establish the importance of TLR7 and guanosine-containing self-ligands for human lupus pathogenesis, which paves the way for therapeutic TLR7 or MyD88 inhibition.

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.

Alchemical Binding Free Energy Calculations in AMBER20: Advances and Best Practices for Drug Discovery

Predicting protein-ligand binding affinities and the associated thermodynamics of biomolecular recognition is a primary objective of structure-based drug design. Alchemical free energy simulations offer a highly accurate and computationally efficient route to achieving this goal. While the AMBER molecular dynamics package has successfully been used for alchemical free energy simulations in academic research groups for decades, widespread impact in industrial drug discovery settings has been minimal because of the previous limitations within the AMBER alchemical code, coupled with challenges in system setup and postprocessing workflows. Through a close academia-industry collaboration we have addressed many of the previous limitations with an aim to improve accuracy, efficiency, and robustness of alchemical binding free energy simulations in industrial drug discovery applications. Here, we highlight some of the recent advances in AMBER20 with a focus on alchemical binding free energy (BFE) calculations, which are less computationally intensive than alternative binding free energy methods where full binding/unbinding paths are explored. In addition to scientific and technical advances in AMBER20, we also describe the essential practical aspects associated with running relative alchemical BFE calculations, along with recommendations for best practices, highlighting the importance not only of the alchemical simulation code but also the auxiliary functionalities and expertise required to obtain accurate and reliable results. This work is intended to provide a contemporary overview of the scientific, technical, and practical issues associated with running relative BFE simulations in AMBER20, with a focus on real-world drug discovery applications.

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.

Rigorous Free Energy Simulations in Virtual Screening

Virtual high throughput screening (vHTS) in drug discovery is a powerful approach to identify hits: when applied successfully, it can be much faster and cheaper than experimental high-throughput screening approaches. However, mainstream vHTS tools have significant limitations: ligand-based methods depend on knowledge of existing chemical matter, while structure-based tools such as docking involve significant approximations that limit their accuracy. Recent advances in scientific methods coupled with dramatic speedups in computational processing with GPUs make this an opportune time to consider the role of more rigorous methods that could improve the predictive power of vHTS workflows. In this Perspective, we assert that alchemical binding free energy methods using all-atom molecular dynamics simulations have matured to the point where they can be applied in virtual screening campaigns as a final scoring stage to prioritize the top molecules for experimental testing. Specifically, we propose that alchemical absolute binding free energy (ABFE) calculations offer the most direct and computationally efficient approach within a rigorous statistical thermodynamic framework for computing binding energies of diverse molecules, as is required for virtual screening. ABFE calculations are particularly attractive for drug discovery at this point in time, where the confluence of large-scale genomics data and insights from chemical biology have unveiled a large number of promising disease targets for which no small molecule binders are known, precluding ligand-based approaches, and where traditional docking approaches have foundered to find progressible chemical matter.

A dual target of Plasmepsin IX and X: Unveiling the atomistic superiority of a core chemical scaffold in malaria therapy

Plasmepsin IX and X, members of the prominent aspartic family of proteases whose function were hitherto unknown have only recently been established as key mediators of erythrocyte invasion and egress of the virulent malarial parasite. Inhibitor 49c, a potent antimalarial peptidomimetic inhibitor initially developed to target Plasmepsin II has lately been proven to exhibit potent inhibitory activity against Plasmepsin IX and X. However, the molecular and structural dynamics supporting its inhibitory activity remain inconclusive. Hindering the motion of the flap and hinge region of an aspartic protease remains essential for disabling the catalytic activity of the enzyme. Integrating molecular dynamic simulations coupled with other advanced biocomputational tools, we reveal the enhanced structural mechanistic competence of 49c in complex with Plasmepsin IX and X relative to Pepstatin. Pepstatin, a known aspartic protease inhibitor which actively hinders the opening and closing of the flap tip and flexible loop and consequently limits access to the catalytic aspartic residues, however, its administration has been related to elevated levels of toxicity. Thermodynamic calculations reveal a higher relative binding free energy associated with Plasmepsin IX and X in complex with 49c as opposed to Pepstatin. A relatively compact and structurally rigid 49c bound complexes sequel into the restriction of the flap and hinge residues by restraining cohesive movement, consequently hindering their "twisting motion" from transpiring. Findings unveil an atomistic perspective into the structural superiority of 49c in complex with Plasmepsin IX and X.

Correlating interfacial water dynamics with protein-protein interaction in complex of GDF-5 and BMPRI receptors

GDF-5 mediated signal transduction regulating chondrogenesis and skeletogenesis involves three different type-I receptors viz. Act-RI, BMPRIA and BMPRIB. BMPRIA and BMPRIB generally shows temporal and spatial co-expression but some spatially different expression pattern has also been observed. BMPRIA receptor is the key receptor implicated in BMP signalling during osteogenesis and is expressed in osteoblasts during the course of bone formation. However, BMPRIB appears to be primarily expressed in mesenchymal pre-cartilage condensations and also found in differentiated osteoblast and chondrocytes. The extracellular pH affects bone cell function and it is experimentally known that mineralization of bone is affected by shift of pH in cultured osteoblast. Here we report the effect of pH on dynamics of water present at the interface of GDF-5:BMPRIA and GDF-5:BMPRIB and binding interaction energy of these complexes. Water dynamics at different pH was analysed using residence time and hydrogen bond relaxation kinetics. pH influences the interaction energy between GDF-5 and BMPRIA and BMPRIB receptors indicating the electrostatic environment modulating the activity of two receptors. This pH dependence of interaction energy is further supported by similar behaviour of hydrogen bond existence of buried water molecules at the interface. In contrast to this the slow and fast exchanging water molecules do not show similar pH dependence of hydrogen bonding relaxation kinetics. Hence; we conclude that only buried water molecule at the interface influences the protein-protein interaction and the electrostatic environment of the extracellular fluid might decide the specificity of the two receptors.

Host-guest complexes of imazalil with cucurbit[8]uril and β-cyclodextrin and their effect on plant pathogenic fungi

We report the control of imazalil (IMZ) antifungal activity utilizing its non-covalent assembly with β-cyclodextrins (β-CD) and cucurbit[8]uril (CB8) macrocycles, as well as its stimuli-responsive disassembly with cadaverine. The NMR results are consistent with inclusion of a single IMZ molecule inside the cavities of either CB8 from its aromatic site or β-CD from its aliphatic end. Efficient complex formation with both host molecules and controlled released upon the addition of cadaverine is supported by NMR measurements. The stimuli-responsiveness of the same host-guest assemblies with cadaverine was validated against seven economically important plant pathogenic fungi which cause agriculturally important plant diseases across the globe. While loading the drug into macrocycles cavities suppressed its activity, subsequent adding of cadaverine efficiently restored it up. The results in the present paper enable researchers working in the area of mycology and plant pathology to inhibit or reduce the fungal growth on demand in order to control these economically important plant pathogenic fungi.

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.

Predicting accurate absolute binding energies in aqueous solution: thermodynamic considerations for electronic structure methods

Recent predictions of absolute binding free energies of host-guest complexes in aqueous solution using electronic structure theory have been encouraging for some systems, while other systems remain problematic. In this paper I summarize some of the many factors that could easily contribute 1-3 kcal mol(-1) errors at 298 K: three-body dispersion effects, molecular symmetry, anharmonicity, spurious imaginary frequencies, insufficient conformational sampling, wrong or changing ionization states, errors in the solvation free energy of ions, and explicit solvent (and ion) effects that are not well-represented by continuum models. While I focus on binding free energies in aqueous solution the approach also applies (with minor adjustments) to any free energy difference such as conformational or reaction free energy differences or activation free energies in any solvent.

Computation of pH-dependent binding free energies

Protein-ligand binding accompanies changes in the surrounding electrostatic environments of the two binding partners and may lead to changes in protonation upon binding. In cases where the complex formation results in a net transfer of protons, the binding process is pH-dependent. However, conventional free energy computations or molecular docking protocols typically employ fixed protonation states for the titratable groups in both binding partners set a priori, which are identical for the free and bound states. In this review, we draw attention to these important yet largely ignored binding-induced protonation changes in protein-ligand association by outlining physical origins and prevalence of the protonation changes upon binding. Following a summary of various theoretical methods for pKa prediction, we discuss the theoretical framework to examine the pH dependence of protein-ligand binding processes.

Conformational Dynamics and Binding Free Energies of Inhibitors of BACE-1: From the Perspective of Protonation Equilibria

BACE-1 is the β-secretase responsible for the initial amyloidogenesis in Alzheimer’s disease, catalyzing hydrolytic cleavage of substrate in a pH-sensitive manner. The catalytic mechanism of BACE-1 requires water-mediated proton transfer from aspartyl dyad to the substrate, as well as structural flexibility in the flap region. Thus, the coupling of protonation and conformational equilibria is essential to a full in silico characterization of BACE-1. In this work, we perform constant pH replica exchange molecular dynamics simulations on both apo BACE-1 and five BACE-1-inhibitor complexes to examine the effect of pH on dynamics and inhibitor binding properties of BACE-1. In our simulations, we find that solution pH controls the conformational flexibility of apo BACE-1, whereas bound inhibitors largely limit the motions of the holo enzyme at all levels of pH. The microscopic pK_a values of titratable residues in BACE-1 including its aspartyl dyad are computed and compared between apo and inhibitor-bound states. Changes in protonation between the apo and holo forms suggest a thermodynamic linkage between binding of inhibitors and protons localized at the dyad. Utilizing our recently developed computational protocol applying the binding polynomial formalism to the constant pH molecular dynamics (CpHMD) framework, we are able to obtain the pH-dependent binding free energy profiles for various BACE-1-inhibitor complexes. Our results highlight the importance of correctly addressing the binding-induced protonation changes in protein-ligand systems where binding accompanies a net proton transfer. This work comprises the first application of our CpHMD-based free energy computational method to protein-ligand complexes and illustrates the value of CpHMD as an all-purpose tool for obtaining pH-dependent dynamics and binding free energies of biological systems.

Formation of insoluble amyloid plaques in the vascular and hippocampal areas of the brain characterizes Alzheimer’s disease, a devastating neurodegenerative disorder causing dementia. Site-specific hydrolytic catalysis of β-secretase, or BACE-1, is responsible for production of oligomerative amyloid β-peptide. As the catalytic activity of BACE-1 is pH-dependent and its structural dynamics are intrinsic to the catalysis, we examine the dependence of dynamics of BACE-1 on solution pH and its implications on the catalytic mechanism of BACE-1. Also, we highlight the importance of accurate description of protonation states of the titratable groups in computer-aided drug discovery targeting BACE-1. We hope the understanding of pH dependence of the dynamics and inhibitor binding properties of BACE-1 will aid the structure-based inhibitor design efforts against Alzheimer’s disease.

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

We present a computational scheme to compute the pH‐dependence of binding free energy with explicit solvent. Despite the importance of pH, the effect of pH has been generally neglected in binding free energy calculations because of a lack of accurate methods to model it. To address this limitation, we use a constant‐pH methodology to obtain a true ensemble of multiple protonation states of a titratable system at a given pH and analyze the ensemble using the Bennett acceptance ratio (BAR) method. The constant pH method is based on the combination of enveloping distribution sampling (EDS) with the Hamiltonian replica exchange method (HREM), which yields an accurate semi‐grand canonical ensemble of a titratable system. By considering the free energy change of constraining multiple protonation states to a single state or releasing a single protonation state to multiple states, the pH dependent binding free energy profile can be obtained. We perform benchmark simulations of a host‐guest system: cucurbit[7]uril (CB[7]) and benzimidazole (BZ). BZ experiences a large pK_a shift upon complex formation. The pH‐dependent binding free energy profiles of the benchmark system are obtained with three different long‐range interaction calculation schemes: a cutoff, the particle mesh Ewald (PME), and the isotropic periodic sum (IPS) method. Our scheme captures the pH‐dependent behavior of binding free energy successfully. Absolute binding free energy values obtained with the PME and IPS methods are consistent, while cutoff method results are off by 2 kcal mol⁻¹. We also discuss the characteristics of three long‐range interaction calculation methods for constant‐pH simulations.

Interpretation of pH-activity profiles for acid-base catalysis from molecular simulations

The measurement of reaction rate as a function of pH provides essential information about mechanism. These rates are sensitive to the pK(a) values of amino acids directly involved in catalysis that are often shifted by the enzyme active site environment. Experimentally observed pH-rate profiles are usually interpreted using simple kinetic models that allow estimation of "apparent pK(a)" values of presumed general acid and base catalysts. One of the underlying assumptions in these models is that the protonation states are uncorrelated. In this work, we introduce the use of constant pH molecular dynamics simulations in explicit solvent (CpHMD) with replica exchange in the pH-dimension (pH-REMD) as a tool to aid in the interpretation of pH-activity data of enzymes and to test the validity of different kinetic models. We apply the methods to RNase A, a prototype acid-base catalyst, to predict the macroscopic and microscopic pK(a) values, as well as the shape of the pH-rate profile. Results for apo and cCMP-bound RNase A agree well with available experimental data and suggest that deprotonation of the general acid and protonation of the general base are not strongly coupled in transphosphorylation and hydrolysis steps. Stronger coupling, however, is predicted for the Lys41 and His119 protonation states in apo RNase A, leading to the requirement for a microscopic kinetic model. This type of analysis may be important for other catalytic systems where the active forms of the implicated general acid and base are oppositely charged and more highly correlated. These results suggest a new way for CpHMD/pH-REMD simulations to bridge the gap with experiments to provide a molecular-level interpretation of pH-activity data in studies of enzyme mechanisms.