Weighted Histogram Analysis Method Research Articles

Investigation of the Driving Force of Replacing Adsorbed Hydrocarbons by CO2 in Organic Matter from an Energy Perspective.

Carbon dioxide (CO2) has been widely used to enhance the recovery of adsorbed hydrocarbons from the organic matter (OM) in shale formations. To reveal the driving force of replacing adsorbed hydrocarbons from OM by CO2, we performed molecular dynamics (MD) simulations of the replacement process and calculated the interaction forces between CO2 and hydrocarbons. In addition, based on the umbrella sampling method, steered MD simulations were performed, and the free energy profiles of hydrocarbons were obtained using the weighted histogram analysis method. Results show that the condition of the hydrocarbon replacement by CO2 requires the hydrocarbon to have sufficient kinetic energy or to have a sufficiently large attractive force exerted to ensure that the hydrocarbon escapes the potential well of the OM. The attractive forces exerted on hydrocarbon molecules by CO2 can significantly decrease the energy barrier associated with hydrocarbon movement away from the OM surface. Furthermore, both CO2 and supercritical CO2 can effectively displace adsorbed hydrocarbon gas (methane) on the OM, while supercritical CO2 is required to enhance the recovery of adsorbed hydrocarbon oil (n-dodecane). The results obtained in this study provide guidance for enhancing the recovery of adsorbed hydrocarbons by CO2 in shale formations.

OGRe: Optimal Grid Refinement Protocol for Accurate Free Energy Surfaces and Its Application in Proton Hopping in Zeolites and 2D COF Stacking.

While free energy surfaces are the crux of our understanding of many chemical and biological processes, their accuracy is generally unknown. Moreover, many developments to improve their accuracy are often complicated, limiting their general use. Luckily, several tools and guidelines are already in place to identify these shortcomings, but they are typically lacking in flexibility or fail to systematically determine how to improve the accuracy of the free energy calculation. To overcome these limitations, this work introduces OGRe, a Python package for optimal grid refinement in an arbitrary number of dimensions. OGRe is based on three metrics that gauge the confinement, consistency, and overlap of each simulation in a series of umbrella sampling (US) simulations, an enhanced sampling technique ubiquitously adopted to construct free energy surfaces for hindered processes. As these three metrics are fundamentally linked to the accuracy of the weighted histogram analysis method adopted to generate free energy surfaces from US simulations, they facilitate the systematic construction of accurate free energy profiles, where each metric is driven by a specific umbrella parameter. This allows for the derivation of a consistent and optimal collection of umbrellas for each simulation, largely independent of the initial values, thereby dramatically increasing the ease-of-use toward accurate free energy surfaces. As such, OGRe is particularly suited to determine complex free energy surfaces with large activation barriers and shallow minima, which underpin many physical and chemical transformations and hence to further our fundamental understanding of these processes.

Application of the weighted histogram method for calculating the thermodynamic parameters of the formation of oligodeoxyribonucleotide duplexes

To date, many derivatives and analogs of nucleic acids (NAs) have been developed. Some of them have found uses in scientific research and biomedical applications. Their effective use is based on the data about their properties. Some of the most important physicochemical properties of oligonucleotides are thermodynamic parameters of the formation of their duplexes with DNA and RNA. These parameters can be calculated only for a few NA derivatives: locked NAs, bridged oligonucleotides, and peptide NAs. Existing predictive approaches are based on an analysis of experimental data and the consequent construction of predictive models. The ongoing pilot studies aimed at devising methods for predicting the properties of NAs by computational modeling techniques are based only on knowledge about the structure of oligonucleotides. In this work, we studied the applicability of the weighted histogram analysis method (WHAM) in combination with umbrella sampling to the calculation of thermodynamic parameters of DNA duplex formation (changes in enthalpy ∆H°, entropy ∆S°, and Gibbs free energy ∆G37° ). A procedure was designed involving WHAM for calculating the hybridization properties of oligodeoxyribonucleotides. Optimal parameters for modeling and calculation of thermodynamic parameters were determined. The feasibility of calculation of ∆H°, ∆S°, and ∆G37° was demonstrated using a representative sample of 21 oligonucleotides 4–16 nucleotides long with a GC content of 14–100 %. Error of the calculation of the thermodynamic parameters was 11.4, 12.9, and 11.8 % for ∆H°, ∆S°, and ∆G37° , respectively, and the melting temperature was predicted with an average error of 5.5 °C. Such high accuracy of computations is comparable with the accuracy of the experimental approach and of other methods for calculating the energy of NA duplex formation. In this paper, the use of WHAM for computation of the energy of DNA duplex formation was systematically investigated for the first time. Our results show that a reliable calculation of the hybridization parameters of new NA derivatives is possible, including derivatives not yet synthesized. This work opens up new horizons for a rational design of constructs based on NAs for solving problems in biomedicine and biotechnology.

Effective Interactions between Double-Stranded DNA Molecules in Aqueous Electrolyte Solutions: Effects of Molecular Architecture and Counterion Valency.

A computational investigation of the effects of molecular topology, namely, linear and circular, as well as counterion valency, on the ensuing pairwise effective interactions between DNA molecules in an unlinked state is presented. Umbrella sampling simulations have been performed through the introduction of bias potential along a reaction coordinate defined as the distance between the centers-of-mass of pairs of DNA molecules, and effective pair interaction potentials have been computed by employing the weighted histogram analysis method. An interesting comparison can be drawn between the different DNA topologies studied here, especially with regard to the contrasting effects of divalent counterions on the effective pair potentials: while DNA-DNA repulsion in short center-of-mass distances decreases significantly in the presence of divalent counterion-ions (as compared to monovalent ions) for linear DNA, the opposite effect occurs for the DNA minicircles. This can be attributed to the fact that linear DNA fragments can easily adopt relative orientations that minimize electrostatic and steric repulsions by rotating relative to one another and by exhibiting more pronounced bending due to the presence of free ends.

A New Approach for Estimating the Free Energy Differences among Multiple Thermodynamic States in Statistical Simulations.

In this letter, a new approach to compute free energy differences (FEDs) between multiple thermodynamics states is introduced. The method directly uses energy probability densities, which can be extracted with high accuracy from equilibrium simulations to obtain FEDs. Methods in current use, such as Bennett acceptance ratio (BAR), its multistate generalization (MBAR), or the weighted histogram analysis method (WHAM), require iterative solution of nonlinear equations which are known to be slowly convergent. The equations providing MBAR FEDs are identical to those derived earlier by Souaille and Roux in a method that has become known informally as "binless WHAM". In contrast, we obtain FEDs by solution of linear equations. For the classic two-state problem, the statistical error of our method, solving linear equations, is shown analytically to match that of BAR under common conditions.

Dynamics of polymer chains confined to a periodic cylinder: molecular dynamics simulation vs. Lifson-Jackson formula.

The dynamics of polymer chains confined to a periodic cylinder is explored using molecular dynamics simulation and theoretical analysis. The cylinder is divided into two cavities in one periodicity: one cavity consists of a channel of length L1 and diameter D1 and another cavity is a channel of length L2 and diameter D2. For L1 = L2 = L/2, the diffusion coefficient D of a single confined polymer chain decreases rapidly with increase in periodicities L. For a fixed periodicity with L = L1 + L2 = constant, the diffusion coefficient D of a single confined polymer chain shows strong dependence on L1 (or L2). Moreover, for a multi-chain system with L1 = L2, the diffusion coefficient D shows strong non-monotonic dependence on the chain monomer density ρ, and the confined polymer chains diffuse fastest for ρ = 0.068, in which there are three polymer chains in two periodicities as well as two opposing effects: one is that the excluded volume effect between polymer chains can reduce the free energy barrier, and another is that when the chain monomer density ρ increases further, the entanglement effect increases, which leads to that the diffusion coefficient D decreases as ρ increases. Finally, we found that the diffusion coefficient D has a similar oscillation relationship with the ratio of R/L for different chain lengths N and different periodicity L, and the oscillation amplitude decreases gradually as R/L increases; here R is the mean end-to-end distance of a single confined polymer chain, i.e., . From the view of free energy potential, both the width of the free energy potential well and the height of the free energy potential barrier govern simultaneously the diffusion behavior of confined polymer chains. According to the mean force potential (PMF) based on the weighted histogram analysis method (WHAM), we found that our results agree very well with the theoretical analysis using the Lifson-Jackson formula. Our investigation may help us understand the dynamics of particles in a periodic medium, which is one of the interesting problems in many different fields of science, such as physics, chemistry and biology.

Absolute binding free energies of the antiviral peptide ATN-161 with protein targets of SARS-CoV-2

The interactions of the antiviral pentapeptide ATN-161 with the closed and open conformations of the α5β1 integrin, the SARS-CoV-2 major protease, and the omicron variant spike protein complexed with hACE2 were studied using molecular docking and molecular dynamics simulation. Molecular docking was performed to obtain ATN-161 binding poses with these studied protein targets. Subsequently, molecular dynamics simulations were performed to verify the ligand stability at the binding site of each protein target. Pulling simulations, umbrella sampling, and weighted histogram analysis method were used to obtain the potential of mean force of each system and calculate the Gibbs free energy of binding for the ATN-161 peptide in each binding site of these protein targets. The results showed that ATN-161 binds to α5β1 integrin in its active and inactive form, binds weakly to the omicron variant spike protein complexed with hACE2, and strongly binds to the main protease target. Communicated by Ramaswamy H. Sarma

Computational Investigation on the Efficiency of Small Molecule Inhibitors Identified from Indian Spices against SARS-CoV-2 Mpro

Recently, small compounds from Indian spices (Carnosol, Arjunglucoside-I, and Rosmanol) have been identified as SARS-CoV-2 main protease (Mpro) inhibitors. The structural dynamics and characteristic features of binding of these small molecules to the SARS-CoV-2 Mpro are not well understood. Here, we have constructed the potential of mean force (PMF) for dissociating Mpro-small molecule inhibitor complexes from the umbrella sampling simulations using the weighted histogram analysis method. Mpro-small molecule inhibitor complexes exhibited relatively higher dissociation energy values than the alpha-ketoamide-Mpro complex (positive control) from the PMF calculations. We found that binding affinity between protein and ligand is higher in Mpro-Arjunglucoside-I complex [∆Gbind = -19.74 kcal mol-1 from MM-GBSA and ∆Gbind = -9.13 kcal mol-1 from MM-PBSA] than in other three SARS-CoV-2 small molecule complexes. The MM-GBSA/MM-PBSA calculations revealed that the small molecule inhibitors studied in this work have substantially higher binding affinity for Mpro. We found the residues present in SARS-CoV-2 Mpro's binding pocket contributed the most binding free energy to SARS-CoV-2 Mpro-small molecule interactions. Our findings emphasize the structural and binding features of the identified small molecule inhibitors with SARS-CoV-2 Mpro, which could be relevant in developing therapeutic candidates to combat SARS-CoV-2.

Modes of adhesion of two Janus nanoparticles on the outer or inner side of lipid vesicles.

Using molecular dynamics simulations of a coarse-grained model, in conjunction with the weighted histogram analysis method, the adhesion modes of two spherical Janus nanoparticles (NPs) on the outer or inner side of lipid vesicles are explored. In particular, the effects of the area fraction, J, of the NPs that interact attractively with lipid head groups, the adhesion strength and the size of the NPs on their adhesion modes are investigated. The NPs are found to exhibit two main modes of adhesion when adhered to the outer side of the vesicle. In the first mode, which occurs at relatively low values of J, the NPs are apart from each other. In the second mode, which occurs at higher values of J, the NPs form an in-plane dimer. Janus NPs, which adhere to the inner side of the vesicle, are always found to be apart from each other, regardless of the value of J and their diameter.

Umbrella Sampling-Based Method to Compute Ligand-Binding Affinity.

Many proteins have a solvent-exposed binding cleft, which permits their inhibitors to bind and unbind without significant protein conformation transforms. The binding/unbinding pathways of these protein-inhibitor complexes can be rather straightforwardly sampled by using umbrella sampling (US) simulation methods. During a US simulation, the Cα atoms of the protein are restrained via a harmonic force. The potential of mean force (PMF) along the binding pathway can be estimated by using the weighted histogram analysis method (WHAM). The binding affinity is then computed as the difference in PMF between the binding and unbinding states.

Empirical optimization of molecular simulation force fields by Bayesian inference

The demands on the accuracy of force fields for classical molecular dynamics simulations are steadily growing as larger and more complex systems are studied over longer times. One way to meet these growing demands is to hand over the learning of force fields and their parameters to machines in a systematic (semi)automatic manner. Doing so, we can take full advantage of exascale computing, the increasing availability of experimental data, and advances in quantum mechanical computations and the calculation of experimental observables from molecular ensembles. Here, we discuss and illustrate the challenges one faces in this endeavor and explore a way forward by adapting the Bayesian inference of ensembles (BioEn) method [Hummer and Köfinger, J. Chem. Phys. (2015)] for force field parameterization. In the Bayesian inference of force fields (BioFF) method developed here, the optimization problem is regularized by a simplified prior on the force field parameters and an entropic prior acting on the ensemble. The latter compensates for the unavoidable over simplifications in the parameter prior. We determine optimal force field parameters using an iterative predictor–corrector approach, in which we run simulations, determine the reference ensemble using the weighted histogram analysis method (WHAM), and update the force field according to the BioFF posterior. We illustrate this approach for a simple polymer model, using the distance between two labeled sites as the experimental observable. By systematically resolving force field issues, instead of just reweighting a structural ensemble, the BioFF corrections extend to observables not included in ensemble reweighting. We envision future force field optimization as a formalized, systematic, and (semi)automatic machine-learning effort that incorporates a wide range of data from experiment and high-level quantum chemical calculations, and takes advantage of exascale computing resources.Graphic abstract

FSATOOL 2.0: An integrated molecular dynamics simulation and trajectory data analysis program.

Molecular dynamics simulation is important in the computational study of the biomolecules. In this paper, we upgrade our previous FSATOOL to version 2.0. It is no longer a plugin as before. Besides the existed enhanced sampling and Markov state model analysis module, FSATOOL 2.0 has three new features now. First, it contains a molecular dynamics simulation engine on both CPU and GPU device. The engine works with an embedded enhanced sampling module. Second, it can do the free energy calculation by various practical methods, including the weighted histogram analysis method and Gaussian mixture model. Third, it has many subroutines to process the trajectory data, such as principal component analysis, time-structure based independent component analysis, contact analysis, and Φ-value analysis. Most importantly, all these calculations are integrated into one package. The trajectory data format is compatible with all the modules. With a proper input parameter file, users can do the molecular dynamics simulation and data analysis work by only a few simplified commands. The capabilities and theoretical backgrounds of FSATOOL 2.0 are introduced in the paper.

Boosting the conformational sampling by combining replica exchange with solute tempering and well-sliced metadynamics.

Methods that combine collective variable (CV) based enhanced sampling and global tempering approaches are used in speeding-up the conformational sampling and free energy calculation of large and soft systems with a plethora of energy minima. In this paper, a new method of this kind is proposed in which the well-sliced metadynamics approach (WSMTD) is united with replica exchange with solute tempering (REST2) method. WSMTD employs a divide-and-conquer strategy wherein high-dimensional slices of a free energy surface are independently sampled and combined. The method enables one to accomplish a controlled exploration of the CV-space with a restraining bias as in umbrella sampling, and enhance-sampling of one or more orthogonal CVs using a metadynamics like bias. The new hybrid method proposed here enables boosting the sampling of more slow degrees of freedom in WSMTD simulations, without the need to specify associated CVs, through a replica exchange scheme within the framework of REST2. The high-dimensional slices of the probability distributions of CVs computed from the united WSMTD and REST2 simulations are subsequently combined using the weighted histogram analysis method to obtain the free energy surface. We show that the new method proposed here is accurate, improves the conformational sampling, and achieves quick convergence in free energy estimates. We demonstrate this by computing the conformational free energy landscapes of solvated alanine tripeptide and Trp-cage mini protein in explicit water.

Mean force based temperature accelerated sliced sampling: Efficient reconstruction of high dimensional free energy landscapes.

Temperature accelerated sliced sampling (TASS) is an efficient method to compute high dimensional free energy landscapes. The original TASS method employs the weighted histogram analysis method (WHAM) which is an iterative post-processing to reweight and stitch high dimensional probability distributions in sliced windows that are obtained in the presence of restraining biases. The WHAM necessitates that TASS windows lie close to each other for proper overlap of distributions and span the collective variable space of interest. On the other hand, increase in number of TASS windows implies more number of simulations, and thus it affects the efficiency of the method. To overcome this problem, we propose herein a new mean-force (MF) based reweighting scheme called TASS-MF, which enables accurate computation with a fewer number of windows devoid of the WHAM post-processing. Application of the technique is demonstrated for alanine di- and tripeptides in vacuo to compute their two- and four-dimensional free energy landscapes, the latter of which is formidable in conventional umbrella sampling and metadynamics. The landscapes are computed within a kcal mol-1 accuracy, ensuring a safe usage for broad applications in computational chemistry.

Spatial arrangements of spherical nanoparticles on lipid vesicles.

We report results of a numerical investigation of the modes of adhesion of two spherical nanoparticles (NPs) on lipid vesicles based on molecular dynamics simulations, in conjunction with the weighted histogram analysis method, of an implicit-solvent model of self-assembled membranes. Our investigation shows that the NPs exhibit a sequence of three modes of adhesion. For low adhesive interactions, the adhering NPs are apart from each other. As the adhesive interaction is increased, the NPs dimerize into in-plane dimers. As the adhesive interaction is further increased for relatively large vesicles, the NPs dimerize into tubular dimers. However, for small vesicles, the tubular dimer state is not observed. For higher values of the adhesive interaction, four endocytosis modes are observed, depending on the initial locations of the NPs on the vesicle and the relative size of the NPs with respect to that of the vesicle. For relatively large vesicles, the NPs are endocytosed individually or as a dimer. For relatively small vesicles, only one NP is endocytosed if the initial distance between the NPs is large, while the second NP remains adhered to the outer leaflet of the vesicle. However, if the initial distance between the NPs is small, one NP is endocytosed, while the other is internalized in the vesicle through a pore.

Finite-temperature effect in the O-acylation of (R,S)-propranolol catalyzed by Candida antarctica lipase B

CalB is a triacylglycerol hydrolase (E.C.3.1.1.3) used in the O-acylation of the beta-adrenergic blocking agent (R,S)-propranolol. The catalytic mechanism involves two steps: enzyme acylation and enzyme deacylation. The enantioselectivity of the O-acylation of (R,S)-propranolol originates from the second step, where the acyl-enzyme transfers the acyl group to the racemic substrate. This step occurs via an initial Michaelis complex (MCC) and a tetrahedral intermediate (TI-2). To gain more insight into the molecular basis of this reaction, we performed an exhaustive conformational sampling along the reaction coordinate of the enantioselective step of the reaction (MCC→TI-2→EPC) applying a QM/MM MD protocol (SCC-DFTB/CHARMM) in combination with umbrella sampling and the weighted histogram analysis method. To identify finite temperature effects we compare the PMF and the potential energy pathway. It is found that the effect of the finite temperature in this reaction is a destabilization of the tetrahedral intermediate and an increase of the barrier height of its formation. This increase is higher for the S-enantiomer.

Phase Transitions and Hysteresis for a Simple Model Liquid Crystal by Replica-Exchange Monte Carlo Simulations.

In this work, the advantages of applying the temperature and pressure replica-exchange method to investigate the phase transitions and the hysteresis for liquid-crystal fluids were demonstrated. In applying this method to the commonly used Hess–Su liquid-crystal model, heat capacity peaks and points of phase co-existence were observed. The absence of a smectic phase at higher densities and a narrow range of the nematic phase were reported. The identity of the crystalline phase of this system was found to a hexagonal close-packed solid. Since the nematic-solid phase transition is strongly first order, care must be taken when using this model not to inadvertently simulate meta-stable nematic states at higher densities. In further analysis, the Weighted Histogram Analysis Method was applied to verify the precise locations of the phase transition points.

Potential of mean force and umbrella sampling simulation for the transport of 5-oxazolidinone in heterotetrameric sarcosine oxidase.

The structure of heterotetrameric sarcosine oxidase (HSO) contains a highly complex system composed of a large cavity and tunnels, which are essential for the reaction and migration of the reactants, products, and intermediates. Previous geometrical analysis using the CAVER program has predicted that there are three possible tunnels, T1, T2, and T3, for the exit pathway of the iminium intermediate, 5-oxazolidinone (5-OXA), of the enzyme reaction. Previous molecular dynamics (MD) simulation of HSO has identified the regions containing the water channels from the density distribution of water. The simulation indicated that tunnel T3 is the most probable exit pathway of 5-OXA. In the present study, the potential of mean force (PMF) for the transport of 5-OXA through tunnels T1, T2, and T3 was calculated using umbrella sampling (US) MD simulations and the weighted histogram analysis method. The PMF profiles for the three tunnels support the notion that tunnel T3 is the exit pathway of 5-OXA, and that 5-OXA tends to stay at the middle of the tunnel. The maximum errors of the calculated PMF for the predicted exit pathway, tunnel T3, were estimated by repeating the US simulations using different sets of initial positions. The PMF profile was also calculated for the transport of glycine within T3. The PMF profiles from the US simulations were in good agreement with the previous predictions that 5-OXA escape through tunnel T3 and how glycine is released to the outside of HSO was discussed.

Molecular dynamics simulation of ion binding in aqueous solution of 18-crown-6 ether

Crown ethers are unique compounds with exceptional ion binding properties. We investigate the 18-crown-6 ether molecule binding selectivity for a set of alkali ions (K+, Rb+, Cs+). We employ combination of umbrella sampling technique with the weighted histogram analysis method to derive free energies of ion–ether association. The obtained results are in a good agreement with experimental values available and represent the ion selectivity feature of the 18-crown-6 ether.

Feature Extraction Technique Using Weighted Histogram Analysis Method (WHAM) for Herbs Discrimination Based on Gas Chromatography Signal

Herbs discrimination by investigating volatile compound using Gas Chromatography Mass Spectrometry (GCMS) is a common method adopted by botanists and scientists. Based on this common method, usually botanists and scientists would only focus on the major volatile compound in order to determine the species of the herbs. However, it is difficult to differentiate the herbs species of the same family group based on the pattern of chromatography signal since they may have almost similar physical features, characteristics, and aroma. In this case, the minor volatile compound needs to be considered in the herbs discrimination analysis. This study proposes the adoption of a Weighted Histogram Analysis Method (WHAM) that utilizes a combination histogram between two single feature histograms of peak area and peak height data in order to extract the new features based on minor and major volatile compound data (chemical properties) derived from chromatography signal patterns. From the results, it is found that WHAM technique results in better discrimination and classification between herbs species in same family group compared to the results without application of WHAM technique for feature extraction. The improvement in reducing the overlap between herbs group clustering can result in better classification as it will increase the classification accuracy.