Proton Binding to Proteins: pKa Calculations with Explicit and Implicit Solvent Models

CHARMM at 45: Enhancements in Accessibility, Functionality, and Speed

Since its inception nearly a half century ago, CHARMM has been playing a central role in computational biochemistry and biophysics. Commensurate with the developments in experimental research and advances in computer hardware, the range of methods and applicability of CHARMM have also grown. This review summarizes major developments that occurred after 2009 when the last review of CHARMM was published. They include the following: new faster simulation engines, accessible user interfaces for convenient workflows, and a vast array of simulation and analysis methods that encompass quantum mechanical, atomistic, and coarse-grained levels, as well as extensive coverage of force fields. In addition to providing the current snapshot of the CHARMM development, this review may serve as a starting point for exploring relevant theories and computational methods for tackling contemporary and emerging problems in biomolecular systems. CHARMM is freely available for academic and nonprofit research at https://academiccharmm.org/program.

A Comprehensive Study of Factors Affecting the Prediction of the pKa Shift of Asp26 in Thioredoxin Protein

The stable protonation state of ionizable amino acids in a protein can be predicted by computing the pKa shift of that residue within the protein environment. Thermodynamic Integration (TI) is an ideal molecular dynamics-based approach for predicting the pKa shift of ionizable protein residues. Here, we probe TI-based simulation protocols for their ability to accurately predict the pKa shift of Asp26 in thioredoxin. While implicit solvent models can predict the pKa shift accurately, explicit solvent models result in substantial errors. To understand the underlying reason for this surprising discrepancy, we investigate the role of various factors such as solvent models, conformational sampling, background charges, and polarization.

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.

Assessment of Nucleobase Protomeric and Tautomeric States in Nucleic Acid Structures for Interaction Analysis and Structure-Based Ligand Design

With increasing interest in RNA as a therapeutic and a potential target, the role of RNA structures has become more important. Even slight changes in nucleobases, such as modifications or protomeric and tautomeric states, can have a large impact on RNA structure and function, while local environments in turn affect protonation and tautomerization. In this work, the application of empirical tools for pKa and tautomer prediction for RNA modifications was elucidated and compared with ab initio quantum mechanics (QM) methods and expanded toward macromolecular RNA structures, where QM is no longer feasible. In this regard, the Protonate3D functionality within the molecular operating environment (MOE) was expanded for nucleobase protomer and tautomer predictions and applied to reported examples of altered protonation states depending on the local environment. Overall, observations of nonstandard protomers and tautomers were well reproduced, including structural C+G:C(A) and A+GG motifs, several mismatches, and protonation of adenosine or cytidine as the general acid in nucleolytic ribozymes. Special cases, such as cobalt hexamine-soaked complexes or the deprotonation of guanosine as the general base in nucleolytic ribozymes, proved to be challenging. The collected set of examples shall serve as a starting point for the development of further RNA protonation prediction tools, while the presented Protonate3D implementation already delivers reasonable protonation predictions for RNA and DNA macromolecules. For cases where higher accuracy is needed, like following catalytic pathways of ribozymes, incorporation of QM-based methods can build upon the Protonate3D-generated starting structures. Likewise, this protonation prediction can be used for structure-based RNA-ligand design approaches.

Understanding Cysteine Reactivity in Protein Environments with Electric Fields

The role cysteine residues play in proteins is mediated by their protonation state, whereby the thiolate form of the side chain is highly reactive while the thiol form is more inert. However, the pKa of cysteine residues is hard to predict as it can differ widely from its reference value in solution, an effect that is accentuated by local effects in the heterogeneous protein environment. Here, we present a new approach to the prediction of cysteine reactivity based on electric field calculations at the thiol/thiolate group. We validated our approach by predicting the protonation state of cysteine residues in different protein environments (in the active site, at the protein surface, and buried within the protein interior), including Cys-25 in papaya protease omega, which was proven problematic for the more traditional constant pH molecular dynamics (MD) technique. We predict pKa shifts consistent with experimental observations, and the decomposition of the electric fields into contributions from molecular fragments provides a direct handle to rationalize local pH and pKa effects in proteins without introducing parameters other than those of the force field used for MD simulations.

Accurately Predicting Protein pKa Values Using Nonequilibrium Alchemy

The stability, solubility, and function of a protein depend on both its net charge and the protonation states of its individual residues. pKa is a measure of the tendency for a given residue to (de)protonate at a specific pH. Although pKa values can be resolved experimentally, theory and computation provide a compelling alternative. To this end, we assess the applicability of a nonequilibrium (NEQ) alchemical free energy method to the problem of pKa prediction. On a data set of 144 residues that span 13 proteins, we report an average unsigned error of 0.77 ± 0.09, 0.69 ± 0.09, and 0.52 ± 0.04 pK for aspartate, glutamate, and lysine, respectively. This is comparable to current state-of-the-art predictors and the accuracy recently reached using free energy perturbation methods (e.g., FEP+). Moreover, we demonstrate that our open-source, pmx-based approach can accurately resolve the pKa values of coupled residues and observe a substantial performance disparity associated with the lysine partial charges in Amber14SB/Amber99SB*-ILDN, for which an underused fix already exists.

Alchemical Calculation of Relative Free Energies for Charge-Changing Mutations at Protein-Protein Interfaces Considering Fixed and Variable Protonation States

The calculation of relative free energies (ΔΔG) for charge-changing mutations at protein-protein interfaces through alchemical methods remains challenging due to variations in the system's net charge during charging steps, the possibility of mutated and contacting ionizable residues occurring in various protonation states, and undersampling issues. In this study, we present a set of strategies, collectively termed TIRST/TIRST-H+, to address some of these challenges. Our approaches combine thermodynamic integration (TI) with the prediction of pKa shifts to calculate ΔΔG values. Moreover, special sets of restraints are employed to keep the alchemically transformed molecules separated. The accuracy of the devised approaches was assessed on a large and diverse data set comprising 164 point mutations of charged residues (Asp, Glu, Lys, and Arg) to Ala at the protein-protein interfaces of complexes with known three-dimensional structures. Mean absolute and root-mean-square errors ranging from 1.38 to 1.66 and 1.89 to 2.44 kcal/mol, respectively, and Pearson correlation coefficients of ∼0.6 were obtained when testing the approaches on the selected data set using the GPU-TI module of Amber18 suite and the ff14SB force field. Furthermore, the inclusion of variable protonation states for the mutated acid residues improved the accuracy of the predicted ΔΔG values. Therefore, our results validate the use of TIRST/TIRST-H+ in prospective studies aimed at evaluating the impact of charge-changing mutations to Ala on the stability of protein-protein complexes.

Combining molecular dynamics simulations and x-ray scattering techniques for the accurate treatment of protonation degree and packing of ionizable lipids in monolayers

The pH-dependent change in protonation of ionizable lipids is crucial for the success of lipid-based nanoparticles as mRNA delivery systems. Despite their widespread application in vaccines, the structural changes upon acidification are not well understood. Molecular dynamics simulations support structure prediction but require an a priori knowledge of the lipid packing and protonation degree. The presetting of the protonation degree is a challenging task in the case of ionizable lipids since it depends on pH and on the local lipid environment and often lacks experimental validation. Here, we introduce a methodology of combining all-atom molecular dynamics simulations with experimental total-reflection x-ray fluorescence and scattering measurements for the ionizable lipid Dlin-MC3-DMA (MC3) in POPC monolayers. This joint approach allows us to simultaneously determine the lipid packing and the protonation degree of MC3. The consistent parameterization is expected to be useful for further predictive modeling of the action of MC3-based lipid nanoparticles.

EGCG-Mediated Protection of Transthyretin Amyloidosis by Stabilizing Transthyretin Tetramers and Disrupting Transthyretin Aggregates

Transthyretin amyloidosis (ATTR) is a progressive and systemic disease caused by the misfolding and amyloid aggregation of transthyretin (TTR). Stabilizing the TTR tetramers and disrupting the formed TTR aggregation are treated as a promising strategy for the treatment of ATTR. Previous studies have reported that epigallocatechin gallate (EGCG) can participate in the whole process of TTR aggregation to prevent ATTR. However, the interaction mechanism of EGCG in this process is still obscure. In this work, we performed molecular dynamics simulations to investigate the interactions between EGCG and TTR tetramers, and between EGCG and TTR aggregates formed by the V30M mutation. The obtained results suggest that EGCG at the binding site of the V30M TTR tetramer can form stable hydrogen bonds with residues in the flexible AB-loop and EF-helix-loop, which reduces the structural mobility of these regions significantly. Additionally, the polyaromatic property of EGCG contributes to the increasement of hydrophobicity at the binding site and thus makes the tetramer difficult to be solvated and dissociated. For V30M-TTR-generated aggregates, EGCG can promote the dissociation of boundary β-strands by destroying key residue interactions of TTR aggregates. Moreover, EGCG is capable of inserting into the side-chain of residues of neighboring β-strands and disrupting the highly structured aggregates. Taken together, this study elucidates the role of EGCG in preventing TTR amyloidosis, which can provide important theoretical support for the future of drug design for ATTR.

Rational engineering of an elevator-type metal transporter ZIP8 reveals a conditional selectivity filter critically involved in determining substrate specificity

Engineering of transporters to alter substrate specificity as desired holds great potential for applications, including metabolic engineering. However, the lack of knowledge on molecular mechanisms of substrate specificity hinders designing effective strategies for transporter engineering. Here, we applied an integrated approach to rationally alter the substrate preference of ZIP8, a Zrt-/Irt-like protein (ZIP) metal transporter with multiple natural substrates, and uncovered the determinants of substrate specificity. By systematically replacing the differentially conserved residues with the counterparts in the zinc transporter ZIP4, we created a zinc-preferring quadruple variant (Q180H/E343H/C310A/N357H), which exhibited largely reduced transport activities towards Cd2+, Fe2+, and Mn2+ whereas increased activity toward Zn2+. Combined mutagenesis, modeling, covariance analysis, and computational studies revealed a conditional selectivity filter which functions only when the transporter adopts the outward-facing conformation. The demonstrated approach for transporter engineering and the gained knowledge about substrate specificity will facilitate engineering and mechanistic studies of other transporters.

Nuclear-Electronic Orbital QM/MM Approach: Geometry Optimizations and Molecular Dynamics

Hybrid quantum mechanical/molecular mechanical (QM/MM) methods allow simulations of chemical reactions in atomistic solvent and heterogeneous environments such as proteins. Herein, the nuclear-electronic orbital (NEO) QM/MM approach is introduced to enable the quantization of specified nuclei, typically protons, in the QM region using a method such as NEO-density functional theory (NEO-DFT). This approach includes proton delocalization, polarization, anharmonicity, and zero-point energy in geometry optimizations and dynamics. Expressions for the energies and analytical gradients associated with the NEO-QM/MM method, as well as the previously developed polarizable continuum model (NEO-PCM), are provided. Geometry optimizations of small organic molecules hydrogen bonded to water in either dielectric continuum solvent or explicit atomistic solvent illustrate that aqueous solvation can strengthen hydrogen-bonding interactions for the systems studied, as indicated by shorter intermolecular distances at the hydrogen-bond interface. We then performed a real-time direct dynamics simulation of a phenol molecule in explicit water using the NEO-QM/MM method. These developments and initial examples provide the foundation for future studies of nuclear-electronic quantum dynamics in complex chemical and biological environments.

Best Practices in Constant pH MD Simulations: Accuracy and Sampling

Various approaches have been proposed to include the effect of pH in molecular dynamics (MD) simulations. Among these, the λ-dynamics approach proposed by Brooks and co-workers [Kong, X.; Brooks III, C. L. J. Chem. Phys.1996, 105, 2414−2423] can be performed with little computational overhead and hfor each typeence be used to routinely perform MD simulations at microsecond time scales, as shown in the accompanying paper [Aho, N. et al. J. Chem. Theory Comput.2022, DOI: 10.1021/acs.jctc.2c00516]. At such time scales, however, the accuracy of the molecular mechanics force field and the parametrization becomes critical. Here, we address these issues and provide the community with guidelines on how to set up and perform long time scale constant pH MD simulations. We found that barriers associated with the torsions of side chains in the CHARMM36m force field are too high for reaching convergence in constant pH MD simulations on microsecond time scales. To avoid the high computational cost of extending the sampling, we propose small modifications to the force field to selectively reduce the torsional barriers. We demonstrate that with such modifications we obtain converged distributions of both protonation and torsional degrees of freedom and hence consistent pK_a estimates, while the sampling of the overall configurational space accessible to proteins is unaffected as compared to normal MD simulations. We also show that the results of constant pH MD depend on the accuracy of the correction potentials. While these potentials are typically obtained by fitting a low-order polynomial to calculated free energy profiles, we find that higher order fits are essential to provide accurate and consistent results. By resolving problems in accuracy and sampling, the work described in this and the accompanying paper paves the way to the widespread application of constant pH MD beyond pK_a prediction.

A Computational Model for the PLP-Dependent Enzyme Methionine γ-Lyase

Pyridoxal-5′-phosphate (PLP) is a cofactor in the reactions of over 160 enzymes, several of which are implicated in diseases. Methionine γ-lyase (MGL) is of interest as a therapeutic protein for cancer treatment. It binds PLP covalently through a Schiff base linkage and digests methionine, whose depletion is damaging for cancer cells but not normal cells. To improve MGL activity, it is important to understand and engineer its PLP binding. We develop a simulation model for MGL, starting with force field parameters for PLP in four main states: two phosphate protonation states and two tautomeric states, keto or enol for the Schiff base moiety. We used the force field to simulate MGL complexes with each form, and showed that those with a fully-deprotonated PLP phosphate, especially keto, led to the best agreement with MGL structures in the PDB. We then confirmed this result through alchemical free energy simulations that compared the keto and enol forms, confirming a moderate keto preference, and the fully-deprotonated and singly-protonated phosphate forms. Extensive simulations were needed to adequately sample conformational space, and care was needed to extrapolate the protonation free energy to the thermodynamic limit of a macroscopic, dilute protein solution. The computed phosphate pK_a was 5.7, confirming that the deprotonated, −2 form is predominant. The PLP force field and the simulation methods can be applied to all PLP enzymes and used, as here, to reveal fine details of structure and dynamics in the active site.

Protein pKa Prediction by Tree-Based Machine Learning

Protonation states of ionizable protein residues modulate many essential biological processes. For correct modeling and understanding of these processes, it is crucial to accurately determine their pKa values. Here, we present four tree-based machine learning models for protein pKa prediction. The four models, Random Forest, Extra Trees, eXtreme Gradient Boosting (XGBoost), and Light Gradient Boosting Machine (LightGBM), were trained on three experimental PDB and pKa datasets, two of which included a notable portion of internal residues. We observed similar performance among the four machine learning algorithms. The best model trained on the largest dataset performs 37% better than the widely used empirical pKa prediction tool PROPKA and 15% better than the published result from the pKa prediction method DelPhiPKa. The overall root-mean-square error (RMSE) for this model is 0.69, with surface and buried RMSE values being 0.56 and 0.78, respectively, considering six residue types (Asp, Glu, His, Lys, Cys, and Tyr), and 0.63 when considering Asp, Glu, His, and Lys only. We provide pKa predictions for proteins in human proteome from the AlphaFold Protein Structure Database and observed that 1% of Asp/Glu/Lys residues have highly shifted pKa values close to the physiological pH.

Electronic Polarization Is Essential for the Stabilization and Dynamics of Buried Ion Pairs in Staphylococcal Nuclease Mutants

Buried charged residues play important roles in the modulation of protein stabilities and conformational dynamics and make crucial contributions to protein functions. Considering the generally nonpolar nature of protein interior, a key question concerns the contribution of electronic polarization to the stabilization and properties of buried charges. We answer this question by conducting free energy simulations using the latest polarizable CHARMM force field based on Drude oscillators for a series of Staphylococcal nuclease mutants that involve a buried Glu-Lys pair in different titration states and orientations. While a nonpolarizable model suggests that the ionized form of the buried Glu-Lys pair is more than 40 kcal/mol less stable than the charge-neutral form, the two titration states are comparable in stability when electronic polarization is included explicitly, a result better reconcilable with available experimental data. Analysis of free energy components suggests that additional stabilization of the ionized Glu-Lys pair has contributions from both the enhanced salt-bridge strength and stronger interaction between the ion-pair and surrounding protein residues and penetrated water. Despite the stronger direct interaction between Glu and Lys, the ion-pair exhibits considerably larger and faster structural fluctuations when polarization is included, due to compensation of interactions in the cavity. Collectively, observations from this work provide compelling evidence that electronic polarization is essential to the stability, hydration, dynamics, and therefore function of buried charges in proteins. Therefore, our study advocates for the explicit consideration of electronic polarization for mechanistic and engineering studies that implicate buried charged residues, such as enzymes and ion transporters.

Prediction of protein pK a with representation learning

The behavior of proteins is closely related to the protonation states of the residues. Therefore, prediction and measurement of pK a are essential to understand the basic functions of proteins. In this work, we develop a new empirical scheme for protein pK a prediction that is based on deep representation learning. It combines machine learning with atomic environment vector (AEV) and learned quantum mechanical representation from ANI-2x neural network potential (J. Chem. Theory Comput. 2020, 16, 4192). The scheme requires only the coordinate information of a protein as the input and separately estimates the pK a for all five titratable amino acid types. The accuracy of the approach was analyzed with both cross-validation and an external test set of proteins. Obtained results were compared with the widely used empirical approach PROPKA. The new empirical model provides accuracy with MAEs below 0.5 for all amino acid types. It surpasses the accuracy of PROPKA and performs significantly better than the null model. Our model is also sensitive to the local conformational changes and molecular interactions.

Why Monoamine Oxidase B Preferably Metabolizes N-Methylhistamine over Histamine: Evidence from the Multiscale Simulation of the Rate-Limiting Step

Histamine levels in the human brain are controlled by rather peculiar metabolic pathways. In the first step, histamine is enzymatically methylated at its imidazole N^τ atom, and the produced N-methylhistamine undergoes an oxidative deamination catalyzed by monoamine oxidase B (MAO-B), as is common with other monoaminergic neurotransmitters and neuromodulators of the central nervous system. The fact that histamine requires such a conversion prior to oxidative deamination is intriguing since MAO-B is known to be relatively promiscuous towards monoaminergic substrates; its in-vitro oxidation of N-methylhistamine is about 10 times faster than that for histamine, yet this rather subtle difference appears to be governing the decomposition pathway. This work clarifies the MAO-B selectivity toward histamine and N-methylhistamine by multiscale simulations of the rate-limiting hydride abstraction step for both compounds in the gas phase, in aqueous solution, and in the enzyme, using the established empirical valence bond methodology, assisted by gas-phase density functional theory (DFT) calculations. The computed barriers are in very good agreement with experimental kinetic data, especially for relative trends among systems, thereby reproducing the observed MAO-B selectivity. Simulations clearly demonstrate that solvation effects govern the reactivity, both in aqueous solution as well as in the enzyme although with an opposing effect on the free energy barrier. In the aqueous solution, the transition-state structure involving histamine is better solvated than its methylated analog, leading to a lower barrier for histamine oxidation. In the enzyme, the higher hydrophobicity of N-methylhistamine results in a decreased number of water molecules at the active side, leading to decreased dielectric shielding of the preorganized catalytic electrostatic environment provided by the enzyme. This renders the catalytic environment more efficient for N-methylhistamine, giving rise to a lower barrier relative to histamine. In addition, the transition state involving N-methylhistamine appears to be stabilized by the surrounding nonpolar residues to a larger extent than with unsubstituted histamine, contributing to a lower barrier with the former.

Oxidation Enhances Binding of Extrahelical 5-Methyl-Cytosines by Thymine DNA Glycosylase

The DNA repair protein thymine DNA glycosylase (TDG) removes mispaired or damaged bases, such as oxidized methyl-cytosine, from DNA by cleavage of the glycosidic bond between the sugar and the target base flipped into the enzyme's active site. The enzyme is active against formyl-cytosine and carboxyl-cytosine, whereas the lower oxidized hydroxymethyl-cytosine and methyl-cytosine itself are not processed by the enzyme. Molecular dynamics simulations with thermodynamic integration of TDG complexed to DNA carrying one of four different (oxidized) methyl-cytosine bases in extrahelcial conformation, methyl-cytosine (mC), hydroxymethyl-cytosine (hmC), formyl-cytosine (fC), or carboxyl-cytosine (caC), show a more favorable binding affinity of the higher oxidized forms, fC and caC, than the nonsubstrate bases hmC and mC. Despite rather comparable, reaction-competent conformations of the flipped bases in the active site of the enzyme, more and stronger interactions with active site residues account for the preferred binding of the higher oxidized bases. Binding of the negatively charged caC and the neutral fC are strengthened by interactions with positively charged His151. Our calculated proton affinities find this protonation state of His151 the preferred one in the presence of caC and conceivable in the presence of fC as well as increasing the binding affinity toward the two bases. Discrimination of the substrate bases is further achieved by the backbone of Tyr152 that forms a strong hydrogen bond to the carboxyl and formyl oxygen atoms of caC and fC, respectively, a contact that is completely lacking in mC and much weaker in hmC. Overall, our computational results indicate that the enzyme discriminates the different oxidation forms of methyl-cytosine already at the formation of the extrahelical complexes.

The electronic temperature and the effective chemical potential parameters of an atom in a molecule. A Fermi-Dirac semi-local variational approach

We developed a numerical procedure to compute the electronic temperature and the effective (local) chemical potential undergone by electrons belonging to a particular molecular species. Our strategy relies on consider atomic basins as open quantum (sub)systems within the context of the quantum theory of atoms in molecules. Each basin is represented by the two parameters, the electronic temperature and the effective chemical potential, which are determined by distributing electrons (fermions) imbedded in each atomic region, through a Fermi-Dirac semi-local variational procedure. The results obtained for 40 different chemical species show that the effective chemical potential is a useful tool to reveal the most acidic/basic atoms in a molecule while the electronic temperature is closely related to the concept of chemical hardness at the local level. Our numerical data also indicate that the electronic temperature values undergone by electrons imbedded in atomic basins are way beyond the room temperature condition, allowing to fractionally occupy several of the one-particle quantum states. In this context, we developed two new indexes useful to reveal outstanding orbitals involved in the chemical reactivity of atoms in molecules.

Quantifying charge state heterogeneity for proteins with multiple ionizable residues

Ionizable residues can release and take up protons and this has an influence on protein structure and function. The extent of protonation is linked to the overall pH of the solution and the local environments of ionizable residues. Binding or unbinding of a single proton generates a distinct charge microstate defined by a specific pattern of charges. Accordingly, the overall partition function is a sum over all charge microstates and Boltzmann weights of all conformations associated with each of the charge microstates. This ensemble-of-ensembles description recast as a q-canonical ensemble allows us to analyze and interpret potentiometric titrations that provide information regarding net charge as a function of pH. In the q-canonical ensemble, charge microstates are grouped into mesostates where each mesostate is a collection of microstates of the same net charge. Here, we show that leveraging the structure of the q-canonical ensemble allows us to decouple contributions of net proton binding and release from proton arrangement and conformational considerations. Through application of the q-canonical formalism to analyze potentiometric measurements of net charge in proteins with repetitive patterns of Lys and Glu residues, we determine the underlying mesostate pKa values and, more importantly, we estimate relative mesostate populations as a function of pH. This is a strength of using the q-canonical approach that cannot be replicated using purely site-specific analyses. Overall, our work shows how measurements of charge equilibria, decoupled from measurements of conformational equilibria, and analyzed using the framework of the q-canonical ensemble, provide protein-specific quantitative descriptions of pH-dependent populations of mesostates. This method is of direct relevance for measuring and understanding how different charge states contribute to conformational, binding, and phase equilibria of proteins.

Investigation of the pKa of the Nucleophilic O2' of the Hairpin Ribozyme

Small ribozymes cleave their RNA phosphodiester backbone by catalyzing a transphosphorylation reaction wherein a specific O2' functions as the nucleophile. While deprotonation of this alcohol through its acidification would increase its nucleophilicity, little is known about the pK_a of this O2' in small ribozymes, in part because high pK_a's are not readily accessible experimentally. Herein, we turn to molecular dynamics to calculate the pK_a of the nucleophilic O2' in the hairpin ribozyme and to study interactions within the active site that may impact its value. We estimate the pK_a of the nucleophilic O2' in the wild-type hairpin ribozyme to be 18.5 ± 0.8, which is higher than the reference compound, and identify a correlation between proper positioning of the O2' for nucleophilic attack and elevation of its pK_a. We find that monovalent ions may play a role in depression of the O2' pK_a, while the exocyclic amine appears to be important for organizing the ribozyme active site. Overall, this study suggests that the pK_a of the O2' is raised in the ground state and lowers during the course of the reaction owing to positioning and metal ion interactions.

Molecular mechanism of sugar transport in plants unveiled by structures of glucose/H+ symporter STP10

Sugars are essential sources of energy and carbon and also function as key signalling molecules in plants. Sugar transport proteins (STP) are proton-coupled symporters responsible for uptake of glucose from the apoplast into plant cells. They are integral to organ development in symplastically isolated tissues such as seed, pollen and fruit. Additionally, STPs play a vital role in plant responses to stressors such as dehydration and prevalent fungal infections like rust and mildew. Here we present a structure of Arabidopsis thaliana STP10 in the inward-open conformation at 2.6 Å resolution and a structure of the outward-occluded conformation at improved 1.8 Å resolution, both with glucose and protons bound. The two structures describe key states in the STP transport cycle. Together with molecular dynamics simulations that establish protonation states and biochemical analysis, they pinpoint structural elements, conserved in all STPs, that clarify the basis of proton-to-glucose coupling. These results advance our understanding of monosaccharide uptake, which is essential for plant organ development, and set the stage for bioengineering strategies in crops.

Reverse Protonation of Buried Ion-Pairs in Staphylococcal Nuclease Mutants

Although buried titratable residues in protein cavities are often of major functional importance, it is generally challenging to understand their properties such as the ionization state and factors of stabilization based on experimental studies alone. A specific set of examples involve buried Glu-Lys pairs in a series of variants of Staphylococcal nuclease, for which recent structural and thermodynamic studies appeared to suggest that both the stability and the ionization state of the buried Glu-Lys pair are sensitive to its orientation (i.e., Glu23-Lys36 vs Lys23-Glu36). To further clarify the situation, especially ionization states of the buried Glu-Lys pairs, we have conducted extensive molecular dynamics simulations and free energy computations. Microsecond molecular dynamics simulations show that the hydration level of the cavity depends on the orientation of the buried ion-pair therein as well as its ionization state; free energy simulations recapitulate the relative stability of Glu23-Lys36 (EK) vs Lys23-Glu36 (KE) mutants measured experimentally, although the difference is similar in magnitude regardless of the ionization state of the Glu-Lys pair. A complementary set of free energy simulations strongly suggests that, in contrast to the original suggestion in the experimental analysis, the Glu and Lys residues prefer to adopt their charge-neutral rather than the ionized states. This result is consistent with the low dielectric constant computed for water in the cavity, which makes it difficult for the protein cavity to stabilize a pair of charged Glu-Lys residues, even with water penetration. The current study highlights the role of free energy simulations in understanding the ionization state of buried titratable residues and the relevant energetic contributions, forming the basis for the rational design of buried charge networks in proteins.

Determination of van der Waals Parameters Using a Double Exponential Potential for Nonbonded Divalent Metal Cations in TIP3P Solvent

A double exponential (DE) functional form for Lennard-Jones (LJ) interactions, proposed in our previous study, has many advantages over LJ potentials including a natural softcore characteristic for the convenience of the pathway-based free-energy calculations, fast convergence, and flexibility in use. In this work, we put the first step on the application of the DE functional form by identifying a DE potential, coined DE-TIP3P, for molecular simulations using the TIP3P water model. The developed DE-TIP3 potential was better than LJ potential in reproducing the experimental water properties. Afterward, we developed the nonbonded models of 15 divalent metal ions, which frequently appear and play vital roles in biological systems, to be consistent with the DE-TIP3P potential and TIP3P water model. Our nonbonded models were as good as the complicated nonbonded dummy cationic models by Jiang et al. and the nonbonded 12-6-4 LJ models by Li and Merz in reproducing the experimental properties of those ions. Moreover, our nonbonded models achieved a better performance than the compromise (CM) LJ models and 12-6-4 LJ models, developed by Li and Merz, in reproducing the properties of MgCl₂ in aqueous solution.

Hybrid MC/MD for protein design

Computational protein design relies on simulations of a protein structure, where selected amino acids can mutate randomly, and mutations are selected to enhance a target property, such as stability. Often, the protein backbone is held fixed and its degrees of freedom are modeled implicitly to reduce the complexity of the conformational space. We present a hybrid method where short molecular dynamics (MD) segments are used to explore conformations and alternate with Monte Carlo (MC) moves that apply mutations to side chains. The backbone is fully flexible during MD. As a test, we computed side chain acid/base constants or pK_a's in five proteins. This problem can be considered a special case of protein design, with protonation/deprotonation playing the role of mutations. The solvent was modeled as a dielectric continuum. Due to cost, in each protein we allowed just one side chain position to change its protonation state and the other position to change its type or mutate. The pK_a's were computed with a standard method that scans a range of pH values and with a new method that uses adaptive landscape flattening (ALF) to sample all protonation states in a single simulation. The hybrid method gave notably better accuracy than standard, fixed-backbone MC. ALF decreased the computational cost a factor of 13.

Biomolecular QM/MM Simulations: What Are Some of the "Burning Issues"?

QM/MM simulations have become an indispensable tool in many chemical and biochemical investigations. Considering the tremendous degree of success, including recognition by a 2013 Nobel Prize in Chemistry, are there still "burning challenges" in QM/MM methods, especially for biomolecular systems? In this short Perspective, we discuss several issues that we believe greatly impact the robustness and quantitative applicability of QM/MM simulations to many, if not all, biomolecules. We highlight these issues with observations and relevant advances from recent studies in our group and others in the field. Despite such limited scope, we hope the discussions are of general interest and will stimulate additional developments that help push the field forward in meaningful directions.

Cyclodipeptide Synthases of the NYH Subfamily Recognize tRNA Using an α-Helix Enriched with Positive Residues

Cyclodipeptide synthases (CDPSs) perform nonribosomal protein synthesis using two aminoacyl-tRNA substrates to produce cyclodipeptides. At present, there are no structural details of the CDPS:tRNA interaction available. Using AlbC, a CDPS that produces cyclo(l-Phe-l-Phe), the interaction between AlbC and its Phe-tRNA substrate was investigated. Simulations of models of the AlbC:tRNA complex, proposed by rigid-body docking or homology modeling, demonstrated that interactions with residues of an AlbC α-helix, α4, significantly contribute to the free energy of binding of AlbC to tRNA. Individual residue contributions to the tRNA binding free energy of the discovered binding mode explain well the available biochemical data, and the results of in vivo assay experiments performed in this work and guided by simulations. In molecular dynamics simulations, the phenylalanyl group predominantly occupied the two positions observed in the experimental structure of AlbC in the dipeptide intermediate state, suggesting that tRNAs of the first and second substrates interact with AlbC in a similar manner. Overall, given the high degree of sequence and structural similarity among the members of the CDPS NYH protein subfamily, the mechanism of the protein:tRNA interaction is expected to be pertinent to a wide range of proteins interacting with tRNA.

Affinity and Selectivity Assessment of Covalent Inhibitors by Free Energy Calculations

Covalent inhibitors have been gaining increased attention in drug discovery due to their beneficial properties such as long residence time, high biochemical efficiency, and specificity. Optimization of covalent inhibitors is a complex task that involves parallel monitoring of the noncovalent recognition elements and the covalent reactivity of the molecules to avoid potential idiosyncratic side effects. This challenge calls for special design protocols, including a variety of computational chemistry methods. Covalent inhibition proceeds through multiple steps, and calculating free energy changes of the subsequent binding events along the overall binding process would help us to better control the design of drug candidates. Inspired by the recent success of free energy calculations on reversible binders, we developed a complex protocol to compute free energies related to the noncovalent and covalent binding steps with thermodynamic integration and hybrid quantum mechanical/molecular mechanical (QM/MM) potential of mean force (PMF) calculations, respectively. In optimization settings, we examined two therapeutically relevant proteins complexed with congeneric sets of irreversible cysteine targeting covalent inhibitors. In the selectivity paradigm, we studied the irreversible binding of covalent inhibitors to phylogenetically close targets by a mutational approach. The results of the calculations are in good agreement with the experimental free energy values derived from the inhibition and kinetic constants (K_i and k_inact) of the enzyme-inhibitor binding. The proposed method might be a powerful tool to predict the potency, selectivity, and binding mechanism of irreversible covalent inhibitors.

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 physics-based energy function allows the computational redesign of a PDZ domain

Computational protein design (CPD) can address the inverse folding problem, exploring a large space of sequences and selecting ones predicted to fold. CPD was used previously to redesign several proteins, employing a knowledge-based energy function for both the folded and unfolded states. We show that a PDZ domain can be entirely redesigned using a "physics-based" energy for the folded state and a knowledge-based energy for the unfolded state. Thousands of sequences were generated by Monte Carlo simulation. Three were chosen for experimental testing, based on their low energies and several empirical criteria. All three could be overexpressed and had native-like circular dichroism spectra and 1D-NMR spectra typical of folded structures. Two had upshifted thermal denaturation curves when a peptide ligand was present, indicating binding and suggesting folding to a correct, PDZ structure. Evidently, the physical principles that govern folded proteins, with a dash of empirical post-filtering, can allow successful whole-protein redesign.

Mechanism of Fully Reversible, pH-Sensitive Inhibition of Human Glutamine Synthetase by Tyrosine Nitration

Glutamine synthetase (GS) catalyzes an ATP-dependent condensation of glutamate and ammonia to form glutamine. This reaction-and therefore GS-are indispensable for the hepatic nitrogen metabolism. Nitration of tyrosine 336 (Y336) inhibits human GS activity. GS nitration and the consequent loss of GS function are associated with a broad range of neurological diseases. The mechanism by which Y336 nitration inhibits GS, however, is not understood. Here, we show by means of unbiased MD simulations, binding, and configurational free energy computations that Y336 nitration hampers ATP binding but only in the deprotonated and negatively charged state of residue 336. By contrast, for the protonated and neutral state, our computations indicate an increased binding affinity for ATP. pK_a computations of nitrated Y336 within GS predict a pK_a of ∼5.3. Thus, at physiological pH, nitrated Y336 exists almost exclusively in the deprotonated and negatively charged state. In vitro experiments confirm these predictions, in that, the catalytic activity of nitrated GS is decreased at pH 7 and 6 but not at pH 4. These results indicate a novel, fully reversible, pH-sensitive mechanism for the regulation of GS activity by tyrosine nitration.

Probing the Accuracy of Explicit Solvent Constant pH Molecular Dynamics Simulations for Peptides

Protonation states of titratable amino acids play a key role in many biomolecular processes. Knowledge of protonatable residue charges at a given pH is essential for a correct understanding of protein catalysis, inter- and intramolecular interactions, substrate binding, and protein dynamics for instance. However, acquiring experimental values for individual amino acid protonation states of complex systems is not straightforward; therefore, several in silico approaches have been developed to tackle this issue. In this work, we assess the accuracy of our previously developed constant pH MD approach by comparing our theoretically obtained pK_a values for titratable residues with experimental values from an equivalent NMR study. We selected a set of four pentapeptides, of adequately small size to ensure comprehensive sampling, but concurrently, due to their charge composition, posing a challenge for protonation state calculation. The comparison of the pK_a values shows good agreement of the experimental and the theoretical approach with a largest difference of 0.25 pK_a units. Further, the corresponding titration curves are in fair agreement, although the shift of the Hill coefficient from a value of 1 was not always reproduced in simulations. The phase space overlap in Cartesian space between trajectories generated in constant pH and standard MD simulations is fair and suggests that our constant pH MD approach reasonably well preserves the dynamics of the system, allowing dynamic protonation MD simulations without introducing structural artifacts.

Thermodynamics of Transition Metal Ion Binding to Proteins

Modeling the thermodynamics of a transition metal (TM) ion assembly be it in proteins or in coordination complexes affords us a better understanding of the assembly and function of metalloclusters in diverse application areas including metal organic framework design, TM-based catalyst design, the trafficking of TM ions in biological systems, and drug design in metalloprotein platforms. While the structural details of TM ions bound to metalloproteins are generally well understood via experimental and computational approaches, accurate studies describing the thermodynamics of TM ion binding are rare. Herein, we demonstrate that we can obtain accurate structural and absolute binding free energies of Co²⁺ and Ni²⁺ to the enzyme glyoxalase I using an optimized 12-6-4 (m12-6-4) potential. Critically, this model simultaneously reproduces the solvation free energy of the individual TM ions and reproduces the thermodynamics of TM ion-ligand coordination as well as the thermodynamics of TM ion binding to a protein active site unlike extant models. We find the incorporation of the thermodynamics associated with protonation state changes for the TM ion (un)binding to be crucial. The high accuracy of m12-6-4 potential in this study presents an accurate route to explore more complicated processes associated with TM cluster assembly and TM ion transport.

Quantum Mechanical and Molecular Mechanics Modeling of Membrane-Embedded Rhodopsins

Computational chemistry provides versatile methods for studying the properties and functioning of biological systems at different levels of precision and at different time scales. The aim of this article is to review the computational methodologies that are applicable to rhodopsins as archetypes for photoactive membrane proteins that are of great importance both in nature and in modern technologies. For each class of computational techniques, from methods that use quantum mechanics for simulating rhodopsin photophysics to less-accurate coarse-grained methodologies used for long-scale protein dynamics, we consider possible applications and the main directions for improvement.

Determination of pKa Values via ab initio Molecular Dynamics and its Application to Transition Metal-Based Water Oxidation Catalysts

The p K a values are important for the in-depth elucidation of catalytic processes, the computational determination of which has been challenging. The first simulation protocols employing ab initio molecular dynamics simulations to calculate p K a values appeared almost two decades ago. Since then several slightly different methods have been proposed. We compare the performance of various evaluation methods in order to determine the most reliable protocol when it comes to simulate p K a values of transition metal-based complexes, such as the here investigated Ru-based water oxidation catalysts. The latter are of high interest for sustainable solar-light driven water splitting, and understanding of the underlying reaction mechanism is crucial for their further development.

Computational structural enzymology methodologies for the study and engineering of fatty acid synthases, polyketide synthases and nonribosomal peptide synthetases

Various computational methodologies can be applied to enzymological studies on enzymes in the fatty acid, polyketide, and non-ribosomal peptide biosynthetic pathways. These multi-domain complexes are called fatty acid synthases, polyketide synthases, and non-ribosomal peptide synthetases. These mega-synthases biosynthesize chemically diverse and complex bioactive molecules, with the intermediates being chauffeured between catalytic partners via a carrier protein. Recent efforts have been made to engineer these systems to expand their product diversity. A major stumbling block is our poor understanding of the transient protein-protein and protein-substrate interactions between the carrier protein and its many catalytic partner domains and product intermediates. The innate reactivity of pathway intermediates in two major classes of polyketide synthases has frustrated our mechanistic understanding of these interactions during the biosynthesis of these natural products, ultimately impeding the engineering of these systems for the generation of engineered natural products. Computational techniques described in this chapter can aid data interpretation or used to generate testable models of these experimentally intractable transient interactions, thereby providing insight into key interactions that are difficult to capture otherwise, with the potential to expand the diversity in these systems.

Conformational sampling and polarization of Asp26 in pKa calculations of thioredoxin

Thioredoxin is a protein that has been used as model system by various computational methods to predict the pK_a of aspartate residue Asp26 which is 3.5 units higher than a solvent exposed one (eg, Asp20). Here, we use extensive atomistic molecular dynamics simulations of two different protonation states of Asp26 in combination with conformational analysis based on RMSD clustering and principle component analysis to identify representative conformations of the protein in solution. For each conformation, the Gibbs free energy of proton transfer between Asp26 and Asp20, which is fully solvated in a loop region of the protein, is calculated with the Amber99sb force field in alchemical transformations. The varying polarization of the two residues in different molecular environments and protonation states is described by Hirshfeld-I (HI) atomic charges obtained from the averaged polarized electron density. Our results show that the Gibbs free energy of proton transfer is dependent on the protein conformation, the proper sampling of the neighboring Lys57 residue orientations and on water molecules entering the hydrophobic cavity upon deprotonating Asp26. The inclusion of the polarization of both aspartate residues in the free energy cycle by HI atomic charges corrects the results from the non-polarizable force field and reproduces the experimental ΔpK_a value of Asp26.

Structural Basis for Isoform Selective Nitric Oxide Synthase Inhibition by Thiophene-2-carboximidamides

The overproduction of nitric oxide in the brain by neuronal nitric oxide synthase (nNOS) is associated with a number of neurodegenerative diseases. Although inhibiting nNOS is an important therapeutic goal, it is important not to inhibit endothelial NOS (eNOS) because of the critical role played by eNOS in maintaining vascular tone. While it has been possible to develop nNOS selective aminopyridine inhibitors, many of the most potent and selective inhibitors exhibit poor bioavailability properties. Our group and others have turned to more biocompatible thiophene-2-carboximidamide (T2C) inhibitors as potential nNOS selective inhibitors. We have used crystallography and computational methods to better understand how and why two commercially developed T2C inhibitors exhibit selectivity for human nNOS over human eNOS. As with many of the aminopyridine inhibitors, a critical active site Asp residue in nNOS versus Asn in eNOS is largely responsible for controlling selectivity. We also present thermodynamic integration results to better understand the change in p Ka and thus the charge of inhibitors once bound to the active site. In addition, relative free energy calculations underscore the importance of enhanced electrostatic stabilization of inhibitors bound to the nNOS active site compared to eNOS.

Enhanced Sampling of Interdomain Motion Using Map-Restrained Langevin Dynamics and NMR: Application to Pin1

Many signaling proteins consist of globular domains connected by flexible linkers that allow for substantial domain motion. Because these domains often serve as complementary functional modules, the possibility of functionally important domain motions arises. To explore this possibility, we require knowledge of the ensemble of protein conformations sampled by interdomain motion. Measurements of NMR residual dipolar couplings (RDCs) of backbone HN bonds offer a per-residue characterization of interdomain dynamics, as the couplings are sensitive to domain orientation. A challenge in reaching this potential is the need to interpret the RDCs as averages over dynamic ensembles of domain conformations. Here, we address this challenge by introducing an efficient protocol for generating conformational ensembles appropriate for flexible, multi-domain proteins. The protocol uses map-restrained self-guided Langevin dynamics simulations to promote collective, interdomain motion while restraining the internal domain motion to near rigidity. Critically, the simulations retain an all-atom description for facile inclusion of site-specific NMR RDC restraints. The result is the rapid generation of conformational ensembles consistent with the RDC data. We illustrate this protocol on human Pin1, a two-domain peptidyl-prolyl isomerase relevant for cancer and Alzheimer's disease. The results include the ensemble of domain orientations sampled by Pin1, as well as those of a dysfunctional variant, I28A-Pin1. The differences between the ensembles corroborate our previous spin relaxation results that showed weakened interdomain contact in the I28A variant relative to wild type. Our protocol extends our abilities to explore the functional significance of protein domain motions.

Physics-based all-atom modeling of RNA energetics and structure

The database of RNA sequences is exploding, but knowledge of energetics, structures, and dynamics lags behind. All-atom computational methods, such as molecular dynamics, hold promise for closing this gap. New algorithms and faster computers have accelerated progress in improving the reliability and accuracy of predictions. Currently, the methods can facilitate refinement of experimentally determined nuclear magnetic resonance and x-ray structures, but are 'unreliable' for predictions based only on sequence. Much remains to be discovered, however, about the many molecular interactions driving RNA folding and the best way to approximate them quantitatively. The large number of parameters required means that a wide variety of experimental results will be required to benchmark force fields and different approaches. As computational methods become more reliable and accessible, they will be used by an increasing number of biologists, much as x-ray crystallography has expanded. Thus, many fundamental physical principles underlying the computational methods are described. This review presents a summary of the current state of molecular dynamics as applied to RNA. It is designed to be helpful to students, postdoctoral fellows, and faculty who are considering or starting computational studies of RNA. WIREs RNA 2017, 8:e1422. doi: 10.1002/wrna.1422.

BAR-based optimum adaptive sampling regime for variance minimization in alchemical transformation: the nonequilibrium stratification

Following the previously proposed equilibrate-state sampling based adaptive sampling regime Optimum Bennett Acceptance Ratio (OBAR), we introduce its nonequilibrium extension, the Optimum Crooks' Equation (OCE) in the current work. The efficiency of the NonEquilibrium Work (NEW) stratification is improved by adaptively manipulating the significance of each nonequilibrium realization followed by importance sampling. As is exhibited in the equilibrium case, the nonequilibrium extension outperforms the simple equal time rule used in nonequilibrium stratification in the sense of minimizing the total variance of the free energy estimate. The speedup of this non-equal time rule is more than 1-fold. The Time Derivative of total Variance (TDV) proposed for the OBAR criterion is extended to determine the importance of each nonequilibrium transformation, which is linearly dependent on the variance. The TDV in the nonequilibrium case gives a totally different importance rank from the standard errors of the free energy differences and OBAR TDV due to the duration of nonequilibrium pulling being added into the OCE equation. The performance of the OCE workflow is demonstrated in the solvation of several small molecules with a series of lambda increments and relaxation times between successive perturbations. To the best of our knowledge, such a nonequilibrium adaptive sampling regime in alchemical transformation is unprecedented.

Entropy and Polarity Control the Partition and Transportation of Drug-like Molecules in Biological Membrane

Partition and transportation of drug in the plasma membrane of a mammalian cell are the prerequisite for its function on target protein. Therefore, comprehensive understanding of the physicochemical properties and mechanism behind these complex phenomena is crucial in pharmaceutical research. By using the state-of-art molecular simulations with polarization effect implicitly or explicitly included, we studied the permeation behavior of 2-aminoethoxydiphenyl borate (2-APB), a broad-spectrum modulator for a number of membrane proteins. We showed that the protonation state and therefore the polarity of the drug is critical for its partition, and that the drug is likely to switch between different protonation states along its permeation pathway. By changing the degrees of freedom, protonation further affects the thermodynamic of the permeation pathway of 2-APB, leading to different entropic contributions. A survey on 54 analog structures with similar backbone to 2-APB showed that delicate balance between entropy and polarity plays an important role in drugs’ potency.

Determination of the Bridging Ligand in the Active Site of Tyrosinase

Tyrosinase is a type-3 copper enzyme that is widely distributed in plants, fungi, insects, and mammals. Developing high potent inhibitors against tyrosinase is of great interest in diverse fields including tobacco curing, food processing, bio-insecticides development, cosmetic development, and human healthcare-related research. In the crystal structure of Agaricus bisporus mushroom tyrosinase, there is an oxygen atom bridging the two copper ions in the active site. It is unclear whether the identity of this bridging oxygen is a water molecule or a hydroxide anion. In the present study, we theoretically determine the identity of this critical bridging oxygen by performing first-principles hybrid quantum mechanics/molecular mechanics/Poisson-Boltzmann-surface area (QM/MM-PBSA) calculations along with a thermodynamic cycle that aim to improve the accuracy. Our results show that the binding with water molecule is energy favored and the QM/MM-optimized structure is very close to the crystal structure, whereas the binding with hydroxide anions causes the increase of energy and significant structural changes of the active site, indicating that the identity of the bridging oxygen must be a water molecule rather than a hydroxide anion. The different binding behavior between water and hydroxide anions may explain why molecules with a carboxyl group or too many negative charges have lower inhibitory activity. In light of this, the design of high potent active inhibitors against tyrosinase should satisfy both the affinity to the copper ions and the charge neutrality of the entire molecule.