Benchmarking pK(a) prediction

Role of Protonation States in the Stability of Molecular Dynamics Simulations of High-Resolution Membrane Protein Structures

Classical molecular dynamics (MD) simulations provide unmatched spatial and time resolution of protein structure and function. However, the accuracy of MD simulations often depends on the quality of force field parameters and the time scale of sampling. Another limitation of conventional MD simulations is that the protonation states of titratable amino acid residues remain fixed during simulations, even though protonation state changes coupled to conformational dynamics are central to protein function. Due to the uncertainty in selecting protonation states, classical MD simulations are sometimes performed with all amino acids modeled in their standard charged states at pH 7. Here, we performed and analyzed classical MD simulations on high-resolution cryo-EM structures of two large membrane proteins that transfer protons by catalyzing protonation/deprotonation reactions. In simulations performed with titratable amino acids modeled in their standard protonation (charged) states, the structure diverges far from its starting conformation. In comparison, MD simulations performed with predetermined protonation states of amino acid residues reproduce the structural conformation, protein hydration, and protein-water and protein-protein interactions of the structure much better. The results support the notion that it is crucial to perform basic protonation state calculations, especially on structures where protonation changes play an important functional role, prior to the launch of any conventional MD simulations. Furthermore, the combined approach of fast protonation state prediction and MD simulations can provide valuable information about the charge states of amino acids in the cryo-EM sample. Even though accurate prediction of protonation states in proteinaceous environments currently remains a challenge, we introduce an approach of combining pKa prediction with cryo-EM density map analysis that helps in improving not only the protonation state predictions but also the atomic modeling of density data.

Computational insights on the hydride and proton transfer mechanisms of L-proline dehydrogenase

L-Proline dehydrogenase (ProDH) is a flavin-dependent oxidoreductase, which catalyzes the oxidation of L-proline to (S)-1-pyrroline-5-carboxylate. Based on the experimental studies, a stepwise proton and hydride transfer mechanism is supported. According to this mechanism, the amino group of L-proline is deprotonated by a nearby Lys residue, which is followed by the hydride transfer process from C5 position of L-proline to N5 position of isoalloxazine ring of FAD. It was concluded that the hydride transfer step is rate limiting in the reductive half-reaction, however, in the overall reaction, the oxidation of FAD is the rate limiting step. In this study, we performed a computational mechanistic investigation based on ONIOM method to elucidate the mechanism of the reductive half-reaction corresponding to the oxidation of L-proline into iminoproline. Our calculations support the stepwise mechanism in which the deprotonation occurs initially as a fast step as result of a proton transfer from L-proline to the Lys residue. Subsequently, a hydride ion transfers from L-proline to FAD with a higher activation barrier. The enzyme-product complex showed a strong interaction between reduced FAD and iminoproline, which might help to explain why a step in the oxidative half-reaction is rate-limiting.

Computational Insights on the Hydride and Proton Transfer Mechanisms of D-Arginine Dehydrogenase

D-Arginine dehydrogenase from Pseudomonas aeruginosa (PaDADH) is an amine oxidase which catalyzes the conversion of D-arginine into iminoarginine. It contains a non-covalent FAD cofactor that is involved in the oxidation mechanism. Based on substrate, solvent, and multiple kinetic isotope effects studies, a stepwise hydride transfer mechanism is proposed. It was shown that D-arginine binds to the active site of enzyme as α-amino group protonated, and it is deprotonated before a hydride ion is transferred from its α-C to FAD. Based on a mutagenesis study, it was concluded that a water molecule is the most likely catalytic base responsible from the deprotonation of α-amino group. In this study, we formulated computational models based on ONIOM method to elucidate the oxidation mechanism of D-arginine into iminoarginine using the crystal structure of enzyme complexed with iminoarginine. The calculations showed that Arg222, Arg305, Tyr249, Glu87, His 48, and two active site water molecules play key roles in binding and catalysis. Model systems showed that the deprotonation step occurs prior to hydride transfer step, and active site water molecule(s) may have participated in the deprotonation process.

Computational Mechanistic Study of l-Aspartate Oxidase by ONIOM Method

l-Aspartate oxidase (Laspo) is responsible for the oxidation of l-aspartate into iminoaspartate using flavin as a cofactor. During this process flavin is reduced, and it can be reoxidized by either molecular oxygen or fumarate. The overall fold and the catalytic residues of Laspo are similar to succinate dehydrogenase and fumarate reductase. On the basis of deuterium kinetic isotope effects as well as other kinetic and structural data, it is proposed that the enzyme can catalyze the oxidation of l-aspartate through a mechanism similar to amino acid oxidases. It is suggested that a proton is removed from the α-amino group, while a hydride is transferred from C2 to flavin. It is also suggested that the hydride transfer is a rate-limiting step. However, there is still an ambiguity about the stepwise or concerted mechanism of hydride- and proton-transfer steps. In this study, we formulated some computational models to study the hydride-transfer mechanism using the crystal structure of Escherichia colil-aspartate oxidase in complexes with succinate. The calculations involved our own N-layered integrated molecular orbital and molecular mechanics method, and we evaluated the geometry and energetics of the hydride/proton-transfer processes while probing the roles of active site residues. Based on the calculations, it is concluded that proton- and hydride-transfer steps are decoupled, and a stepwise mechanism might be operative as opposed to the concerted one.

Fuzzy Drug Targets: Disordered Proteins in the Drug-Discovery Realm

Intrinsically disordered proteins (IDPs) and regions (IDRs) form a large part of the eukaryotic proteome. Contrary to the structure-function paradigm, the disordered proteins perform a myriad of functions in vivo. Consequently, they are involved in various disease pathways and are plausible drug targets. Unlike folded proteins, that have a defined structure and well carved out drug-binding pockets that can guide lead molecule selection, the disordered proteins require alternative drug-development methodologies that are based on an acceptable picture of their conformational ensemble. In this review, we discuss various experimental and computational techniques that contribute toward understanding IDP "structure" and describe representative pursuances toward IDP-targeting drug development. We also discuss ideas on developing rational drug design protocols targeting IDPs.

Sequence-based prediction of pH-dependent protein solubility using CamSol

Solubility is a property of central importance for the use of proteins in research in molecular and cell biology and in applications in biotechnology and medicine. Since experimental methods for measuring protein solubility are material intensive and time consuming, computational methods have recently emerged to enable the rapid and inexpensive screening of solubility for large libraries of proteins, as it is routinely required in development pipelines. Here, we describe the development of one such method to include in the predictions the effect of the pH on solubility. We illustrate the resulting pH-dependent predictions on a variety of antibodies and other proteins to demonstrate that these predictions achieve an accuracy comparable with that of experimental methods. We make this method publicly available at https://www-cohsoftware.ch.cam.ac.uk/index.php/camsolph, as the version 3.0 of CamSol.

Bioinspired Polyethylene Glycol Coatings for Reduced Nanoparticle-Protein Interactions

Nanoparticles (NPs) and other engineered nanomaterials have great potential as nanodrugs or nanomedical devices for biomedical applications. However, the adsorption of proteins in blood circulation or similar physiological fluids can significantly alter the surface properties and therapeutic response induced by most nanomaterials. For example, interaction with proteins can change the bloodstream circulation time and availability of therapeutic NPs or hinder the accumulation in their desired target organs. Proteins can also trigger or prevent agglomeration. By combining experimental and computational approaches, we have developed NPs carrying polyethylene glycol (PEG) polymeric coatings that mimic the surface charge distribution of proteins typically found in blood, which are known to show low aggregation under normal blood conditions. Here, we show that NPs with coatings based on apoferritin or human serum albumin display better antifouling properties and weaker protein interaction compared to similar NPs carrying conventional PEG polymeric coatings.

Combined MD and QM/MM Calculations Reveal Allostery-Driven Promiscuity in Dipeptide Epimerases of Enolase Family

The adaptability of the active site to amplify the secondary function is supposed to be the fundamental cause of the promiscuity and the evolution of new functions in enzymes. In most cases, mutations occur close to the active site and/or in the catalytic site to change the active site plasticity to accommodate the non-native substrate. In the present study, using MD simulations and hybrid QM/MM calculations, we have shown a way to enhance the promiscuity, i. e., the allostery-driven promiscuity. Using a case study of the AEE enzyme where the capping loop recognizes the substrate, herein, we show that a single site mutation (D321G) far from the capping loop can induce a large conformational change in the capping loop to recognize different substrates for different functions. The QM/MM calculations for the WT and mutated enzyme provide a first validation of the mechanism of 1,1-proton transfer and dehydration by the AEE enzyme. Since AEE epimerase possesses a highly conserved TIM-barrel fold, we believe that our study provides a crucial lead to understanding the mechanism of emergence of secondary function which can be useful to repurpose ancient enzymes for modern usage.

Computational Analysis of the Inhibition Mechanism of NOTUM by the ONIOM Method

Notum is a member of serine hydrolyses that cleaves the palmitoleate moiety from Wingless-related integration site (Wnt) ligands. This enzyme plays crucial functions through modulating the Wnt signaling pathway. Inhibition of Notum carries therapeutic effects against a number of maladies including osteoporosis, cancer, and Alzheimer's disease. Recently, a class of irreversible inhibitors based on esters of 4-(indolin-1-yl)-4-oxobutanoic acid have been reported. Using the crystal structures of enzyme-4-(indolin-1-yl)-4-oxobutanoate adduct and 4-(indolin-1-yl)-4-oxobutanoic acid-enzyme complex, we studied computationally the proposed inhibition mechanism using model systems based on the own n-layered integrated molecular orbital and molecular mechanics (ONIOM) method. In the first place, model systems were formulated to investigate the transesterification between the catalytic serine residue, Ser-232, and the methyl ester of 4-(indolin-1-yl)-4-oxobutanoate. In the second place, the hydrolysis mechanism of the resultant enzyme-inhibitor adduct was studied. The energetics of these steps were analyzed using a density functional theory functional in the ONIOM method. In addition, the roles of active-site residues during these steps were highlighted. It was found that the hydrolysis of the covalent adduct is highly endergonic corroborating the irreversible inhibition mechanism. These results will shed light not only on the inhibition mechanism but also on the catalytic mechanism.

Modeling a unit cell: crystallographic refinement procedure using the biomolecular MD simulation platform Amber

A procedure has been developed for the refinement of crystallographic protein structures based on the biomolecular simulation program Amber. The procedure constructs a model representing a crystal unit cell, which generally contains multiple protein molecules and is fully hydrated with TIP3P water. Periodic boundary conditions are applied to the cell in order to emulate the crystal lattice. The refinement is conducted in the form of a specially designed short molecular-dynamics run controlled by the Amber ff14SB force field and the maximum-likelihood potential that encodes the structure-factor-based restraints. The new Amber-based refinement procedure has been tested on a set of 84 protein structures. In most cases, the new procedure led to appreciably lower R free values compared with those reported in the original PDB depositions or obtained by means of the industry-standard phenix.refine program. In particular, the new method has the edge in refining low-accuracy scrambled models. It has also been successful in refining a number of molecular-replacement models, including one with an r.m.s.d. of 2.15 Å. In addition, Amber-refined structures consistently show superior MolProbity scores. The new approach offers a highly realistic representation of protein-protein interactions in the crystal, as well as of protein-water interactions. It also offers a realistic representation of protein crystal dynamics (akin to ensemble-refinement schemes). Importantly, the method fully utilizes the information from the available diffraction data, while relying on state-of-the-art molecular-dynamics modeling to assist with those elements of the structure that do not diffract well (for example mobile loops or side chains). Finally, it should be noted that the protocol employs no tunable parameters, and the calculations can be conducted in a matter of several hours on desktop computers equipped with graphical processing units or using a designated web service.

How does α1Histidine102 affect the binding of modulators to α1β2γ2 GABAA receptors? molecular insights from in silico experiments

The activation of GABAA receptors by the neurotransmitter gamma-aminobutyric acid mediates the rapid inhibition response in the central nervous system of mammals. Many neurological and mental health disorders arise from alterations in the structure or function of these pentameric ion channels. GABAA receptors are targets for numerous drugs, including benzodiazepines, which bind to α1β2γ2 GABAA receptors with high affinity to a site in the extracellular domain, between subunits α1 and γ2. It has been established experimentally that the binding of these drugs depends on the presence of one particular amino acid in the α1 subunit: histidine 102. However, the specific role it plays in the intermolecular interaction has not been elucidated. In this work, we applied in silico methods to understand whether certain protonation and rotamer states of α1His102 facilitate the binding of modulators. We analysed diazepam binding, a benzodiazepine, and the antagonist flumazenil to the GABAA receptor using molecular dynamics simulations and adaptive biasing force simulations. The binding free energy follows changes in the protonation state for both ligands, and rotameric states of α1His102 were specific for the different compounds, suggesting distinct preferences for positive allosteric modulators and antagonists. Moreover, in the presence of diazepam and favoured by a neutral tautomer, we identified a water molecule that links loops A, B, and C and may be relevant to the modulation mechanism.

SAMPL7 TrimerTrip host-guest binding poses and binding affinities from spherical-coordinates-biased simulations

Host-guest binding remains a major challenge in modern computational modelling. The newest 7th statistical assessment of the modeling of proteins and ligands (SAMPL) challenge contains a new series of host-guest systems. The TrimerTrip host binds to 16 structurally diverse guests. Previously, we have successfully employed the spherical coordinates as the collective variables coupled with the enhanced sampling technique metadynamics to enhance the sampling of the binding/unbinding event, search for possible binding poses and calculate the binding affinities in all three host-guest binding cases of the 6th SAMPL challenge. In this work, we report a retrospective study on the TrimerTrip host-guest systems by employing the same protocol to investigate the TrimerTrip host in the SAMPL7 challenge. As no binding pose is provided by the SAMPL7 host, our simulations initiate from randomly selected configurations and are proceeded long enough to obtain converged free energy estimates and search for possible binding poses. The calculated binding affinities are in good agreement with the experimental reference, and the obtained binding poses serve as a nice starting point for end-point or alchemical free energy calculations. Note that as the work is performed after the close of the SAMPL7 challenge, we do not participate in the challenge and the results are not formally submitted to the SAMPL7 challenge.

Investigation of bound and unbound phosphoserine phosphatase conformations through elastic network models and molecular dynamics simulations

The human phosphoserine phosphatase (hPSP) catalyses the last step in the biosynthesis of L-serine. It involves conformational changes of the enzyme lid once the substrate, phosphoserine (PSer), is bound in the active site. Here, Elastic Network Model (ENM) is applied to the crystal structure of hPSP to probe the transition between open and closed conformations of hPSP. Molecular Dynamics (MD) simulations are carried out on several PSer-hPSP systems to characterise the intermolecular interactions and their effect on the dynamics of the enzyme lid. Systems involving either Ca⁺⁺ or Mg⁺⁺ are considered. The first ENM normal mode shows that an open-closed transition can be explained from a simple description of the enzyme in terms of harmonic potentials. Principal Component Analyses applied to the MD trajectories also highlight a trend for a closing/opening motion. Different PSer orientations inside the enzyme cavity are identified, i.e. either the carboxylate, the phosphate group of PSer, or both, are oriented towards the cation. The interaction patterns are analysed in terms of hydrogen bonds, electrostatics, and bond critical points of the electron density distributions. The latter approach yields a global description of the bonding intermolecular interactions. The PSer orientation determines the content of the cation coordination shell and the mobility of the substrate, while Lys158 and Thr182, involved in the reaction mechanism, are always in interaction with the substrate. Closed enzyme conformations involve Met52-Gln204, Arg49-Glu29, and Arg50-Glu29 interactions. Met52, as well as Arg49 and Arg50, also stabilize PSer inside the cavity. Communicated by Ramaswamy H. Sarma.

Theoretical understanding of the thermodynamics and interactions in transcriptional regulator TtgR-ligand binding

The transcriptional regulator TtgR belongs to the TetR family of transcriptional repressors. It depresses the transcription of the TtgABC operon and itself and thus regulates the extrusion of noxious chemicals with efflux pumps in bacterial cells. As the ligand-binding domain of TtgR is rather flexible, it can bind with a number of structurally diverse ligands, such as antibiotics, flavonoids and aromatic solvents. In the current work, we perform equilibrium and nonequilibrium alchemical free energy simulation to predict the binding affinities of a series of ligands targeting the TtgR protein and an agreement between the theoretical prediction and the experimental result is observed. End-point methods MM/PBSA and MM/GBSA are also employed for comparison. We further study the interaction maps and contacts between the protein and the ligand and identify important interactions in the protein-ligand binding cases. The dynamics fluctuation and secondary structures are also investigated. The current work sheds light on atomic and thermodynamic understanding of the TtgR-ligand interactions.

a-ARM: Automatic Rhodopsin Modeling with Chromophore Cavity Generation, Ionization State Selection, and External Counterion Placement

The Automatic Rhodopsin Modeling (ARM) protocol has recently been proposed as a tool for the fast and parallel generation of basic hybrid quantum mechanics/molecular mechanics (QM/MM) models of wild type and mutant rhodopsins. However, in its present version, input preparation requires a few hours long user's manipulation of the template protein structure, which also impairs the reproducibility of the generated models. This limitation, which makes model building semiautomatic rather than fully automatic, comprises four tasks: definition of the retinal chromophore cavity, assignment of protonation states of the ionizable residues, neutralization of the protein with external counterions, and finally congruous generation of single or multiple mutations. In this work, we show that the automation of the original ARM protocol can be extended to a level suitable for performing the above tasks without user's manipulation and with an input preparation time of minutes. The new protocol, called a-ARM, delivers fully reproducible (i.e., user independent) rhodopsin QM/MM models as well as an improved model quality. More specifically, we show that the trend in vertical excitation energies observed for a set of 25 wild type and 14 mutant rhodopsins is predicted by the new protocol better than when using the original. Such an agreement is reflected by an estimated (relative to the probed set) trend deviation of 0.7 ± 0.5 kcal mol-1 (0.03 ± 0.02 eV) and mean absolute error of 1.0 kcal mol-1 (0.04 eV).

Catalytic mechanism of type C sialidase from Streptococcus pneumoniae: from covalent intermediate to final product

Streptococcus pneumoniae is a Gram-positive human pathogenic bacterium, which is the main cause of pneumonia and meningitis in children and the elderly. Three sialidases (or neuraminidases) encoded from Streptococcus pneumoniae could catalyze the cleavage of sialic acid linkages. This mechanism is directly connected with infection, apoptosis, and signaling, and usually considered to be one of the critical virulence factors. Type C neuraminidase (NanC) is unique because its primary product of Neu5Ac2en is considered to be an inhibitor to the other two sialidases. Experimentally, there are two different pathways for the formation mechanism of Neu5Ac2en catalyzed by NanC. In this work, a combined quantum mechanical and molecular mechanical approach was employed in all calculations. Starting from the covalent sialylated intermediate, we first examined the reaction to Neu5Ac2en and found the reaction prefers a direct proton abstraction mechanism rather than the water mediated proton abstraction mechanism. Free energy profiles can confirm that Neu5Ac2en is the major product of NanC. Functional roles of some important residues were also investigated, e.g., D315 acts as the proton acceptor during the formation of Neu5Ac2en, while the general base for the hydrolytic reaction to Neu5Ac. This study can facilitate the understanding of the catalytic mechanism of NanC and has the potential to aid in future inhibitor design studies.

Determination of the protonation preferences of bilin pigments in cryptophyte antenna complexes

The light-harvesting mechanisms of cryptophyte antenna complexes have attracted considerable attention due to their ability to exhibit maximal photosynthetic activity under very low-light conditions and to display several colors, as well as the observation of vibronic coherent features in their two-dimensional electronic spectra. However, detailed investigations on the interplay between the protein environment and their light-harvesting properties are hampered by the uncertainty related to the protonation state of the underlying bilin pigments. Here we study the protonation preferences of four types of bilin pigments including 15,16-dihydrobiliverdin (DBV), phycoerythrobilin (PEB), phycocyanobilin (PCB) and mesobiliverdin (MBV), which are found in phycoerythrin PE545 and phycocyanin PC577, PC612, PC630 and PC645 complexes. We apply quantum chemical calculations coupled to continuum solvation calculations to predict the intrinsic acidity of bilins in aqueous solution, and then combine molecular dynamics simulations with empirical pKa estimates to investigate the impact of the local protein environment on the acidity of the pigments. We also report measurements of the absorption spectra of the five complexes in a wide range of pH in order to validate our simulations and investigate possible changes in the light harvesting properties of the complexes in the range of physiological pH found in the lumen (pH ∼ 5-7). The results suggest a pKa > 7 for DBV and MBV pigments in the α polypeptide chains of PE545 and PC630/PC645 complexes, which are not coordinated to a negatively charged amino acid. For the other PEB, DBV and PCB pigments, which interact with a Glu or Asp side chain, higher pKa values (pKa > 8) are estimated. Overall, the results support a preferential population of the fully protonated state for bilins in cryptophyte complexes under physiological conditions regardless of the specific type of pigment and local protein environment.

Coupled molecular dynamics and continuum electrostatic method to compute the ionization pKa's of proteins as a function of pH. Test on a large set of proteins

A computational method, to predict the pKa values of the ionizable residues Asp, Glu, His, Tyr, and Lys of proteins, is presented here. Calculation of the electrostatic free-energy of the proteins is based on an efficient version of a continuum dielectric electrostatic model. The conformational flexibility of the protein is taken into account by carrying out molecular dynamics simulations of 10 ns in implicit water. The accuracy of the proposed method of calculation of pKa values is estimated from a test set of experimental pKa data for 297 ionizable residues from 34 proteins. The pKa-prediction test shows that, on average, 57, 86, and 95% of all predictions have an error lower than 0.5, 1.0, and 1.5 pKa units, respectively. This work contributes to our general understanding of the importance of protein flexibility for an accurate computation of pKa, providing critical insight about the significance of the multiple neutral states of acid and histidine residues for pKa-prediction, and may spur significant progress in our effort to develop a fast and accurate electrostatic-based method for pKa-predictions of proteins as a function of pH.

Membrane Position Dependency of the pKa and Conductivity of the Protein Ion Channel

The dependency of current-voltage characteristics of the α-hemolysin channel on the channel position within the membrane was studied using Poisson-Nernst-Planck theory of ion conductivity with soft repulsion between mobile ions and protein atoms (SP-PNP). The presence of the membrane environment also influences the protonation state of the residues at the boundary of the water-lipid interface. In this work, we predict that Asp and Lys residues at the protein rim change their protonation state upon penetration to the lipid environment. Free energies of protein insertion in the membrane for different penetration depths were estimated using the Poisson-Boltzmann/solvent-accessible surface area (PB/SASA) model. The results show that rectification and reversal potentials are very sensitive to the relative position of channel in the membrane, which in turn contributes to alternative protonation states of lipid-penetrating ionizable groups. The prediction of channel position based on the matching of calculated rectification with experimentally determined rectification is in good agreement with recent neutron reflection experiments. Based on the results, we conclude that α-hemolysin membrane position is determined by a combination of factors and not only by the pattern of the surface hydrophobicity as is typically assumed.

In Silico Evaluation of Variable pH on the Binding of Epidermal Growth Factor Receptor Ectodomain to its Ligand Through Molecular Dynamics Simulation in Tumors

Many aggressive and metastatic cancer cell types show Warburg Effect; therefore, it is a possible adaptation helping cancer cells to rapidly divide and utilize the glycolytic intermediates for biosynthesis of ribose sugars (for nucleotide biosynthesis), fatty acid synthesis (lipids for membrane synthesis), NADPH (cellular currency for reductive biosynthesis) and lactate. This in due course results in decrease of extracellular pH, leading to acidic tumor micro-environment. EGFR is a crucial cell surface signaling receptor implicated in cancer cell survival and progression. This warrants studying the effect of the acidic micro-environmental conditions on the binding of the EGFR cell surface receptor to one of its natural extracellular ligand EGF. We exploited in silico approaches: molecular dynamics simulation at variable pH and MM-GBSA free energy of binding calculation method to evaluate the effect of this change in microenvironmental pH. Through the present study it is reported that at pH 6.6 the EGFR binds to EGF with decreased free energy of binding as compared to pH 7.2.

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.

Extensive Assessment of Various Computational Methods for Aspartate's pKa Shift

A series of computational methods for pK_a shift prediction are extensively tested on a set of benchmark protein systems, aiming at identifying pitfalls and evaluating their performance on high variants. Including 19 ASP residues in 10 protein systems, the benchmark set consists of both residues with highly shifted pK_a values as well as those varying little from the reference value, with an experimental RMS free energy differences of 2.49 kcal/mol with respect to blocked amino acid, namely the RMS pK_a shift being 1.82 pK_a units. The constant pH molecular dynamics (MD), alchemical methods, PROPKA3.1, and multiconformation continuum electrostatics give RMSDs of 1.52, 2.58, 1.37, and 3.52 pK_a units, respectively, on the benchmark set. The empirical scoring method is the most accurate one with extremely low computational cost, and the pH-dependent model is also able to provide accurate results, while the accuracy of MD sampling incorporating alchemical free energy simulation is prohibited by convergence achievement and the performance of conformational search incorporating multiconformation continuum electrostatics is bad. Former research works did not define statistical uncertainty with care and yielded the questionable conclusion that alchemical methods perform well in most benchmarks. In this work the traditional alchemical methods are thoroughly tested for high variants. We also performed the first application of nonequilibrium alchemical methods to the pK_a cases.

PlayMolecule ProteinPrepare: A Web Application for Protein Preparation for Molecular Dynamics Simulations

Protein preparation is a critical step in molecular simulations that consists of refining a Protein Data Bank (PDB) structure by assigning titration states and optimizing the hydrogen-bonding network. In this application note, we describe ProteinPrepare, a web application designed to interactively support the preparation of protein structures. Users can upload a PDB file, choose the solvent pH value, and inspect the resulting protonated residues and hydrogen-bonding network within a 3D web interface. Protonation states are suggested automatically but can be manually changed using the visual aid of the hydrogen-bonding network. Tables and diagrams provide estimated pK_a values and charge states, with visual indication for cases where review is required. We expect the graphical interface to be a useful instrument to assess the validity of the preparation, but nevertheless, a script to execute the preparation offline with the High-Throughput Molecular Dynamics (HTMD) environment is also provided for noninteractive operations.

Effect of a pH Gradient on the Protonation States of Cytochrome c Oxidase: A Continuum Electrostatics Study

Cytochrome c oxidase (CcO) couples the reduction of dioxygen to water with transmembrane proton pumping, which leads to the generation of an electrochemical gradient. In this study we analyze how one of the components of the electrochemical gradient, the difference in pH across the membrane, or ΔpH, influences the protonation states of residues in CcO. We modified our continuum electrostatics/Monte Carlo (CE/MC) method in order to include the ΔpH and applied it to the study of CcO, in what is, to our best knowledge, the first CE/MC study of CcO in the presence of a pH gradient. The inclusion of a transmembrane pH gradient allows for the identification of residues whose titration behavior depends on the pH on both sides of the membrane. Among the several residues with unusual titration profiles, three are well-known key residues in the proton transfer process of CcO: E286_I, Y288_I, and K362_I. All three residues have been previously identified as being critical for the catalytic or proton pumping functions of CcO. Our results suggest that when the pH gradient increases, these residues may be part of a regulatory mechanism to stem the proton flow.

Predicting stabilizing mutations in proteins using Poisson-Boltzmann based models: study of unfolded state ensemble models and development of a successful binary classifier based on residue interaction energies

In many cases the stability of a protein has to be increased to permit its biotechnological use. Rational methods of protein stabilization based on optimizing electrostatic interactions have provided some fine successful predictions. However, the precise calculation of stabilization energies remains challenging, one reason being that the electrostatic effects on the unfolded state are often neglected. We have explored here the feasibility of incorporating Poisson-Boltzmann model electrostatic calculations performed on representations of the unfolded state as large ensembles of geometrically optimized conformations calculated using the ProtSA server. Using a data set of 80 electrostatic mutations experimentally tested in two-state proteins, the predictive performance of several such models has been compared to that of a simple one that considers an unfolded structure of non-interacting residues. The unfolded ensemble models, while showing correlation between the predicted stabilization values and the experimental ones, are worse than the simple model, suggesting that the ensembles do not capture well the energetics of the unfolded state. A more attainable goal is classifying potential mutations as either stabilizing or non-stabilizing, rather than accurately calculating their stabilization energies. To implement a fast classification method that can assist in selecting stabilizing mutations, we have used a much simpler electrostatic model based only on the native structure and have determined its precision using different stabilizing energy thresholds. The binary classifier developed finds 7 true stabilizing mutants out of every 10 proposed candidates and can be used as a robust tool to propose stabilizing mutations.

Benchmarking pKa Prediction Methods for Residues in Proteins

Methods for estimation of pKa values of residues in proteins were tested on a set of benchmark proteins with experimentally known pKa values. The benchmark set includes 80 different residues (20 each for Asp, Glu, Lys, and His), half of which consists of significantly variant cases (ΔpKa ≥ 1 pKa unit from the amino acid in solution). The method introduced by Case and co-workers [J. Am. Chem. Soc. 2004, 126, 4167-4180], referred to as the molecular dynamics/generalized-Born/thermodynamic integration (MD/GB/TI) technique, gives a root-mean-square deviation (rmsd) of 1.4 pKa units on the benchmark set. The use of explicit waters in the immediate region surrounding the residue was shown to generally reduce high errors for this method. Longer simulation time was also shown to increase the accuracy of this method. The empirical approach developed by Jensen and co-workers [Proteins 2005, 61, 704-721], PROPKA, also gives an overall rmsd of 1.4 pKa units and is more or less accurate based on residue type-the method does very well for Lys and Glu, but less so for Asp and His. Likewise, the absolute deviation is quite similar for the two methods-5.2 for PROPKA and 5.1 for MD/GB/TI. A comparison of these results with several prediction methods from the literature is presented. The error in pKa prediction is analyzed as a function of variation of the pKa from that in water and the solvent accessible surface area (SASA) of the residue. A case study of the catalytic lysine residue in 2-deoxyribose-5-phosphate aldolase (DERA) is also presented.

Expanding the Catalytic Triad in Epoxide Hydrolases and Related Enzymes

Potato epoxide hydrolase 1 exhibits rich enantio- and regioselectivity in the hydrolysis of a broad range of substrates. The enzyme can be engineered to increase the yield of optically pure products as a result of changes in both enantio- and regioselectivity. It is thus highly attractive in biocatalysis, particularly for the generation of enantiopure fine chemicals and pharmaceuticals. The present work aims to establish the principles underlying the activity and selectivity of the enzyme through a combined computational, structural, and kinetic study using the substrate trans-stilbene oxide as a model system. Extensive empirical valence bond simulations have been performed on the wild-type enzyme together with several experimentally characterized mutants. We are able to computationally reproduce the differences between the activities of different stereoisomers of the substrate and the effects of mutations of several active-site residues. In addition, our results indicate the involvement of a previously neglected residue, H104, which is electrostatically linked to the general base H300. We find that this residue, which is highly conserved in epoxide hydrolases and related hydrolytic enzymes, needs to be in its protonated form in order to provide charge balance in an otherwise negatively charged active site. Our data show that unless the active-site charge balance is correctly treated in simulations, it is not possible to generate a physically meaningful model for the enzyme that can accurately reproduce activity and selectivity trends. We also expand our understanding of other catalytic residues, demonstrating in particular the role of a noncanonical residue, E35, as a "backup base" in the absence of H300. Our results provide a detailed view of the main factors driving catalysis and regioselectivity in this enzyme and identify targets for subsequent enzyme design efforts.

Effect of mutations on the thermostability of Aspergillus aculeatus β-1,4-galactanase

New variants of β-1,4-galactanase from the mesophilic organism Aspergillus aculeatus were designed using the structure of β-1,4-galactanase from the thermophile organism Myceliophthora thermophila as a template. Some of the variants were generated using PROPKA 3.0, a validated pKa prediction tool, to test its usefulness as an enzyme design tool. The PROPKA designed variants were D182N and S185D/Q188T, G104D/A156R. Variants Y295F and G306A were designed by a consensus approach, as a complementary and validated design method. D58N was a stabilizing mutation predicted by both methods. The predictions were experimentally validated by measurements of the melting temperature (Tm ) by differential scanning calorimetry. We found that the Tm is elevated by 1.1 °C for G306A, slightly increased (in the range of 0.34 to 0.65 °C) for D182N, D58N, Y295F and unchanged or decreased for S185D/Q188T and G104D/A156R. The Tm changes were in the range predicted by PROPKA. Given the experimental errors, only the D58N and G306A show significant increase in thermodynamic stability. Given the practical importance of kinetic stability, the kinetics of the irreversible enzyme inactivation process were also investigated for the wild-type and three variants and found to be biphasic. The half-lives of thermal inactivation were approximately doubled in G306A, unchanged for D182N and, disappointingly, a lot lower for D58N. In conclusion, this study tests a new method for estimating Tm changes for mutants, adds to the available data on the effect of substitutions on protein thermostability and identifies an interesting thermostabilizing mutation, which may be beneficial also in other galactanases.

Protein flexibility and cysteine reactivity: influence of mobility on the H-bond network and effects on pKa prediction

Thanks to its chemical plasticity, cysteine (Cys) is a very versatile player in proteins. A major determinant of Cys reactivity is pKa: the ability to predict it is deemed critical in redox bioinformatics. I considered different computational methods for pKa predictions and ultimately applied one (propka, ppka1) to various datasets; for all residues I assessed the effect of (1) hydrogen bonding, electrostatics and solvation on predictions and (2) protein mobility on pKa variability. Particularly for Cys, exposure and H-bond contributions heavily dictated propka predictions. The prominence of H-bond contributions was previously reported: this may explain the effectiveness of ppka1 (with Cys, tested in a benchmark). However ppka1 was also very sensitive to protein mobility; I assessed the effects of mobility on particularly large (compared to previous studies) datasets of structural ensembles; I found that exposed Cys presented the highest pKa variability, ascribable to correspondingly high H-bond fluctuations associated with protein flexibility. The benefit of including protein dynamics in pKa predictions was previously proposed, but empirical methods were never tested in this sense; instead, giving their outstanding speed, they could lend particularly well to this purpose. I devised a strategy combining short range molecular dynamics with ppka1; the protocol aimed to mitigate high ppka1 variability by including a "statistical view" of fast conformational changes. Tested in a benchmark, the strategy lead to improved performances. These results provide new insights on Cys bioinformatics (pKa prediction protocols) and Cys biology (effect of mobility on exposed Cys properties).

The alteration of the C-terminal region of human frataxin distorts its structural dynamics and function

Friedreich's ataxia (FRDA) is linked to a deficiency of frataxin (FXN), a mitochondrial protein involved in iron-sulfur cluster synthesis. FXN is a small protein with an α/β fold followed by the C-terminal region (CTR) with a nonperiodic structure that packs against the protein core. In the present study, we explored the impact of the alteration of the CTR on the stability and dynamics of FXN. We analyzed several pathological and rationally designed CTR mutants using complementary spectroscopic and biophysical approaches. The pathological mutation L198R yields a global destabilization of the structure correlating with a significant and highly localized alteration of dynamics, mainly involving residues that are in contact with L198 in wild-type FXN. Variant FXN 90-195, which is closely related to the FRDA-associated mutant FXN 81-193, conserves a globular shape with a native-like structure. However, the truncation of the CTR results in an extreme alteration of global stability and protein dynamics over a vast range of timescales and encompassing regions far from the CTR, as shown by proton-water exchange rates and (15) N-relaxation measurements. Increased sensitivity to proteolysis, observed in vitro for both mutants, suggests a faster degradation rate in vivo, whereas the enhanced tendency to aggregate exhibited by the truncated variant may account for the loss of functional FXN, with both phenomena providing an explanation as to why the alteration of the CTR causes FRDA. These results contribute to understanding how stability and activity are linked to protein motions and they might be useful for the design of target-specific ligands to control local protein motions for stability enhancement.

Quantum mechanics/molecular mechanics restrained electrostatic potential fitting

We present a quantum mechanics/molecular mechanics (QM/MM) method to evaluate the partial charges of amino acid residues for use in MM potentials based on their protein environment. For each residue of interest, the nearby residues are included in the QM system while the rest of the protein is treated at the MM level of theory. After a short structural optimization, the partial charges of the central residue are fit to the electrostatic potential using the restrained electrostatic potential (RESP) method. The resulting charges and electrostatic potential account for the individual environment of the residue, although they lack the transferable nature of library partial charges. To evaluate the quality of the QM/MM RESP charges, thermodynamic integration is used to measure the pKa shift of the aspartic acid residues in three different proteins, turkey egg lysozyme, beta-cryptogein, and Thioredoxin. Compared to the AMBER ff99SB library values, the QM/MM RESP charges show better agreement between the calculated and experimental pK(a) values for almost all of the residues considered.

Predicting pKa Values in Aqueous Solution for the Guanidine Functional Group from Gas Phase Ab Initio Bond Lengths

Here we applied a novel method1a to predict pKa values of the guanidine functional group, which is a notoriously difficult. This method, which was developed in our lab, uses only one ab initio bond length obtained at a low level of theory. The method is shown to work for drug molecules, delivers prediction errors of less than 0.5 log units, successfully treats tautomerisation in close relation with experiment, and demonstrates strong correlations with only a few data points. The high structural content of the ab initio bond length makes a given data set essentially divide itself into high correlation subsets. One then observes that molecules within a subset possess a common substructure. Each high correlation subset exists in its own region of chemical space. The high correlation subset method is explored with respect to this position in chemical space, in particular tautomerisation. The proposed method is able to distinguish between different tautomeric forms and the preferred tautomeric form emerges naturally, in agreement with experiment.

Benchmark Study for the Cysteine-Histidine Proton Transfer Reaction in a Protein Environment: Gas Phase, COSMO, QM/MM Approaches

Proton transfer reactions are of crucial interest for the investigation of proteins. We have investigated the accuracy of commonly used quantum chemical methods for the description of proton transfer reactions in different environments (gas phase, COSMO, QM/MM) using the proton transfer between the catalytic dyad residues cysteine 145 and histidine 41 of SARS coronavirus main protease as a case study. The test includes thermodynamic, kinetic, and structural properties. The study comprises computationally demanding ab initio approaches (HF, CC2, MP2, SCS-CC2, SCS-MP2, CCSD(T)), popular density functional theories (BLYP, B3LYP, M06-2X), and semiempirical methods (MNDO/d, AM1, RM1, PM3, PM6). The approximated coupled cluster approach LCCSD(T) is taken as a reference method. We find that the robustness of the tested methods with respect to the environment correlates well with the level of theory. As an example HF, CC2, MP2, and their SCS variants show similar errors for gas phase, COSMO, or QM/MM computations. In contrast for semiempirical methods, the errors strongly diversify if one goes from gas phase to COSMO or QM/MM. Particular problems are observed for the recent semiempirical methods PM6 and RM1, which show the best performance for gas phase calculations but possess larger errors in conjunction with COSMO. Finally, a combination of SCS-MP2 and B3LYP or M06-2X allows reliable estimates about remaining errors.

Understanding the pK(a) of redox cysteines: the key role of hydrogen bonding

Many cellular functions involve cysteine chemistry via thiol-disulfide exchange pathways. The nucleophilic cysteines of the enzymes involved are activated as thiolate. A thiolate is much more reactive than a neutral thiol. Therefore, determining and understanding the pK(a)s of functional cysteines are important aspects of biochemistry and molecular biology with direct implications for redox signaling. Here, we describe the experimental and theoretical methods to determine cysteine pK(a) values, and we examine the factors that control these pK(a)s. Drawing largely on experience gained with the thioredoxin superfamily, we examine the roles of solvation, charge-charge, helix macrodipole, and hydrogen bonding interactions as pK(a)-modulating factors. The contributions of these factors in influencing cysteine pK(a)s and the associated chemistry, including the relevance for the reaction kinetics and thermodynamics, are discussed. This analysis highlights the critical role of direct hydrogen bonding to the cysteine sulfur as a key factor modulating the equilibrium between thiol S-H and thiolate S(-). This role is easily understood intuitively and provides a framework for biochemical functional insights.

Rigorous incorporation of tautomers, ionization species, and different binding modes into ligand-based and receptor-based 3D-QSAR methods

Speciation of drug candidates and receptors caused by ionization, tautomerism, and/or covalent hydration complicates ligandand receptor-based predictions of binding affinities by 3-dimensional structure-activity relationships (3D-QSAR). The speciation problem is exacerbated by tendency of tautomers to bind in multiple conformations or orientations (modes) in the same binding site. New forms of the 3D-QSAR correlation equations, capable of capturing this complexity, can be developed using the time hierarchy of all steps that lie behind the monitored biological process - binding, enzyme inhibition or receptor activity. In most cases, reversible interconversions of individual ligand and receptor species can be treated as quickly established equilibria because they are finished in a small fraction of the exposure time that is used to determine biological effects. The speciation equilibria are satisfactorily approximated by invariant fractions of individual ligand and receptor species for buffered experimental or in vivo conditions. For such situations, the observed drug-receptor association constant of a ligand is expressed as the sum of products, for each ligand and receptor species pair, of the association microconstant and the fractions of involved species. For multiple binding modes, each microconstant is expressed as the sum of microconstants of individual modes. This master equation leads to new 3D-QSAR correlation equations integrating the results of all molecular simulations or calculations, which are run for each ligand-receptor species pair separately. The multispecies, multimode 3D-QSAR approach is illustrated by a ligand-based correlation of transthyretin binding of thyroxine analogs and by a receptor-based correlation of inhibition of MK2 by benzothiophenes and pyrrolopyrimidines.

pKa modulation of the acid/base catalyst within GH32 and GH68: a role in substrate/inhibitor specificity?

Glycoside hydrolases of families 32 (GH32) and 68 (GH68) belong to clan GH-J, containing hydrolytic enzymes (sucrose/fructans as donor substrates) and fructosyltransferases (sucrose/fructans as donor and acceptor substrates). In GH32 members, some of the sugar substrates can also function as inhibitors, this regulatory aspect further adding to the complexity in enzyme functionalities within this family. Although 3D structural information becomes increasingly available within this clan and huge progress has been made on structure-function relationships, it is not clear why some sugars bind as inhibitors without being catalyzed. Conserved aspartate and glutamate residues are well known to act as nucleophile and acid/bases within this clan. Based on the available 3D structures of enzymes and enzyme-ligand complexes as well as docking simulations, we calculated the pKa of the acid-base before and after substrate binding. The obtained results strongly suggest that most GH-J members show an acid-base catalyst that is not sufficiently protonated before ligand entrance, while the acid-base can be fully protonated when a substrate, but not an inhibitor, enters the catalytic pocket. This provides a new mechanistic insight aiming at understanding the complex substrate and inhibitor specificities observed within the GH-J clan. Moreover, besides the effect of substrate entrance on its own, we strongly suggest that a highly conserved arginine residue (in the RDP motif) rather than the previously proposed Tyr motif (not conserved) provides the proton to increase the pKa of the acid-base catalyst.

Who cares for the protons?

It is tempting to use standard protonation states for the analysis of protein-ligand interactions. Two different pK(a) calculation methods, PROPKA (protein pK(a)) and MCCE (multi conformation continuum electrostatics), were applied to challenge this convenient behavior. As data basis, we selected five recently approved drugs for which structural information of the protein-drug complex is available. We analyzed the pK(a) calculations in terms of a measure termed BIPS (binary protonation states) recently introduced by us. Both methods agree in detecting the majority of the sites with atypical BIPS values. However, when using only one method, some of the atypcial BIPS value would have been missed. Therefore, we recommend using both methods to set such an interpretation on a solid basis.