Free Energy Difference Research Articles

Optimizing the Calculation of Free Energy Differences in Nonequilibrium Work SQM/MM Switching Simulations.

A key step during indirect alchemical free energy simulations using quantum mechanical/molecular mechanical (QM/MM) hybrid potential energy functions is the calculation of the free energy difference ΔAlow→high between the low level (e.g., pure MM) and the high level of theory (QM/MM). A reliable approach uses nonequilibrium work (NEW) switching simulations in combination with Jarzynski’s equation; however, it is computationally expensive. In this study, we investigate whether it is more efficient to use more shorter switches or fewer but longer switches. We compare results obtained with various protocols to reference free energy differences calculated with Crooks’ equation. The central finding is that fewer longer switches give better converged results. As few as 200 sufficiently long switches lead to ΔAlow→high values in good agreement with the reference results. This optimized protocol reduces the computational cost by a factor of 40 compared to earlier work. We also describe two tools/ways of analyzing the raw data to detect sources of poor convergence. Specifically, we find it helpful to analyze the raw data (work values from the NEW switching simulations) in a quasi-time series-like manner. Principal component analysis helps to detect cases where one or more conformational degrees of freedom are different at the low and high level of theory.

Is styrene competitive for dopamine receptor binding?

The potential role of styrene oxide in altering the dopaminergic pathway in the ear is investigated by means of molecular docking and molecular dynamics simulations. We estimate the binding affinity of both styrene oxide and dopamine to the dopaminergic receptor DrD2 by computing the free-energy difference, ∆G, between the configuration where the ligand is bound to the receptor and the situation in which it is "infinitely" far away from it. The results show that the styrene oxide has a somewhat lower affinity for binding with respect to dopamine, which, however, may not be enough to prevent exogenous high concentration styrene oxide to compete with endogenous dopamine for DrD2 binding.

Peptide Permeation across a Phosphocholine Membrane: An Atomically Detailed Mechanism Determined through Simulations and Supported by Experimentation.

Cell-penetrating peptides (CPPs) facilitate translocation across biological membranes and are of significant biological and medical interest. Several CPPs can permeate into specific cells and organelles. We examine the incorporation and translocation of a novel anticancer CPP in a dioleoylphosphatidylcholine (DOPC) lipid bilayer membrane. The peptide, NAF-144-67, is a short fragment of a transmembrane protein, consisting of hydrophobic N-terminal and charged C-terminal segments. Experiments using fluorescently labeled NAF-144-67 in ∼100 nm DOPC vesicles and atomically detailed simulations conducted with Milestoning support a model in which a significant barrier for peptide-membrane entry is found at the interface between the aqueous solution and membrane. The initial step is the insertion of the N-terminal segment and the hydrophobic helix into the membrane, passing the hydrophilic head groups. Both experiments and simulations suggest that the free energy difference in the first step of the permeation mechanism in which the hydrophobic helix crosses the phospholipid head groups is -0.4 kcal mol-1 slightly favoring motion into the membrane. Milestoning calculations of the mean first passage time and the committor function underscore the existence of an early polar barrier followed by a diffusive barrierless motion in the lipid tail region. Permeation events are coupled to membrane fluctuations that are examined in detail. Our study opens the way to investigate in atomistic resolution the molecular mechanism, kinetics, and thermodynamics of CPP permeation to diverse membranes.

Effect of Al content on the δ→γ transformation and M23C6 precipitation behavior in 18Cr-xAl-Si ferritic heat-resistant stainless steel during solidification process

Aimed at clarifying δ→γ transformation and M23C6 precipitation behavior in18Cr-xAl-Si (x = 0, 0.43wt.%, 0.91wt.% and 1.82wt.%) ferritic heat-resistant stainless steel (FHSS), high temperature confocal scanning laser microscopy (CSLM) was used for in-situ observation during solidification. The results demonstrated that increasing of Al content reduced the nucleation temperature of δ-ferrite, and at the same time reduced the free energy difference between δ-ferrite and γ-austenite, which increased the difficulty of δ→γ transformation, in turn reduced the temperature of δ→γ transformation. The γ-austenite preferentially precipitated at the δ-ferrite phase boundaries and gradually along the sides of the δ-ferrite phase boundaries as the temperature decreased, and part of the martensite transformed inside the γ-austenite phases. M23C6 carbides mainly precipitated at the phase boundaries between δ-ferrite and γ-austenite/martensite. With the increase of Al content, the amount of γ-austenite/martensite and M23C6 carbides decreased. In addition, the solidification mechanism of the complex nucleus was also discussed.

Predicting Melt Curves of Energetic Materials Using Molecular Models

AbstractIn this work, the solid–liquid coexistence curves of classical fully flexible atomistic models of α‐RDX and β‐HMX were calculated using thermodynamically rigorous methodologies that identify where the free energy difference between the phases is zero. The free energy difference between each phase at a given state point was computed using the pseudosupercritical path (PSCP) method, and Gibbs–Helmholtz integration was used to evaluate the solid–liquid free energy difference as a function of temperature. This procedure was repeated for several pressures to determine points along the coexistence curve, which were then fit to the Simon–Glatzel functional form. While effective, this method is computationally expensive. An alternative approach is to compute the melting point at a single pressure via the PSCP method, and then use the Gibbs–Duhem integration technique to trace out the coexistence curve in a more computationally economical manner. Both approaches were used to determine the coexistence curve of α‐RDX. The Gibbs–Duhem integration method was shown to generate a melt curve that is in good agreement with the PSCP‐derived melt curve, while only costing ∼10 % of the computational resources used for the PSCP method. For α‐RDX, the predicted melting temperature increases significantly more for a given increase in pressure when compared to available experimental data.

Obtaining QM/MM binding free energies in the SAMPL8 drugs of abuse challenge: indirect approaches

Accurately predicting free energy differences is essential in realizing the full potential of rational drug design. Unfortunately, high levels of accuracy often require computationally expensive QM/MM Hamiltonians. Fortuitously, the cost of employing QM/MM approaches in rigorous free energy simulation can be reduced through the use of the so-called “indirect” approach to QM/MM free energies, in which the need for QM/MM simulations is avoided via a QM/MM “correction” at the classical endpoints of interest. Herein, we focus on the computation of QM/MM binding free energies in the context of the SAMPL8 Drugs of Abuse host–guest challenge. Of the 5 QM/MM correction coupled with force-matching submissions, PM6-D3H4/MM ranked submission proved the best overall QM/MM entry, with an RMSE from experimental results of 2.43 kcal/mol (best in ranked submissions), a Pearson’s correlation of 0.78 (second-best in ranked submissions), and a Kendall tau correlation of 0.52 (best in ranked submissions).

Insight into anisotropic magnetocaloric effect of CrI[formula omitted

CrI3 is considered to be a promising candidate for spintronic devices and data storage. We derived the Heisenberg Hamiltonian for CrI3 from density functional calculations using the Liechtenstein formula. Moreover, the Monte–Carlo simulations with the Sucksmith–Thompson method were performed to analyze the effect of magnetic anisotropy energy on the thermodynamic properties. Our method successfully reproduced the negative sign of isothermal magnetic entropy changes when a magnetic field was applied along the hard plane. We found that the temperature dependence of the magnetocrystalline anisotropy energy is not negligible at temperatures slightly above the Curie temperature. We clarified that the origin of this phenomenon is attributed to anisotropic magnetic susceptibility and magnetization anisotropy. The difference between the entropy change of the easy axis and the hard plane is proportional to the temperature dependence of the magnetic anisotropy energy, implying that the anisotropic entropy term is the main source of the temperature dependence of the free energy difference when magnetizing in a specific direction other than the easy axis. We also investigated the magnetic susceptibility that can be used for the characterization of the negative sign of the entropy change in the case of a hard plane. The competition of the magnetocrystalline anisotropy energy and external magnetic field at a low temperature and low magnetic field causes a high magnetic susceptibility as the magnetization fluctuates. Meanwhile, the anisotropy energy is suppressed as a sufficient magnetic field is applied along the hard axis, and the magnetization is fully rotated to the direction of the external magnetic field.

Stress Regulation of the Oxygen Reduction Reaction on a Pt (100) Surface Using First Principles Calculations

The oxygen reduction reaction (ORR) on a strained Pt (100) surface has been studied using periodic density functional theory (DFT) calculations. The results showed a linear relationship between strain and the d-band center. The adsorption energy displayed greater sensitivity to compressive strain and the adsorption energy results indicated that the bridge site was more stable after the adsorption of OH, O and OOH. Concomitantly, the calculated free energy diagram results confirmed that compression deformation could reduce the over potential of the ORR reaction and increase the catalyst activity on the Pt (100) surface. The evolution of H2O exhibited the largest free energy difference, suggesting that it was a potentially rate limiting step. Notably, the limiting potentials for the individual steps were fitted by the OOH vs. OH or O scaling relationship and the maximum overpotential value for the ORR reaction was 0.820 V when the compressive strain reached −10%. This study provided mechanistic and quantitative guidance for the rational selection of stress regulation strategies to fabricate Pt catalysts with desirable electrocatalytic activity in emerging energy applications.

Influence of Air Cold Plasma Modification on the Surface Properties of Paper Used for Packaging Production

In this study, the effect of air plasma on the surface properties of printed and coated cardboard was investigated. The material was activated by low-pressure cold plasma for 1, 10, 20, and 30 s. Wettability changes on the surface were examined by contact angle measurements using the sessile droplet technique. The differences in the surface free energy were calculated using the Lifshitz–van der Waals/acid–base and Contact Angle Hysteresis approaches. Optical profilometry was used for the surface roughness evaluation and an X-ray photoelectron spectroscopy analysis was performed to find changes in surface chemistry. Adhesive strength tests were carried out to estimate the adhesion changes after the material’s modification. It was found that the water and formamide contact angles increased after the plasma treatment while the diiodomethane contact angle did not change. As a result of the modification, the surface free energy also increased significantly and the surface roughness increased. The pull-off tests confirmed the improvement in the material’s surface properties. Moreover, it was demonstrated that the optimal effect can be obtained after just 10 s of the plasma process.

Improved Binding Affinity of Omicron's Spike Protein for the Human Angiotensin-Converting Enzyme 2 Receptor Is the Key behind Its Increased Virulence.

The new variant of severe acute respiratory syndrome coronavirus type 2 (SARS-CoV-2), Omicron, has been quickly spreading in many countries worldwide. Compared to the original virus, Omicron is characterized by several mutations in its genomic region, including the spike protein’s receptor-binding domain (RBD). We have computationally investigated the interaction between the RBD of both the wild type and Omicron variant of SARS-CoV-2 with the human angiotensin-converting enzyme 2 (hACE2) receptor using molecular dynamics and molecular mechanics-generalized Born surface area (MM-GBSA)-based binding free energy calculations. The mode of the interaction between Omicron’s RBD with the hACE2 receptor is similar to the original SARS-CoV-2 RBD except for a few key differences. The binding free energy difference shows that the spike protein of Omicron has an increased affinity for the hACE2 receptor. The mutated residues in the RBD showed strong interactions with a few amino acid residues of hACE2. More specifically, strong electrostatic interactions (salt bridges) and hydrogen bonding were observed between R493 and R498 residues of the Omicron RBD with D30/E35 and D38 residues of the hACE2, respectively. Other mutated amino acids in the Omicron RBD, e.g., S496 and H505, also exhibited hydrogen bonding with the hACE2 receptor. A pi-stacking interaction was also observed between tyrosine residues (RBD-Tyr501: hACE2-Tyr41) in the complex, which contributes majorly to the binding free energies and suggests that this is one of the key interactions stabilizing the formation of the complex. The resulting structural insights into the RBD:hACE2 complex, the binding mode information within it, and residue-wise contributions to the free energy provide insight into the increased transmissibility of Omicron and pave the way to design and optimize novel antiviral agents.

Inferring work by quantum superposing forward and time-reversal evolutions

The study of thermodynamic fluctuations allows one to relate the free energy difference between two equilibrium states with the work done on a system through processes far from equilibrium. This finding plays a crucial role in the quantum regime, where the definition of work becomes non-trivial. Based on these relations, here we develop a simple interferometric method allowing a direct estimation of the work distribution and the average dissipative work during a driven thermodynamic process by superposing the forward and time-reversal evolutions of the process. We show that our scheme provides useful upper bounds on the average dissipative work even without full control over the thermodynamic process, and we propose methodological variations depending on the possible experimental limitations encountered. Finally, we exemplify its applicability by an experimental proposal for implementing our method on a quantum photonics system, on which the thermodynamic process is performed through polarization rotations induced by liquid crystals acting in a discrete temporal regime.

Exploring the Effect of Enhanced Sampling on Protein Stability Prediction.

Changes in protein stability due to side-chain mutations are evaluated by alchemical free-energy calculations based on classical molecular dynamics (MD) simulations in explicit solvent using the GROMOS force field. Three proteins of different complexity with a total number of 93 single-point mutations are analyzed, and the relative free-energy differences are discussed with respect to configurational sampling and (dis)agreement with experimental data. For the smallest protein studied, a 34-residue WW domain, the starting structure dependence of the alchemical free-energy changes, is discussed in detail. Deviations from previous simulations for the two other proteins are shown to result from insufficient sampling in the earlier studies. Hamiltonian replica exchange in combination with multiple starting structures and sufficient sampling time of more than 100 ns per intermediate alchemical state is required in some cases to reach convergence.

A density-functional benchmark of vibrational free-energy corrections for molecular crystal polymorphism.

Many crystal structure prediction protocols only concern themselves with the electronic energy of molecular crystals. However, vibrational contributions to the free energy (Fvib) can be significant in determining accurate stability rankings for crystal candidates. While force-field studies have been conducted to gauge the magnitude of these free-energy corrections, highly accurate results from quantum mechanical methods, such as density-functional theory (DFT), are desirable. Here, we introduce the PV17 set of 17 polymorphic pairs of organic molecular crystals, for which plane wave DFT is used to calculate the vibrational free energies and free-energy differences (ΔFvib) between each pair. Our DFT results confirm that the vibrational free-energy corrections are small, having a mean value of 1.0 kJ/mol and a maximum value of 2.3 kJ/mol for the PV17 set. Furthermore, we assess the accuracy of a series of lower-cost DFT, semi-empirical, and force-field models for computing ΔFvib that have been proposed in the literature. It is found that calculating Fvib using the Γ-point frequencies does not provide ΔFvib values of sufficiently high quality. In addition, ΔFvib values calculated using various approximate methods have mean absolute errors relative to our converged DFT results of equivalent or larger magnitude than the vibrational free-energy corrections themselves. Thus, we conclude that, in a crystal structure prediction protocol, it is preferable to forego the inclusion of vibrational free-energy corrections than to estimate them with any of the approximate methods considered here.

Searching for a mechanistic description of pairwise epistasis in protein systems.

When two or more amino acid mutations occur in protein systems, they can interact in a nonadditive fashion termed epistasis. One way to quantify epistasis between mutation pairs in protein systems is by using free energy differences: ϵ=ΔΔG1,2 - (ΔΔG1 +ΔΔG2 ) where ΔΔG refers to the change in the Gibbs free energy, subscripts 1 and 2 refer to single mutations in arbitrary order and 1,2 refers to the double mutant. In this study, we explore possible biophysical mechanisms that drive pairwise epistasis in both protein-protein binding affinity and protein folding stability. Using the largest available datasets containing experimental protein structures and free energy data, we derived statistical models for both binding and folding epistasis (ϵ) with similar explanatory power (R2 ) of .299 and .258, respectively. These models contain terms and interactions that are consistent with intuition. For example, increasing the Cartesian separation between mutation sites leads to a decrease in observed epistasis for both folding and binding. Our results provide insight into factors that contribute to pairwise epistasis in protein systems and their importance in explaining epistasis. However, the low explanatory power indicates that more study is needed to fully understand this phenomenon.

SARS-CoV-2 Membrane Protein: From Genomic Data to Structural New Insights.

Severe Acute Respiratory Syndrome CoronaVirus-2 (SARS-CoV-2) is composed of four structural proteins and several accessory non-structural proteins. SARS-CoV-2’s most abundant structural protein, Membrane (M) protein, has a pivotal role both during viral infection cycle and host interferon antagonism. This is a highly conserved viral protein, thus an interesting and suitable target for drug discovery. In this paper, we explain the structural nature of M protein homodimer. To do so, we developed and applied a detailed and robust in silico workflow to predict M protein dimeric structure, membrane orientation, and interface characterization. Single Nucleotide Polymorphisms (SNPs) in M protein were retrieved from over 1.2 M SARS-CoV-2 genomes and proteins from the Global Initiative on Sharing All Influenza Data (GISAID) database, 91 of which were located at the predicted dimer interface. Among those, we identified SNPs in Variants of Concern (VOC) and Variants of Interest (VOI). Binding free energy differences were evaluated for dimer interfacial SNPs to infer mutant protein stabilities. A few high-prevalent mutated residues were found to be especially relevant in VOC and VOI. This realization may be a game-changer to structure-driven formulation of new therapeutics for SARS-CoV-2.

Determination of Equilibrium Adsorbed Morphologies of Surfactants at Metal-Water Interfaces Using a Modified Umbrella Sampling-Based Methodology.

Surfactants adsorb to metal-water interfaces in various morphologies, including self-assembled monolayers (SAMs), cylindrical and spherical micelles, or hemimicelles. Current molecular simulation methods are unable to efficiently sample the formation of these morphologies because of the large diffusive/energetic barriers. We introduce a modified umbrella sampling-based methodology that allows sampling of these morphologies from any initial configuration and provides free energy differences between them. Using this methodology, we have studied adsorption behavior of cationic [quaternary ammonium (quat) of 4 and 12 carbon long alkyl tails], uncharged [decanethiol], and anionic [phosphate monoester] surfactants and their mixtures at a gold-water interface. We find that while Coulombic repulsion between the charged head groups of quat-4 limits their adsorption to a sparse layer, stronger hydrophobic interactions between the alkyl tails of quat-12 promote adsorption resulting in a morphology with adsorbed hemispherical micelles sitting atop a monolayer. Decanethiol molecules adsorb in a densely packed bilayer with the molecules standing-up on the surface in the first layer and lying parallel to the surface in the second layer. Cationic and anionic surfactant mixtures display a synergistic adsorption behavior. These results elucidate the role of molecular characteristics in dictating the nature of adsorbed morphologies of surfactants at metal-water interfaces.

Mechanistic understanding of the aspect ratio-dependent adjuvanticity of engineered aluminum oxyhydroxide nanorods in prophylactic vaccines

Aluminum oxyhydroxide (AlOOH) adjuvants are widely used in human vaccines. However, the interaction mechanisms at the material-bio interface, and further understandings on physicochemical property-dependent modulation of the immune responses still remain uncertain. Herein, a library of AlOOH nanorods with well-defined aspect ratios is designed to explore the mechanisms of adjuvanticity. The aspect ratios of AlOOH nanorods were demonstrated to be intrinsically modulated by the hydroxide supersaturation level during crystal growth, leading to the differences in surface free energy (SFE). As a result, higher aspect ratio AlOOH nanoadjuvants with lower SFE exhibited more hydrophobic surface, resulting in more membrane depolarization, cellular uptake and dendritic cell (DC) activation. By using hepatitis B surface antigen (HBsAg) virus-like particles (VLPs) or SARS-CoV-2 spike protein receptor-binding domain (RBD) as model antigens, AlOOH nanorods with higher aspect ratio were determined to elicit more potent humoral immune responses, which could be attributed to the enhanced DC activation and the efficient antigen trafficking to the draining lymph nodes. Our findings highlight the critical role of aspect ratio of AlOOH nanorods in modulating adjuvanticity, and further provide a design strategy for engineered nanoadjuvants for prophylactic vaccines.

Multiscale Study of the Charge Transport Properties of Silicone Rubber Oligomers.

Polymer dielectrics are widely used as insulation in electronic and electrical equipment. The charge transport process in polymer dielectrics under a high electric field, which is considered to result in electrical aging and even breakdown of equipment insulation, has attracted considerable attention. Based on a localized charge transfer model, the charge transport mechanism for typical silicone rubber (SR) insulating materials was studied by multiscale methods. A model of silicone rubber oligomers under standard conditions was generated by classical molecular dynamics. The frontier molecular orbitals and projected density of states for the SR oligomer were obtained via quantum chemistry methods. The electronic coupling, reorganization energy, and free energy difference for both electron and hole transfer processes between adjacent SR molecules in the molecular dynamics model were calculated. Both hole transfer and electron transfer in SR conform to an intermolecular hopping mechanism due to their high intramolecular reorganization energy and low intermolecular electronic coupling. The results of normal mode analysis for reorganization energy indicate that the high-temperature approximation holds for charge transport in SR around or above room temperature. The charge transfer trajectory and charge mobility in SR were simulated based on kinetic Monte Carlo simulations. The hole and electron mobilities at room temperature were calculated to be 7.24 × 10-11 and 2.76 × 10-9 cm2/Vs, respectively, which agrees with the experimental data. Both electron transport and hole transport in SR show thermal activation characteristics, with corresponding activation energies of 358 and 314 meV, respectively. This work suggests a physical model to describe the charge transport process in SR polymer dielectrics.

Computational prediction of the effect of mutations in the receptor-binding domain on the interaction between SARS-CoV-2 and human ACE2.

The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) causing COVID-19 continues to mutate. Numerous studies have indicated that this viral mutation, particularly in the receptor-binding domain area, may increase the viral affinity for human angiotensin-converting enzyme 2 (hACE2), the receptor for viral entry into host cells, thereby increasing viral virulence and transmission. In this study, we investigated the binding affinity of SARS-CoV-2 variants (Delta plus, Iota, Kappa, Mu, Lambda, and C.1.2) on hACE2 using computational modeling with a protein-protein docking approach. The simulation results indicated that there were differences in the interactions between the RBD and hACE2, including hydrogen bonding, salt bridge interactions, non-bonded interactions, and binding free energy differences among these variants. Molecular dynamics simulations revealed that mutations in the RBD increase the stability of the hACE2-spike protein complex relative to the wild type, following the global stability trend and increasing the binding affinity. The value of binding-free energy calculated using molecular mechanics/Poisson-Boltzmann surface area (MM/PBSA) indicated that all mutations in the spike protein increased the contagiousness of SARS-CoV-2 variants. The findings of this study provide a foundation for developing effective interventions against these variants. Computational modeling elucidates that the spike protein of SARS-CoV-2 variants binds considerably stronger than the wild-type to hACE2.

A complete description of thermodynamic stabilities of molecular crystals

Predictions of relative stabilities of (competing) molecular crystals are of great technological relevance, most notably for the pharmaceutical industry. However, they present a long-standing challenge for modeling, as often minuscule free energy differences are sensitively affected by the description of electronic structure, the statistical mechanics of the nuclei and the cell, and thermal expansion. The importance of these effects has been individually established, but rigorous free energy calculations for general molecular compounds, which simultaneously account for all effects, have hitherto not been computationally viable. Here we present an efficient "end to end" framework that seamlessly combines state-of-the art electronic structure calculations, machine-learning potentials, and advanced free energy methods to calculate ab initio Gibbs free energies for general organic molecular materials. The facile generation of machine-learning potentials for a diverse set of polymorphic compounds-benzene, glycine, and succinic acid-and predictions of thermodynamic stabilities in qualitative and quantitative agreement with experiments highlight that predictive thermodynamic studies of industrially relevant molecular materials are no longer a daunting task.