Free Energy Approaches Research Articles

CL-FEP: An End-State Free Energy Perturbation Approach.

The Central Limit Free Energy Perturbation (CL-FEP) approach, based on the Free Energy Perturbation (FEP) theory and the Central Limit (CL) theorem, allows evaluating the FEP identity directly from the energy samples of the end states of a system transformation. In CL-FEP, energies obtained from explicit solvent simulations are used to estimate the absolute free energy change. No fitted parameters are introduced in our implementation, and no stratification is needed to obtain accurate free energy evaluations. We first illustrate the applicability of CL-FEP in four dissimilar benchmark systems from host-guest to proteins-peptide complexes, for which the deviation from the experimental values was below 1 kcal/mol. Next, we extended the validation to three sets of complexes comprising 25 host-guest and protein-ligand systems. There, the mean absolute error of the free energy changes calculated with CL-FEP is between 1.1 and 1.4 kcal/mol among the different sets, which is in the range of accuracy of the most reliable free energy approaches applied to these systems. CL-FEP is an unbiased free energy change estimator that also profits from a bootstrapping algorithm, which makes the evaluation of convergence and confidence an inherent component of the procedure. Finally, we developed CLFEP.pl, an automatic tool to compute free energy changes directly from the energy output of standard molecular dynamics simulations.

A review on neural network models of schizophrenia and autism spectrum disorder

This survey presents the most relevant neural network models of autism spectrum disorder and schizophrenia, from the first connectionist models to recent deep neural network architectures. We analyzed and compared the most representative symptoms with its neural model counterpart, detailing the alteration introduced in the network that generates each of the symptoms, and identifying their strengths and weaknesses. We additionally cross-compared Bayesian and free-energy approaches, as they are widely applied to model psychiatric disorders and share basic mechanisms with neural networks. Models of schizophrenia mainly focused on hallucinations and delusional thoughts using neural dysconnections or inhibitory imbalance as the predominating alteration. Models of autism rather focused on perceptual difficulties, mainly excessive attention to environment details, implemented as excessive inhibitory connections or increased sensory precision. We found an excessively tight view of the psychopathologies around one specific and simplified effect, usually constrained to the technical idiosyncrasy of the used network architecture. Recent theories and evidence on sensorimotor integration and body perception combined with modern neural network architectures could offer a broader and novel spectrum to approach these psychopathologies. This review emphasizes the power of artificial neural networks for modeling some symptoms of neurological disorders but also calls for further developing of these techniques in the field of computational psychiatry.

Fast anharmonic free energy method with an application to vacancies in ZrC

We propose an approach to calculate the anharmonic part of the volumetric-strain and temperature dependent free energy of a crystal. The method strikes an effective balance between accuracy and computational efficiency, showing a $\times10$ speed-up on comparable free energy approaches at the level of density functional theory, with average errors less than 1 meV/atom. As a demonstration we make new predictions on the thermodynamics of substoichiometric ZrC$_x$, including vacancy concentration and heat capacity.

Impact of short-term aging temperatures of asphalt binder and aggregate roughness levels on bond strength

In the field, roughness of aggregates varies during production process due to abrasion, stripping and ravelling, which influences its bonding potential with asphalt binder. However, no study relating aggregate roughness levels and bond strength of aggregate-binder system was conducted. In addition, change in short-term aging temperature of asphalt binder during production processes may affect aggregate-binder bond. Therefore, present study evaluates effects of aggregate roughness and short-term aging temperature of binder on bond strength based on surface free energy (SFE) and asphalt bond strength (ABS) test approaches. For the purpose of the study, basalt aggregate was polished using 220, 400, 600 and 800 grade silicon carbide (SiC) powders and their roughness values were measured using zeta-profiler. Also, AC30 binder was short-term aged using rolling thin-film oven (RTFO) at temperatures 143 °C, 163 °C and 173 °C. Bond strength was evaluated using energy ratio (ER) and bond strength ratio (BSR) parameters, obtained from SFE and ABS approaches, respectively. Bond strength of aggregates having different roughness values was evaluated with standard short-term aged binder (163 °C for 85 min) whereas, bond strength of binder aged at different temperatures were evaluated with 220SiC polished aggregate. The present study also validates use of 280SiC polished aggregates as recommended by AASHTO T-361, by examining relationship between aggregate roughness and BSR. Results showed that bond strength increases with increase in aggregate roughness and short-term aging temperature of binder. Further, one-way ANOVA and Post-Hoc test showed significant effect of aggregate roughness and binder aging levels on aggregate-binder bond strength. Markedly, linear correlation was observed between ER and BSR with aggregate roughness and binder short-term aging conditions.

Gō model revisited.

This review discusses Gō models broadly used in biomolecular simulations. I start with a brief description of the original lattice model study by Nobuhiro Gō. Then, the theory of protein folding behind Gō model, free energy approaches, and off-lattice Gō models are reviewed. I also mention a stringent test for the assumption in Gō models given from all-atom molecular dynamics simulations. Subsequently, I move to application of Gō models to protein dynamical functions. Various extension of Gō models is also reviewed. Finally, some publicly available tools to use Gō models are listed.

Configurational entropy of polydisperse supercooled liquids.

We propose a computational method to measure the configurational entropy in generic polydisperse glass-formers. In particular, our method resolves issues related to the diverging mixing entropy term due to a continuous polydispersity. The configurational entropy is measured as the difference between the well-defined fluid entropy and a more problematic glass entropy. We show that the glass entropy can be computed by a simple generalisation of the Frenkel-Ladd thermodynamic integration method, which takes into account permutations of the particle diameters. This approach automatically provides a physically meaningful mixing entropy for the glass entropy and includes contributions that are not purely vibrational. The proposed configurational entropy is thus devoid of conceptual and technical difficulties due to continuous polydispersity, while being conceptually closer, but technically simpler, than alternative free energy approaches.

Molecular simulations and free-energy calculations suggest conformation-dependent anion binding to a cytoplasmic site as a mechanism for Na+/K+-ATPase ion selectivity

Na+/K+-ATPase transports Na+ and K+ ions across the cell membrane via an ion-binding site becoming alternatively accessible to the intra- and extracellular milieu by conformational transitions that confer marked changes in ion-binding stoichiometry and selectivity. To probe the mechanism of these changes, we used molecular simulation and free-energy perturbation approaches to identify probable protonation states of Na+- and K+-coordinating residues in E1P and E2P conformations of Na+/K+-ATPase. Analysis of these simulations revealed a molecular mechanism responsible for the change in protonation state: the conformation-dependent binding of an anion (a chloride ion in our simulations) to a previously unrecognized cytoplasmic site in the loop between transmembrane helices 8 and 9, which influences the electrostatic potential of the crucial Na+-coordinating residue Asp926 This mechanistic model is consistent with experimental observations and provides a molecular-level picture of how E1P to E2P enzyme conformational transitions are coupled to changes in ion-binding stoichiometry and selectivity.

Wall depletion length of a channel-confined polymer.

Numerous experiments have taken advantage of DNA as a model system to test theories for a channel-confined polymer. A tacit assumption in analyzing these data is the existence of a well-defined depletion length characterizing DNA-wall interactions such that the experimental system (a polyelectrolyte in a channel with charged walls) can be mapped to the theoretical model (a neutral polymer with hard walls). We test this assumption using pruned-enriched Rosenbluth method (PERM) simulations of a DNA-like semiflexible polymer confined in a tube. The polymer-wall interactions are modeled by augmenting a hard wall interaction with an exponentially decaying, repulsive soft potential. The free energy, mean span, and variance in the mean span obtained in the presence of a soft wall potential are compared to equivalent simulations in the absence of the soft wall potential to determine the depletion length. We find that the mean span and variance about the mean span have the same depletion length for all soft potentials we tested. In contrast, the depletion length for the confinement free energy approaches that for the mean span only when depletion length no longer depends on channel size. The results have implications for the interpretation of DNA confinement experiments under low ionic strengths.

Modulating O2 Permeability of the Central Pore of RH50 by in Silico Site Directed Mutagenesis

The Rhesus (Rh) family of integral membrane proteins exhibit ubiquitous expression among animals and their structural homology to ammonium transport (Amt) proteins suggests their function as ammonia/ammonium transporters. They have also been implicated in facilitating permeation of gaseous species, such as dioxygen and carbon dioxide, across the cell membrane. Rh50 proteins are expressed in erythrocytes, as well as in epithelial tissues of organs such as the kidneys where critical nitrogen processing occurs. Using the solved Rh50 structure from bacterial homologue NeRh50 found in Nitrosomonas europaea, we employed a combination of molecular dynamics and free energy approaches including implicit ligand sampling, explicit ligand sampling, umbrella sampling, and in silico mutagenesis to characterize O2 permeation pathway, mechanism, and energetics. Free-energy calculations reveal high-energy barriers both in the monomeric pore as well as the pore formed in the trimerization domain similar to those previously reported for CO2 and NH3. Based on the the free-energy calculation results, we identified key amino acids potentially responsible for the high barrier. Bioinformatic analysis of these sites shows that some of the residues lining the high-energy region of the pore exhibit a high degree of conservation only for members of the Rh50 family, with the exception of NeRh50. Residues Thr199 in Arg202, which are localized near the pore opening, exist as hydrophobic amino acids in the protein's eukaryotic counterparts and Leu206 in NeRh50 presents itself primarily as an isoleucine. Mutation of these residues to the conserved amino acids results in a reduction in the energy barrier. This suggests a permeation pathway for gaseous species not observable in the native NeRh50 structure which can be verified by functional characterization of the mutants.

Role of R292K mutation in influenza H7N9 neuraminidase toward oseltamivir susceptibility: MD and MM/PB(GB)SA study.

The H7N9 avian influenza virus is a novel re-assortment from at least four different strains of virus. Neuraminidase, which is a glycoprotein on the surface membrane, has been the target for drug treatment. However, some H7N9 strains that have been isolated from patient after drug treatment have a R292K mutation in neuraminidase. This substitution was found to facilitate drug resistance using protein- and virus- assays, in particular it gave a high resistance to the most commonly used drug, oseltamivir. The aim of this research is to understand the source of oseltamivir resistance using MD simulations and the MM/PB(GB)SA binding free energy approaches. Both methods can predict the reduced susceptibility of oseltamivir in good agreement to the IC 50 binding energy, although MM/GBSA underestimates this prediction compared to the MM/PBSA calculation. Electrostatic interaction is the main contribution for oseltamivir binding in terms of both interaction and solvation. We found that the source of the drug resistance is a decrease in the binding interaction combined with the reduction of the dehydration penalty. The smaller K292 mutated residue has a larger binding pocket cavity compared to the wild-type resulting in the loss of drug carboxylate-K292 hydrogen bonding and an increased accessibility for water molecules around the K292 mutated residue. In addition, oseltamivir does not bind well to the R292K mutant complex as shown by the high degree of fluctuation in ligand RMSD during the simulation and the change in angular distribution of bulky side chain groups.

Identification of possible siRNA molecules for TDP43 mutants causing amyotrophic lateral sclerosis: In silico design and molecular dynamics study

The DNA binding protein, TDP43 is a major protein involved in amyotrophic lateral sclerosis and other neurological disorders such as frontotemporal dementia, Alzheimer disease, etc. In the present study, we have designed possible siRNAs for the glycine rich region of tardbp mutants causing ALS disorder based on a systematic theoretical approach including (i) identification of respective codons for all mutants (reported at the protein level) based on both minimum free energy and probabilistic approaches, (ii) rational design of siRNA, (iii) secondary structure analysis for the target accessibility of siRNA, (iii) determination of the ability of siRNA to interact with mRNA and the formation/stability of duplex via molecular dynamics study for a period of 15ns and (iv) characterization of mRNA–siRNA duplex stability based on thermo-physical analysis. The stable GC-rich siRNA expressed strong binding affinity towards mRNA and forms stable duplex in A-form. The linear dependence between the thermo-physical parameters such as Tm, GC content and binding free energy revealed the ability of the identified siRNAs to interact with mRNA in comparable to that of the experimentally reported siRNAs. Hence, this present study proposes few siRNAs as the possible gene silencing agents in RNAi therapy based on the in silico approach.

Pressure and surface tension of an active simple liquid: a comparison between kinetic, mechanical and free-energy based approaches.

We discuss different definitions of pressure for a system of active spherical particles driven by a non-thermal coloured noise. We show that mechanical, kinetic and free-energy based approaches lead to the same result up to first order in the non-equilibrium expansion parameter. The first prescription is based on a generalisation of the kinetic mesoscopic virial equation and expresses the pressure exerted on the walls in terms of the average of the virial of the inter-particle forces. In the second approach, the pressure and the surface tension are identified with the volume and area derivatives, respectively, of the partition function associated with the known stationary non-equilibrium distribution of the model. The third method is a mechanical approach and is related to the work necessary to deform the system. The pressure is obtained by comparing the expression of the work in terms of local stress and strain with the corresponding expression in terms of microscopic distribution. This is determined from the force balance encoded in the Born-Green-Yvon equation. Such a method has the advantage of giving a formula for the local pressure tensor and the surface tension even in inhomogeneous situations. By direct inspection, we show that the three procedures lead to the same values of the pressure, and give support to the idea that the partition function, obtained via the unified coloured noise approximation, is more than a formal property of the system, but determines the stationary non-equilibrium thermodynamics of the model.

Exploring the mechanochemical cycle of dynein motor proteins: structural evidence of crucial intermediates.

Dyneins, a class of motor proteins consisting of six AAA+ modules (AAA1-AAA6), convert chemical energy derived from the hydrolysis of ATP into mechanical energy to walk along the microtubule track towards its minus end while accomplishing various cellular tasks including the transportation of various intracellular cargos. In a full mechanochemical cycle, dynein goes through ATP binding induced open to closed state transition of AAA1, hydrolysis of that ATP and closed to open state transition induced by the release of hydrolysed products along with linker remodelling in different nucleotide states. Here we built structure based models (SBMs) to explore the sequence of events of this mechanochemical cycle from structural aspects. Free energy and kinetic simulation approaches on a multi-basin SBM of dynein reveal the following pathways: (1) in the closing pathway, the AAA1 domain first converts to a closed state followed by the movement of the linker and (2) in the opening transition, initially the AAA1 domain partially opens up and then the complete linker movement takes place followed by the complete opening of the AAA1 domain. In the opening transition, we have observed two intermediate states from our simulations where the AAA1 domain is partially opened. However, in one state the linker is at a closed position and in the other the linker is at an open position. The existence of such intermediates (Pi released, ADP bound state) of dynein has been suggested by numerous experimental studies earlier. Finally, we discuss the biological relevance of this sequence of events in terms of processivity and efficiency of the cycle. The current study also shows how the basic principle of protein folding can be extended to understand complex phenomena like the stepping mechanism of motor proteins.

FEW: A workflow tool for free energy calculations of ligand binding

In the later stages of drug design projects, accurately predicting relative binding affinities of chemically similar compounds to a biomolecular target is of utmost importance for making decisions based on the ranking of such compounds. So far, the extensive application of binding free energy approaches has been hampered by the complex and time-consuming setup of such calculations. We introduce the free energy workflow (FEW) tool that facilitates setup and execution of binding free energy calculations with the AMBER suite for multiple ligands. FEW allows performing free energy calculations according to the implicit solvent molecular mechanics (MM-PB(GB)SA), the linear interaction energy, and the thermodynamic integration approaches. We describe the tool's architecture and functionality and demonstrate in a show case study on Factor Xa inhibitors that the time needed for the preparation and analysis of free energy calculations is considerably reduced with FEW compared to a fully manual procedure.

Free Energy Approaches for Modeling Atomic Force Microscopy in Liquids.

High resolution atomic force microscopy (AFM) in liquids offers atomic scale insight into the structure at water/solid interfaces and is perhaps the only tool capable of resolving the nature of formed hydration layers. However, convolution between the imaging signal and the tip/surface interactions and hydration layers means that interpretation is far from straightforward. Modeling the complex imaging mechanism of atomic force microscopy in liquids requires calculation of the free energy profile as a function of the distance between AFM tip and surface. Its derivative is the best approximation for the force acting on the AFM tip, including entropic contributions from interactions with water molecules in hydration layers over the surface and around the tip apex. In order to establish a reliable approach for these simulations, we compare two methods of calculating free energy profiles from atomistic molecular dynamics simulations, umbrella sampling and free energy perturbation, on two model surfaces, calcium fluoride and calcium carbonate. Our results demonstrate that both methods effectively provide equivalent free energy profiles but offer different possibilities in terms of efficiency, constraints, and analysis of the free energy components.

The thermal Casimir–Polder interaction of an atom with a spherical plasma shell

The van der Waals and Casimir–Polder interaction energy of an atom with an infinitely thin sphere with finite conductivity is investigated in the framework of the hydrodynamic approach at finite temperature. This configuration models the real interaction of an atom with fullerene. The Lifshitz approach is used to find the free energy. We find the explicit expression for the free energy and perform its analysis for (i) high and low temperatures, (ii) large radii of the sphere and (iii) short separation between an atom and sphere. At low temperatures the thermal part of the free energy approaches zero as the fourth power of the temperature, while for high temperatures it is proportional to the first degree of the temperature. The entropy of this system is positive for small radii of the sphere and it becomes negative at low temperatures and for large radii of the sphere.

Disease-Associated Mutations That Alter the RNA Structural Ensemble

Genome-wide association studies (GWAS) often identify disease-associated mutations in intergenic and non-coding regions of the genome. Given the high percentage of the human genome that is transcribed, we postulate that for some observed associations the disease phenotype is caused by a structural rearrangement in a regulatory region of the RNA transcript. To identify such mutations, we have performed a genome-wide analysis of all known disease-associated Single Nucleotide Polymorphisms (SNPs) from the Human Gene Mutation Database (HGMD) that map to the untranslated regions (UTRs) of a gene. Rather than using minimum free energy approaches (e.g. mFold), we use a partition function calculation that takes into consideration the ensemble of possible RNA conformations for a given sequence. We identified in the human genome disease-associated SNPs that significantly alter the global conformation of the UTR to which they map. For six disease-states (Hyperferritinemia Cataract Syndrome, β-Thalassemia, Cartilage-Hair Hypoplasia, Retinoblastoma, Chronic Obstructive Pulmonary Disease (COPD), and Hypertension), we identified multiple SNPs in UTRs that alter the mRNA structural ensemble of the associated genes. Using a Boltzmann sampling procedure for sub-optimal RNA structures, we are able to characterize and visualize the nature of the conformational changes induced by the disease-associated mutations in the structural ensemble. We observe in several cases (specifically the 5′ UTRs of FTL and RB1) SNP–induced conformational changes analogous to those observed in bacterial regulatory Riboswitches when specific ligands bind. We propose that the UTR and SNP combinations we identify constitute a “RiboSNitch,” that is a regulatory RNA in which a specific SNP has a structural consequence that results in a disease phenotype. Our SNPfold algorithm can help identify RiboSNitches by leveraging GWAS data and an analysis of the mRNA structural ensemble.

Secondary Structure Models for the Internal Transcribed Spacer (ITS) Region 1 from Symbiotic Dinoflagellates

Ribosomal genes and their spacers have been extensively utilized to examine the biodiversity and phylogenetics of protists. Among these, the internal transcribed spacer regions 1 and 2 (ITS1 and ITS2) are known to form secondary structures that are critically important for proper processing of the pre-rRNA into mature ribosomes. Although the secondary structure of ITS2 has been widely investigated, considerably less is known about ITS1 and its secondary structure. Here, secondary structures of the ITS1 were modeled for 46 ITS "types" from Symbiodinium, a diverse dinoflagellate genus that forms symbioses with many protists and metazoans, using comparative phylogenetic and minimum free energy approaches. The predicted ITS1 secondary structures for each Symbiodinium "type" were highly stable (DeltaG=-46.40 to -85.30 kcal mol(-1) at 37 degrees C) and consisted of an open loop with five helices separated by single-stranded regions. Several structural characteristics were conserved within monophyletic sub-groups, providing additional support for the predicted structures and the relationships within this genus. Finally, the structures were applied to identify potential pseudogenes from five Symbiodinium ITS1 datasets. Consequently, ITS1 secondary structures are useful in understanding the biology and phylogenetics, as well as recognizing and excluding questionable sequences from datasets, of protists such as Symbiodinium.

Free energy profiles of amino acid side chain analogs near water‐vapor interface obtained via MD simulations

Solvation of 13 neutral amino acid side chain analogs at water-vapor interface was studied by computing high precision free energy profiles of the molecules across the interface using molecular dynamics (MD) simulations. The SPC water model (Berendsen, H. J. C., Postma, J. P. M., van Gunsteren, W. F., Hermans, P. A. K. J., Dixon, R., Cornell, W., Fox, T., Chipot, C., Pohorille, A. In: Wilkinson, A., Weiner, P. and van Gunsteren, W. F. editors. Intermolecular Forces, 1981, 3, 331) and OPLS-AA (Jorgensen et al., J Am Chem Soc, 1996, 118, 11225) potential parameter sets were used. A rigorous approach for the computation of high precision free energy profiles at water-vapor interface using constraint force technique is implemented. Methodology of obtaining high precision potential of mean force (PMF) profiles free of simulation artifacts in MD simulations is outlined and discussed. The accuracy of the calculations is examined by comparing the hydration free energies of studied solutes obtained from PMF calculations and separately using Bennett acceptance ratio technique by decoupling solvent-solute interactions in the bulk. All molecules exhibit a free energy minimum at interface. No significant desolvation barrier is observed for any of the studied species. Adsorption energies for studied molecules at water-vapor interface are estimated and compared with experimental observations. We find that for modeled neutral compounds pronounced surface influence on solute solvation vanishes already at 6-7 A behind the water surface as the solvation free energy approaches the bulk value. The possibility of force field refinement using adsorption free energies is outlined.

Evaluating the Molecular Mechanics Poisson−Boltzmann Surface Area Free Energy Method Using a Congeneric Series of Ligands to p38 MAP Kinase

The recently described molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method for calculating free energies is applied to a congeneric series of 16 ligands to p38 MAP kinase whose binding constants span approximately 2 orders of magnitude. These compounds have previously been used to test and compare other free energy calculation methods, including thermodynamic integration (TI), OWFEG, ChemScore, PLPScore, and Dock Energy Score. We find that the MM-PBSA performs relatively poorly for this set of ligands, yielding results much inferior to those from TI or OWFEG, inferior to Dock Energy Score, and not appreciably better than ChemScore or PLPScore but at an appreciably larger computational cost than any of these other methods. This suggests that one should be selective in applying the MM-PBSA method and that for systems that are amenable to other free energy approaches, these other approaches may be preferred. We also examine the single simulation approximation for MM-PBSA, whereby the required ligand and protein trajectories are extracted from a single MD simulation rather than two separate MD runs. This assumption, sometimes used to speed the MM-PBSA calculation, is found to yield significantly inferior results with only a moderate net percentage reduction in total simulation time.