Alchemical Free Energy Perturbation Research Articles

Absolute binding affinities of phospholipids to Erwinia ligand-gated ion channel (ELIC) by streamlined alchemical free energy perturbations (SAFEP).

An increasing number of cryo-EM structures of membrane proteins contain putative lipid binding sites with partially resolved lipid densities. In principle, these densities might be identified by computing the binding free energies of candidate ligands by alchemical free energy perturbation (FEP). In practice, however, these systems have several challenging features including lipid flexibility, system heterogeneity, and low binding affinity. Each of these confounds the distinction between the bound and unbound state: flexible molecules have large accessible volumes in conformation space; heterogeneous environments complicate the unbound state; and open, low-affinity sites mean that a “bound” ligand may still have considerably higher entropy than more in a more traditional, high-affinity binding pocket. Our SAFEP method seeks to address these problems by defining a single collective variable to define the bound state. This distance-from-bound configuration (DBC) metric reliably classifies ligand-protein conformations as “bound” or “unbound.” SAFEP has been successfully applied to the computation of absolute binding affinities of cholesterol to three GPCRs and relative binding affinities of phospholipids to two mutants of ELIC based on structural data. We have extended SAFEP to the calculation of absolute binding affinities of phospholipids. This has required additional sampling of each state and refinement to the collective variables used, especially in the unbound state. While this is somewhat more computationally expensive than computing relative free energies, absolute affinities make functional questions more accessible and mitigate the error propagation inherent to relative binding affinities. It is our hope that these methods will aid in the interpretation of similar systems.

Oncogenic mutations in the DNA-binding domain of FOXO1 that disrupt folding: Quantitative insights from experiments and molecular simulation.

FOXO1, a member of the family of winged-helix motif Forkhead box (FOX) transcription factors, is the most abundantly expressed FOXO member in mature B cells. Sequencing of diffuse large B-cell lymphoma (DLBCL) tumors and cell lines identified specific mutations in the forkhead domain linked to loss of function. Differential scanning calorimetry and thermal shift assays were used to characterize how eight of these mutations affect the stability of the FOX domain. Mutations L183P and L183R were found to be particularly destabilizing. Electrophoresis mobility shift assays show these same mutations also disrupt FOXO1 binding to their canonical DNA sequences, suggesting that the loss of function is due to destabilization of the folded structure. Computational modeling of the effect of mutations on FOXO1 folding was performed using alchemical free energy perturbation (FEP), and a Markov model of the entire folding reaction was constructed from massively parallel molecular simulations, which predicts folding pathways involving the late folding of helix α3. Although FEP can qualitatively predict the destabilization from L183 mutations, we find that a simple hydrophobic transfer model, combined with estimates of unfolded-state solvent-accessible surface areas from molecular simulations, is able to more accurately predict changes in folding free energies due to mutations. These results suggest that the atomic detail provided by simulations is important for the accurate prediction of mutational effects on folding stability. Corresponding disease-associated mutations in other FOX family members support further experimental and computational studies of the folding mechanism of FOX domains.

Comparison of cyclic AMP binding activity within HCN channel subfamily using atomistic molecular dynamics simulations of isolated cyclic nucleotide binding domains.

Hyperpolarization-activated cyclic nucleotide-gated (HCN) channels are a subfamily of homotetrameric, voltage gated potassium/sodium channels whose activity is modulated by cyclic AMP (cAMP). The binding of cAMP to an intracellular cyclic nucleotide binding domain (CNBD) promotes an increase in channel activity through a depolarizing shift in the voltage threshold necessary for channel activation, as well as kinetically favoring the active channel conformation, leading to improved conductance. While the precise mechanisms of this change are still being investigated, the shift in the CNDB region which is induced by the binding of cAMP is able to be propagated to the remainder of the channel through an adjacent C-linker, causing an alteration in the channel that affects conductance. Although various static structures of the four HCN isoforms with and without bound cyclic nucleotide have been solved and deposited on the Protein Data Bank, there is still a need to understand the precise binding interactions which lead to different sensitivities and selectivities for cyclic nucleotides within the subfamily. Using atomistic molecular dynamics, we have simulated the isolated CNBD of HCN1-4 bound to cAMP and have used both equilibrium simulations and alchemical free energy perturbation to characterize and compare the binding activities of the different isoforms. These simulations highlight the differential cyclic nucleotide binding behavior between HCN channels.

Regioselective metabolism of cannabinoids by cytochrome P450 2J2.

Cannabichromene (CBC) is a major non-psychoactive phytocannabinoid extracted from the cannabis plant. CBC has been shown to have anti-inflammatory effects; hence it can serve as a non-narcotic alternative to tetrahydrocannabinol (THC). Nevertheless, the metabolic pathways and metabolites of CBC have not been fully characterized. Metabolism of phytocannabinoids typically occurs via a superfamily of hemoproteins called cytochrome P450s (CYP450s). Our liquid chromatography-mass spectrometry (LC/MS) experiments show that CBC is metabolized to three different products by various human CYP450s. We identified that CYP450 2J2 (CYP2J2) predominantly produces 8'-hydroxy-CBC, an allylic hydroxide, out of the three metabolites. To investigate the CBC metabolic regioselectivity in CYP2J2, we leveraged atomistic molecular dynamics (MD) and modeling of the system. Using 200 unique conformations of CYP2J2 obtained from MD, we performed ensemble molecular docking of CBC to reveal three putative binding modes. We then assessed the stability of the obtained CBC binding modes by simulating them individually for 1 μs. Our simulations captured an average distance of 4.8 Å between the CYP2J2 heme and the 8’ carbon atom on CBC for one of the predominant binding modes, at which the allylic hydroxide metabolite occurs experimentally. At the regioselective binding mode, we observe conserved hydrogen bonding networks, mainly between the CBC hydroxyl group and the backbone oxygen of I487, as well as N379, stabilized by a water molecule. We further corroborated the stability of the regioselective binding modes by calculating relative binding free energies using alchemical free energy perturbation calculations. In conjunction, our computational and experimental methodologies elucidate a likely mode for regioselective binding and metabolism of CBC by CYP2J2, furthering the way towards accessible plant-based non-psychoactive anti-inflammatory therapeutics.

Open-channel structure of a pentameric ligand-gated ion channel reveals a mechanism of leaflet-specific phospholipid modulation

Pentameric ligand-gated ion channels (pLGICs) mediate synaptic transmission and are sensitive to their lipid environment. The mechanism of phospholipid modulation of any pLGIC is not well understood. We demonstrate that the model pLGIC, ELIC (Erwinia ligand-gated ion channel), is positively modulated by the anionic phospholipid, phosphatidylglycerol, from the outer leaflet of the membrane. To explore the mechanism of phosphatidylglycerol modulation, we determine a structure of ELIC in an open-channel conformation. The structure shows a bound phospholipid in an outer leaflet site, and structural changes in the phospholipid binding site unique to the open-channel. In combination with streamlined alchemical free energy perturbation calculations and functional measurements in asymmetric liposomes, the data support a mechanism by which an anionic phospholipid stabilizes the activated, open-channel state of a pLGIC by specific, state-dependent binding to this site.

Free energy perturbation calculations of tetrahydroquinolines complexed to the first bromodomain of BRD4

Alchemical free energy perturbation (FEP) theory is widely used nowadays to calculate protein–ligand binding energies, often in support of drug discovery endeavours. We assess the accuracy and sensitivity of absolute FEP binding energies with respect to the CHARMM/CGenFF and the AMBER/GAFF force field parameterisations for a set of tetrahydroquinoline inhibitors of the first bromodomain of BRD4, a target of keen interest for the development of anti-cancer drugs. We find that AMBER/GAFF is better able than CHARMM/CGenFF to cover the range of and to distinguish between the relative binding energies of the 16 ligands.

Oncogenic Mutations in the DNA-Binding Domain of FOXO1 that Disrupt Folding: Quantitative Insights from Experiments and Molecular Simulations.

FOXO1, a member of the family of winged-helix motif Forkhead box (FOX) transcription factors, is the most abundantly expressed FOXO member in mature B cells. Sequencing of diffuse large B-cell lymphoma (DLBCL) tumors and cell lines identified specific mutations in the forkhead domain linked to loss of function. Differential scanning calorimetry and thermal shift assays were used to characterize how eight of these mutations affect the stability of the FOX domain. Mutations L183P and L183R were found to be particularly destabilizing. Electrophoresis mobility shift assays show these same mutations also disrupt FOXO1 binding to their canonical DNA sequences, suggesting that the loss of function is due to destabilization of the folded structure. Computational modeling of the effect of mutations on FOXO1 folding was performed using alchemical free energy perturbation (FEP), and a Markov model of the entire folding reaction was constructed from massively parallel molecular simulations, which predicts folding pathways involving the late folding of helix α3. Although FEP can qualitatively predict the destabilization from L183 mutations, we find that a simple hydrophobic transfer model, combined with estimates of unfolded-state solvent-accessible surface areas from molecular simulations, is able to more accurately predict changes in folding free energies due to mutations. These results suggest that the atomic detail provided by simulations is important for the accurate prediction of mutational effects on folding stability. Corresponding disease-associated mutations in other FOX family members support further experimental and computational studies of the folding mechanism of FOX domains.

Large-scale application of free energy perturbation calculations for antibody design

Alchemical free energy perturbation (FEP) is a rigorous and powerful technique to calculate the free energy difference between distinct chemical systems. Here we report our implementation of automated large-scale FEP calculations, using the Amber software package, to facilitate antibody design and evaluation. In combination with Hamiltonian replica exchange, our FEP simulations aim to predict the effect of mutations on both the binding affinity and the structural stability. Importantly, we incorporate multiple strategies to faithfully estimate the statistical uncertainties in the FEP results. As a case study, we apply our protocols to systematically evaluate variants of the m396 antibody for their conformational stability and their binding affinity to the spike proteins of SARS-CoV-1 and SARS-CoV-2. By properly adjusting relevant parameters, the particle collapse problems in the FEP simulations are avoided. Furthermore, large statistical errors in a small fraction of the FEP calculations are effectively reduced by extending the sampling, such that acceptable statistical uncertainties are achieved for the vast majority of the cases with a modest total computational cost. Finally, our predicted conformational stability for the m396 variants is qualitatively consistent with the experimentally measured melting temperatures. Our work thus demonstrates the applicability of FEP in computational antibody design.

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.

Computational Assessment of Protein-Protein Binding Specificity within a Family of Synaptic Surface Receptors.

Atomic-level information is essential to explain the formation of specific protein complexes in terms of structure and dynamics. The set of Dpr and DIP proteins, which play a key role in the neuromorphogenesis in the nervous system of Drosophila melanogaster, offer a rich paradigm to learn about protein-protein recognition. Many members of the DIP subfamily cross-react with several members of the Dpr family and vice versa. While there exists a total of 231 possible Dpr-DIP heterodimer complexes from the 21 Dpr and 11 DIP proteins, only 57 "cognate" pairs have been detected by surface plasmon resonance (SPR) experiments, suggesting that the remaining 174 pairs have low or unreliable binding affinity. Our goal is to assess the performance of computational approaches to characterize the global set of interactions between Dpr and DIP proteins and identify the specificity of binding between each DIP with their corresponding Dpr binding partners. In addition, we aim to characterize how mutations influence the specificity of the binding interaction. In this work, a wide range of knowledge-based and physics-based approaches are utilized, including mutual information, linear discriminant analysis, homology modeling, molecular dynamics simulations, Poisson-Boltzmann continuum electrostatics calculations, and alchemical free energy perturbation to decipher the origin of binding specificity of the Dpr-DIP complexes examined. Ultimately, the results show that those two broad strategies are complementary, with different strengths and limitations. Biological inter-relations are more clearly revealed through knowledge-based approaches combining evolutionary and structural features, the molecular determinants controlling binding specificity can be predicted accurately with physics-based approaches based on atomic models.

Symmetry-Adapted Restraints for Binding Free Energy Calculations.

Binding free energy calculations rely critically on a precise definition of the bound state and well-designed ligand restraints to ensure that binding free energy calculations converge rapidly and yield estimates of well-defined thermodynamic quantities. The distance-to-bound-configuration (DBC) is a single variable that can precisely delineate the bound state of a ligand including translational, rotational and conformational degrees of freedom and has been successfully used to capture binding modes with complex geometries. DBC is defined as the root-mean-square deviation (RMSD) of ligand coordinates in the frame of reference of the binding site. In the special case where the ligand features symmetry-equivalent atoms, a standard RMSD arbitrarily distinguishes equivalent poses, mixing equivalent and nonequivalent degrees of freedom, and preventing the precise delineation of the bound state ensemble, which negates the benefits of defining a flat-bottom binding restraint. To remedy this, we introduce a symmetry-adapted DBC coordinate where the RMSD is minimized over permutations of equivalent ligand atoms. This coordinate is implemented in a portable software library, the Collective Variables Module. We tested the approach by computing the absolute binding free energy of benzene to the engineered site of a mutant lysozyme (L99A/M102H) using alchemical free energy perturbation. We found that the symmetry-adapted restraint leads to well-behaved convergence of both the decoupling free energy in the binding site and the restrained free energy in the gas phase, recovering the affinity computed using a classic center-of-mass restraint. Thus, symmetry-adapted DBC seamlessly generalizes the benefits of DBC restraints to the case of symmetric ligands. The underlying symmetric RMSD coordinate can also be used for analyzing or biasing simulations in other contexts than affinity predictions.

Modeling the Binding and Conformational Energetics of a Targeted Covalent Inhibitor to Bruton's Tyrosine Kinase.

Targeted covalent inhibitors (TCIs) bind to their targets in both covalent and noncovalent modes, providing exceptionally high affinity and selectivity. These inhibitors have been effectively employed as inhibitors of protein kinases, with Taunton and coworkers (Nat. Chem. Biol. 2015, 11, 525-531) reporting a notable example of a TCI with a cyanoacrylamide warhead that forms a covalent thioether linkage to an active-site cysteine (Cys481) of Bruton's tyrosine kinase (BTK). The specific mechanism of the binding and the relative importance of the covalent and noncovalent interactions is difficult to determine experimentally, and established simulation methods for calculating the absolute binding affinity of an inhibitor cannot describe the covalent bond-forming steps. Here, an integrated approach using alchemical free-energy perturbation and QM/MM molecular dynamics methods was employed to model the complete Gibbs energy profile for the covalent inhibition of BTK by a cyanoacrylamide TCI. These calculations provide a rigorous and complete absolute Gibbs energy profile of the covalent modification binding process. Following a classic thiol-Michael addition mechanism, the target cysteine is deprotonated to form a nucleophilic thiolate, which then undergoes a facile conjugate addition to the electrophilic functional group to form a bond with the noncovalently bound ligand. This model predicts that the formation of the covalent linkage is highly exergonic relative to the noncovalent binding alone. Nevertheless, noncovalent interactions between the ligand and individual amino acid residues in the binding pocket of the enzyme are also essential for ligand binding, particularly van der Waals dispersion forces, which have a larger contribution to the binding energy than the covalent component in absolute terms. This model also shows that the mechanism of covalent modification of a protein occurs through a complex series of steps and that entropy, conformational flexibility, noncovalent interactions, and the formation of covalent linkage are all significant factors in the ultimate binding affinity of a covalent drug to its target.

CHARMM-GUI Free Energy Calculator for Absolute and Relative Ligand Solvation and Binding Free Energy Simulations.

Alchemical free energy simulations have long been utilized to predict free energy changes for binding affinity and solubility of small molecules. However, while the theoretical foundation of these methods is well established, seamlessly handling many of the practical aspects regarding the preparation of the different thermodynamic end states of complex molecular systems and the numerous processing scripts often remains a burden for successful applications. In this work, we present CHARMM-GUI Free Energy Calculator (http://www.charmm-gui.org/input/fec) that provides various alchemical free energy perturbation molecular dynamics (FEP/MD) systems with input and post-processing scripts for NAMD and GENESIS. Four submodules are available: Absolute Ligand Binder (for absolute ligand binding FEP/MD), Relative Ligand Binder (for relative ligand binding FEP/MD), Absolute Ligand Solvator (for absolute ligand solvation FEP/MD), and Relative Ligand Solvator (for relative ligand solvation FEP/MD). Each module is designed to build multiple systems of a set of selected ligands at once for high-throughput FEP/MD simulations. The capability of Free Energy Calculator is illustrated by absolute and relative solvation FEP/MD of a set of ligands and absolute and relative binding FEP/MD of a set of ligands for T4-lysozyme in solution and the adenosine A2A receptor in a membrane. The calculated free energy values are overall consistent with the experimental and published free energy results (within ∼1 kcal/mol). We hope that Free Energy Calculator is useful to carry out high-throughput FEP/MD simulations in the field of biomolecular sciences and drug discovery.

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.

Open and Closed Structures of a Barium-Blocked Potassium Channel

Barium (Ba2+) is a classic permeant blocker of potassium (K+) channels. The “external lock-in effect” in barium block experiments, whereby the binding of external K+ impedes the forward translocation of the blocker, provides a powerful avenue to investigate the selectivity of the binding sites along the pore of potassium channels. Barium block experiments show that the external lock-in site is highly selective for K+ over Na+. Wild-type KcsA was crystallized in low K+ conditions, and the crystals were soaked in solutions containing various concentrations of barium. Structural analysis reveals open and closed gate conformations of the KcsA channel. Anomalous diffraction experiments show that Ba2+ primarily binds to the innermost site S4 of the selectivity filter of the open-gate conformation and also the site S2, but no binding is detected with the closed-gate conformation. Alchemical free-energy perturbation calculations indicate that the presence of a Ba2+ ion in the selectivity filter boosts the specificity of K+ binding relative to Na+ in the external sites S0–S2.

On Calculating Free Energy Differences Using Ensembles of Transition Paths.

The free energy of a process is the fundamental quantity that determines its spontaneity or propensity at a given temperature. In particular, the binding free energy of a drug candidate to its biomolecular target is used as an objective quantity in drug design. Recently, binding kinetics—rates of association (kon) and dissociation (koff)—have also demonstrated utility for their ability to predict efficacy and in some cases have been shown to be more predictive than the binding free energy alone. Some methods exist to calculate binding kinetics from molecular simulations, although these are typically more difficult to calculate than the binding affinity as they depend on details of the transition path ensemble. Assessing these rate constants can be difficult, due to uncertainty in the definition of the bound and unbound states, large error bars and the lack of experimental data. As an additional consistency check, rate constants from simulation can be used to calculate free energies (using the log of their ratio) which can then be compared to free energies obtained experimentally or using alchemical free energy perturbation. However, in this calculation it is not straightforward to account for common, practical details such as the finite simulation volume or the particular definition of the “bound” and “unbound” states. Here we derive a set of correction terms that can be applied to calculations of binding free energies using full reactive trajectories. We apply these correction terms to revisit the calculation of binding free energies from rate constants for a host-guest system that was part of a blind prediction challenge, where significant deviations were observed between free energies calculated with rate ratios and those calculated from alchemical perturbation. The correction terms combine to significantly decrease the error with respect to computational benchmarks, from 3.4 to 0.76 kcal/mol. Although these terms were derived with weighted ensemble simulations in mind, some of the correction terms are generally applicable to free energies calculated using physical pathways via methods such as Markov state modeling, metadynamics, milestoning, or umbrella sampling.

Menthol Binding to the Human α4β2 Nicotinic Acetylcholine Receptor Facilitated by Its Strong Partitioning in the Membrane.

We utilize various computational methodologies to study menthol's interaction with multiple organic phases, a lipid bilayer, and the human α4β2 nicotinic acetylcholine receptor (nAChR), the most abundant nAChR in the brain. First, force field parameters developed for menthol are validated in alchemical free energy perturbation simulations to calculate solvation free energies of menthol in water, dodecane, and octanol and compare the results against experimental data. Next, umbrella sampling is used to construct the free energy profile of menthol permeation across a 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) bilayer. The results from a flooding simulation designed to study the water-membrane partitioning of menthol in a POPC lipid bilayer are used to determine the penetration depth and the preferred orientation of menthol in the bilayer. Finally, employing both docking and flooding simulations, menthol is shown to bind to different sites on the human α4β2 nAChR. The most likely binding mode of menthol to a desensitized membrane-embedded α4β2 nAChR is identified to be via a membrane-mediated pathway in which menthol binds to the sites at the lipid-protein interface after partitioning in the membrane. A rare but distinct binding mode in which menthol binds to the extracellular opening of receptor's ion permeation pore is also reported.

Relative Affinities of General Anesthetics for Experimentally Identified Binding Sites in Ryanodine Receptors (RyR1)

Malignant hyperthermia (MH) is a pharmacogenetic disorder that manifests as a hypermetabolic response to certain anesthetic agents (halothane, isoflurane, sevoflurane, desflurane, succinylcholine). A defective or disordered Calcium ion channel, ryanodine receptor 1 (RyR1), located in the sarcoplasmic reticulum (SR) membrane, underlies MH susceptibility. RyR1, the largest known ion channel, is responsible for the rapid release of Ca2+ ions from SR into the cytoplasm, a key event that triggers skeletal muscle contraction. It has a homotetrameric organization and approximately 5000 residues in each protomer. There are a large number (>100) of known mutations in RyR1 or its binding partners which are associated with MH sensitivity. In the presence of a triggering agent and an MH-causative mutation, the channel may be biased open, allowing uncontrolled Ca2+ release, generating an extreme metabolic load. In order to understand how MH is triggered mechanistically by anesthetics, we use molecular dynamics simulations to study experimentally-identified binding sites in specific domains of the channel and alchemical free energy perturbation calculations to evaluate binding configurations and affinities. We computationally refine residue-level binding sites identified by photoaffinity labeling for selected anesthetic ligands on RyR1. We further evaluate MH-causative mutations that may precipitate an increase or decrease in the binding affinity of anesthetic molecules to RyR1. Evaluating and comparing binding sites of MH-causative anesthetics, in wild-type and mutated RyR1, helps us understand combinations of an RYR1 mutation and a bound anesthetic configuration that bias an open conformation of the channel and trigger MH.

Prediction of octanol-water partition coefficients for the SAMPL6-$$\log P$$ molecules using molecular dynamics simulations with OPLS-AA, AMBER and CHARMM force fields

All-atom molecular dynamics simulations with stratified alchemical free energy calculations were used to predict the octanol-water partition coefficient [Formula: see text] of eleven small molecules as part of the SAMPL6-[Formula: see text] blind prediction challenge using four different force field parametrizations: standard OPLS-AA with transferable charges, OPLS-AA with non-transferable CM1A charges, AMBER/GAFF, and CHARMM/CGenFF. Octanol parameters for OPLS-AA, GAFF and CHARMM were validated by comparing the density as a function of temperature, the chemical potential, and the hydration free energy to experimental values. The partition coefficients were calculated from the solvation free energy for the compounds in water and pure ("dry") octanol or "wet" octanol with 27 mol% water dissolved. Absolute solvation free energies were computed by thermodynamic integration (TI) and the multistate Bennett acceptance ratio with uncorrelated samples from data generated by an established protocol using 5-ns windowed alchemical free energy perturbation (FEP) calculations with the Gromacs molecular dynamics package. Equilibration of sets of FEP simulations was quantified by a new measure of convergence based on the analysis of forward and time-reversed trajectories. The accuracy of the [Formula: see text] predictions was assessed by descriptive statistical measures such as the root mean square error (RMSE) of the data set compared to the experimental values. Discarding the first 1 ns of each 5-ns window as an equilibration phase had a large effect on the GAFF data, where it improved the RMSE by up to 0.8 log units, while the effect for other data sets was smaller or marginally worsened the agreement. Overall, CGenFF gave the best prediction with RMSE 1.2 log units, although for only eight molecules because the current CGenFF workflow for Gromacs does not generate files for certain halogen-containing compounds. Over all eleven compounds, GAFF gave an RMSE of 1.5. The effect of using a mixed water/octanol solvent slightly decreased the accuracy for CGenFF and GAFF and slightly increased it for OPLS-AA. The GAFF and OPLS-AA results displayed a systematic error where molecules were too hydrophobic whereas CGenFF appeared to be more balanced, at least on this small data set.

Molecular Mechanism of T-2 Toxin-Induced Cerebral Edema by Aquaporin-4 Blocking and Permeation.

The present study aimed to reveal the molecular mechanism of T-2 toxin-induced cerebral edema by aquaporin-4 (AQP4) blocking and permeation. AQP4 is a class of aquaporin channels that is mainly expressed in the brain, and its structural changes lead to life-threatening complications such as cardio-respiratory arrest, nephritis, and irreversible brain damage. We employed molecular dynamics simulation, text mining, and in vitro and in vivo analysis to study the structural and functional changes induced by the T-2 toxin on AQP4. The action of the toxin leads to disrupted permeation of water and permeation coefficients are found to be affected, from the native (2.49 ± 0.02 × 10-14 cm3/s) to toxin-treated AQP4 (7.68 ± 0.15 × 10-14 cm3/s) channels. Furthermore, the T-2 toxin forms strong electrostatic interactions at the binding site and pushes the key residues (Ala210, Phe77, Arg216, and His201) outward at the selectivity filter. Also, the role of a histidine residue in the AQP4 channel was identified by alchemical transformation and umbrella sampling methods. Alchemical free-energy perturbation energy for H201A ↔ A201H, which was found to be 3.07 ± 0.18 kJ/mol, indicates the structural importance of the histidine residue at 201. In addition, histopathology and expression of AQP4 in the Mus musculus brain tissues show the damaged and altered expression of the protein. Text mining reveals the co-occurrence of genes/proteins associated with the AQP4 expression and T-2 toxin-induced cell apoptosis, which leads to cerebral edema.