Free-Energy Profiles of Confined Reactions: Influence of Confinement Type and Challenges for Metadynamics Methods

A powerful tool to study the mechanism of reactions in solutions or enzymes is to perform the ab initio quantum mechanical/molecular mechanical (QM/MM) molecular dynamics (MD) simulations. However, the computational cost is too high due to the explicit electronic structure calculations at every time step of the simulation. A neural network (NN) method can accelerate the QM/MM-MD simulations, but it has long been a problem to accurately describe the QM/MM electrostatic coupling by NN in the electrostatic embedding (EE) scheme. In this work, we developed a new method to accelerate QM/MM calculations in the mechanic embedding (ME) scheme. The potentials and partial point charges of QM atoms are first learned in vacuo by the embedded atom neural networks (EANN) approach. MD simulations are then performed on this EANN/MM potential energy surface (PES) to obtain free energy (FE) profiles for reactions, in which the QM/MM electrostatic coupling is treated in the mechanic embedding (ME) scheme. Finally, a weighted thermodynamic perturbation (wTP) corrects the FE profiles in the ME scheme to the EE scheme. For two reactions in water and one in methanol, our simulations reproduced the B3LYP/MM free energy profiles within 0.5 kcal/mol with a speed-up of 30-60-fold. The results show that the strategy of combining EANN potential in the ME scheme with the wTP correction is efficient and reliable for chemical reaction simulations in liquid. Another advantage of our method is that the QM PES is independent of the MM subsystem, so it can be applied to various MM environments as demonstrated by an SN2 reaction studied in water and methanol individually, which used the same EANN PES. The free energy profiles are in excellent accordance with the results obtained from B3LYP/MM-MD simulations. In future, this method will be applied to the reactions of enzymes and their variants.

Determination of membrane permeability to small molecules from first-principles represents a promising approach for screening lead compounds according to their permeation properties upstream in the drug discovery process and prior to their synthesis. Theoretical investigation of permeation events requires, at its core, a molecular model of the membrane, and the choice of this model impacts not only the predicted permeability but also its relation to the experimental measurements commonly performed in pharmaceutical settings with a variety of cell lines capable of mimicking intestinal passive permeation. Homogeneous single-lipid bilayers have traditionally been utilized in computer simulations of membrane permeability predictions due to the ease of sampling all the relevant configurations, as well as the availability of parameters for a range of components of the biological membrane. To assess the influence of the membrane heterogeneity on the permeability to small molecules, we have examined the permeation of ethanol in six different single-lipid bilayers and compared the computed free-energy and diffusivity profiles with those obtained using a mammalian cell membrane model consisting of 26 components. Our results suggest that the membrane permeability only mildly depends on the lipid composition, spanning only 1 order of magnitude between the small phosphoethanolamine and the large phosphocholine head groups, or the short, saturated lauryl and the long, unsaturated oleyl acyl chains, that is, nearly as close as current theoretical estimates can get to experiment. The staggering computer time required to obtain an accurate free-energy profile, devoid of hysteresis between the upper and the lower leaflets of the lipid bilayer, in excess of several microseconds, provides an impetus for the development of approximate routes for membrane permeability predictions. Here, we have modeled the free-energy profile underlying permeation by means of a series of free-energy perturbation calculations, whereby the substrate is reversibly coupled to its environment at fixed values in the direction normal to the lipid bilayer. The diffusivity profile is modeled based on the bulk self-diffusion of the permeant, and the membrane permeability is recovered without significant loss of accuracy. The proposed numerical approach can be seamlessly extended to the determination of the relative membrane permeability to alternate substrates, thereby allowing large sets of permeants to be screened at a fraction of the computational cost of a rigorous determination of their respective free-energy profile.

A novel variational method for construction of free energy profiles from molecular simulation data is presented. The variational free energy profile (VFEP) method uses the maximum likelihood principle applied to the global free energy profile based on the entire set of simulation data (e.g from multiple biased simulations) that spans the free energy surface. The new method addresses common obstacles in two major problems usually observed in traditional methods for estimating free energy surfaces: the need for overlap in the re-weighting procedure and the problem of data representation. Test cases demonstrate that VFEP outperforms other methods in terms of the amount and sparsity of the data needed to construct the overall free energy profiles. For typical chemical reactions, only ~5 windows and ~20-35 independent data points per window are sufficient to obtain an overall qualitatively correct free energy profile with sampling errors an order of magnitude smaller than the free energy barrier. The proposed approach thus provides a feasible mechanism to quickly construct the global free energy profile and identify free energy barriers and basins in free energy simulations via a robust, variational procedure that determines an analytic representation of the free energy profile without the requirement of numerically unstable histograms or binning procedures. It can serve as a new framework for biased simulations and is suitable to be used together with other methods to tackle with the free energy estimation problem.

By employing mechanical work analogies, we derive a convenient computational approach for evaluation of the free energy profile (FEP) along some discretized path defined as a sequence of hyperplanes. A hyperplane is fully specified by any of its point and a tangent vector. The FEP is obtained as an integral of two components. The translational component of the free energy is computed by integrating the hyperplane constraint force. The rotational component is evaluated via the hyperplane torque. Both ingredients--the constraint force and the hyperplane torque-are evaluated on each hyperplane independently. The integration procedure utilizes a set of reference points defining a point of rotation on each hyperplane, and these points can be chosen before or after the sampling takes place. A shift in the reference points redistributes the FEP contributions between the translational and rotational components. For systems where the FEP is dominated by the potential energy differences, reference points residing on the minimum energy path present a natural choice. We demonstrate the validity of our approach on two examples, a simple two-dimensional (2D) potential, and a seven-atom Lennard-Jones cluster. In each case, we compare the numerical FEP with the harmonic approximation estimates. Our results for the 2D potential are also verified by the data available in the literature. In both cases, the rotational component of the FEP represents a sizable contribution to the total FEP, so ignoring it would yield clearly incorrect results.

The elucidation of the chemical mechanisms whereby biological molecules control, regulate and catalyze phosphoryl transfer reactions has profound implications for processes such as transcription, energy storage and transfer, cell signalling and gene regulation.[1, 2] The catalytic properties of RNA, in particular, have application in the design of new biotechnology and implications into the evolutionary origins of life itself.[3] Of primary importance to the understanding of mechanism is the characterization of the transition state for these reactions. Kinetic isotope effects (KIEs) offer one of the most powerful and sensitive experimental probes to interrogate the chemical environment of the transition state.[4-6] However, for complex reactions, theoretical methods are required to interpret the experimental measurements in terms of a detailed mechanistic model that traces the pathway from the reactant state through the transition state and into the product state.[7, 8] This paper presents experimental and computational results to characterize the mechanism of model phosphoryl transfer reactions that mimic RNA cleavage transesterification catalyzed by enzymes such as RNase A[1] as well as endonucleolytic ribozymes such as the hammerhead, hairpin, hepatitis delta virus (HDV), VS and glmS ribozymes.[9-11] Herein, secondary kinetic isotope effects are reported for the cleavage transesterification of a dinucleotide system, which, together with previously reported primary isotope effect measurements,[12, 13] represent a comprehensive characterization of isotope effects for a native (unmodified) RNA system. Scheme 1 illustrates the general mechanism for the reverse dianionic in-line methanolysis of ethylene phosphate, a model for base-catalysed RNA phosphate transesterification, with phosphoryl oxygen positions labelled in accord with their RNA counterparts. In this study, the free energy profiles for Scheme 1 were determined with density-functional quantum mechanical/molecular mechanical (QM/MM) simulations in explicit solvent.[14-17] These simulations are state-of-the-art, and take into account the dynamical fluctuations of the solute and solvent degrees of freedom in determination of the free energy profiles. In addition, adiabatic reaction energy profiles were determined with solvation effects treated implicitly with a polarizable continuum model (PCM)[18] specifically calibrated for the native and 3′ and 5′ thio-substituted compounds (Figure 1). The 3′ and 5′ thio-substituted compounds model the corresponding chemically modified RNAs that serve as valuable experimental probes of phosphoryl transfer mechanisms catalyzed by ribozymes.[19] The S5′ substitution, for example, in the HDV ribozyme serves as an enhanced leaving group that suppresses the deleterious effect of mutation of a critical cytosine residue, which has been interpreted to support its role as a general acid catalyst.[20] The energy values for stationary points of the native and thio-substituted reactions are in Table 1. Using our recently developed ab initio path-integral method based on Kleinert’s variational perturbation theory,[7, 21-23] we also calculated kinetic isotope effects, which are shown along with the most relevant experimental values for comparison in Table 2. The agreement that is achieved between the theoretical and experimental results allows a detailed mechanistic interpretation based on the theoretical models.[7, 8] Figure 1 (Color online) Comparison of density-functional QM/MM free-energy and adiabatic PCM profiles for the native reaction (top), and density-functional adiabatic PCM profiles for native and thio-substituted reactions (bottom) as a function of the difference ... Scheme 1 General reaction scheme for the (associative) reverse of dianionic in-line methanolysis of ethylene phosphate: a model for RNA phosphate transesterification under alkaline conditions. In the present work, the native reaction shown in the scheme is studied ... Table 1 Relative free energy (kcal/mol) and reaction coordinate (Δbond) values (A) computed for stationary points along the coordinate of the native and thio-substituted RNA phosphate transesterification reaction models shown in Scheme 1.[[a] ... Table 2 Primary kinetic isotope effects (KIEs) on 2′ nucleophile (18kNu) and 5′ leaving (18kLg) oxygens, and secondary KIEs on O1P (18kO1P) and O2P (18kO2P) oxygens in aqueous solution for the TS1 and TS2 transition states, along with the most ... All of the profiles calculated in this work correspond to associative (or concerted) mechanisms characterized by initial nucleophilic attack, as is typical of phosphate diesters.[6] The departure of the leaving group can be concerted with nucleophilic attack (as in the S5′ substituted reaction) or can occur in a stepwise fashion resulting in the formation of a stable pentavalent phosphorane intermediate. In the stepwise mechanism, two transition states occur: one in which nucleophilic attack occurs (TS1) and another one in which leaving group departure takes place (TS2) as indicated in Scheme 1. These transition states themselves can be characterized as either “early” or “late”, depending on the degree of P-O2′ and P-O5′ bond formation/cleavage.

The concept of a protein diffusing in its free-energy folding landscape has been fruitful for both theory and experiment. Yet the choice of the reaction coordinate (RC) introduces an undesirable degree of arbitrariness into the problem. We analyze extensive simulation data of an alpha-helix in explicit water solvent as it stochastically folds and unfolds. The free-energy profiles for different RCs exhibit significant variations, some having an activation barrier, while others not. We show that this variation has little effect on the predicted folding kinetics if the diffusivity profiles are properly taken into account. This kinetic quasi-universality is rationalized by an RC rescaling, which, due to the reparameterization invariance of the Fokker-Planck equation, allows the combination of free-energy and diffusivity effects into a single function, the rescaled free-energy profile. This rescaled free energy indeed shows less variation among different RCs than the bare free energy and diffusivity profiles separately do, if we properly distinguish between RCs that contain knowledge of the native state and those that are purely geometric in nature. Our method for extracting diffusivity profiles is easily applied to experimental single molecule time series data and might help to reconcile conflicts that arise when comparing results from different experimental probes for the same protein.

The problem of defining efficient strategies for partitioning the cavity surface in QM solvation procedures based on boundary elements methods is considered here. The GEPOL procedure to get the cavity surface, and its partition into tesserae is adopted as a starting point: a version with variable tesselation is presented. The procedure to build the new sphere tesselations is described and several different options to select the surface partition have been implemented. The effects of the variation of the surface partition on the free energy of solvation of several solutes are also presented. Two free energy of solvation profiles evaluated with several different cavity partitions are analysed. We find that a radius-driven tesselation for every sphere reduces the number and extension of the cavity artefacts.

Membrane permeability of drug molecules plays a significant role in the development of new therapeutic agents. Accordingly, methods to predict the passive permeability of drug candidates during a medicinal chemistry campaign offer the potential to accelerate the drug design process. In this work, we combine the physics-based site identification by ligand competitive saturation (SILCS) method and data-driven artificial intelligence (AI) to create a high-throughput predictive model for the passive permeability of druglike molecules. In this study, we present a comparative analysis of four regression models to predict membrane permeabilities of small druglike molecules; of the tested models, Random Forest was the most predictive yielding an R2 of 0.81 for the independent data set. The input feature vector used to train the developed prediction model includes absolute free energy profiles of ligands through a POPC-cholesterol bilayer based on ligand grid free energy (LGFE) profiles obtained from the SILCS approach. The use of the membrane free energy profiles from SILCS offers information on the physical forces contributing to ligand permeability, while the use of AI yields a more predictive model trained on experimental PAMPA permeability data for a collection of 229 molecules. This combination allows for rapid estimations of ligand permeability at a level of accuracy beyond currently available predictive models while offering insights into the contributions of the functional groups in the ligands to the permeability barrier, thereby offering quantitative information to facilitate rational ligand design.

Free energy profiles form the cornerstone in the study of protein folding and function. In this study, the free energy profile of SUMO1 protein is directly reconstructed using an extension of the Jarzynski equality from atomic force microscope (AFM) based single-molecule force spectroscopy (SMFS) experiments. SUMO1 is a ubiquitin-like posttranslational modifier protein having a β clamp motif in its structure, imparting it with mechanical stability. We use the Jarzynski equality to obtain the equilibrium free energy profile from repeated nonequilibrium single-molecule pulling experiments. Indeed, the free energy values determined by the Jarzynski equality are lesser than the normal work average at all extensions. The free energy profiles constructed for the two velocities (100 and 400 nm/s) overlap with each other. The unfolding free energy barrier is estimated to be ∼7.5 kcal/mol. We anticipate that the Jarzynski equality can be applied in a similar manner to other ubiquitin-like proteins to extract their differences in the free energy profile, and hence, the effect of sequence diversity of structurally homologous proteins on the free energy landscape can be studied.

Density functional theory (DFT) together with Car-Parrinello ab initio molecular dynamics (CP-AIMD) simulation has been used to investigate the free energy profiles of two representative SN2 reactions: (A) Cl- + CH3Cl → ClCH3 + Cl-; (B) NH3 + H3BNH3 → H3NBH3 + NH3. The free energy profiles along the reaction coordinates at 300 K and 600 K were determined directly by a pointwise thermodynamic integration (PTI) technique. Comparison between the well-known double-well potential energy profile (PEP) and the free energy profiles (FEP) has been made. The results show that, for reaction A, the double-well profile is maintained for the FEP at 300 K due to the stronger ion−dipole interaction between chloromethane and the chloride anion. In comparison with the PEP, the FEP has a higher central barrier and a more shallow well depth. However, at 600 K the double wells almost disappear on the FEP, whereas the central barrier increases further. For reaction B, the 300 K FEP also presents a higher central barrier peak and a more shallow well depth compared to the PEP. However, when the temperature increases to 600 K, a saddle shape FEP is obtained, which indicates that the reaction has changed mechanism from an associative SN2 reaction to a dissociative SN1 reaction. This change is driven by entropy.

First-principles determination of free energy profiles for condensed-phase chemical reactions is hampered by the daunting costs associated with configurational sampling on ab initio quantum mechanical/molecular mechanical (AI/MM) potential energy surfaces. Here, we report a new method that enables efficient AI/MM free energy simulations through mean force fitting. In this method, a free energy path in collective variables (CVs) is first determined on an efficient reactive aiding potential. Based on the configurations sampled along the free energy path, correcting forces to reproduce the AI/MM forces on the CVs are determined through force matching. The AI/MM free energy profile is then predicted from simulations on the aiding potential in conjunction with the correcting forces. Such cycles of correction-prediction are repeated until convergence is established. As the instantaneous forces on the CVs sampled in equilibrium ensembles along the free energy path are fitted, this procedure faithfully restores the target free energy profile by reproducing the free energy mean forces. Due to its close connection with the reaction path-force matching (RP-FM) framework recently introduced by us, we designate the new method as RP-FM in collective variables (RP-FM-CV). We demonstrate the effectiveness of this method on a type-II solution-phase SN2 reaction, NH3 + CH3Cl (the Menshutkin reaction), simulated with an explicit water solvent. To obtain the AI/MM free energy profiles, we employed the semiempirical AM1/MM Hamiltonian as the base level for determining the string minimum free energy pathway, along which the free energy mean forces are fitted to various target AI/MM levels using the Hartree-Fock (HF) theory, density functional theory (DFT), and the second-order Møller-Plesset perturbation (MP2) theory as the AI method. The forces on the bond-breaking and bond-forming CVs at both the base and target levels are obtained by force transformation from Cartesian to redundant internal coordinates under the Wilson B-matrix formalism, where the linearized FM is facilitated by the use of spline functions. For the Menshutkin reaction tested, our FM treatment greatly reduces the deviations on the CV forces, originally in the range of 12-33 to ∼2 kcal/mol/Å. Comparisons with the experimental and benchmark AI/MM results, tests of the new method under a variety of simulation protocols, and analyses of the solute-solvent radial distribution functions suggest that RP-FM-CV can be used as an efficient, accurate, and robust method for simulating solution-phase chemical reactions.

Quantitative predictions of biomembrane/water partition coefficients are important, as they are a key property in pharmaceutical applications and toxicological studies. Molecular dynamics (MD) simulations are used to calculate free energy profiles for different solutes in lipid bilayers. How to calculate partition coefficients from these profiles is discussed in detail and different definitions of partition coefficients are compared. Importantly, it is shown that the calculated coefficients are in quantitative agreement with experimental results. Furthermore, we compare free energy profiles from MD simulations to profiles obtained by the recent method COSMOmic, which is an extension of the conductor-like screening model for realistic solvation to micelles and biomembranes. The free energy profiles from these molecular methods are in good agreement. Additionally, solute orientations calculated with MD and COSMOmic are compared and again a good agreement is found. Four different solutes are investigated in detail: 4-ethylphenol, propanol, 5-phenylvaleric acid, and dibenz[a,h]anthracene, whereby the latter belongs to the class of polycyclic aromatic hydrocarbons. The convergence of the free energy profiles from biased MD simulations is discussed and the results are shown to be comparable to equilibrium MD simulations. For 5-phenylvaleric acid the influence of the carboxyl group dihedral angle on free energy profiles is analyzed with MD simulations.

The free energy of pore formation in lipid bilayers has been previously calculated using a variety of reaction coordinates. Here, we use free energy perturbation of a cylindrical lipid exclusion restraint to compute the free energy profile as a function of pore radius in dimyristoylphosphatidylcholine (DMPC) and dioleoylphosphatidylcholine (DOPC) bilayers. Additionally restraining the headgroups to lie on the membrane surface allows us to also calculate the free energy profile of hydrophobic pores, i.e., cylindrical pores lined by acyl chains. For certain pore radii, the free energy of wetting of hydrophobic pores is calculated using the density bias method. It is found that wetting of hydrophobic pores becomes thermodynamically favorable at 5.0 Å for DMPC and 6.5 Å for DOPC, although significant barriers prevent spontaneous wetting of the latter on a nanosecond time scale. The free energy of transformation of hydrophilic pores to hydrophobic ones is also calculated using free energy perturbation of headgroup restraints along the bilayer normal. This quantity, along with wetting and pore growth free energies, provides complete free energy profiles as a function of radius. Pore line tension values for the hydrophilic pores obtained from the slope of the free energy profiles are 37.6 pN for DMPC and 53.7 pN for DOPC. The free energy profiles for the hydrophobic pores are analyzed in terms of elementary interfacial tensions. It is found that a positive three-phase line tension is required to explain the results. The estimated value for this three-phase line tension (51.2 pN) lies within the expected range.

Liposomes and micelles find various applications as potential solubilizers in extraction processes or in drug delivery systems. Thermodynamic and transport processes governing the interactions of different kinds of solutes in liposomes or micelles can be analyzed regarding the free energy profiles of the solutes in the system. However, free energy profiles in heterogeneous systems such as micelles are experimentally almost not accessible. Therefore, the development of predictive methods is desirable. Molecular dynamics (MD) simulations reliably simulate the structure and dynamics of lipid membranes and micelles, whereas COSMO-RS accurately reproduces solvation free energies in different solvents. For the first time, free energy profiles in micellar systems, as well as mixed lipid bilayers, are investigated, taking advantage of both methods: MD simulations and COSMO-RS, referred to as COSMOmic (Klamt, A.; Huniar, U.; Spycher, S.; Keldenich, J. COSMOmic: A Mechanistic Approach to the Calculation of Membrane-Water Partition Coefficients and Internal Distributions within Membranes and Micelles. J. Phys. Chem. B 2008, 112, 12148-12157). All-atom molecular dynamics simulations of the system SDS/water and CTAB/water have been applied in order to retrieve representative micelle structures for further analysis with COSMOmic. For the system CTAB/water, different surfactant concentrations were considered, which results in different micelle sizes. Free energy profiles of more than 200 solutes were predicted and validated by means of experimental partition coefficients. To our knowledge, these are the first quantitative predictions of micelle/water partition coefficients, which are based on whole free energy profiles from molecular methods. Further, the partitioning in lipid bilayer systems containing different hydrophobic tail groups (DOPC (1,2-dioleoyl-sn-glycero-3-phosphocholine), SOPC (stearoyl-oleoylphosphatidylcholine), DMPC (1,2-dimyristoyl-sn-glycero-3-phosphocholine), and POPC (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine)) as well as mixed bilayers was calculated. Experimental partition coefficients (log P) were reproduced with a root-mean-square error (RMSE) of 0.62. To determine the influence of cholesterol as an important component of cellular membranes, free energy profiles in the presence of cholesterol were calculated and shown to be in good agreement with experimental data.

We develop a new Deep Potential - Range Correction (DPRc) machine learning potential for combined quantum mechanical/molecular mechanical (QM/MM) simulations of chemical reactions in the condensed phase. The new range correction enables short-ranged QM/MM interactions to be tuned for higher accuracy, and the correction smoothly vanishes within a specified cutoff. We further develop an active learning procedure for robust neural network training. We test the DPRc model and training procedure against a series of 6 non-enzymatic phosphoryl transfer reactions in solution that are important in mechanistic studies of RNA-cleaving enzymes. Specifically, we apply DPRc corrections to a base QM model and test its ability to reproduce free energy profiles generated from a target QM model. We perform comparisons using the MNDO/d and DFTB2 semiempirical models because they produce free energy profiles which differ significantly from each other, thereby providing us a rigorous stress test for the DPRc model and training procedure. The comparisons show that accurate reproduction of the free energy profiles requires correction of the QM/MM interactions out to 6 Å. We further find that the model's initial training benefits from generating data from temperature replica exchange simulations and including high-temperature configurations into the fitting procedure so the resulting models are trained to properly avoid high-energy regions. A single DPRc model was trained to reproduce 4 different reactions and yielded good agreement with the free energy profiles made from the target QM/MM simulations. The DPRc model was further demonstrated to be transferable to 2D free energy surfaces and 1D free energy profiles that were not explicitly considered in the training. Examination of the computational performance of the DPRc model showed that it was fairly slow when run on CPUs, but was sped up almost 100-fold when using an NVIDIA V100 GPUs, resulting in almost negligible overhead. The new DPRc model and training procedure provide a potentially powerful new tool for the creation of next-generation QM/MM potentials for a wide spectrum of free energy applications ranging from drug discovery to enzyme design.