Free Energy Perturbation Techniques Research Articles

Integrating molecular dynamics simulations and experimental data for azeotrope predictions in binary mixtures.

An azeotrope is a constant boiling point mixture, and its behavior is important for fluid separation processes. Predicting azeotropes from atomistic simulations is difficult due to the complexities and convergence problems of Monte Carlo and free-energy perturbation techniques. Here, we present a methodology for predicting the azeotropes of binary mixtures, which computes the compositional dependence of chemical potentials from molecular dynamics simulations using the S0 method and employs experimental boiling point and vaporization enthalpy data. Using this methodology, we reproduce the azeotropes, or lack thereof, in five case studies, including ethanol/water, ethanol/isooctane, methanol/water, hydrazine/water, and acetone/chloroform mixtures. We find that it is crucial to use the experimental boiling point and vaporization enthalpy for reliable azeotrope predictions, as empirical force fields are not accurate enough for these quantities. Finally, we use regular solution models to rationalize the azeotropes and reveal that they tend to form when the mixture components have similar boiling points and strong interactions.

A comparative study on the solubility of CH3SH and H2S in imidazolium-based ionic liquids: A connectionist approach between free volume and thermodynamics

Sulfur contaminants in industrial gas streams, such as hydrogen sulfide and mercaptans, must be eliminated due to safety and environmental considerations. Ionic liquids, as a new class of chemical compounds, are an alternative to conventional aqueous alkanolamine solutions currently used in absorption gas treatment processes. In this work, the solubilities of methyl mercaptan (CH3SH) and hydrogen sulfide gases in 1-butyl-3-methylimidazolium-based ionic liquids with tetrafluoroborate ([C4mim][BF4]), hexafluorophosphate ([C4mim][PF6]), trifluoromethanesulfonate ([C4mim][OTf]), and bis(trifluoromethylsulfonyl)imide ([C4mim][Tf2N]) anions have been studied using molecular dynamics (MD) simulation and a recently developed OPLS all-atom (OPLS-AA-0.8) force field. A new concept, the induced free volume, has been introduced in this work, which provides a realistic and precise connection between the thermodynamic entropic contributions to the dissolution process of CH3SH and H2S in the ionic liquids, on one hand, and the microscopic free/void volume induced as a result of the entrance of a gas solute molecule into the liquid structure, on the other hand. The thermodynamic energetic contributions were correlated to the CH3SH/H2S interaction with the ionic species of the ionic liquids. Quantum chemical energy decomposition analysis (EDA) was employed to calculate and analyze the binding interaction and nature of CH3SH/H2S interactions with the ionic species of ionic liquids. It was found that the van der Waals contributions to the interaction of CH3SH with ionic liquids are more significant compared to that of H2S, while the opposite conclusion could be deduced for the corresponding electrostatic contributions. Both CH3SH and H2S established equivalent hydrogen bonding interactions with the ionic liquids studied in this work. The diffusion coefficients of the ionic species, as well as the excess chemical potentials and Henry’s constants for the dissolution of CH3SH and H2S in the ionic liquids, were calculated by free energy perturbation technique (FEP) and compared with the MD simulated and experimental data, reported in the literature. It was concluded that the OPLS-AA-0.8 force field yields Henry’s constants and diffusivities that are in better agreement with the reliable experimental data corresponding to other OPLS versions reported to date.

Does Hamiltonian Replica Exchange via Lambda-Hopping Enhance the Sampling in Alchemical Free Energy Calculations?

In the context of computational drug design, we examine the effectiveness of the enhanced sampling techniques in state-of-the-art free energy calculations based on alchemical molecular dynamics simulations. In a paradigmatic molecule with competition between conformationally restrained E and Z isomers whose probability ratio is strongly affected by the coupling with the environment, we compare the so-called -hopping technique to the Hamiltonian replica exchange methods assessing their convergence behavior as a function of the enhanced sampling protocols (number of replicas, scaling factors, simulation times). We found that the pure -hopping, commonly used in solvation and binding free energy calculations via alchemical free energy perturbation techniques, is ineffective in enhancing the sampling of the isomeric states, exhibiting a pathological dependence on the initial conditions. Correct sampling can be restored in -hopping simulation by the addition of a “hot-zone” scaling factor to the -stratification (FEP+ approach), provided that the additive hot-zone scaling factors are tuned and optimized using preliminary ordinary replica-exchange simulation of the end-states.

Molecular dynamic insight into solubility of H2S in ionic liquids [emim][BF4], [emim][OTf] and [emim][Tf2N

Natural gas sweetening is one of the important processes in petroleum industry. Ionic liquids (ILs) as a new class of compounds are an alternative to conventional solvents to achieve this goal. In this study, the solubility of H2S gas in ionic liquids 1-Ethyl-3-methylimidazolium tetrafluoroborate ([emim][BF4]), 1-Ethyl-3-methylimidazolium trifluoromethanesulfonate ([emim][OTf]) and 1-Ethyl-3-methylimidazolium bis (trifluoromethylsulfonyl)imide was calculated based on the excess chemical potential energy by molecular dynamics free energy perturbation (FEP) technique. In this regard, validation was carried out for force fields of gas and ILs. Two force fields, All-Atom (A-A) and United-Atom (U-A), were compared to calculate the H2S solubility in the ILs. The performance of A-A force field was better for [emim][BF4] and [emim][OTf]; whereas, that of U-A force field for [emim][Tf2N]. The calculated Henry’s constant and chemical potential energy were in a good agreement with experimental data. The results showed that higher negative potential energy resulted in higher molar fraction and the ionic liquid [emim][BF4] had the best performance for H2S removal from a gas stream.

Insights on small molecule binding to the Hv1 proton channel from free energy calculations with molecular dynamics simulations

Hv1 is a voltage-gated proton channel whose main function is to facilitate extrusion of protons from the cell. The development of effective channel blockers for Hv1 can lead to new therapeutics for the treatment of maladies related to Hv1 dysfunction. Although the mechanism of proton permeation in Hv1 remains to be elucidated, a series of small molecules have been discovered to inhibit Hv1. Here, we computed relative binding free energies of a prototypical Hv1 blocker on a model of human Hv1 in an open state. We used alchemical free energy perturbation techniques based on atomistic molecular dynamics simulations. The results support our proposed open state model and shed light on the preferred tautomeric state of the channel blocker. This work lays the groundwork for future studies on adapting the blocker molecule for more effective inhibition of Hv1.

Toward Automated Free Energy Calculation with Accelerated Enveloping Distribution Sampling (A-EDS).

Free-energy perturbation (FEP) methods are commonly used in drug design to calculate relative binding free energies of different ligands to a common host protein. Alchemical ligand transformations are usually performed in multiple steps which need to be chosen carefully to ensure sufficient phase-space overlap between neighboring states. With one-step or single-step FEP techniques, a single reference state is designed that samples phase-space not only representative of a full transformation but also ideally resembles multiple ligand end states and hence allows for efficient multistate perturbations. Enveloping distribution sampling (EDS) is one example for such a method in which the reference state is created by a mathematical combination of the different ligand end states based on solid statistical mechanics. We have recently proposed a novel approach to EDS which enables efficient barrier crossing between the different end states, termed accelerated EDS (A-EDS). In this work, we further simplify the parametrization of the A-EDS reference state and demonstrate the automated calculation of multiple free-energy differences between different ligands from a single simulation in three different well-described drug design model systems.

Molecular Dynamics Study of the Hybridization between RNA and Modified Oligonucleotides.

MicroRNAs (miRNAs) are attractive drug candidates for many diseases as they can modulate the expression of gene networks. Recently, we discovered that DNAs targeting microRNA-22-3p (miR-22-3p) hold the potential for treating obesity and related metabolic disorders (type 2 diabetes mellitus, hyperlipidemia, and nonalcoholic fatty liver disease (NAFLD)) by turning fat-storing white adipocytes into fat-burning adipocytes. In this work, we explored the effects of chemical modifications, including phosphorothioate (PS), locked nucleic acid (LNA), and peptide nucleic acid (PNA), on the structure and energy of DNA analogs by using molecular dynamics (MD) simulations. To achieve a reliable prediction of the hybridization free energy, the AMOEBA polarizable force field and the free energy perturbation technique were employed. The calculated hybridization free energies are generally compatible with previous experiments. For LNA and PNA, the enhanced duplex stability can be explained by the preorganization mechanism, i.e., the single strands adopt stable helical structures similar to those in the duplex. For PS, the S and R isomers (Sp and Rp) have preferences for C2'-endo and C3'-endo sugar puckering conformations, respectively, and therefore Sp is less stable than Rp in DNA/RNA hybrids. In addition, the solvation penalty of Rp accounts for its destabilization effect. PS-LNA is similar to LNA as the sugar puckering is dominated by the locked sugar ring. This work demonstrated that MD simulations with polarizable force fields are useful for the understanding and design of modified nucleic acids.

A Streamlined, General Approach for Computing Ligand Binding Free Energies and Its Application to GPCR-Bound Cholesterol.

The theory of receptor-ligand binding equilibria has long been well-established in biochemistry, and was primarily constructed to describe dilute aqueous solutions. Accordingly, few computational approaches have been developed for making quantitative predictions of binding probabilities in environments other than dilute isotropic solution. Existing techniques, ranging from simple automated docking procedures to sophisticated thermodynamics-based methods, have been developed with soluble proteins in mind. Biologically and pharmacologically relevant protein-ligand interactions often occur in complex environments, including lamellar phases like membranes and crowded, nondilute solutions. Here, we revisit the theoretical bases of ligand binding equilibria, avoiding overly specific assumptions that are nearly always made when describing receptor-ligand binding. Building on this formalism, we extend the asymptotically exact Alchemical Free Energy Perturbation technique to quantifying occupancies of sites on proteins in a complex bulk, including phase-separated, anisotropic, or nondilute solutions, using a thermodynamically consistent and easily generalized approach that resolves several ambiguities of current frameworks. To incorporate the complex bulk without overcomplicating the overall thermodynamic cycle, we simplify the common approach for ligand restraints by using a single distance-from-bound-configuration (DBC) ligand restraint during AFEP decoupling from protein. DBC restraints should be generalizable to binding modes of most small molecules, even those with strong orientational dependence. We apply this approach to compute the likelihood that membrane cholesterol binds to known crystallographic sites on three GPCRs (β2-adrenergic, 5HT-2B, and μ-opioid) at a range of concentrations. Nonideality of cholesterol in a binary cholesterol:phosphatidylcholine (POPC) bilayer is characterized and consistently incorporated into the interpretation. We find that the three sites exhibit very different affinities for cholesterol: The site on the adrenergic receptor is predicted to be high affinity, with 50% occupancy for 1:109 CHOL:POPC mixtures. The sites on the 5HT-2B and μ-opioid receptor are predicted to be lower affinity, with 50% occupancy for 1:103 CHOL:POPC and 1:102 CHOL:POPC, respectively. These results could not have been predicted from the crystal structures alone.

Lanosterol Disrupts Aggregation of Human γD-Crystallin by Binding to the Hydrophobic Dimerization Interface.

Cataracts are a leading cause of vision impairment, which stem from the misfolding and aggregation of crystallins in the eye lens. Despite its prevalence and severity, the detailed mechanism by which misfolded crystallins aggregate into cataracts remains elusive. Recently, in vitro and in vivo experiments demonstrated that lanosterol, a steroid-type compound found in human and animal eyes, can not onlyprevent cataract formation but also reverse the formation. Inspired by these experimental observations, we investigate the preventive activity of lanosterol in the aggregate formation of human γD-crystallins (HγD-Crys) using all atom molecular dynamics (MD) simulation and free energy perturbation (FEP) techniques. Our results reveal that lanosterol preferentially binds to the HγD-Crys hydrophobic dimerization interface, in particular, to the structured C-terminus (near residues 135-165) with a stronger binding affinity than the unfolded N-terminus. Furthermore, we observe that the C-terminal binding is more favorable than lanosterol self-aggregation, further attesting to lanosterol's efficacy. Finally, we compare the binding free energy of lanosterol with cholesterol using alchemical transformation and discuss the possible correlation of the molecular geometry of steroids with binding affinity.

Relative Affinities of General Anesthetics for Pseudo-Symmetric Intersubunit Binding Sites of Heteromeric GABA(A) Receptors

GABA(A), a pentameric ligand gated ion channel is critical for regulating neuronal excitability. These inhibitory receptors, gated by ɣ-amino butyric acid (GABA), can be potentiated and also directly activated by intravenous and inhalational anesthetics. Although this receptor is a widely-studied target for general anesthetics, the mechanism of receptor modulation remains unclear. These receptors are predominantly found in 2α:2β:1γ stoichiometry, with four unique inter-subunit interfaces. Here we use thermodynamically rigorous free energy perturbation (AFEP) techniques and Molecular Dynamics simulations to rank the different intersubunit sites by affinity.

Free energy barrier for dissociation of the guanosine monophosphate anion in water

We report free energy barriers for the ground-state dissociation of the guanosine nucleotide anion in solution, employing implicit and explicit solvation models. The latter was based on the Free Energy Perturbation technique and Monte Carlo simulations. For the lowest-energy structure, both solvation models indicate a solvent-induced transition from a dipole-bound state in the gas phase to a π∗ valence state in solution. The free barrier estimates obtained from explicit and implicit solvation are in fair agreement with each other, although significantly overestimated in comparison to a previously reported explicit solvation model based on ab initio molecular dynamics simulations. Accounting for corrections related to the different DFT functionals used in the present and previous studies significantly improves the agreement.

Relative Affinities of Positive and Negative Modulators of Heteromeric GABA(A) Receptors for Pseudo-Symmetric Intersubunit Binding Sites

Type-A γ-aminobutyric acid receptors (GABAARs) are pentameric ligand-gated ion channels that are ubiquitous to the central nervous system (CNS) and critical for regulating neuronal excitability. These inhibitory receptors, gated by γ-aminobutyric acid (GABA), can be potentiated and also directly activated by certain exogenous and endogenous lipophilic small molecules. Among the various endogenous modulators are derivatives of cholesterol called neurosteroids; neurosteroids that positively modulate GABAARs can function as natural sedatives, anesthetics, anxiolytics and anti-convulsants. The mechanisms underlying modulation of GABAARs by neurosteroids and exogenous modulators like general anesthetics are still poorly understood, in part due to the location of likely quasi-symmetric interaction sites in the transmembrane domain of the receptor, as well as inconclusive results from mutagenesis experiments. In particular, the extent to which binding modes of positively and negatively modulating neurosteroids overlap is unknown. Here we use molecular dynamics simulations and the thermodynamically rigorous alchemical free energy perturbation technique to rank pseudo-symmetric intersubunit binding modes by affinity for the positively modulating neurosteroid allopregnanolone, the negatively modulating neurosteroid pregnenolone sulfate, and a novel analog with complex dose-response: triiodothyronine (thyroid hormone). Differential rankings of intersubunit sites for positive and negative modulating neurosteroids, triiodothyronine, and the general anesthetics propofol and sevoflurane, are interpreted in the context of their distinctive functional effects.

Solvation Thermodynamic Properties of Hydrogen Sulfide in [C4mim][PF6], [C4mim][BF4], and [C4mim][Cl] Ionic Liquids, Determined by Molecular Simulations.

Removal of hydrogen sulfide (H2S) and acid gases from natural gas is accomplished by absorption processes using a solvent. The gas solubility in a liquid can be used to measure the degree of removal of the gas and is quantified by the Henry's constant, the free energy of solvation at infinite dilution, or the excess chemical potential. In this work, Henry's constants and thermodynamic properties of solvation of H2S were calculated in three ionic liquids: [C4mim][PF6], [C4mim][BF4], and [C4mim][Cl] ([C4mim], 1-butyl-3-methyl imidazolium). The first step in this work was the evaluation of the force fields for the gas and condensed phases in order to obtain accurate values for the excess chemical potential for H2S on each ionic liquid using free energy perturbation techniques. In the H2S-[C4mim][PF6] and H2S-[C4mim][BF4] systems, the results obtained by molecular simulation agree with the experimental values reported in the literature. However, the solvation free energy calculated for the H2S-[C4mim][Cl] system can be considered predictive because of the lack of experimental data at the simulated conditions. Based on these results, the best solvent for removing H2S is [C4mim][Cl] because it has the highest affinity for this species (lowest value of the Henry's constant). Also, solvation thermodynamic properties such as enthalpy and entropy were calculated in order to evaluate their contribution to the free energy of solvation.

QM/MM modeling of the hydroxylation of the androstenedione substrate catalyzed by cytochrome P450 aromatase (CYP19A1).

CYP19A1 aromatase is a member of the Cytochrome P450 family of hemeproteins, and is the enzyme responsible for the final step of the androgens conversion into the corresponding estrogens, via a three-step oxidative process. For this reason, the inhibition of this enzyme plays an important role in the treatment of hormone-dependent breast cancer. The first catalytic subcycle, corresponding to the hydroxilation of androstenedione, has been proposed to occur through a first hydrogen abstraction and a subsequent oxygen rebound step. In present work, we have studied the mechanism of the first catalytic subcycle by means of hybrid quantum mechanics/molecular mechanics methods. The inclusion of the protein flexibility has been achieved by means of Free Energy Perturbation techniques, giving rise to a free energy of activation for the hydrogen abstraction step of 13.5 kcal/mol. The subsequent oxygen rebound step, characterized by a small free energy barrier (1.5 kcal/mol), leads to the hydroxylated products through a highly exergonic reaction. In addition, an analysis of the primary deuterium kinetic isotopic effects, calculated for the hydrogen abstraction step, reveals values (∼10) overpassing the semiclassical limit for the CH, indicating the presence of a substantial tunnel effect. Finally, a decomposition analysis of the interaction energy for the substrate and cofactor in the active site is also discussed. According to our results, the role of the enzymatic environment consists of a transition state stabilization by means of dispersive and polarization effects.

Qgui: A high-throughput interface for automated setup and analysis of free energy calculations and empirical valence bond simulations in biological systems

Structural information and activity data has increased rapidly for many protein targets during the last decades. In this paper, we present a high-throughput interface (Qgui) for automated free energy and empirical valence bond (EVB) calculations that use molecular dynamics (MD) simulations for conformational sampling. Applications to ligand binding using both the linear interaction energy (LIE) method and the free energy perturbation (FEP) technique are given using the estrogen receptor (ERα) as a model system. Examples of free energy profiles obtained using the EVB method for the rate-limiting step of the enzymatic reaction catalyzed by trypsin are also shown. In addition, we present calculation of high-precision Arrhenius plots to obtain the thermodynamic activation enthalpy and entropy with Qgui from running a large number of EVB simulations.

Efficient nonequilibrium method for binding free energy calculations in molecular dynamics simulations.

We introduce an effective technique for the calculation of the binding free energy in drug-receptor systems using nonequilibrium molecular dynamics and application of the Jarzynski theorem. In essence, this novel methodology constitutes the nonequilibrium adaptation of an ancient free energy perturbation technique called Double Annihilation Method, invented more than 25 years ago [J. Chem. Phys. 1988, 89, 3742-3746] and upon which modern approaches of binding free energy computation in drug-receptor systems are heavily based. The proposed computational approach, termed Fast Switching Double Annihilation Methods (FS-DAM) in honor of its ancient ancestor, is applied to a prototypical example system with multiple binding sites, proving its computational potential and versatility in unraveling multiple site or allosteric binding processes.

Mechanism for the decomposition of 5-aza-2′-deoxycytidine: a theoretical study using Monte Carlo simulation plus local microhydration model

The decomposition mechanism of 5-aza-2′-deoxycytidine (decitabine) has been investigated using the Monte Carlo simulation with free energy perturbation technique based on the local microhydration model. The irreversible deformylation process (process 2) was firstly evaluated. Four possible pathways, two stepwise (paths 1 and 3) and two concerted (paths 2 and 4) were examined. In two stepwise pathways, a tetrahedral intermediate is given by the addition of H2O to the C(1)=O(2) bond, followed by the proton shift to the N(3) atom leading to the final products. In two concerted pathways, the nucleophilic attack of H2O at the C(1) atom and the simultaneous movement of hydroxyl H to the N(3) atom induce the cleavage of the C(1)–N(3) single bond. The potential energy profiles along the minimum energy path for four pathways in the gas phase and in water were obtained. Our computed results exhibit that two stepwise mechanisms (paths 1 and 3) are more favorable than two concerted mechanisms (paths 2 and 4) when the solvent effects of bulk water are taken into account. Then, in collaboration with the previous study (Gao et al. in Theor Chem Acc 131:1108–1123, 2012), the whole decomposition processes of decitabine were discussed. Two possible pathways were considered: pathway t1 be made up of B-path 2 (in process 1) and path 1 as well as pathway t2 be composed of B-path 3 (in process 1) and path 3. The results clearly manifest that pathway t2 is more favorable than pathway t1 in the decomposition reaction of decitabine, which leads to IM-dR2 and GuaUre-dR2 as ring-opening-formylated intermediate and decomposition product, respectively. The third step is the rate-controlling step with the free energy barrier height of 29.8 kcal/mol, which is in agreement with the experimentally determined activation free energy (25.4 kcal/mol).

Irreversible Thermodynamics

We present a maximum-entropy theory of mesoscopic kinetics. The theory gives fully nonlinear nonequilibrium thermodynamic relationships and has no explicit requirement for either microscopic bath variables, an equilibrium energy, or an equilibrium partition function. The entropy maximization process is instead carried out over transition probability distributions with constraints on particle position and velocity updates. The Lagrange multipliers for these constraints express the instantaneous temperature and pressure of external (or microscopic) thermostatic driving systems, with which the distinguished system may or may not eventually reach equilibrium. We show that the analogues of the Gibbs-Maxwell relations and free energy perturbation techniques carry over to fluctuation-dissipation theorems and nonequilibrium ensemble reweighting techniques as should be expected. The result is a fully time-dependent, non-local description of a nonequilibrium ensemble coupled to reservoirs at possibly time-varying thermostatic or mechanical states. We also show that the thermodynamic entropy production extends the generalized fluctuation theorem through the addition of an instantaneous information entropy term for the end-points, leading to a concise statement of the second law of thermodynamics.

Theoretical studies on the chemical decomposition of 5-aza-2′-deoxycytidine: DFT study and Monte Carlo simulation

The decomposition mechanism of 5-Aza-2′-deoxycytidine has been studied by the use of computational techniques. Optimized structures for all of the stationary points in the gas phase were investigated at B3LYP/6-31+G(d,p) level of theory. Single-point energies were determined employing the ab initio MP2 method in conjunction with the 6-311++G(d,p) basis set. Five possible pathways, paths 1–5, were evaluated. In each pathway, the direct (A-paths 1–5) and water-assisted (B-paths 1–5) processes were considered. Meanwhile, the local microhydration model with the direct participation of three water molecules around the reaction centers was adopted to mimic the system for the water-assisted decomposition mechanisms above, where one water molecule is the nucleophilic reactant and the other two are the auxiliary molecules located on each side of the nucleophilic water. The results in the gas phase exhibit that the energy barriers of the water-assisted pathways based on the local microhydration model decrease dramatically by about 15–20 kcal/mol as compared with those of the direct pathways because of the contribution of the auxiliary water molecules. In addition, bulk solvent effects of water were determined by means of the self-consistent reaction field based on the conductor-like polarized continuum model and Monte Carlo simulation with free energy perturbation (MC-FEP) technique, respectively. Our computational results indicate that B-path 3 in the decomposition reaction of 5-azadC is the most favorable, where the calculated rate constant (1.68 × 10−3 min−1) using the MC-FEP method is within the range of the experimentally determined values [(5.89 ± 0.54) × 10−3 min−1 by UV and (1.46 ± 0.08) × 10−3 min−1 by NMR].

Biophysical Probes Reveal a “Compromise” Nature of the Methyl-lysine Binding Pocket in L3MBTL1

Histone lysine methylation (Kme) encodes essential information modulating many biological processes including gene expression and transcriptional regulation. However, the atomic-level recognition mechanisms of methylated histones by their respective adaptor proteins are still elusive. For instance, it is unclear how L3MBTL1, a methyl-lysine histone code reader, recognizes equally well both mono- and dimethyl marks but ignores unmodified and trimethylated lysine residues. We made use of molecular dynamics (MD) and free energy perturbation (FEP) techniques in order to investigate the energetics and dynamics of the methyl-lysine recognition. Isothermal titration calorimetry (ITC) was employed to experimentally validate the computational findings. Both computational and experimental methods were applied to a set of designed "biophysical" probes that mimic the shape of a single lysine residue and reproduce the binding affinities of cognate histone peptides. Our results suggest that, besides forming favorable interactions, the L3MBTL1 binding pocket energetically penalizes both methylation states and has most probably evolved as a "compromise" that nonoptimally fits to both mono- and dimethyl-lysine marks.