pH dependence of binding reactions from free energy simulations and macroscopic continuum electrostatic calculations: application to 2'GMP/3'GMP binding to ribonuclease T1 and implications for catalysis

Estimating the Roles of Protonation and Electronic Polarization in Absolute Binding Affinity Simulations

Accurate prediction of binding free energies is critical to streamlining the drug development and protein design process. With the advent of GPU acceleration, absolute alchemical methods, which simulate the removal of ligand electrostatics and van der Waals interactions with the protein, have become routinely accessible and provide a physically rigorous approach that enables full consideration of flexibility and solvent interaction. However, standard explicit solvent simulations are unable to model protonation or electronic polarization changes upon ligand transfer from water to the protein interior, leading to inaccurate prediction of binding affinities for charged molecules. Here, we perform extensive simulation totaling ∼540 μs to benchmark the impact of modeling conditions on predictive accuracy for absolute alchemical simulations. Binding to urokinase plasminogen activator (UPA), a protein frequently overexpressed in metastatic tumors, is evaluated for a set of 10 inhibitors with extended flexibility, highly charged character, and titratable properties. We demonstrate that the alchemical simulations can be adapted to utilize the MBAR/PBSA method to improve the accuracy upon incorporating electronic polarization, highlighting the importance of polarization in alchemical simulations of binding affinities. Comparison of binding energy prediction at various protonation states indicates that proper electrostatic setup is also crucial in binding affinity prediction of charged systems, prompting us to propose an alternative binding mode with protonated ligand phenol and Hid-46 at the binding site, a testable hypothesis for future experimental validation.

pKa Calculations with the Polarizable Drude Force Field and Poisson-Boltzmann Solvation Model

Electronic polarization effects have been suggested to play an important role in proton binding to titratable residues in proteins. In this work, we describe a new computational method for pK_a calculations, using Monte Carlo (MC) simulations to sample protein protonation states with the Drude polarizable force field and Poisson-Boltzmann (PB) continuum electrostatic solvent model. While the most populated protonation states at the selected pH, corresponding to residues that are half-protonated at that pH, are sampled using the exact relative free energies computed with Drude particles optimized in the field of the PB implicit solvation model, we introduce an approximation for the protein polarization of low-populated protonation states to reduce the computational cost. The highly populated protonation states used to compute the polarization and pK_a's are then iteratively improved until convergence. It is shown that for lysozyme, when considering 9 of the 18 titratable residues, the new method converged within two iterations with computed pK_a's differing only by 0.02 pH units from pK_a's estimated with the exact approach. Application of the method to predict pK_a's of 94 titratable side chains in 8 proteins shows the Drude-PB model to produce physically more correct results as compared to the additive CHARMM36 (C36) force field (FF). With a dielectric constant of two assigned to the protein interior the Root Mean Square (RMS) deviation between computed and experimental pK_a's is 2.07 and 3.19 pH units with the Drude and C36 models, respectively, and the RMS deviation using the Drude-PB model is relatively insensitive to the choice of the internal dielectric constant in contrast to the additive C36 model. At the higher internal dielectric constant of 20, pK_a's computed with the additive C36 model converge to the results obtained with the Drude polarizable force field, indicating the need to artificially overestimate electrostatic screening in a nonphysical way with the additive FF. In addition, inclusion of both syn and anti orientations of the proton in the neutral state of acidic groups is shown to yield improved agreement with experiment. The present work, which is the first example of the use of a polarizable model for the prediction of pK_a's in proteins, shows that the use of a polarizable model represents a more physically correct model for the treatment of electrostatic contributions to pK_a shifts in proteins.

A Role for the 2' OH of peptidyl-tRNA substrate in peptide release on the ribosome revealed through RF-mediated rescue

The 2' OH of the peptidyl-tRNA substrate is thought to be important for catalysis of both peptide bond formation and peptide release in the ribosomal active site. The release reaction also specifically depends on a release factor protein (RF) to hydrolyze the ester linkage of the peptidyl-tRNA upon recognition of stop codons in the A site. Here, we demonstrate that certain amino acid substitutions (in particular those containing hydroxyl or thiol groups) in the conserved GGQ glutamine of release factor RF1 can rescue defects in the release reaction associated with peptidyl-tRNA substrates lacking a 2' OH. We explored this rescue effect through biochemical and computational approaches that support a model where the 2' OH of the P-site substrate is critical for orienting the nucleophile in a hydrogen-bonding network productive for catalysis.

In silico modeling of pH-optimum of protein-protein binding

Protein-protein association is a pH-dependent process and thus the binding affinity depends on the local pH. In vivo the association occurs in a particular cellular compartment, where the individual monomers are supposed to meet and form a complex. Since the monomers and the complex exist in the same micro environment, it is plausible that they coevolved toward its properties, in particular, toward the characteristic subcellular pH. Here we show that the pH at which the monomers are most stable (pH-optimum) or the pH at which stability is almost pH-independent (pH-flat) of monomers are correlated with the pH-optimum of maximal affinity (pH-optimum of binding) or pH interval at which affinity is almost pH-independent (pH-flat of binding) of the complexes made of the corresponding monomers. The analysis of interfacial properties of protein complexes demonstrates that pH-dependent properties can be roughly estimated using the interface charge alone. In addition, we introduce a parameter beta, proportional to the square root of the absolute product of the net charges of monomers, and show that protein complexes characterized with small or very large beta tend to have neutral pH-optimum. Further more, protein complexes made of monomers carrying the same polarity net charge at neutral pH have either very low or very high pH-optimum of binding. These findings are used to propose empirical rule for predicting pH-optimum of binding provided that the amino acid compositions of the corresponding monomers are available.

Phosphorylation facilitates the integrin binding of filamin under force

Filamins are actin binding proteins that contribute to cytoskeletal integrity and biochemical scaffolds during mechanochemical signal transductions. Structurally, human filamins are dimers composed of an actin-binding domain with 24 immunoglobulin (Ig)-like repeats. In this study, we focus on the recently solved high-resolution crystal structure of Ig-like repeats 19-21 of filamin-A (IgFLNa-R19-R21). IgFLNa-R19-21 is of marked importance because it contains the binding site for integrins and facilitates the dynamic ability of filamin-A to communicate with the extracellular environment. However, the structure of filamin-A shows an interesting domain arrangement where the integrin binding site on IgFLNa-R21 is hindered sterically by IgFLNa-R20. Thus, a number of hypotheses on the regulation of filamin-A exist. Using molecular dynamics simulations we evaluated the effects of two primary regulators of filamin-A, force and phosphorylation. We find that a tensile force of 40 pN is sufficient to initiate the partial removal of the autoinhibition on the integrin binding site of IgFLNa-R21. Force coupled to phosphorylation at Ser(2152), however, affords complete dissociation of autoinhibition with a decreased force requirement. Phosphorylation seems to decrease the threshold for removing the IgFLNa-R20 beta-strand inhibitor within 300 ps with 40 pN tensile force. Furthermore, the molecular dynamic trajectories illustrate phosphorylation of Ser(2152) without force is insufficient to remove autoinhibition. We believe the results of this study implicate filamin-A as a tunable mechanosensor, where its sensitivity can be modulated by the degree of phosphorylation.

Molecular mechanics of filamin's rod domain

Rearrangement of the actin cytoskeleton is integral to cell shape and function. Actin-binding proteins, e.g., filamin, can naturally contribute to the mechanics and function of the actin cytoskeleton. The molecular mechanical bases for filamin's function in actin cytoskeletal reorganization are examined here using molecular dynamics simulations. Simulations are performed by applying forces ranging from 25 pN to 125 pN for 2.5 ns to the rod domain of filamin. Applying small loads ( approximately 25 pN) to filamin's rod domain supplies sufficient energy to alter the conformation of the N-terminal regions of the rod. These forces break local hydrogen bond coordination often enough to allow side chains to find new coordination partners, in turn leading to drastic changes in the conformation of filamin, for example, increasing the hydrophobic character of the N-terminal rod region and, alternatively, activating the C-terminal region to become increasingly stiff. These changes in conformation can lead to changes in the affinity of filamin for its binding partners. Therefore, filamin can function to transduce mechanical signals as well as preserve topology of the actin cytoskeleton throughout the rod domain.

Proton binding to proteins: a free-energy component analysis using a dielectric continuum model

Proton binding plays a critical role in protein structure and function. We report pK(a) calculations for three aspartates in two proteins, using a linear response approach, as well as a "standard" Poisson-Boltzmann approach. Averaging over conformations from the two endpoints of the proton-binding reaction, the protein's atomic degrees of freedom are explicitly modeled. Treating macroscopically the protein's electronic polarizability and the solvent, a meaningful model is obtained, without adjustable parameters. It reproduces qualitatively the electrostatic potentials, proton-binding free energies, Marcus reorganization free energies, and pK(a) shifts from explicit solvent molecular dynamics simulations, and the pK(a) shifts from experiment. For thioredoxin Asp-26, which has a large pK(a) upshift, we correctly capture the balance between unfavorable carboxylate desolvation and favorable interactions with a nearby lysine; similarly for RNase A Asp-14, which has a large pK(a) downshift. For the unshifted thioredoxin Asp-20, desolvation by the protein cavity is overestimated by 2.9 pK(a) units; several effects could explain this. "Standard" Poisson-Boltzmann methods sidestep this problem by using a large, ad hoc protein dielectric; but protein charge-charge interactions are then incorrectly downscaled, giving an unbalanced description of the reaction and a large error for the shifted pK(a) values of Asp-26 and Asp-14.