Modulation of human transthyretin stability by the mutations at histidine 88 studied by free energy simulation.

Human transthyretin (TTR) is a homotetrameric plasma protein associated with a high percentage of ß-sheet, which forms amyloid fibrils and accumulates in tissues or extracellular matrix to cause amyloid diseases. Free energy simulations based on all-atom molecular dynamics simulations were carried out to analyze the effects of the His88 → Arg, Phe, and Tyr mutations on the stability of human TTR. The calculated free energy change differences (ΔΔG) caused by the His → Arg, Phe, and Tyr mutations at position 88 are 6.48 ± 0.45, -9.99 ± 0.54, and 2.66 ± 0.33 kcal/mol, respectively. These calculated free energy change differences between wild type and the mutants are in excellent agreement with prior experimental values. Our simulation results show that the wild type of the TTR is more stable than H88R and H88Y mutants, whereas it is less stable than the H88F mutant. The free energy component analysis shows that the primary contribution to the free energy change difference (ΔΔG) for the His → Arg mutation arises from electrostatic interaction; the ΔΔG for the His → Phe mutation is from van der Waals and electrostatic interactions and that for the His → Tyr mutation from covalent interaction. The simulation results show that the free energy calculation with thermodynamic integration is beneficial for understanding the detailed microscopic mechanism of protein stability. The implications of the results for understanding stabilizing and destabilizing effect of the mutation and the contribution to protein stability are discussed.

Prediction of protein-protein complexes using replica exchange with repulsive scaling.

The realistic prediction of protein-protein complex structures is import to ultimately model the interaction of all proteins in a cell and for the design of new protein-protein interactions. In principle, molecular dynamics (MD) simulations allow one to follow the association process under realistic conditions including full partner flexibility and surrounding solvent. However, due to the many local binding energy minima at the surface of protein partners, MD simulations are frequently trapped for long times in transient association states. We have designed a replica-exchange based scheme employing different levels of a repulsive biasing between partners in each replica simulation. The bias acts only on intermolecular interactions based on an increase in effective pairwise van der Waals radii (repulsive scaling (RS)-REMD) without affecting interactions within each protein or with the solvent. For a set of five protein test cases (out of six) the RS-REMD technique allowed the sampling of near-native complex structures even when starting from the opposide site with respect to the native binding site for one partner. Using the same start structures and same computational demand regular MD simulations sampled near native complex structures only for one case. The method showed also improved results for the refinement of docked structures in the vicinity of the native binding geometry compared to regular MD refinement.

BAR-based optimum adaptive steered MD for configurational sampling.

The equilibrium and nonequilibrium adaptive alchemical free energy simulation methods optimum Bennett's acceptance ratio and optimum crooks' equation (OCE), based on the statistically optimal bidirectional reweighting estimator named Bennett's Acceptance Ratio or Crooks' equation, perform initial sampling in the staging alchemical transformation and then determine the importance rank of different states via the time-derivative of the variance. The method is proven to give speedups compared with the equal time rule. In the current work, we extend the time derivative of variance guided adaptive sampling method to the configurational space, falling in the term of steered MD (SMD). The SMD approach biasing physically meaningful collective variable (CV) such as one dihedral or one distance to pulling the system from one conformational state to another. By minimizing the variance of the free energy differences along the pathway in an optimized way, a new type of adaptive SMD (ASMD) is introduced. As exhibits in the alchemical case, this adaptive sampling method outperforms the traditional equal-time SMD in nonequilibrium stratification. Also, the method gives much more efficient calculation of potential of mean force than the selection criterion-based ASMD scheme, which is proven to be more efficient than traditional SMD. The OCE workflow is periodicity-of-CV dependent while ASMD is not. The performance is demonstrated in a dihedral flipping case and two distance pulling cases, accounting for periodic and nonperiodic CVs, respectively. © 2019 Wiley Periodicals, Inc.

New QM/MM implementation of the DFTB3 method in the gromacs package.

The approximate density-functional tight-binding theory method DFTB3 has been implemented in the quantum mechanics/molecular mechanics (QM/MM) framework of the Gromacs molecular simulation package. We show that the efficient smooth particle-mesh Ewald implementation of Gromacs extends to the calculation of QM/MM electrostatic interactions. Further, we make use of the various free-energy functionalities provided by Gromacs and the PLUMED plugin. We exploit the versatility and performance of the current framework in three typical applications of QM/MM methods to solve biophysical problems: (i) ultrafast proton transfer in malonaldehyde, (ii) conformation of the alanine dipeptide, and (iii) electron-induced repair of a DNA lesion. Also discussed is the further development of the framework, regarding mostly the options for parallelization.

Monte Carlo simulations of proteins at constant pH with generalized Born solvent, flexible sidechains, and an effective dielectric boundary.

Titratable residues determine the acid/base behavior of proteins, strongly influencing their function; in addition, proton binding is a valuable reporter on electrostatic interactions. We describe a method for pK(a) calculations, using constant-pH Monte Carlo (MC) simulations to explore the space of sidechain conformations and protonation states, with an efficient and accurate generalized Born model (GB) for the solvent effects. To overcome the many-body dependency of the GB model, we use a "Native Environment" approximation, whose accuracy is shown to be good. It allows the precalculation and storage of interactions between all sidechain pairs, a strategy borrowed from computational protein design, which makes the MC simulations themselves very fast. The method is tested for 12 proteins and 167 titratable sidechains. It gives an rms error of 1.1 pH units, similar to the trivial "Null" model. The only adjustable parameter is the protein dielectric constant. The best accuracy is achieved for values between 4 and 8, a range that is physically plausible for a protein interior. For sidechains with large pKa shifts, ≥2, the rms error is 1.6, compared to 2.5 with the Null model and 1.5 with the empirical PROPKA method.

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.

Effect of alanine versus serine at position 88 of human transthyretin mutants on the protein stability.

Human transthyretin (TTR) is a homo-tetrameric plasma protein associated with a high percentage of ß-sheet forming amyloid fibrils. It accumulates in tissues or extracellular matrices to cause amyloid diseases. Free energy simulations with thermodynamic integration based on all-atom molecular dynamics simulations have been carried out to analyze the effects of the His88 â†’ Ala and Ser mutations on the stability of human TTR. The calculated free energy change differences (ΔΔG) caused by the His88 â†’ Ala and His88 â†’ Ser mutations are -1.84 ± 0.86 and 7.56 ± 0.55 kcal/mol, respectively, which are in excellent agreement with prior reported experimental values. The simulation results show that the H88A mutant is more stable than the wild type, whereas the H88S mutant is less stable than the wild type. The free energy component analysis shows that the contribution to the free energy change difference (ΔΔG) for the His88 â†’ Ala and His88 â†’ Ser mutations mainly arise from electrostatic and van der Waals interactions, respectively. The electrostatic term stabilizes the H88A mutant more than the wild type, but the van der Waals interaction destabilizes the H88S mutant relative to the wild type. Individual residue contributions to the free energy change show neighboring residues exert stabilizing and destabilizing influence on the mutants. The implications of the simulation results for understanding the stabilizing and destabilizing effect and its contribution to protein stability are discussed.

Collective Variable for Metadynamics Derived From AlphaFold Output.

AlphaFold is a neural network-based tool for the prediction of 3D structures of proteins. In CASP14, a blind structure prediction challenge, it performed significantly better than other competitors, making it the best available structure prediction tool. One of the outputs of AlphaFold is the probability profile of residue-residue distances. This makes it possible to score any conformation of the studied protein to express its compliance with the AlphaFold model. Here, we show how this score can be used to drive protein folding simulation by metadynamics and parallel tempering metadynamics. Using parallel tempering metadynamics, we simulated the folding of a mini-protein Trp-cage and ß hairpin and predicted their folding equilibria. We observe the potential of the AlphaFold-based collective variable in applications beyond structure prediction, such as in structure refinement or prediction of the outcome of a mutation.

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.

Free energy simulations to understand the effect of Met → Ala mutations at positions 205, 206 and 213 on stability of human prion protein.

Prion diseases are a family of infectious amyloid diseases affecting human and animals. Prion propagation in transmissible spongiform encephalopathies is associated with the unfolding and conversion of normal cellular prion protein into its pathogenic scrapie form. Understanding the fundamentals of prion protein aggregation caused by mutations is crucial to unravel the pathology of prion diseases. To help understand the contributions of individual residues to the stability of the human prion protein, we have carried out free energy simulations based on atomistic molecular dynamics trajectories. We focus on Met → Ala mutations at positions 205, 206 and 213, which are mostly buried residues located on helix 3 of the protein. The simulations predicted that all three mutations destabilize the prion protein. Changes in unfolding free energy upon mutation, ∆∆G, are 3.10 ± 0.79, 2.00 ± 0.26 and 3.06 ± 0.66 kcal/mol for M205A, M206A and M213A, respectively, in excellent agreement with the corresponding experimental values of 3.09 ± 0.28, 1.50 ± 0.34 and 3.12 ± 0.27 kcal/mol [T. Hart et al. (2009) PNAS 106, 5651-5656]. Component analysis indicates that the major contributions to the loss of protein stability arise from van der Waals interactions for the M205A and M206A mutations, and from van der Waals and covalent energy terms for M213A. Interestingly, while free energy contributions from a majority of residues neighboring the mutation sites tend to stabilize the wild type, there are a few residues stabilizing the mutant side chains. Our results show that this approach to free energy calculation can be very useful for understanding the detailed mechanism of human prion protein stability.

Molecular Dynamics Free Energy Simulations Reveal the Mechanism for the Antiviral Resistance of the M66I HIV-1 Capsid Mutation.

While drug resistance mutations can often be attributed to the loss of direct or solvent-mediated protein-ligand interactions in the drug-mutant complex, in this study we show that a resistance mutation for the picomolar HIV-1 capsid (CA)-targeting antiviral (GS-6207) is mainly due to the free energy cost of the drug-induced protein side chain reorganization in the mutant protein. Among several mutations, M66I causes the most suppression of the GS-6207 antiviral activity (up to ~84,000-fold), and only 83- and 68-fold reductions for PF74 and ZW-1261, respectively. To understand the molecular basis of this drug resistance, we conducted molecular dynamics free energy simulations to study the structures, energetics, and conformational free energy landscapes involved in the inhibitors binding at the interface of two CA monomers. To minimize the protein-ligand steric clash, the I66 side chain in the M66I-GS-6207 complex switches to a higher free energy conformation from the one adopted in the apo M66I. In contrast, the binding of GS-6207 to the wild-type CA does not lead to any significant M66 conformational change. Based on an analysis that decomposes the absolute binding free energy into contributions from two receptor conformational states, it appears that it is the free energy cost of side chain reorganization rather than the reduced protein-ligand interaction that is largely responsible for the drug resistance against GS-6207.

Mechanism of recognition of parallel G-quadruplexes by DEAH/RHAU helicase DHX36 explored by molecular dynamics simulations.

Because of high stability and slow unfolding rates of G-quadruplexes (G4), cells have evolved specialized helicases that disrupt these non-canonical DNA and RNA structures in an ATP-dependent manner. One example is DHX36, a DEAH-box helicase, which participates in gene expression and replication by recognizing and unwinding parallel G4s. Here, we studied the molecular basis for the high affinity and specificity of DHX36 for parallel-type G4s using all-atom molecular dynamics simulations. By computing binding free energies, we found that the two main G4-interacting subdomains of DHX36, DSM and OB, separately exhibit high G4 affinity but they act cooperatively to recognize two distinctive features of parallel G4s: the exposed planar face of a guanine tetrad and the unique backbone conformation of a continuous guanine tract, respectively. Our results also show that DSM-mediated interactions are the main contributor to the binding free energy and rely on making extensive van der Waals contacts between the GXXXG motifs and hydrophobic residues of DSM and a flat guanine plane. Accordingly, the sterically more accessible 5'-G-tetrad allows for more favorable van der Waals and hydrophobic interactions which leads to the preferential binding of DSM to the 5'-side. In contrast to DSM, OB binds to G4 mostly through polar interactions by flexibly adapting to the 5'-terminal guanine tract to form a number of strong hydrogen bonds with the backbone phosphate groups. We also identified a third DHX36/G4 interaction site formed by the flexible loop missing in the crystal structure.

Recent Developments in Free Energy Calculations for Drug Discovery.

The grand challenge in structure-based drug design is achieving accurate prediction of binding free energies. Molecular dynamics (MD) simulations enable modeling of conformational changes critical to the binding process, leading to calculation of thermodynamic quantities involved in estimation of binding affinities. With recent advancements in computing capability and predictive accuracy, MD based virtual screening has progressed from the domain of theoretical attempts to real application in drug development. Approaches including the Molecular Mechanics Poisson Boltzmann Surface Area (MM-PBSA), Linear Interaction Energy (LIE), and alchemical methods have been broadly applied to model molecular recognition for drug discovery and lead optimization. Here we review the varied methodology of these approaches, developments enhancing simulation efficiency and reliability, remaining challenges hindering predictive performance, and applications to problems in the fields of medicine and biochemistry.

Computational design of small molecular modulators of protein-protein interactions with a novel thermodynamic cycle: Allosteric inhibitors of HIV-1 integrase.

Targeting protein-protein interactions for therapeutic development involves designing small molecules to either disrupt or enhance a known PPI. For this purpose, it is necessary to compute reliably the effect of chemical modifications of small molecules on the protein-protein association free energy. Here we present results obtained using a novel thermodynamic free energy cycle, for the rational design of allosteric inhibitors of HIV-1 integrase (ALLINI) that act specifically in the early stage of the infection cycle. The new compounds can serve as molecular probes to dissect the multifunctional mechanisms of ALLINIs to inform the discovery of new allosteric inhibitors. The free energy protocol developed here can be more broadly applied to study quantitatively the effects of small molecules on modulating the strengths of protein-protein interactions.

In vitro studies on non-canonical DNA binding specificities of KAP6 and HMO1 and mechanistic insights into DNA bound and unbinding dynamics of KAP6.

High mobility group box (HMGB) members are DNA binding proteins with varied functions present across kingdoms. The mechanism by which HMGBs with varying number of HMG boxes are able to carry out similar functions, are poorly understood. Moreover, how non-canonical DNAs are recognized by HMGB proteins is not clear. To address these, we carried out detailed biochemical and computational studies to characterize two HMGB members- Kinetoplast associated protein (KAP6) of Trypanosoma and High mobility group protein 1 (HMO1) from yeast. Here, we report that KAP6 binds non-canonical DNAs tighter than B-form DNA. Among non-canonical DNAs, KAP6 has the highest affinity for splayed and flap structures, but least for Holliday Junction (HJ). In contrast, HMO1 binds tighter to HJ. Computational analysis show that the secondary structural elements involved in DNA interaction are conserved in HMGB members KAP6 and mitochondrial transcription factor A. Simulation analyses revealed that the ~90° bend in DNA induced by KAP6 HMG box is a result of two ~45° bends, by Helix 1 and Helix 2 of the protein. Our data also suggests that the orthologs of HMO1 and KAP6 are oligomers in solution, which could be necessary for their functioning such as DNA bending and looping.

Analysis of Phosphoryl-Transfer Enzymes with QM/MM Free Energy Simulations.

We discuss the application of quantum mechanics/molecular mechanics (QM/MM) free energy simulations to the analysis of phosphoryl transfers catalyzed by two enzymes: alkaline phosphatase and myosin. We focus on the nature of the transition state and the issue of mechanochemical coupling, respectively, in the two enzymes. The results illustrate unique insights that emerged from the QM/MM simulations, especially concerning the interpretation of experimental data regarding the nature of enzymatic transition states and coupling between global structural transition and catalysis in the active site. We also highlight a number of technical issues worthy of attention when applying QM/MM free energy simulations, and comment on a number of technical and mechanistic issues that require further studies.

Recent Developments and Applications of the MMPBSA Method.

The Molecular Mechanics Poisson-Boltzmann Surface Area (MMPBSA) approach has been widely applied as an efficient and reliable free energy simulation method to model molecular recognition, such as for protein-ligand binding interactions. In this review, we focus on recent developments and applications of the MMPBSA method. The methodology review covers solvation terms, the entropy term, extensions to membrane proteins and high-speed screening, and new automation toolkits. Recent applications in various important biomedical and chemical fields are also reviewed. We conclude with a few future directions aimed at making MMPBSA a more robust and efficient method.

Toward Determining ATPase Mechanism in ABC Transporters: Development of the Reaction Path-Force Matching QM/MM Method.

Adenosine triphosphate (ATP)-binding cassette (ABC) transporters are ubiquitous ATP-dependent membrane proteins involved in translocations of a wide variety of substrates across cellular membranes. To understand the chemomechanical coupling mechanism as well as functional asymmetry in these systems, a quantitative description of how ABC transporters hydrolyze ATP is needed. Complementary to experimental approaches, computer simulations based on combined quantum mechanical and molecular mechanical (QM/MM) potentials have provided new insights into the catalytic mechanism in ABC transporters. Quantitatively reliable determination of the free energy requirement for enzymatic ATP hydrolysis, however, requires substantial statistical sampling on QM/MM potential. A case study shows that brute force sampling of ab initio QM/MM (AI/MM) potential energy surfaces is computationally impractical for enzyme simulations of ABC transporters. On the other hand, existing semiempirical QM/MM (SE/MM) methods, although affordable for free energy sampling, are unreliable for studying ATP hydrolysis. To close this gap, a multiscale QM/MM approach named reaction path-force matching (RP-FM) has been developed. In RP-FM, specific reaction parameters for a selected SE method are optimized against AI reference data along reaction paths by employing the force matching technique. The feasibility of the method is demonstrated for a proton transfer reaction in the gas phase and in solution. The RP-FM method may offer a general tool for simulating complex enzyme systems such as ABC transporters.