Molecular Mechanism of Labelling Functional Cysteines by Heterocyclic Thiones

Heterocyclic thiones have recently been identified as reversible covalent warheads, consistent with their mild electrophilic nature. Little is known so far about their mechanism of action in labelling nucleophilic sidechains, especially cysteines. The vast number of tractable cysteines promotes a wide range of target proteins to examine; however, our focus was put on functional cysteines. We chose the main protease of SARS-CoV-2 harboring Cys145 at the active site that is a structurally characterized and clinically validated target of covalent inhibitors. We screened an in-house, cysteine-targeting covalent inhibitor library which resulted in several covalent fragment hits with benzoxazole, benzothiazole and benzimidazole cores. Thione derivatives and Michael acceptors were selected for further investigations with the objective of exploring the mechanism of inhibition of the thiones and using the thoroughly characterized Michael acceptors for benchmarking our studies. Classical and hybrid quantum mechanical/molecular mechanical (QM/MM) molecular dynamics simulations were carried out that revealed a new mechanism of covalent cysteine labelling by thione derivatives, which was supported by QM and free energy calculations and by a wide range of experimental results. Our study shows that the molecular recognition step plays a crucial role in the overall binding of both sets of molecules.

Assessing the Performance of Non-Equilibrium Thermodynamic Integration in Flavodoxin Redox Potential Estimation

Flavodoxins are enzymes that contain the redox-active flavin mononucleotide (FMN) cofactor and play a crucial role in numerous biological processes, including energy conversion and electron transfer. Since the redox characteristics of flavodoxins are significantly impacted by the molecular environment of the FMN cofactor, the evaluation of the interplay between the redox properties of the flavin cofactor and its molecular surroundings in flavoproteins is a critical area of investigation for both fundamental research and technological advancements, as the electrochemical tuning of flavoproteins is necessary for optimal interaction with redox acceptor or donor molecules. In order to facilitate the rational design of biomolecular devices, it is imperative to have access to computational tools that can accurately predict the redox potential of both natural and artificial flavoproteins. In this study, we have investigated the feasibility of using non-equilibrium thermodynamic integration protocols to reliably predict the redox potential of flavodoxins. Using as a test set the wild-type flavodoxin from Clostridium Beijerinckii and eight experimentally characterized single-point mutants, we have computed their redox potential. Our results show that 75% (6 out of 8) of the calculated reaction free energies are within 1 kcal/mol of the experimental values, and none exceed an error of 2 kcal/mol, confirming that non-equilibrium thermodynamic integration is a trustworthy tool for the quantitative estimation of the redox potential of this biologically and technologically significant class of enzymes.

The effect of Mg2+ on Ca2+ binding to cardiac troponin C in hypertrophic cardiomyopathy associated TNNC1 variants

Cardiac troponin C (cTnC) is the critical Ca2+ -sensing component of the troponin complex. Binding of Ca2+ to cTnC triggers a cascade of conformational changes within the myofilament that culminate in force production. Hypertrophic cardiomyopathy (HCM)-associated TNNC1 variants generally induce a greater degree and duration of Ca2+ binding, which may underly the hypertrophic phenotype. Regulation of contraction has long been thought to occur exclusively through Ca2+ binding to site II of cTnC. However, work by several groups including ours suggest that Mg2+ , which is several orders of magnitude more abundant in the cell than Ca2+ , may compete for binding to the same cTnC regulatory site. We previously used isothermal titration calorimetry (ITC) to demonstrate that physiological concentrations of Mg2+ may decrease site II Ca2+ -binding in both N-terminal and full-length cTnC. Here, we explore the binding of Ca2+ and Mg2+ to cTnC harbouring a series of TNNC1 variants thought to be causal in HCM. ITC and thermodynamic integration (TI) simulations show that A8V, L29Q and A31S elevate the affinity for both Ca2+ and Mg2+ . Further, L48Q, Q50R and C84Y that are adjacent to the EF hand binding motif of site II have a more significant effect on affinity and the thermodynamics of the binding interaction. To the best of our knowledge, this work is the first to explore the role of Mg2+ in modifying the Ca2+ affinity of cTnC mutations linked to HCM. Our results indicate a physiologically significant role for cellular Mg2+ both at baseline and when elevated on modifying the Ca2+ binding properties of cTnC and the subsequent conformational changes which precede cardiac contraction.

Free energy simulations to study mutational effect of a conserved residue, Trp24, on stability of human serum retinol-binding protein

Human serum retinol-binding protein (RBP) is a plasma transport protein for vitamin A. RBP is a prime subclass of lipocalins, which bind nonpolar ligands within a β-barrel. To understand the role of Trp 24, one of the highly conserved residues in RBP, free energy simulations have been carried out to understand the effects of the mutations from Trp at position 24 to Leu, Phe, and Tyr in the apo-RBP on its thermal stability. We examine various unfolded systems to study the dependence of the free energy differences on the denatured structure. Our calculated free energy difference values for the three mutations are in excellent agreement with the experimental values when the initial coordinates of the seven-residue peptide segments truncated from the crystal structure are used for the denatured systems. Our free energy change differences for the Trp→Leu, Trp→Phe, and Trp→Tyr mutations are 2.50 ± 0.69, 2.58 ± 0.50, and 2.49 ± 0.48 kcal/mol, respectively, when the native-like seven-residue peptides are used as models for the denatured systems. The main contributions to the free energy change differences for the Trp24→Leu and Trp24→Phe mutations are mainly from van der Waals and covalent interactions, respectively. Electrostatic, van der Waals and covalent terms equally contribute to the free energy change difference for the Trp24→Tyr mutation. The free energy simulation helps understand the detailed microscopic mechanism of the stability of the RBP mutants relative to the wild type and the role of the highly conserved residue, Trp24, of the human RBP.Communicated by Ramaswamy H. Sarma.

The phase diagrams of beryllium and magnesium oxide at megabar pressures

We performab initiosimulations of beryllium (Be) and magnesium oxide (MgO) at megabar pressures and compare their structural and thermodynamic properties. We make a detailed comparison of our two recently derived phase diagrams of Be (Wuet al2021Phys. Rev.B104014103) and MgO (Soubiran and Militzer 2020Phys. Rev. Lett.125175701) using the thermodynamic integration technique, as they exhibit striking similarities regarding their shape. We explore whether the Lindemann criterion can explain the melting temperatures of these materials through the calculation of the Debye temperature at high pressure. From our free energy calculations, we find that the melting line of both materials is well represented by the Simon-Glazel fitT_m(P) =T₀(1 +P/a)^1/c, whereT₀= 1564 K,a= 15.8037 GPa andc= 2.4154 for Be, whileT₀= 3010 K,a= 10.5797 GPa andc= 2.8683 for the MgO in the B1. For the B2 phase, we use the valuesa= 26.1163 GPa andc= 2.2426. Both materials exhibit negative Clapeyron slopes on the boundaries between the two solid phases that are strongly affected by anharmonic effects, which also influence the location of the solid-solid-liquid triple point. We find that the quasi-harmonic approximation underestimates the stability range of the low-pressure phases, namely hcp for Be and B1 for MgO. We also compute the phonon dispersion relations at low and high pressure for each of the phases of these materials, and also explore how the phonon density of states is modified by temperature. Finally, we derive secondary shock Hugoniot curves in addition to the principal Hugoniot curve for both materials, and study their offsets in pressure between solid and liquid branches.

Tissue-Nonspecific Alkaline Phosphatase (TNAP) as the Enzyme Involved in the Degradation of Nucleotide Analogues in the Ligand Docking and Molecular Dynamics Approaches

Tissue-nonspecific alkaline phosphatase (TNAP) is known to be involved in the degradation of extracellular ATP via the hydrolysis of pyrophosphate (PPi). We investigated, using three different computational methods, namely molecular docking, thermodynamic integration (TI) and conventional molecular dynamics (MD), whether TNAP may also be involved in the utilization of β,γ-modified ATP analogues. For that, we analyzed the interaction of bisphosphonates with this enzyme and evaluated the obtained structures using in silico studies. Complexes formed between pyrophosphate, hypophosphate, imidodiphosphate, methylenediphosphonic acid monothiopyrophosphate, alendronate, pamidronate and zoledronate with TNAP were generated and analyzed based on ligand docking, molecular dynamics and thermodynamic integration. The obtained results indicate that all selected ligands show high affinity toward this enzyme. The forming complexes are stabilized through hydrogen bonds, electrostatic interactions and van der Waals forces. Short- and middle-term molecular dynamics simulations yielded very similar affinity results and confirmed the stability of the protein and its complexes. The results suggest that certain effectors may have a significant impact on the enzyme, changing its properties.

Binding of calcium and magnesium to human cardiac troponin C

Cardiac muscle thin filaments are composed of actin, tropomyosin, and troponin that change conformation in response to Ca²⁺ binding, triggering muscle contraction. Human cardiac troponin C (cTnC) is the Ca²⁺-sensing component of the thin filament. It contains structural sites (III/IV) that bind both Ca²⁺ and Mg²⁺ and a regulatory site (II) that has been thought to bind only Ca²⁺. Binding of Ca²⁺ at this site initiates a series of conformational changes that culminate in force production. However, the mechanisms that underpin the regulation of binding at site II remain unclear. Here, we have quantified the interaction between site II and Ca²⁺/Mg²⁺ through isothermal titration calorimetry and thermodynamic integration simulations. Direct and competitive binding titrations with WT N-terminal cTnC and full-length cTnC indicate that physiologically relevant concentrations of both Ca²⁺/Mg²⁺ interacted with the same locus. Moreover, the D67A/D73A N-terminal cTnC construct in which two coordinating residues within site II were removed was found to have significantly reduced affinity for both cations. In addition, 1 mM Mg²⁺ caused a 1.4-fold lower affinity for Ca²⁺. These experiments strongly suggest that cytosolic-free Mg²⁺ occupies a significant population of the available site II. Interaction of Mg²⁺ with site II of cTnC likely has important functional consequences for the heart both at baseline as well as in diseased states that decrease or increase the availability of Mg²⁺, such as secondary hyperparathyroidism or ischemia, respectively.

Generating Excess Protons in Microsolvated Acid Clusters under Ambient Conditions: An Issue of Configurational Entropy versus Internal Energy

Acid dissociation, and thus liberation of excess protons in small water droplets, impacts on diverse fields such as interstellar, atmospheric or environmental chemistry. At cryogenic temperatures below 1 K, it is now well established that as few as four water molecules suffice to dissociate the generic strong acid HCl, yet temperature-driven recombination sets in simply upon heating that cluster. Here, the fundamental question is posed of how many more water molecules are required to stabilize a hydrated excess proton at room temperature. Ab initio path integral simulations disclose that not five, but six water molecules are needed at 300 K to allow for HCl dissociation independently from nuclear quantum effects. In order to provide the molecular underpinnings of these observations, the classical and quantum free energy profiles were decomposed along the dissociation coordinate in terms of the corresponding internal energy and entropy profiles. What decides in the end about acid dissociation, and thus ion pair formation, in a specific microsolvated water cluster at room temperature is found to be a fierce competition between classical configurational entropy and internal energy, where the former stabilizes the undissociated state whereas the latter favors dissociation. It is expected that these are generic findings with broad implications on acid-base chemistry depending on temperature in small water assemblies.

Thermodynamic Integration in 3n Dimensions without Biases or Alchemy for Protein Interactions

Thermodynamic integration (TI), a powerful formalism for computing Gibbs free energy, has been implemented for many biophysical processes with alchemical schemes that require delicate human efforts to choose/design biasing potentials for sampling the desired biophysical events and to remove their artifactitious consequences afterwards. Theoretically, an alchemical scheme is exact but practically, an unsophisticated implementation of this exact formula can cause error amplifications. Small relative errors in the input parameters can be amplified many times in their propagation into the computed free energy [due to subtraction of similar numbers such as (105 ± 5)‒(100 ± 5) = 5 ± 7]. In this paper, we present an unsophisticated implementation of TI in 3n dimensions (3nD) (n=1,2,3…) for the potential of mean force along a 3nD path connecting one state in the bound state ensemble to one state in the unbound state ensemble. Fluctuations in these 3nD are integrated in the bound and unbound state ensembles but not along the 3nD path. Using TI3nD, we computed the standard binding free energies of three protein complexes: trometamol in Salmonella effector SpvD (n=1), biotin-avidin (n=2), and Colicin E9 endonuclease with cognate immunity protein Im9 (n=3). We employed three different protocols in three independent computations of E9-Im9 to show TI3nD's robustness. We also computed the hydration energies of ten biologically relevant compounds (n=1 for water, acetamide, urea, glycerol, trometamol, ammonium and n=2 for erythritol, 1,3-propanediol, xylitol, biotin). Each of the 15 computations is accomplishable within one (for hydration) to ten (for E9-Im9) days on an inexpensive GPU workstation. The computed results all agree with the available experimental data.

Molecular dynamics study of solubilization of cyclohexane, benzene, and phenol into mixed micelles composed of sodium dodecyl sulfate and octaethylene glycol monododecyl ether

Molecular dynamics calculations of a mixed micelle composed of sodium dodecyl sulfate (SDS) and octaethylene glycol monododecyl ether (C₁₂ E₈ ) were performed for six compositions (SDS/C₁₂ E₈ = 100/0, 80/20, 60/40, 40/60, 20/80, and 0/100) to investigate the composition dependence of the mixed micelle structure and solubilization of cyclohexane, benzene, and phenol molecules by the micelle. The radial density distribution of the hydrophilic polyoxyethylene (POE) group of C₁₂ E₈ as a function of distance from the micelle center is very sharp for micelles with high SDS content because the POE group captures a Na⁺ ion in solution and wraps around it to form a compact crown-ether-like complex. The hydrophobic dodecyl groups of SDS and C₁₂ E₈ were separately distributed in the mixed micelle core. ΔG(r) evaluated for each solute showed that despite the structural changes of the micelle the binding strength of the solute molecules to the micelle did not change significantly. © 2019 Wiley Periodicals, Inc.

Studying effects of different protonation states of His11 and His102 in ribose-5-phosphate isomerase of Trypanosoma cruzi: an example of cooperative behavior

The Trypanosoma cruzi ribose-5-phosphate isomerase B (TcRpiB) is a crucial piece in the pentose phosphate pathway and thus is a potential drug target for treatment of Chagas' disease. TcRpiB residues, such as Cys69, Asp45, Glu149 and Pro47, have confirmed their roles in substrate recognition, catalytic reaction and binding site conformation. However, the joint performance of His11 and His102, in the D-ribose-5-phosphate (R5P) in the catalysis is not well understood. In this work, we probed the influence of different protonation states of His11 and His102 on the behavior of the ligand R5P using molecular dynamics simulations, network analysis and thermodynamic integration. Simulations revealed that a protonated His11 combined with a neutral His102 (His11⁺‒His102) was able to stabilize the ligand R5P in the binding site. Moreover, calculated relative free energy differences showed that when protonated His11 was coupled to a neutral His102 an exergonic process takes place. On the other hand, neutral His11 combined with a protonated His102 (His11‒His102⁺), sampled conformations that resembled the catalyzed product D-ribulose-5-phosphate (Ru5P). Network analysis also demonstrated some peculiarities for these systems with some negatively correlated nodes in the binding site for His11‒His102⁺, and exclusive suboptimal paths for His11⁺‒His102. Therefore, the combined approach presented in this paper proposes two suitable protonation states for the TcRpiB catalytic mechanism, where an extra proton in either histidines might favor R5P binding or influence isomerization reaction to Ru5P. Our results may guide further in silico drug discovery studies. Communicated by Ramaswamy H. Sarma.

Configurational Entropy in Multicomponent Alloys: Matrix Formulation from Ab Initio Based Hamiltonian and Application to the FCC Cr-Fe-Mn-Ni System

Configuration entropy is believed to stabilize disordered solid solution phases in multicomponent systems at elevated temperatures over intermetallic compounds by lowering the Gibbs free energy. Traditionally, the increment of configuration entropy with temperature was computed by time-consuming thermodynamic integration methods. In this work, a new formalism based on a hybrid combination of the Cluster Expansion (CE) Hamiltonian and Monte Carlo simulations is developed to predict the configuration entropy as a function of temperature from multi-body cluster probability in a multi-component system with arbitrary average composition. The multi-body probabilities are worked out by explicit inversion and direct product of a matrix formulation within orthonomal sets of point functions in the clusters obtained from symmetry independent correlation functions. The matrix quantities are determined from semi canonical Monte Carlo simulations with Effective Cluster Interactions (ECIs) derived from Density Functional Theory (DFT) calculations. The formalism is applied to analyze the 4-body cluster probabilities for the quaternary system Cr-Fe-Mn-Ni as a function of temperature and alloy concentration. It is shown that, for two specific compositions (Cr 25Fe 25Mn 25Ni 25 and Cr 18Fe 27Mn 27Ni 28), the high value of probabilities for Cr-Fe-Fe-Fe and Mn-Mn-Ni-Ni are strongly correlated with the presence of the ordered phases L1 2 -CrFe 3 and L1 0-MnNi, respectively. These results are in an excellent agreement with predictions of these ground state structures by ab initio calculations. The general formalism is used to investigate the configuration entropy as a function of temperature and for 285 different alloy compositions. It is found that our matrix formulation of cluster probabilities provides an efficient tool to compute configuration entropy in multi-component alloys in a comparison with the result obtained by the thermodynamic integration method. At high temperatures, it is shown that many-body cluster correlations still play an important role in understanding the configuration entropy before reaching the solid solution limit of high-entroy alloys (HEAs).

Rapid and accurate structure-based therapeutic peptide design using GPU accelerated thermodynamic integration

Peptide-based therapeutics are an alternative to small molecule drugs as they offer superior specificity, lower toxicity, and easy synthesis. Here we present an approach that leverages the dramatic performance increase afforded by the recent arrival of GPU accelerated thermodynamic integration (TI). GPU TI facilitates very fast, highly accurate binding affinity optimization of peptides against therapeutic targets. We benchmarked TI predictions using published peptide binding optimization studies. Prediction of mutations involving charged side-chains was found to be less accurate than for non-charged, and use of a more complex 3-step TI protocol was found to boost accuracy in these cases. Using the 3-step protocol for non-charged side-chains either had no effect or was detrimental. We use the benchmarked pipeline to optimize a peptide binding to our recently discovered cancer target: EME1. TI calculations predict beneficial mutations using both canonical and non-canonical amino acids. We validate these predictions using fluorescence polarization and confirm that binding affinity is increased. We further demonstrate that this increase translates to a significant reduction in pancreatic cancer cell viability.

Structural consequences of hereditary spastic paraplegia disease-related mutations in kinesin

A wide range of mutations in the kinesin motor Kif5A have been linked to a neuronal disorder called hereditary spastic paraplegia (HSP). The position of these mutations can vary, and a range of different motile behaviors have been observed, indicating that the HSP mutants can alter distinct aspects of kinesin mechanochemistry. While focusing on four key HSP-associated mutants, this study examined the structural and dynamic perturbations that arise from these mutations using a series of different computational methods, ranging from bioinformatics analyses to all-atom simulations, that account for solvent effects explicitly. We show that two catalytic domain mutations (R280S and K253N) reduce the microtubule (MT) binding affinity of the kinesin head domains appreciably, while N256S has a much smaller impact. Bioinformatics analysis suggests that the stalk mutation A361V perturbs motor dimerization. Subsequent integration of these effects into a coarse-grained structure-based model of dimeric kinesin revealed that the order-disorder transition of the neck linker is substantially affected, indicating a hampered directionality and processivity of kinesin. The present analyses therefore suggest that, in addition to kinesin-MT binding and coiled-coil dimerization, HSP mutations affecting motor stepping transitions and processivity can lead to disease.

Revised mechanism of carboxylation of ribulose-1,5-biphosphate by rubisco from large scale quantum chemical calculations

Here, we describe a computational approach for studying enzymes that catalyze complex multi-step reactions and apply it to Ribulose 1,5-bisphosphate carboxylase-oxygenase (Rubisco), the enzyme that fixes atmospheric carbon dioxide within photosynthesis. In the 5-step carboxylase reaction, the substrate Ribulose-1,5-bisphosphate (RuBP) first binds Rubisco and undergoes enolization before binding the second substrate, CO₂ . Hydration of the RuBP.CO₂ complex is followed by CC bond scission and stereospecific protonation. However, details of the roles and protonation states of active-site residues, and sources of protons and water, remain highly speculative. Large-scale computations on active-site models provide a means to better understand this complex chemical mechanism. The computational protocol comprises a combination of hybrid semi-empirical quantum mechanics and molecular mechanics within constrained molecular dynamics simulations, together with constrained gradient minimization calculations using density functional theory. Alternative pathways for hydration of the RuBP.CO₂ complex and associated active-site protonation networks and proton and water sources were investigated. The main findings from analysis of the resulting energetics advocate major revision to existing mechanisms such that: hydration takes place anti to the CO₂ ; both hydration and CC bond scission require early protonation of CO₂ in the RuBP.CO₂ complex; CC bond scission and stereospecific protonation reactions are concerted and, effectively, there is only one stable intermediate, the C3-gemdiolate complex. Our main conclusions for interpreting enzyme kinetic results are that the gemdiolate may represent the elusive Michaelis-Menten-like complex corresponding to the empirical K_m (=K_c ) with turnover to product via bond scission concerted with stereospecific protonation consistent with the observed catalytic rate. © 2018 Wiley Periodicals, Inc.

An Accurate Estimate of the Free Energy and Phase Diagram of All-DNA Bulk Fluids

We present a numerical study in which large-scale bulk simulations of self-assembled DNA constructs have been carried out with a realistic coarse-grained model. The investigation aims at obtaining a precise, albeit numerically demanding, estimate of the free energy for such systems. We then, in turn, use these accurate results to validate a recently proposed theoretical approach that builds on a liquid-state theory, the Wertheim theory, to compute the phase diagram of all-DNA fluids. This hybrid theoretical/numerical approach, based on the lowest-order virial expansion and on a nearest-neighbor DNA model, can provide, in an undemanding way, a parameter-free thermodynamic description of DNA associating fluids that is in semi-quantitative agreement with experiments. We show that the predictions of the scheme are as accurate as those obtained with more sophisticated methods. We also demonstrate the flexibility of the approach by incorporating non-trivial additional contributions that go beyond the nearest-neighbor model to compute the DNA hybridization free energy.

Bayesian model evidence as a practical alternative to deviance information criterion

While model evidence is considered by Bayesian statisticians as a gold standard for model selection (the ratio in model evidence between two models giving the Bayes factor), its calculation is often viewed as too computationally demanding for many applications. By contrast, the widely used deviance information criterion (DIC), a different measure that balances model accuracy against complexity, is commonly considered a much faster alternative. However, recent advances in computational tools for efficient multi-temperature Markov chain Monte Carlo algorithms, such as steppingstone sampling (SS) and thermodynamic integration schemes, enable efficient calculation of the Bayesian model evidence. This paper compares both the capability (i.e. ability to select the true model) and speed (i.e. CPU time to achieve a given accuracy) of DIC with model evidence calculated using SS. Three important model classes are considered: linear regression models, mixed models and compartmental models widely used in epidemiology. While DIC was found to correctly identify the true model when applied to linear regression models, it led to incorrect model choice in the other two cases. On the other hand, model evidence led to correct model choice in all cases considered. Importantly, and perhaps surprisingly, DIC and model evidence were found to run at similar computational speeds, a result reinforced by analytically derived expressions.

Effect of Atomic Charges on Octanol-Water Partition Coefficient Using Alchemical Free Energy Calculation

The octanol-water partition coefficient (logPow) is an important index for measuring solubility, membrane permeability, and bioavailability in the drug discovery field. In this paper, the logPow values of 58 compounds were predicted by alchemical free energy calculation using molecular dynamics simulation. In free energy calculations, the atomic charges of the compounds are always fixed. However, they must be recalculated for each solvent. Therefore, three different sets of atomic charges were tested using quantum chemical calculations, taking into account vacuum, octanol, and water environments. The calculated atomic charges in the different environments do not necessarily influence the correlation between calculated and experimentally measured ∆Gwater values. The largest correlation coefficient values of the solvation free energy in water and octanol were 0.93 and 0.90, respectively. On the other hand, the correlation coefficient of logPow values calculated from free energies, the largest of which was 0.92, was sensitive to the combination of the solvation free energies calculated from the calculated atomic charges. These results reveal that the solvent assumed in the atomic charge calculation is an important factor determining the accuracy of predicted logPow values.

Investigation of the bindings of a class of inhibitors with GSK3β kinase using thermodynamic integration MD simulation and kinase assay

GSK3β kinase is a noteworthy target for discovery of the drugs that will be used to treat several diseases. In the effort to identify a new inhibitor lead compound, we utilized thermodynamic integration (TI)-molecular dynamics (MD) simulation and kinase assay to investigate the bindings between GSK3β kinase and five compounds that were analogous to a known inhibitor with an available crystal structure. TI-MD simulations of the first two compounds (analogs 1 and 2) were used for calibration. The computed binding affinities of analogs 1 and 2 agreed well with the experimental results. The rest three compounds (analogs 3-5) were newly obtained from a database search, and their affinity data were newly measured in our labs. TI-MD simulations predicted the binding modes and the computed ΔΔG values have a reasonably good correlation with the experimental affinity data. These newly identified inhibitors appear to be new leads according to our survey of GSK3β inhibitors listed in recent review articles. The predicted binding modes of these compounds should aid in designing new derivatives of these compounds in the future.

Experimental and Computational Analysis of Protein Stabilization by Gly-to-d-Ala Substitution: A Convolution of Native State and Unfolded State Effects

The rational and predictable enhancement of protein stability is an important goal in protein design. Most efforts target the folded state, however stability is the free energy difference between the folded and unfolded states thus both are suitable targets. Strategies directed at the unfolded state usually seek to decrease chain entropy by introducing cross-links or by replacing glycines. Cross-linking has led to mixed results. Replacement of glycine with an l-amino acid, while reducing the entropy of the unfolded state, can introduce unfavorable steric interactions in the folded state, since glycine is often found in conformations that require a positive φ angle such as helical C-capping motifs or type I' and II″ β-turns. l-Amino acids are strongly disfavored in these conformations, but d-amino acids are not. However, there are few reported examples and conflicting results have been obtained when glycines are replaced with d-Ala. We critically examine the effect of Gly-to-d-Ala substitutions on protein stability using experimental approaches together with molecular dynamics simulations and free energy calculations. The data, together with a survey of high resolution structures, show that the vast majority of proteins can be stabilized by substitution of C-capping glycines with d-Ala. Sites suitable for substitutions can be identified via sequence alignment with a high degree of success. Steric clashes in the native state due to the new side chain are rarely observed, but are likely responsible for the destabilizing or null effect observed for the small subset of Gly-to-d-Ala substitutions which are not stabilizing. Changes in backbone solvation play less of a role. Favorable candidates for d-Ala substitution can be identified using a rapid algorithm based on molecular mechanics.

Combined Simulation and Mutation Studies to Elucidate Selectivity of Unsubstituted Amphetamine-like Cathinones at the Dopamine Transporter

The dopamine and serotonin transporter proteins (DAT, SERT) play a vital role in behavior and mental illness. Although their substrate transport has been studied extensively, the molecular basis of their selectivity is not completely understood yet. In this study, we exploit molecular dynamics simulations combined with mutagenesis studies to shed light on the driving factors for DAT-over-SERT selectivity of a set of cathinones. Results indicate that these compounds can adopt two binding modes of which one is more favorable. In addition, free energy calculations indicated the substrate binding site (S1) as the primary recognition site for these ligands. By simulating DAT with SERT-like mutations, we hypothesize unsubstituted cathinones to bind more favorably to DAT, due to a Val152 offering more space, as compared to the bulkier Ile172 in SERT. This was supported by uptake inhibition measurements, which showed an increase in activity in SERT-I172V.

Free-energy calculations of residue mutations in a tripeptide using various methods to overcome inefficient sampling

Previous free-energy calculations have shown that the seemingly simple transformation of the tripeptide KXK to KGK in water holds some unobvious challenges concerning the convergence of the forward and backward thermodynamic integration processes (i.e., hysteresis). In the current study, the central residue X was either alanine, serine, glutamic acid, lysine, phenylalanine, or tyrosine. Interestingly, the transformation from alanine to glycine yielded the highest hysteresis in relation to the extent of the chemical change of the side chain. The reason for that could be attributed to poor sampling of φ2 /ψ2 dihedral angles along the transformation. Altering the nature of alanine's Cβ atom drastically improved the sampling and at the same time led to the identification of high energy barriers as cause for it. Consequently, simple strategies to overcome these barriers are to increase simulation time (computationally expensive) or to use enhanced sampling techniques such as Hamiltonian replica exchange molecular dynamics and one-step perturbation. © 2016 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.

Estimating the Number of Subpopulations (K) in Structured Populations

A key quantity in the analysis of structured populations is the parameter K, which describes the number of subpopulations that make up the total population. Inference of K ideally proceeds via the model evidence, which is equivalent to the likelihood of the model. However, the evidence in favor of a particular value of K cannot usually be computed exactly, and instead programs such as Structure make use of heuristic estimators to approximate this quantity. We show-using simulated data sets small enough that the true evidence can be computed exactly-that these heuristics often fail to estimate the true evidence and that this can lead to incorrect conclusions about K Our proposed solution is to use thermodynamic integration (TI) to estimate the model evidence. After outlining the TI methodology we demonstrate the effectiveness of this approach, using a range of simulated data sets. We find that TI can be used to obtain estimates of the model evidence that are more accurate and precise than those based on heuristics. Furthermore, estimates of K based on these values are found to be more reliable than those based on a suite of model comparison statistics. Finally, we test our solution in a reanalysis of a white-footed mouse data set. The TI methodology is implemented for models both with and without admixture in the software MavericK1.0.

Calculation of Relative Binding Free Energy in the Water-Filled Active Site of Oligopeptide-Binding Protein A

The periplasmic oligopeptide binding protein A (OppA) represents a well-known example of water-mediated protein-ligand interactions. Here, we perform free-energy calculations for three different ligands binding to OppA, using a thermodynamic integration approach. The tripeptide ligands share a high structural similarity (all have the sequence KXK), but their experimentally-determined binding free energies differ remarkably. Thermodynamic cycles were constructed for the ligands, and simulations conducted in the bound and (freely solvated) unbound states. In the unbound state, it was observed that the difference in conformational freedom between alanine and glycine leads to a surprisingly slow convergence, despite their chemical similarity. This could be overcome by increasing the softness parameter during alchemical transformations. Discrepancies remained in the bound state however, when comparing independent simulations of the three ligands. These difficulties could be traced to a slow relaxation of the water network within the active site. Fluctuations in the number of water molecules residing in the binding cavity occur mostly on a timescale larger than the simulation time along the alchemical path. After extensive simulations, relative binding free energies that were converged to within thermal noise could be obtained, which agree well with available experimental data.

Evaluating thermodynamic integration performance of the new amber molecular dynamics package and assess potential halogen bonds of enoyl-ACP reductase (FabI) benzimidazole inhibitors

Thermodynamic integration (TI) can provide accurate binding free energy insights in a lead optimization program, but its high computational expense has limited its usage. In the effort of developing an efficient and accurate TI protocol for FabI inhibitors lead optimization program, we carefully compared TI with different Amber molecular dynamics (MD) engines (sander and pmemd), MD simulation lengths, the number of intermediate states and transformation steps, and the Lennard-Jones and Coulomb Softcore potentials parameters in the one-step TI, using eleven benzimidazole inhibitors in complex with Francisella tularensis enoyl acyl reductase (FtFabI). To our knowledge, this is the first study to extensively test the new AMBER MD engine, pmemd, on TI and compare the parameters of the Softcore potentials in the one-step TI in a protein-ligand binding system. The best performing model, the one-step pmemd TI, using 6 intermediate states and 1 ns MD simulations, provides better agreement with experimental results (RMSD = 0.52 kcal/mol) than the best performing implicit solvent method, QM/MM-GBSA from our previous study (RMSD = 3.00 kcal/mol), while maintaining similar efficiency. Briefly, we show the optimized TI protocol to be highly accurate and affordable for the FtFabI system. This approach can be implemented in a larger scale benzimidazole scaffold lead optimization against FtFabI. Lastly, the TI results here also provide structure-activity relationship insights, and suggest the parahalogen in benzimidazole compounds might form a weak halogen bond with FabI, which is a well-known halogen bond favoring enzyme.

Insights into regioselective metabolism of mefenamic acid by cytochrome P450 BM3 mutants through crystallography, docking, molecular dynamics, and free energy calculations

Cytochrome P450 BM3 (CYP102A1) mutant M11 is able to metabolize a wide range of drugs and drug-like compounds. Among these, M11 was recently found to be able to catalyze formation of human metabolites of mefenamic acid and other nonsteroidal anti-inflammatory drugs (NSAIDs). Interestingly, single active-site mutations such as V87I were reported to invert regioselectivity in NSAID hydroxylation. In this work, we combine crystallography and molecular simulation to study the effect of single mutations on binding and regioselective metabolism of mefenamic acid by M11 mutants. The heme domain of the protein mutant M11 was expressed, purified, and crystallized, and its X-ray structure was used as template for modeling. A multistep approach was used that combines molecular docking, molecular dynamics (MD) simulation, and binding free-energy calculations to address protein flexibility. In this way, preferred binding modes that are consistent with oxidation at the experimentally observed sites of metabolism (SOMs) were identified. Whereas docking could not be used to retrospectively predict experimental trends in regioselectivity, we were able to rank binding modes in line with the preferred SOMs of mefenamic acid by M11 and its mutants by including protein flexibility and dynamics in free-energy computation. In addition, we could obtain structural insights into the change in regioselectivity of mefenamic acid hydroxylation due to single active-site mutations. Our findings confirm that use of MD and binding free-energy calculation is useful for studying biocatalysis in those cases in which enzyme binding is a critical event in determining the selective metabolism of a substrate.

Multiscale computational models in physical systems biology of intracellular trafficking

In intracellular trafficking, a definitive understanding of the interplay between protein binding and membrane morphology remains incomplete. The authors describe a computational approach by integrating coarse-grained molecular dynamics (CGMD) simulations with continuum Monte Carlo (CM) simulations of the membrane to study protein-membrane interactions and the ensuing membrane curvature. They relate the curvature field strength discerned from the molecular level to its effect at the cellular length-scale. They perform thermodynamic integration on the CM model to describe the free energy landscape of vesiculation in clathrin-mediated endocytosis. The method presented here delineates membrane morphologies and maps out the free energy changes associated with membrane remodeling due to varying coat sizes, coat curvature strengths, membrane bending rigidities, and tensions; furthermore several constraints on mechanisms underlying clathrin-mediated endocytosis have also been identified, Their CGMD simulations have revealed the importance of PIP2 for stable binding of proteins essential for curvature induction in the bilayer and have provided a molecular basis for the positive curvature induction by the epsin N-terminal homology (EIMTH) domain. Calculation of the free energy landscape for vesicle budding has identified the critical size and curvature strength of a clathrin coat required for nucleation and stabilisation of a mature vesicle.

Importance of polarization effect in the study of metalloproteins: application of polarized protein specific charge scheme in predicting the reduction potential of azurin

Molecular dynamics (MD) simulation is commonly used in the study of protein dynamics, and in recent years, the extension of MD simulation to the study of metalloproteins is gaining much interest. Choice of force field is crucial in MD studies, and the inclusion of metal centers complicates the process of accurately describing the electrostatic environment that surrounds the redox centre. Herein, we would like to explore the importance of including electrostatic contribution from both protein and solvent in the study of metalloproteins. MD simulations with the implementation of thermodynamic integration will be conducted to model the reduction process of azurin from Pseudomonas aeruginosa. Three charge schemes will be used to derive the partial charges of azurin. These charge schemes differ in terms of the amount of immediate environment, respective to copper, considered during charge fitting, which ranges from the inclusion of copper and residues in the first coordination sphere during density functional theory charge fitting to the comprehensive inclusion of protein and solvent effect surrounding the metal centre using polarized protein-specific charge scheme. From the simulations conducted, the relative reduction potential of the mutated azurins respective to that of wild-type azurin (ΔEcal) were calculated and compared with experimental values. The ΔEcal approached experimental value with increasing consideration of environmental effect hence substantiating the importance of polarization effect in the study of metalloproteins. This study also attests the practicality of polarized protein-specific charge as a computational tool capable of incorporating both protein environment and solvent effect into MD simulations.

Protonation states of active-site lysines of penicillin-binding protein 6 from Escherichia coli and the mechanistic implications

The protonation states of the two active-site lysines (Lys69 and Lys235) of PBP 6 of Escherichia coli were explored to understand the active site chemistry of this enzyme. Each lysine was individually mutated to cysteine, and the resultant two mutant proteins were purified to homogeneity. Each protein was denatured, and its cysteine was chemically modified to produce an S-aminoethylated cysteine (γ-thialysine) residue. Following renaturation, the evaluation of the kinetics of the dd-carboxypeptidase activity of PBP 6 as a function of pH was found consistent with one lysine in its free-base (Lys69) and the other in the protonated state (Lys235) for optimal catalysis. The experimental estimates for their pKa values were compared with the pKa values calculated computationally, using molecular-dynamics simulations and a thermodynamic cycle. Study of the γ-thialysine69 showed that lysine at position 69 influenced the basic limb of catalysis, consistent with the fact that the two lysine side chains are in proximity to each other in the active site. Based on these observations, a reaction sequence for PBP 6 is proposed, wherein protonated Lys235 serves as the electrostatic substrate anchor and Lys69 as the conduit for protons in the course of the acylation and deacylation half-reactions.

pH-triggered conformational switching of the diphtheria toxin T-domain: the roles of N-terminal histidines

pH-induced conformational switching is essential for functioning of diphtheria toxin, which undergoes a membrane insertion/translocation transition triggered by endosomal acidification as a key step of cellular entry. In order to establish the sequence of molecular rearrangements and side-chain protonation accompanying the formation of the membrane-competent state of the toxin's translocation (T) domain, we have developed and applied an integrated approach that combines multiple techniques of computational chemistry [e.g., long-microsecond-range, all-atom molecular dynamics (MD) simulations; continuum electrostatics calculations; and thermodynamic integration (TI)] with several experimental techniques of fluorescence spectroscopy. TI calculations indicate that protonation of H257 causes the greatest destabilization of the native structure (6.9 kcal/mol), which is consistent with our early mutagenesis results. Extensive equilibrium MD simulations with a combined length of over 8 μs demonstrate that histidine protonation, while not accompanied by the loss of structural compactness of the T-domain, nevertheless results in substantial molecular rearrangements characterized by the partial loss of secondary structure due to unfolding of helices TH1 and TH2 and the loss of close contact between the C- and N-terminal segments. The structural changes accompanying the formation of the membrane-competent state ensure an easier exposure of the internal hydrophobic hairpin formed by helices TH8 and TH9, in preparation for its subsequent transmembrane insertion.

Thermodynamic Cycle Without Turning Off Self-Interactions: Formal Discussion and a Numerical Example

The efficiency and accuracy of thermodynamic cycle calculations are considered. It is rigorously shown that the energy of the mutated part (MP) need not be scaled in a thermodynamic cycle computed with dual topology. Hence, there is no need to scale to zero any of the self-interactions (i.e. the interactions involving only particles of the same MP) regardless of whether the MP is bound or not to the main system. This observation carries a promise to lower computational resources and increase accuracy. A numerical test of a complete thermodynamic cycle illustrates cost and accuracy considerations.

Absolute free energy of binding of avidin/biotin, revisited

The binding of biotin to avidin is one of the strongest in nature with absolute free energy of binding, ΔA(0) = -20.4 kcal/mol. Therefore, this complex became a target for a large number of computational studies, which all, however, are based on approximate techniques or simplified models and have led to a wide range of results Therefore, ΔA(0) is calculated here by rigorous statistical mechanical methods and models that consider long-range electrostatics. (1) We apply our method, "hypothetical scanning molecular dynamics with thermodynamic integration" (HSMD-TI) to avidin-biotin modeled by periodic boundary conditions with particle mesh ewald (PME). (2) We apply the double decoupling method (DDM) to this system modeled by the spherical solvent boundary potential (SSBP) and the generalized solvent boundary potential (GSBP). The corresponding results for neutral biotin, ΔA(0) = -29.1 ± 0.8 and -25.2 ± 0.5 kcal/mol are significantly lower than the experimental value; we also provide the result for a charged biotin, ΔA(0) = -33.3 ± 0.8 kcal/mol. It is plausible to suggest that this disagreement with the experiment may stem from ignoring the (positive) contribution of a mobile loop that changes its structure upon ligand binding.

Fungal bis-Naphthopyrones as Inhibitors of Botulinum Neurotoxin Serotype A

An in silico screen of the NIH Molecular Library Small Molecule Repository (MLSMR) of ∼350000 compounds and confirmatory bioassays led to identification of chaetochromin A (1) as an inhibitor of botulinum neurotoxin serotype A (BoNT A). Subsequent acquisition and testing of analogues of 1 uncovered two compounds, talaroderxines A (2) and B (3), with improved activity. These are the first fungal metabolites reported to exhibit BoNT/A inhibitory activity.

Comparison of free energy surfaces calculations from ab initio molecular dynamic simulations at the example of two transition metal catalyzed reactions

We carried out ab initio molecular dynamic simulations in order to determine the free energy surfaces of two selected reactions including solvents, namely a rearrangement of a ruthenium oxoester in water and a carbon dioxide addition to a palladium complex in carbon dioxide. For the latter reaction we also investigated the gas phase reaction in order to take solvent effects into account. We used two techniques to reconstruct the free energy surfaces: thermodynamic integration and metadynamics. Furthermore, we gave a reasonable error estimation of the computed free energy surface. We calculated a reaction barrier of ΔF = 59.5 ± 8.5 kJ mol(-1) for the rearrangement of a ruthenium oxoester in water from thermodynamic integration. For the carbon dioxide addition to the palladium complex in carbon dioxide we found a ΔF = 44.9 ± 3.3 kJ mol(-1) from metadynamics simulations with one collective variable. The investigation of the same reactions in the gas phase resulted in ΔF = 24.9 ± 6.7 kJ mol(-1) from thermodynamic integration, in ΔF = 26.7 ± 2.3 kJ mol(-1) from metadynamics simulations with one collective variable, and in ΔF = 27.1 ± 5.9 kJ mol(-1) from metadynamics simulations with two collective variables.

Molecular dynamics guided study of salt bridge length dependence in both fluorinated and non-fluorinated parallel dimeric coiled-coils

The alpha-helical coiled-coil is one of the most common oligomerization motifs found in both native and engineered proteins. To better understand the stability and dynamics of the coiled-coil motifs, including those modified by fluorination, several fluorinated and nonfluorinated parallel dimeric coiled-coil protein structures were designed and modeled. We also attempt to investigate how changing the length and geometry of the important stabilizing salt bridges influences the coiled-coil protein structure. Molecular dynamics (MD) and free energy simulations with AMBER used a particle mesh Ewald treatment of the electrostatics in explicit TIP3P solvent with balanced force field treatments. Preliminary studies with legacy force fields (ff94, ff96, and ff99) show a profound instability of the coiled-coil structures in short MD simulation. Significantly, better behavior is evident with the more balanced ff99SB and ff03 protein force fields. Overall, the results suggest that the coiled-coil structures can readily accommodate the larger acidic arginine or S-2,7-diaminoheptanedoic acid mutants in the salt bridge, whereas substitution of the smaller L-ornithine residue leads to rapid disruption of the coiled-coil structure on the MD simulation time scale. This structural distortion of the secondary structure allows both the formation of large hydration pockets proximal to the charged groups and within the hydrophobic core. Moreover, the increased structural fluctuations and movement lead to a decrease in the water occupancy lifetimes in the hydration pockets. In contrast, analysis of the hydration in the stable dimeric coiled-coils shows high occupancy water sites along the backbone residues with no water occupancy in the hydrophobic core, although transitory water interactions with the salt bridge residues are evident. The simulations of the fluorinated coiled-coils suggest that in some cases fluorination electrostatically stabilizes the intermolecular coiled-coil salt bridges. Structural analyses also reveal different side chain rotamer preferences for leucine when compared with 5,5,5,5',5',5'-hexafluoroleucine mutants. These observed differences in the side chain rotamer populations suggest differential changes in the side chain conformational entropy upon coiled-coil formation when the protein is fluorinated. The free energy of hydration of the isolated 5,5,5,5',5',5'-hexafluoroleucine amino acid is calculated to be 1.1 kcal/mol less stable than leucine; this hydrophobic penalty in the monomer may provide a driving force for coiled-coil dimer formation. Estimation of the ellipticity at 222 nm from a series of snapshots from the MD simulations with DicroCalc shows distinct increases in the ellipticity when the coiled-coil is fluorinated, which suggests that the helicity in the folded coiled-coils is greater when fluorinated.

Involvement of water in carbohydrate-protein binding: concanavalin A revisited

Ordered water molecules bound to protein surfaces, or in protein-ligand interfaces, are frequently observed by crystallography. The investigation of the impact of such conserved water molecules on protein stability and ligand affinity requires detailed structural, dynamic, and thermodynamic analyses. Several crystal structures of the legume lectin concanavalin A (Con A) bound to closely related carbohydrate ligands show the presence of a conserved water molecule that mediates ligand binding. Experimental thermodynamic and theoretical studies have examined the role of this conserved water in the complexation of Con A with a synthetic analog of the natural trisaccharide, in which a hydroxyethyl side chain replaces the hydroxyl group at the C-2 position in the central mannosyl residue. Molecular modeling earlier indicated (Clarke, C.; Woods, R. J.; Glushka, J.; Cooper, A.; Nutley, M. A.; Boons, G.-J. J. Am. Chem. Soc. 2001, 123, 12238-12247) that the hydroxyl group in this synthetic side chain could occupy a position equivalent to that of the conserved water, and thus might displace it. An interpretation of the experimental thermodynamic data, which was consistent with the displacement of the conserved water, was also presented. The current work reports the crystal structure of Con A with this synthetic ligand and shows that even though the position and interactions of the conserved water are distorted, this key water is not displaced by the hydroxyethyl moiety. This new structural data provides a firm basis for molecular dynamics simulations and thermodynamic integration calculations whose results indicate that differences in van der Waals contacts (insertion energy), rather than electrostatic interactions (charging energy) are fundamentally responsible for the lower affinity of the synthetic ligand. When combined with the new crystallographic data, this study provides a straightforward interpretation for the lower affinity of the synthetic analog; specifically, that it arises primarily from weaker interactions with the protein via the positionally perturbed conserved water. This interpretation is fully consistent with the experimental observations that the free energy of binding is enthalpy driven, that there is both less enthalpic gain and less entropic penalty for binding the synthetic ligand, relative to the natural trisaccharide, and that the entropic component does not arise from releasing an ordered water molecule from the protein surface to the bulk solvent.

Secondary structure bias in generalized Born solvent models: comparison of conformational ensembles and free energy of solvent polarization from explicit and implicit solvation

The effects of the use of three generalized Born (GB) implicit solvent models on the thermodynamics of a simple polyalanine peptide are studied via comparing several hundred nanoseconds of well-converged replica exchange molecular dynamics (REMD) simulations using explicit TIP3P solvent to REMD simulations with the GB solvent models. It is found that when compared to REMD simulations using TIP3P the GB REMD simulations contain significant differences in secondary structure populations, most notably an overabundance of alpha-helical secondary structure. This discrepancy is explored via comparison of the differences in the electrostatic component of the free energy of solvation (DeltaDeltaG(pol)) between TIP3P (via thermodynamic Integration calculations), the GB models, and an implicit solvent model based on the Poisson equation (PE). The electrostatic components of the solvation free energies are calculated using each solvent model for four representative conformations of Ala10. Since the PE model is found to have the best performance with respect to reproducing TIP3P DeltaDeltaG(pol) values, effective Born radii from the GB models are compared to effective Born radii calculated with PE (so-called perfect radii), and significant and numerous deviations in GB radii from perfect radii are found in all GB models. The effect of these deviations on the solvation free energy is discussed, and it is shown that even when perfect radii are used the agreement of GB with TIP3P DeltaDeltaG(pol) values does not improve. This suggests a limit to the optimization of the effective Born radius calculation and that future efforts to improve the accuracy of GB models must extend beyond such optimizations.