CHARMM-GUI QM/MM Interfacer for a Quantum Mechanical and Molecular Mechanical (QM/MM) Simulation Setup: 1. Semiempirical Methods.

Quantum mechanical (QM) treatments, when combined with molecular mechanical (MM) force fields, can effectively handle enzyme-catalyzed reactions without significantly increasing the computational cost. In this context, we present CHARMM-GUI QM/MM Interfacer, a web-based cyberinfrastructure designed to streamline the preparation of various QM/MM simulation inputs with ligand modification. The development of QM/MM Interfacer has been achieved through integration with existing CHARMM-GUI modules, such as PDB Reader and Manipulator, Solution Builder, and Membrane Builder. In addition, new functionalities have been developed to facilitate the one-stop preparation of QM/MM systems and enable interactive and intuitive ligand modifications and QM atom selections. QM/MM Interfacer offers support for a range of semiempirical QM methods, including AM1(+/d), PM3(+/PDDG), MNDO(+/d, +/PDDG), PM6, RM1, and SCC-DFTB, tailored for both AMBER and CHARMM. A nontrivial setup related to ligand modification, link-atom insertion, and charge distribution is automatized through intuitive user interfaces. To illustrate the robustness of QM/MM Interfacer, we conducted QM/MM simulations of three enzyme-substrate systems: dihydrofolate reductase, insulin receptor kinase, and oligosaccharyltransferase. In addition, we have created three tutorial videos about building these systems, which can be found at https://www.charmm-gui.org/demo/qmi. QM/MM Interfacer is expected to be a valuable and accessible web-based tool that simplifies and accelerates the setup process for hybrid QM/MM simulations.

Perspectives on Computational Enzyme Modeling: From Mechanisms to Design and Drug Development.

Understanding enzyme mechanisms is essential for unraveling the complex molecular machinery of life. In this review, we survey the field of computational enzymology, highlighting key principles governing enzyme mechanisms and discussing ongoing challenges and promising advances. Over the years, computer simulations have become indispensable in the study of enzyme mechanisms, with the integration of experimental and computational exploration now established as a holistic approach to gain deep insights into enzymatic catalysis. Numerous studies have demonstrated the power of computer simulations in characterizing reaction pathways, transition states, substrate selectivity, product distribution, and dynamic conformational changes for various enzymes. Nevertheless, significant challenges remain in investigating the mechanisms of complex multistep reactions, large-scale conformational changes, and allosteric regulation. Beyond mechanistic studies, computational enzyme modeling has emerged as an essential tool for computer-aided enzyme design and the rational discovery of covalent drugs for targeted therapies. Overall, enzyme design/engineering and covalent drug development can greatly benefit from our understanding of the detailed mechanisms of enzymes, such as protein dynamics, entropy contributions, and allostery, as revealed by computational studies. Such a convergence of different research approaches is expected to continue, creating synergies in enzyme research. This review, by outlining the ever-expanding field of enzyme research, aims to provide guidance for future research directions and facilitate new developments in this important and evolving field.

Elucidating Dynamics of Adenylate Kinase from Enzyme Opening to Ligand Release.

This study explores ligand-driven conformational changes in adenylate kinase (AK), which is known for its open-to-close conformational transitions upon ligand binding and release. By utilizing string free energy simulations, we determine the free energy profiles for both enzyme opening and ligand release and compare them with profiles from the apoenzyme. Results reveal a three-step ligand release process, which initiates with the opening of the adenosine triphosphate-binding subdomain (ATP lid), followed by ligand release and concomitant opening of the adenosine monophosphate-binding subdomain (AMP lid). The ligands then transition to nonspecific positions before complete dissociation. In these processes, the first step is energetically driven by ATP lid opening, whereas the second step is driven by ATP release. In contrast, the AMP lid opening and its ligand release make minor contributions to the total free energy for enzyme opening. Regarding the ligand binding mechanism, our results suggest that AMP lid closure occurs via an induced-fit mechanism triggered by AMP binding, whereas ATP lid closure follows conformational selection. This difference in the closure mechanisms provides an explanation with implications for the debate on ligand-driven conformational changes of AK. Additionally, we determine an X-ray structure of an AK variant that exhibits significant rearrangements in the stacking of catalytic arginines, explaining its reduced catalytic activity. In the context of apoenzyme opening, the sequence of events is different. Here, the AMP lid opens first while the ATP lid remains closed, and the free energy associated with ATP lid opening varies with orientation, aligning with the reported AK opening and closing rate heterogeneity. Finally, this study, in conjunction with our previous research, provides a comprehensive view of the intricate interplay between various structural elements, ligands, and catalytic residues that collectively contribute to the robust catalytic power of the enzyme.

Free Energy Profile Decomposition Analysis for QM/MM Simulations of Enzymatic Reactions.

In enzyme mechanistic studies and mutant design, it is highly desirable to know the individual residue contributions to the reaction free energy and barrier. In this work, we show that such free energy contributions from each residue can be readily obtained by postprocessing ab initio quantum mechanical molecular mechanical (ai-QM/MM) free energy simulation trajectories. Specifically, through a mean force integration along the minimum free energy pathway, one can obtain the electrostatic, polarization, and van der Waals contributions from each residue to the free energy barrier. Separately, a similar analysis procedure allows us to assess the contribution from different collective variables along the reaction coordinate. The chorismate mutase reaction is used to demonstrate the utilization of these two trajectory analysis tools.

Molecular Simulations of Conformational Transitions within the Insulin Receptor Kinase Reveal Consensus Features in a Multistep Activation Pathway.

Modulating the transitions between active and inactive conformations of protein kinases is the primary means of regulating their catalytic activity, achieved by phosphorylation of the activation loop (A-loop). To elucidate the mechanism of this conformational activation, we applied the string method to determine the conformational transition path of insulin receptor kinase between the active and inactive conformations and the corresponding free-energy profiles with and without A-loop phosphorylation. The conformational change was found to proceed in three sequential steps: first, the flipping of the DFG motif of the active site; second, rotation of the A-loop; finally, the inward movement of the αC helix. The main energetic bottleneck corresponds to the conformational change in the A-loop, while changes in the DFG motif and αC helix occur before and after A-loop conformational change, respectively. In accordance with this, two intermediate states are identified, the first state just after the DFG flipping and the second state after the A-loop rotation. These intermediates exhibit structural features characteristic of the corresponding inactive and active conformations of other protein kinases. To understand the impact of A-loop phosphorylation on kinase conformation, the free energies of A-loop phosphorylation were determined at several states along the conformational transition path using the free-energy perturbation simulations. The calculated free energies reveal that while the unphosphorylated kinase interconverts between the inactive and active conformations, A-loop phosphorylation restricts access to the inactive conformation, thereby increasing the active conformation population. Overall, this study suggests a consensus mechanism of conformational activation between different protein kinases.

Protein dynamics: The future is bright and complicated!

Biological life depends on motion, and this manifests itself in proteins that display motion over a formidable range of time scales spanning from femtoseconds vibrations of atoms at enzymatic transition states, all the way to slow domain motions occurring on micro to milliseconds. An outstanding challenge in contemporary biophysics and structural biology is a quantitative understanding of the linkages among protein structure, dynamics, and function. These linkages are becoming increasingly explorable due to conceptual and methodological advances. In this Perspective article, we will point toward future directions of the field of protein dynamics with an emphasis on enzymes. Research questions in the field are becoming increasingly complex such as the mechanistic understanding of high-order interaction networks in allosteric signal propagation through a protein matrix, or the connection between local and collective motions. In analogy to the solution to the "protein folding problem," we argue that the way forward to understanding these and other important questions lies in the successful integration of experiment and computation, while utilizing the present rapid expansion of sequence and structure space. Looking forward, the future is bright, and we are in a period where we are on the doorstep to, at least in part, comprehend the importance of dynamics for biological function.

Mechanistic Basis for a Connection between the Catalytic Step and Slow Opening Dynamics of Adenylate Kinase.

Escherichia coli adenylate kinase (AdK) is a small, monomeric enzyme that synchronizes the catalytic step with the enzyme's conformational dynamics to optimize a phosphoryl transfer reaction and the subsequent release of the product. Guided by experimental measurements of low catalytic activity in seven single-point mutation AdK variants (K13Q, R36A, R88A, R123A, R156K, R167A, and D158A), we utilized classical mechanical simulations to probe mutant dynamics linked to product release, and quantum mechanical and molecular mechanical calculations to compute a free energy barrier for the catalytic event. The goal was to establish a mechanistic connection between the two activities. Our calculations of the free energy barriers in AdK variants were in line with those from experiments, and conformational dynamics consistently demonstrated an enhanced tendency toward enzyme opening. This indicates that the catalytic residues in the wild-type AdK serve a dual role in this enzyme's function─one to lower the energy barrier for the phosphoryl transfer reaction and another to delay enzyme opening, maintaining it in a catalytically active, closed conformation for long enough to enable the subsequent chemical step. Our study also discovers that while each catalytic residue individually contributes to facilitating the catalysis, R36, R123, R156, R167, and D158 are organized in a tightly coordinated interaction network and collectively modulate AdK's conformational transitions. Unlike the existing notion of product release being rate-limiting, our results suggest a mechanistic interconnection between the chemical step and the enzyme's conformational dynamics acting as the bottleneck of the catalytic process. Our results also suggest that the enzyme's active site has evolved to optimize the chemical reaction step while slowing down the overall opening dynamics of the enzyme.

Insights into the evolution of enzymatic specificity and catalysis: From Asgard archaea to human adenylate kinases.

Enzymatic catalysis is critically dependent on selectivity, active site architecture, and dynamics. To contribute insights into the interplay of these properties, we established an approach with NMR, crystallography, and MD simulations focused on the ubiquitous phosphotransferase adenylate kinase (AK) isolated from Odinarchaeota (OdinAK). Odinarchaeota belongs to the Asgard archaeal phylum that is believed to be the closest known ancestor to eukaryotes. We show that OdinAK is a hyperthermophilic trimer that, contrary to other AK family members, can use all NTPs for its phosphorylation reaction. Crystallographic structures of OdinAK-NTP complexes revealed a universal NTP-binding motif, while 19F NMR experiments uncovered a conserved and rate-limiting dynamic signature. As a consequence of trimerization, the active site of OdinAK was found to be lacking a critical catalytic residue and is therefore considered to be "atypical." On the basis of discovered relationships with human monomeric homologs, our findings are discussed in terms of evolution of enzymatic substrate specificity and cold adaptation.

Self-assembling Peptides with Internal Ionizable Unnatural Amino Acids: A General Approach to pH-responsive Peptide Materials.

Self-assembled peptides are an emerging family of biomaterials that show great promise for a range of biomedical and biotechnological applications. Introducing and tuning the pH-responsiveness of the assembly is highly desirable for improving their biological activities. Inspired by proteins with internal ionizable residues, we report a simple but effective approach to constructing pH-responsive peptide assembly containing unnatural ionic amino acids with an aliphatic tertiary amine side chain. Through a combined experimental and computational investigation, we demonstrate that these residues can be accommodated and stabilized within the internal hydrophobic compartment of the peptide assembly. The hydrophobic microenvironment shifts their pKa significantly from a basic pH typically found for free amines to a more biologically relevant pH in the weakly acidic range. The pH-induced ionization and ionization-dependent self-assembly and disassembly are thoroughly investigated and correlated with the biological activity of the assembly. This new approach has unique advantages in tuning the pH-responsiveness of self-assembled peptides across a large pH range in a complex biological environment. We anticipate the ionizable amino acids developed here can be widely applicable to the synthesis and self-assembly of many amphiphilic peptides with endowed pH-responsive properties to enhance their biological activities toward applications ranging from targeted therapeutic delivery to proton transport.

Accelerating Ab Initio Quantum Mechanical and Molecular Mechanical (QM/MM) Molecular Dynamics Simulations with Multiple Time Step Integration and a Recalibrated Semiempirical QM/MM Hamiltonian.

Molecular dynamics (MD) simulations employing ab initio quantum mechanical and molecular mechanical (ai-QM/MM) potentials are considered to be the state of the art, but the high computational cost associated with the ai-QM calculations remains a theoretical challenge for their routine application. Here, we present a modified protocol of the multiple time step (MTS) method for accelerating ai-QM/MM MD simulations of condensed-phase reactions. Within a previous MTS protocol [Nam J. Chem. Theory Comput. 2014, 10, 4175], reference forces are evaluated using a low-level (semiempirical QM/MM) Hamiltonian and employed at inner time steps to propagate the nuclear motions. Correction forces, which arise from the force differences between high-level (ai-QM/MM) and low-level Hamiltonians, are applied at outer time steps, where the MTS algorithm allows the time-reversible integration of the correction forces. To increase the outer step size, which is bound by the highest-frequency component in the correction forces, the semiempirical QM Hamiltonian is recalibrated in this work to minimize the magnitude of the correction forces. The remaining high-frequency modes, which are mainly bond stretches involving hydrogen atoms, are then removed from the correction forces. When combined with a Langevin or SIN(R) thermostat, the modified MTS-QM/MM scheme remains robust with an up to 8 (with Langevin) or 10 fs (with SIN(R)) outer time step (with 1 fs inner time steps) for the chorismate mutase system. This leads to an over 5-fold speedup over standard ai-QM/MM simulations, without sacrificing the accuracy in the predicted free energy profile of the reaction.

Machine-Learning-Assisted Free Energy Simulation of Solution-Phase and Enzyme Reactions.

Despite recent advances in the development of machine learning potentials (MLPs) for biomolecular simulations, there has been limited effort on developing stable and accurate MLPs for enzymatic reactions. Here we report a protocol for performing machine-learning-assisted free energy simulation of solution-phase and enzyme reactions at the ab initio quantum-mechanical/molecular-mechanical (ai-QM/MM) level of accuracy. Within our protocol, the MLP is built to reproduce the ai-QM/MM energy and forces on both QM (reactive) and MM (solvent/enzyme) atoms. As an alternative strategy, a delta machine learning potential (ΔMLP) is trained to reproduce the differences between the ai-QM/MM and semiempirical (se) QM/MM energies and forces. To account for the effect of the condensed-phase environment in both MLP and ΔMLP, the DeePMD representation of a molecular system is extended to incorporate the external electrostatic potential and field on each QM atom. Using the Menshutkin and chorismate mutase reactions as examples, we show that the developed MLP and ΔMLP reproduce the ai-QM/MM energy and forces with errors that on average are less than 1.0 kcal/mol and 1.0 kcal mol-1 Å-1, respectively, for representative configurations along the reaction pathway. For both reactions, MLP/ΔMLP-based simulations yielded free energy profiles that differed by less than 1.0 kcal/mol from the reference ai-QM/MM results at only a fraction of the computational cost.

Dynamic Connection between Enzymatic Catalysis and Collective Protein Motions.

Enzymes employ a wide range of protein motions to achieve efficient catalysis of chemical reactions. While the role of collective protein motions in substrate binding, product release, and regulation of enzymatic activity is generally understood, their roles in catalytic steps per se remain uncertain. Here, molecular dynamics simulations, enzyme kinetics, X-ray crystallography, and nuclear magnetic resonance spectroscopy are combined to elucidate the catalytic mechanism of adenylate kinase and to delineate the roles of catalytic residues in catalysis and the conformational change in the enzyme. This study reveals that the motions in the active site, which occur on a time scale of picoseconds to nanoseconds, link the catalytic reaction to the slow conformational dynamics of the enzyme by modulating the free energy landscapes of subdomain motions. In particular, substantial conformational rearrangement occurs in the active site following the catalytic reaction. This rearrangement not only affects the reaction barrier but also promotes a more open conformation of the enzyme after the reaction, which then results in an accelerated opening of the enzyme compared to that of the reactant state. The results illustrate a linkage between enzymatic catalysis and collective protein motions, whereby the disparate time scales between the two processes are bridged by a cascade of intermediate-scale motion of catalytic residues modulating the free energy landscapes of the catalytic and conformational change processes.

Copper(I) and silver(I) chemistry of vinyltrifluoroborate supported by a bis(pyrazolyl)methane ligand.

Although unsaturated organotrifluoroborates are common synthons in metal-organic chemistry, their transition metal complexes have received little attention. [CH2(3,5-(CH3)2Pz)2]Cu(CH2[double bond, length as m-dash]CHBF3), (SIPr)Cu(MeCN)(CH2[double bond, length as m-dash]CHBF3) and [CH2(3,5-(CH3)2Pz)2]Ag(CH2[double bond, length as m-dash]CHBF3) represent rare, isolable molecules featuring a vinyltrifluoroborate ligand on coinage metals. The X-ray crystal structures show the presence of three-coordinate metal sites in these complexes. The vinyltrifluoroborate group binds asymmetrically to the metal site in [CH2(3,5-(CH3)2Pz)2]M(CH2[double bond, length as m-dash]CHBF3) (M = Cu, Ag) with relatively closer M-C(H)2 distances. The computed structures of [CH2(3,5-(CH3)2Pz)2]M(CH2[double bond, length as m-dash]CHBF3) and M(CH2[double bond, length as m-dash]CHBF3), however, have shorter M-C(H)BF3 distances than M-C(H)2. These molecules feature various inter- or intra-molecular contacts involving fluorine of the BF3 group, possibly affecting these M-C distances. The binding energies of [CH2[double bond, length as m-dash]CHBF3]- to Cu+, Ag+ and Au+ have been calculated at the wB97XD/def2-TZVP level of theory, in the presence and absence of the supporting ligand CH2(3,5-(CH3)2Pz)2. The calculation shows that Au+ has the strongest binding to the [CH2[double bond, length as m-dash]CHBF3]- ligand, followed by Cu+ and Ag+, irrespective of the presence of the supporting ligand. However, in all three metals, the supporting ligand weakens the binding of olefin to the metal. The same trends were also found from the analysis of the σ-donation and π-backbonding interactions between the metal fragment and the π and π* orbitals of [CH2[double bond, length as m-dash]CHBF3]-.

A simplified charge projection scheme for long-range electrostatics in ab initio QM/MM calculations.

In a previous work [Pan et al., Molecules 23, 2500 (2018)], a charge projection scheme was reported, where outer molecular mechanical (MM) charges [>10 Å from the quantum mechanical (QM) region] were projected onto the electrostatic potential (ESP) grid of the QM region to accurately and efficiently capture long-range electrostatics in ab initio QM/MM calculations. Here, a further simplification to the model is proposed, where the outer MM charges are projected onto inner MM atom positions (instead of ESP grid positions). This enables a representation of the long-range MM electrostatic potential via augmentary charges (AC) on inner MM atoms. Combined with the long-range electrostatic correction function from Cisneros et al. [J. Chem. Phys. 143, 044103 (2015)] to smoothly switch between inner and outer MM regions, this new QM/MM-AC electrostatic model yields accurate and continuous ab initio QM/MM electrostatic energies with a 10 Å cutoff between inner and outer MM regions. This model enables efficient QM/MM cluster calculations with a large number of MM atoms as well as QM/MM calculations with periodic boundary conditions.

Demonstration of Selective Single-Barium Ion Detection with Dry Diazacrown Ether Naphthalimide Turn-on Chemosensors.

Single-molecule fluorescence imaging (SMFI) of gas-phase ions has been proposed for "barium tagging," a burgeoning area of research in particle physics to detect individual barium daughter ions. This has potential to significantly enhance the sensitivity of searches for neutrinoless double-beta decay (0νßß) that is obscured by background radiation events. The chemistry required to make such sensitive detection of Ba2+ by SMFI in dry Xe gas at solid interfaces has implications for solid-phase detection methods but has not been demonstrated. Here, we synthesized simple, robust, and effective Ba2+-selective chemosensors capable of function within ultrapure high-pressure 136Xe gas. Turn-on fluorescent naphthalimide-(di)azacrown ether chemosensors were Ba2+-selective and achieved SMFI in a polyacrylamide matrix. Fluorescence and NMR experiments supported a photoinduced electron transfer mechanism for turn-on sensing. Ba2+ selectivity was achieved with computational calculations correctly predicting the fluorescence responses of sensors to barium, mercury, and potassium ions. With these molecules, dry-phase single-Ba2+ ion imaging with turn-on fluorescence was realized using an oil-free microscopy technique for the first time-a significant advance toward single-Ba2+ ion detection within large volumes of 136Xe, plausibly enabling a background-independent technique to search for the hypothetical process of 0νßß.

The trajectory of intrahelical lesion recognition and extrusion by the human 8-oxoguanine DNA glycosylase.

Efficient search for DNA damage embedded in vast expanses of the DNA genome presents one of the greatest challenges to DNA repair enzymes. We report here crystal structures of human 8-oxoguanine (oxoG) DNA glycosylase, hOGG1, that interact with the DNA containing the damaged base oxoG and the normal base G while they are nested in the DNA helical stack. The structures reveal that hOGG1 engages the DNA using different protein-DNA contacts from those observed in the previously determined lesion recognition complex and other hOGG1-DNA complexes. By applying molecular dynamics simulations, we have determined the pathways taken by the lesion and normal bases when extruded from the DNA helix and their associated free energy profiles. These results reveal how the human oxoG DNA glycosylase hOGG1 locates the lesions inside the DNA helix and facilitates their extrusion for repair.

Rhizolutin, a Novel 7/10/6-Tricyclic Dilactone, Dissociates Misfolded Protein Aggregates and Reduces Apoptosis/Inflammation Associated with Alzheimer's Disease.

Rhizolutin (1) was discovered as a natural product of ginseng-rhizospheric Streptomyces sp. WON17. Its structure features an unprecedented 7/10/6-tricyclic dilactone carbon skeleton composed of dimethylcyclodecatriene flanked by a 7-membered and a 6-membered lactone ring based on spectroscopic analysis. During an unbiased screening of natural product libraries, this novel compound was found to dissociate amyloid-ß (Aß) plaques and tau tangles, which are key pathological hallmarks of Alzheimer's disease (AD). Rhizolutin treatment of APP/PS1 double transgenic mice with AD significantly dissociated hippocampal plaques. In vitro, rhizolutin substantially decreased Aß-induced apoptosis and inflammation in neuronal and glial cells. Our findings introduce a unique chemical entity that targets Aß and tau concurrently by mimicking misfolded protein clearance mechanisms of immunotherapy, which is prominently investigated in clinical trials.

Structural Basis for GTP versus ATP Selectivity in the NMP Kinase AK3.

ATP and GTP are exceptionally important molecules in biology with multiple, and often discrete, functions. Therefore, enzymes that bind to either of them must develop robust mechanisms to selectively utilize one or the other. Here, this specific problem is addressed by molecular studies of the human NMP kinase AK3, which uses GTP to phosphorylate AMP. AK3 plays an important role in the citric acid cycle, where it is responsible for GTP/GDP recycling. By combining a structural biology approach with functional experiments, we present a comprehensive structural and mechanistic understanding of the enzyme. We discovered that AK3 functions by recruitment of GTP to the active site, while ATP is rejected and nonproductively bound to the AMP binding site. Consequently, ATP acts as an inhibitor with respect to GTP and AMP. The overall features with specific recognition of the correct substrate and nonproductive binding by the incorrect substrate bear a strong similarity to previous findings for the ATP specific NMP kinase adenylate kinase. Taken together, we are now able to provide the fundamental principles for GTP and ATP selectivity in the large NMP kinase family. As a side-result originating from nonlinearity of chemical shifts in GTP and ATP titrations, we find that protein surfaces offer a general and weak binding affinity for both GTP and ATP. These nonspecific interactions likely act to lower the available intracellular GTP and ATP concentrations and may have driven evolution of the Michaelis constants of NMP kinases accordingly.

Development of Noonan syndrome by deregulation of allosteric SOS autoactivation.

Ras family proteins play an essential role in several cellular functions, including growth, differentiation, and survival. The mechanism of action of Ras mutants in Costello syndrome and cancers has been identified, but the contribution of Ras mutants to Noonan syndrome, a genetic disorder that prevents normal development in various parts of the body, is unknown. Son of Sevenless (SOS) is a Ras guanine nucleotide exchange factor. In response to Ras-activating cell signaling, SOS autoinhibition is released and is followed by accelerative allosteric feedback autoactivation. Here, using mutagenesis-based kinetic and pulldown analyses, we show that Noonan syndrome Ras mutants I24N, T50I, V152G, and D153V deregulate the autoactivation of SOS to populate their active form. This previously unknown process has been linked so far only to the development of Noonan syndrome. In contrast, other Noonan syndrome Ras mutants-V14I, T58I, and G60E-populate their active form by deregulation of the previously documented Ras GTPase activities. We propose a novel mechanism responsible for the deregulation of SOS autoactivation, where I24N, T50I, V152G, and D153V Ras mutants evade SOS autoinhibition. Consequently, they are capable of forming a complex with the SOS allosteric site, thus aberrantly promoting SOS autoactivation, resulting in the population of active Ras mutants in cells. The results of this study elucidate the molecular mechanism of the Ras mutant-mediated development of Noonan syndrome.

Repulsive Soft-Core Potentials for Efficient Alchemical Free Energy Calculations.

In alchemical free energy (FE) simulations, annihilation and creation of atoms are generally achieved with the soft-core potential that shifts the interparticle separations. While this soft-core potential eliminates the numerical instability occurring near the two end states of the transformation, it makes the hybrid Hamiltonian vary nonlinearly with respect to the parameter λ, which interpolates between the Hamiltonians representing the two end states. This complicates FE estimation by Bennett acceptance ratio (BAR), free energy perturbation (FEP), and thermodynamic integration (TI) methods, thus reducing their calculation efficiency. In this work, we develop a new type of repulsive soft-core potential, called Gaussian soft-core (GSC) potential, with two parameters controlling its maximum and width. The main advantage of this potential is the linearity of the hybrid Hamiltonian with respect to λ, thus permitting the direct application of BAR, FEP, TI, and other variant FE methods. The accuracy and efficiency of the GSC potential are demonstrated by comparing the free energies of annihilation determined for 13 small molecules and an alchemical mutation of a protein side chain. In addition, in combination with a TI integrand (∂H/∂λ) estimation strategy, we show that GSC can considerably reduce the number of λ simulations compared to the commonly used separation-shifted soft-core potential.