Improved convergence in dual-topology free energy calculations through use of harmonic restraints

GENESIS 2.1: High-Performance Molecular Dynamics Software for Enhanced Sampling and Free-Energy Calculations for Atomistic, Coarse-Grained, and Quantum Mechanics/Molecular Mechanics Models

GENeralized-Ensemble SImulation System (GENESIS) is a molecular dynamics (MD) software developed to simulate the conformational dynamics of a single biomolecule, as well as molecular interactions in large biomolecular assemblies and between multiple biomolecules in cellular environments. To achieve the latter purpose, the earlier versions of GENESIS emphasized high performance in atomistic MD simulations on massively parallel supercomputers, with or without graphics processing units (GPUs). Here, we implemented multiscale MD simulations that include atomistic, coarse-grained, and hybrid quantum mechanics/molecular mechanics (QM/MM) calculations. They demonstrate high performance and are integrated with enhanced conformational sampling algorithms and free-energy calculations without using external programs except for the QM programs. In this article, we review new functions, molecular models, and other essential features in GENESIS version 2.1 and discuss ongoing developments for future releases.

A computationally designed ACE2 decoy has broad efficacy against SARS-CoV-2 omicron variants and related viruses in vitro and in vivo

SARS-CoV-2, especially B.1.1.529/omicron and its sublineages, continues to mutate to evade monoclonal antibodies and antibodies elicited by vaccination. Affinity-enhanced soluble ACE2 (sACE2) is an alternative strategy that works by binding the SARS-CoV-2 S protein, acting as a 'decoy' to block the interaction between the S and human ACE2. Using a computational design strategy, we designed an affinity-enhanced ACE2 decoy, FLIF, that exhibited tight binding to SARS-CoV-2 delta and omicron variants. Our computationally calculated absolute binding free energies (ABFE) between sACE2:SARS-CoV-2 S proteins and their variants showed excellent agreement to binding experiments. FLIF displayed robust therapeutic utility against a broad range of SARS-CoV-2 variants and sarbecoviruses, and neutralized omicron BA.5 in vitro and in vivo. Furthermore, we directly compared the in vivo therapeutic efficacy of wild-type ACE2 (non-affinity enhanced ACE2) against FLIF. A few wild-type sACE2 decoys have shown to be effective against early circulating variants such as Wuhan in vivo. Our data suggest that moving forward, affinity-enhanced ACE2 decoys like FLIF may be required to combat evolving SARS-CoV-2 variants. The approach described herein emphasizes how computational methods have become sufficiently accurate for the design of therapeutics against viral protein targets. Affinity-enhanced ACE2 decoys remain highly effective at neutralizing omicron subvariants.

Enhancing Paraoxon Binding to Organophosphorus Hydrolase Active Site

Organophosphorus hydrolase (OPH) is a metalloenzyme that can hydrolyze organophosphorus agents resulting in products that are generally of reduced toxicity. The best OPH substrate found to date is diethyl p-nitrophenyl phosphate (paraoxon). Most structural and kinetic studies assume that the binding orientation of paraoxon is identical to that of diethyl 4-methylbenzylphosphonate, which is the only substrate analog co-crystallized with OPH. In the current work, we used a combined docking and molecular dynamics (MD) approach to predict the likely binding mode of paraoxon. Then, we used the predicted binding mode to run MD simulations on the wild type (WT) OPH complexed with paraoxon, and OPH mutants complexed with paraoxon. Additionally, we identified three hot-spot residues (D253, H254, and I255) involved in the stability of the OPH active site. We then experimentally assayed single and double mutants involving these residues for paraoxon binding affinity. The binding free energy calculations and the experimental kinetics of the reactions between each OPH mutant and paraoxon show that mutated forms D253E, D253E-H254R, and D253E-I255G exhibit enhanced substrate binding affinity over WT OPH. Interestingly, our experimental results show that the substrate binding affinity of the double mutant D253E-H254R increased by 19-fold compared to WT OPH.

Computationally Designed ACE2 Decoy Receptor Binds SARS-CoV-2 Spike (S) Protein with Tight Nanomolar Affinity

Even with the availability of vaccines, therapeutic options for COVID-19 still remain highly desirable, especially in hospitalized patients with moderate or severe disease. Soluble ACE2 (sACE2) is a promising therapeutic candidate that neutralizes SARS CoV-2 infection by acting as a decoy. Using computational mutagenesis, we designed a number of sACE2 derivatives carrying three to four mutations. The top-predicted sACE2 decoy based on the in silico mutagenesis scan was subjected to molecular dynamics and free-energy calculations for further validation. After illuminating the mechanism of increased binding for our designed sACE2 derivative, the design was verified experimentally by flow cytometry and BLI-binding experiments. The computationally designed sACE2 decoy (ACE2-FFWF) bound the receptor-binding domain of SARS-CoV-2 tightly with low nanomolar affinity and ninefold affinity enhancement over the wild type. Furthermore, cell surface expression was slightly greater than wild-type ACE2, suggesting that the design is well-folded and stable. Having an arsenal of high-affinity sACE2 derivatives will help to buffer against the emergence of SARS CoV-2 variants. Here, we show that computational methods have become sufficiently accurate for the design of therapeutics for current and future viral pandemics.

An in silico approach for identification of novel inhibitors as potential therapeutics targeting COVID-19 main protease

Respiratory disease caused by a novel coronavirus, COVID-19, has been labeled a pandemic by the World Health Organization. Very little is known about the infection mechanism for this virus. More importantly, there are no drugs or vaccines that can cure or prevent a person from getting COVID-19. In this study, the binding affinity of 2692 protease inhibitor compounds that are known in the protein data bank, are calculated against the main protease of the novel coronavirus with docking and molecular dynamics (MD). Both the docking and MD methods predict the macrocyclic tissue factor-factor VIIa (PubChem ID: 118098670) inhibitor to bind strongly with the main protease with a binding affinity of -10.6 and -10.0 kcal/mol, respectively. The TF-FVIIa inhibitors are known to prevent the coagulation of blood and have antiviral activity as shown in the case of SARS coronavirus. Two more inhibitors, phenyltriazolinones (PubChem ID: 104161460) and allosteric HCV NS5B polymerase thumb pocket 2 (PubChem ID: 163632044) have shown antiviral activity and also have high affinity towards the main protease of COVID-19. Furthermore, these inhibitors interact with the catalytic dyad in the active site of the COVID-19 main protease that is especially important in viral replication. The calculated theoretical dissociation constants of the proposed COVID-19 inhibitors are found to be very similar to the experimental dissociation constant values of similar protease-inhibitor systems.Communicated by Ramaswamy H. Sarma.

In-silico screening and microsecond molecular dynamics simulations to identify single point mutations that destabilize β-hexosaminidase A causing Tay-Sachs disease

β-hexosaminidase A (HexA) protein is responsible for the degradation of GM2 gangliosides in the central and peripheral nervous systems. Tay-Sachs disease occurs when HexA within Hexosaminidase does not properly function and harmful GM2 gangliosides begin to build up within the neurons. In this study, in silico methods such as SIFT, PolyPhen-2, PhD-SNP, and MutPred were utilized to analyze the effects of nonsynonymous single nucleotide polymorphisms (nsSNPs) on HexA in order to identify possible pathogenetic and deleterious variants. Molecular dynamics (MD) simulations showed that two mutants, P25S and W485R, experienced an increase in structural flexibility compared to the native protein. Particularly, there was a decrease in the overall number and frequencies of hydrogen bonds for the mutants compared to the wildtype. MM/GBSA calculations were performed to help assess the change in binding affinity between the wildtype and mutant structures and a mechanism-based inhibitor, NGT, which is known to help increase the residual activity of HexA. Both of the mutants experienced a decrease in the binding affinity from -23.8 kcal/mol in wildtype to -20.9 and -18.7 kcal/mol for the P25S and W485R variants of HexA, respectively.

Dummy Atoms in Alchemical Free Energy Calculations

In calculations of relative free energy differences, the number of atoms of the initial and final states is rarely the same. This necessitates the introduction of dummy atoms. These placeholders interact with the physical system only by bonded energy terms. We investigate the conditions necessary so that the presence of dummy atoms does not influence the result of a relative free energy calculation. On the one hand, one has to ensure that dummy atoms only give a multiplicative contribution to the partition function so that their contribution cancels from double-free energy differences. On the other hand, the bonded terms used to attach a dummy atom (or group of dummy atoms) to the physical system have to maintain it in a well-defined position and orientation relative to the physical system. A detailed theoretical analysis of both aspects is provided, illustrated by 24 calculations of relative solvation free energy differences, for which all four legs of the underlying thermodynamic cycle were computed. Cycle closure (or lack thereof) was used as a sensitive indicator to probing the effects of dummy atom treatment on the resulting free energy differences. We find that a naive (but often practiced) treatment of dummy atoms results in errors of up to kBT when calculating the relative solvation free energy difference between two small solutes, such as methane and ammonia. While our analysis focuses on the so-called single topology approach to set up alchemical transformations, similar considerations apply to dual topology, at least many widely used variants thereof.

Prediction and evaluation of deleterious and disease causing non-synonymous SNPs (nsSNPs) in human NF2 gene responsible for neurofibromatosis type 2 (NF2)

The majority of genetic variations in the human genome that lead to variety of different diseases are caused by non-synonymous single nucleotide polymorphisms (nsSNPs). Neurofibromatosis type 2 (NF2) is a deadly disease caused by nsSNPs in the NF2 gene that encodes for a protein called merlin. This study used various in silico methods, SIFT, Polyphen-2, PhD-SNP and MutPred, to investigate the pathogenic effect of 14 nsSNPs in the merlin FERM domain. The G197C and L234R mutations were found to be two deleterious and disease mutations associated with the mild and severe forms of NF2, respectively. Molecular dynamics (MD) simulations were conducted to understand the stability, structure and dynamics of these mutations. Both mutant structures experienced larger flexibility compared to the wildtype. The L234R mutant suffered from more prominent structural instability, which may help to explain why it is associated with the more severe form of NF2. The intramolecular hydrogen bonding in L234R mutation decreased from the wildtype, while intermolecular hydrogen bonding of L234R mutation with solvent greatly increased. The native contacts were also found to be important. Protein-protein docking revealed that L234R mutation decreased the binding complementarity and binding affinity of LATS2 to merlin, which may have an impact on merlin's ability to regulate the Hippo signaling pathway. The calculated binding affinity of the LATS2 to L234R mutant and wildtype merlin protein is found to be 21.73 and -11 kcal/mol, respectively. The binding affinity of the wildtype merlin agreed very well with the experimental value, -8 kcal/mol.Communicated by Ramaswamy H. Sarma.

Reproducibility of Free Energy Calculations across Different Molecular Simulation Software Packages

Alchemical free energy calculations are an increasingly important modern simulation technique to calculate free energy changes on binding or solvation. Contemporary molecular simulation software such as AMBER, CHARMM, GROMACS, and SOMD include support for the method. Implementation details vary among those codes, but users expect reliability and reproducibility, i.e., for a given molecular model and set of force field parameters, comparable free energy differences should be obtained within statistical bounds regardless of the code used. Relative alchemical free energy (RAFE) simulation is increasingly used to support molecule discovery projects, yet the reproducibility of the methodology has been less well tested than its absolute counterpart. Here we present RAFE calculations of hydration free energies for a set of small organic molecules and demonstrate that free energies can be reproduced to within about 0.2 kcal/mol with the aforementioned codes. Absolute alchemical free energy simulations have been carried out as a reference. Achieving this level of reproducibility requires considerable attention to detail and package-specific simulation protocols, and no universally applicable protocol emerges. The benchmarks and protocols reported here should be useful for the community to validate new and future versions of software for free energy calculations.

Interactions of the DNA Repair Enzyme Human Thymine DNA Glycosylase with Cognate and Noncognate DNA

Glycosylases specifically recognize and flip their target base out of the DNA helix into the enzyme's active site. Our simulations show that a partially flipped state, already present in free DNA carrying a T:G mispair, becomes the more probable state compared to the closed state after binding of thymine DNA glycosylase (TDG). Paired thymine (T:A) or methyl-cytosine (mC:G) does not exhibit a partially flipped state in free or complexed DNA. Important enzyme-DNA interactions exhibit significant strength in the intrahelical and extrahelical TDG-DNA complexes. The computed binding free energy differences suggest these interactions account for the stabilization of the partially flipped state, thereby driving the T:G mispair toward base flip. In the fully flipped state, the cognate base thymine is significantly better accommodated in the enzyme's active site than noncognate bases are, suggesting the hydrolysis step as the last of several stages at which base recognition can be achieved.

A generalized linear response framework for expanded ensemble and replica exchange simulations

Expanded ensemble simulation is a powerful technique for enhancing sampling over a range of thermodynamic parameters. However, although the premise is relatively simple, running successful simulations in practice still presents something of an ad hoc challenge. Three main difficulties exist: (1) the selection of the thermodynamic states, (2) the selection of the sampling weights, and (3) efficient sampling of the expanded parameter space. Here we consider these problems in the context of a pairwise linear response approach to the work fluctuation theorem. The approach offers comprehensive tactics for addressing the three difficulties and can be used as either an alternative or a complement to replica exchange simulations. Importantly, the results are trivially implemented for multi-dimensional parameter spaces and they recover results from the literature aimed at the special cases of simulated/parallel tempering and replica exchange umbrella sampling. Illustrative examples are shown using the NAMD simulation engine.

Constant-pH Molecular Dynamics Simulations for Large Biomolecular Systems

An increasingly important endeavor is to develop computational strategies that enable molecular dynamics (MD) simulations of biomolecular systems with spontaneous changes in protonation states under conditions of constant pH. The present work describes our efforts to implement the powerful constant-pH MD simulation method, based on a hybrid nonequilibrium MD/Monte Carlo (neMD/MC) technique within the highly scalable program NAMD. The constant-pH hybrid neMD/MC method has several appealing features; it samples the correct semigrand canonical ensemble rigorously, the computational cost increases linearly with the number of titratable sites, and it is applicable to explicit solvent simulations. The present implementation of the constant-pH hybrid neMD/MC in NAMD is designed to handle a wide range of biomolecular systems with no constraints on the choice of force field. Furthermore, the sampling efficiency can be adaptively improved on-the-fly by adjusting algorithmic parameters during the simulation. Illustrative examples emphasizing medium- and large-scale applications on next-generation supercomputing architectures are provided.

Efficiency in nonequilibrium molecular dynamics Monte Carlo simulations

Hybrid algorithms combining nonequilibrium molecular dynamics and Monte Carlo (neMD/MC) offer a powerful avenue for improving the sampling efficiency of computer simulations of complex systems. These neMD/MC algorithms are also increasingly finding use in applications where conventional approaches are impractical, such as constant-pH simulations with explicit solvent. However, selecting an optimal nonequilibrium protocol for maximum efficiency often represents a non-trivial challenge. This work evaluates the efficiency of a broad class of neMD/MC algorithms and protocols within the theoretical framework of linear response theory. The approximations are validated against constant pH-MD simulations and shown to provide accurate predictions of neMD/MC performance. An assessment of a large set of protocols confirms (both theoretically and empirically) that a linear work protocol gives the best neMD/MC performance. Finally, a well-defined criterion for optimizing the time parameters of the protocol is proposed and demonstrated with an adaptive algorithm that improves the performance on-the-fly with minimal cost.

Changes in active site histidine hydrogen bonding trigger cryptochrome activation

Cryptochrome (CRY) is the principal light sensor of the insect circadian clock. Photoreduction of the Drosophila CRY (dCRY) flavin cofactor to the anionic semiquinone (ASQ) restructures a C-terminal tail helix (CTT) that otherwise inhibits interactions with targets that include the clock protein Timeless (TIM). All-atom molecular dynamics (MD) simulations indicate that flavin reduction destabilizes the CTT, which undergoes large-scale conformational changes (the CTT release) on short (25 ns) timescales. The CTT release correlates with the conformation and protonation state of conserved His378, which resides between the CTT and the flavin cofactor. Poisson-Boltzmann calculations indicate that flavin reduction substantially increases the His378 pKa Consistent with coupling between ASQ formation and His378 protonation, dCRY displays reduced photoreduction rates with increasing pH; however, His378Asn/Arg variants show no such pH dependence. Replica-exchange MD simulations also support CTT release mediated by changes in His378 hydrogen bonding and verify other responsive regions of the protein previously identified by proteolytic sensitivity assays. His378 dCRY variants show varying abilities to light-activate TIM and undergo self-degradation in cellular assays. Surprisingly, His378Arg/Lys variants do not degrade in light despite maintaining reactivity toward TIM, thereby implicating different conformational responses in these two functions. Thus, the dCRY photosensory mechanism involves flavin photoreduction coupled to protonation of His378, whose perturbed hydrogen-bonding pattern alters the CTT and surrounding regions.

Redox potential tuning through differential quinone binding in the photosynthetic reaction center of Rhodobacter sphaeroides

Ubiquinone forms an integral part of the electron transport chain in cellular respiration and photosynthesis across a vast number of organisms. Prior experimental results have shown that the photosynthetic reaction center (RC) from Rhodobacter sphaeroides is only fully functional with a limited set of methoxy-bearing quinones, suggesting that specific interactions with this substituent are required to drive electron transport and the formation of quinol. The nature of these interactions has yet to be determined. Through parameterization of a CHARMM-compatible quinone force field and subsequent molecular dynamics simulations of the quinone-bound RC, we have investigated and characterized the interactions of the protein with the quinones in the Q(A) and Q(B) sites using both equilibrium simulation and thermodynamic integration. In particular, we identify a specific interaction between the 2-methoxy group of ubiquinone in the Q(B) site and the amide nitrogen of GlyL225 that we implicate in locking the orientation of the 2-methoxy group, thereby tuning the redox potential difference between the quinones occupying the Q(A) and Q(B) sites. Disruption of this interaction leads to weaker binding in a ubiquinone analogue that lacks a 2-methoxy group, a finding supported by reverse electron transfer electron paramagnetic resonance experiments of the Q(A)⁻Q(B)⁻ biradical and competitive binding assays.

Computations of standard binding free energies with molecular dynamics simulations

An increasing number of studies have reported computations of the standard (absolute) binding free energy of small ligands to proteins using molecular dynamics (MD) simulations and explicit solvent molecules that are in good agreement with experiments. This encouraging progress suggests that physics-based approaches hold the promise of making important contributions to the process of drug discovery and optimization in the near future. Two types of approaches are principally used to compute binding free energies with MD simulations. The most widely known is the alchemical double decoupling method, in which the interaction of the ligand with its surroundings are progressively switched off. It is also possible to use a potential of mean force (PMF) method, in which the ligand is physically separated from the protein receptor. For both of these computational approaches, restraining potentials may be activated and released during the simulation for sampling efficiently the changes in translational, rotational, and conformational freedom of the ligand and protein upon binding. Because such restraining potentials add bias to the simulations, it is important that their effects be rigorously removed to yield a binding free energy that is properly unbiased with respect to the standard state. A review of recent results is presented, and differences in computational methods are discussed. Examples of computations with T4-lysozyme mutants, FKBP12, SH2 domain, and cytochrome P450 are discussed and compared. Remaining difficulties and challenges are highlighted.

Configurational entropy and cooperativity between ligand binding and dimerization in glycopeptide antibiotics

Oligomerization and ligand binding are thermodynamically cooperative processes in many biochemical systems, and the mechanisms giving rise to cooperative behavior are generally attributed to changes in structure. In glycopeptide antibiotics, however, these cooperative processes are not accompanied by significant structural changes. To investigate the mechanism by which cooperativity arises in these compounds, fully solvated molecular dynamics simulations and quasiharmonic normal-mode analysis were performed on chloroeremomycin, vancomycin, and dechlorovancomycin. Configurational entropies were derived from the vibrational modes recovered from ligand-free and ligand-bound forms of the monomeric and dimeric species. Results indicate that both ligand binding and dimerization incur an entropic cost as vibrational activity in the central core of the antibiotic is shifted to higher frequencies with lower amplitudes. Nevertheless, ligand binding and dimerization are cooperative because the entropic cost of both processes occurring together is less than the cost of these processes occurring separately. These reductions in configurational entropy are more than sufficient in magnitude to account for the experimentally observed cooperativity between dimerization and ligand binding. We conclude that biochemical cooperativity can be mediated through changes in vibrational activity, irrespective of the presence or absence of concomitant structural change. This may represent a general mechanism of allostery underlying cooperative phenomena in diverse macromolecular systems.

The role of configurational entropy in biochemical cooperativity

Cooperativity is a common biochemical phenomenon in which two or more otherwise independent processes are thermodynamically coupled. Because cooperative processes are usually attended by changes in molecular conformation, thermodynamic coupling is usually attributed to an enthalpy-driven mechanism. In the family of glycopeptide antibiotics that includes vancomycin, however, cooperative phenomena occur that cannot be explained by conformational change. In this communication, we demonstrate that cooperativity in these systems can arise solely from changes in vibrational activity.

A rational strategy for enhancing the affinity of vancomycin towards depsipeptide ligands

Glycopeptide antibiotics with enhanced affinity for model depsipeptide ligands may also exhibit enhanced efficacy against bacteria exhibiting the vanA resistance phenotype. To design modified agents with enhanced affinity for these ligands, and better understand why traditional agents have low affinity for depsipeptide ligands, free energy perturbation studies were performed on vancomycin derivatives by means of molecular dynamics simulation. The results suggest that modifications of the asparagine side chain on residue 3 of the antibiotic which enhance its hydrophobicity will enhance the affinity of glycopeptide antibiotics for depsipeptide ligands, and act synergistically with other modifications that enhance the efficacy of these agents against vanA-positive bacteria.