Quantification of Solvent Contribution to the Stability of Noncovalent Complexes

Computational Insights into Cyclodextrin Inclusion Complexes with the Organophosphorus Flame Retardant DOPO

Cyclodextrins (CDs) were used as green char promoters in the formulation of organophosphorus flame retardants (OPFRs) for polymeric materials, and they could reduce the amount of usage of OPFRs and their release into the environment by forming [host:guest] inclusion complexes with them. Here, we report a systematic study on the inclusion complexes of natural CDs (α-, β-, and γ-CD) with a representative OPFR of DOPO using computational methods of molecular docking, molecular dynamics (MD) simulations, and quantum mechanical (QM) calculations. The binding modes and energetics of [host:guest] inclusion complexes were analyzed in details. α-CD was not able to form a complete inclusion complex with DOPO, and the center of mass distance [host:guest] distance amounted to 4-5 Å. β-CD and γ-CD allowed for a deep insertion of DOPO into their hydrophobic cavities, and DOPO was able to frequently change its orientation within the γ-CD cavity. The energy decomposition analysis based on the dispersion-corrected density functional theory (sobEDAw) indicated that electrostatic, orbital, and dispersion contributions favored [host:guest] complexation, while the exchange-repulsion term showed the opposite. This work provides an in-depth understanding of using CD inclusion complexes in OPFRs formulations.

Repurposing Drugs for Inhibition against ALDH2 via a 2D/3D Ligand-Based Similarity Search and Molecular Simulation

Aldehyde dehydrogenase-2 (ALDH2) is a crucial enzyme participating in intracellular aldehyde metabolism and is acknowledged as a potential therapeutic target for the treatment of alcohol use disorder and other addictive behaviors. Using previously reported ALDH2 inhibitors of Daidzin, CVT-10216, and CHEMBL114083 as reference molecules, here we perform a ligand-based virtual screening of world-approved drugs via 2D/3D similarity search methods, followed by the assessments of molecular docking, toxicity prediction, molecular simulation, and the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) analysis. The 2D molecular fingerprinting of ECFP4 and FCFP4 and 3D molecule-shape-based USRCAT methods show good performances in selecting compounds with a strong binding behavior with ALDH2. Three compounds of Zeaxanthin (q = 0), Troglitazone (q = 0), and Sequinavir (q = +1 e) are singled out as potential inhibitors; Zeaxanthin can only be hit via USRCAT. These drugs displayed a stronger binding strength compared to the reported potent inhibitor CVT-10216. Sarizotan (q = +1 e) and Netarsudil (q = 0/+1 e) displayed a strong binding strength with ALDH2 as well, whereas they displayed a shallow penetration into the substrate-binding tunnel of ALDH2 and could not fully occupy it. This likely left a space for substrate binding, and thus they were not ideal inhibitors. The MM-PBSA results indicate that the selected negatively charged compounds from the similarity search and Vina scoring are thermodynamically unfavorable, mainly due to electrostatic repulsion with the receptor (q = -6 e for ALDH2). The electrostatic attraction with positively charged compounds, however, yielded very strong binding results with ALDH2. These findings reveal a deficiency in the modeling of electrostatic interactions (in particular, between charged moieties) in the virtual screening via the 2D/3D similarity search and molecular docking with the Vina scoring system.

Virtual Screening of FDA-Approved Drugs for Enhanced Binding with Mitochondrial Aldehyde Dehydrogenase

Mitochondrial aldehyde dehydrogenase (ALDH2) is a potential target for the treatment of substance use disorders such as alcohol addiction. Here, we adopted computational methods of molecular dynamics (MD) simulation, docking, and molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) analysis to perform a virtual screening of FDA-approved drugs, hitting potent inhibitors against ALDH2. Using MD-derived conformations as receptors, butenafine (net charge q = +1 e) and olaparib (q = 0) were selected as promising compounds with a low toxicity and a binding strength equal to or stronger than previously reported potent inhibitors of daidzin and CVT-10216. A few negatively charged compounds were also hit from the docking with the Autodock Vina software, while the MM-PBSA analysis yielded positive binding energies (unfavorable binding) for these compounds, mainly owing to electrostatic repulsion in association with a negatively charged receptor (q = -6 e for ALDH2 plus the cofactor NAD⁺). This revealed a deficiency of the Vina scoring in dealing with strong charge-charge interactions between binding partners, due to its built-in protocol of not using atomic charges for electrostatic interactions. These observations indicated a requirement of further verification using MD and/or MM-PBSA after docking prediction. The identification of key residues for the binding implied that the receptor residues at the bottom and entrance of the substrate-binding hydrophobic tunnel were able to offer additional interactions with different inhibitors such as π-π, π-alkyl, van der Waals contacts, and polar interactions, and that the rational use of these interactions is beneficial to the design of potent inhibitors against ALDH2.

Computational Investigation of Structural Basis for Enhanced Binding of Isoflavone Analogues with Mitochondrial Aldehyde Dehydrogenase

Isoflavone compounds are potent inhibitors against mitochondrial aldehyde dehydrogenase (ALDH2) for the treatment of alcoholism and drug addiction, and an in-depth understanding of the underlying structural basis helps design new inhibitors for enhanced binding. Here, we investigated the binding poses and strengths of eight isoflavone analogues (including CVT-10216 and daidzin) with ALDH2 via computational methods of molecular docking, molecular dynamics (MD) simulation, molecular mechanics Poisson-Boltzmann surface area (MM-PBSA), steered MD, and umbrella sampling. Neither the Vina scoring of docked and MD-sampled complexes nor the nonbonded protein-inhibitor interaction energy from MD simulations is able to reproduce the relative binding strength of the inhibitors compared to experimental IC₅₀ values. Considering the solvation contribution, MM-PBSA and relatively expensive umbrella sampling yield good performance for the relative binding (free) energies. The isoflavone skeleton prefers to form π-π stacking, π-sulfur, and π-alkyl interactions with planar (Phe and Trp) or sulfur-containing (Cys and Met) residues. The enhanced inhibition of CVT-10216 originates from both end groups of the isoflavone skeleton offering strong van der Waals contacts and from the methylsulfonamide group at the 4' position by hydrogen bonding (HB) with neighboring receptor residues. These results indicate that the hydrophobic binding tunnel of ALDH2 is larger than the isoflavone skeleton in length and thus an extended hydrophobic core is likely a premise for potent inhibitors.

Thermodynamic driving forces of guest confinement in a photoswitchable cage

Photoswitchable cages that confine small guest molecules inside their cavities offer a way to control the binding/unbinding process through irradiation with light of different wavelengths. However, detailed characterization of the structural and thermodynamic consequences of photoswitching is very challenging to achieve by experiments alone. Thus, all-atom molecular dynamics (MD) simulations were carried out to gain insight into the relationship between the structure and binding affinity. Binding free energies of the B₁₂F₁₂^2- guest were obtained for all photochemically accessible forms of a photoswitchable dithienylethene (DTE) based coordination cage. The MD simulations show that successive photo-induced closure of the four individual DTE ligands that form the cage gradually decreases the binding affinity. Closure of the first ligand significantly lowers the unbinding barrier and the binding free energy, and therefore favours guest unbinding both kinetically and thermodynamically. The analysis of different enthalpy contributions to the free energy shows that binding is enthalpically unfavourable and thus is an entropy-driven process, in agreement with the experimental data. Separating the enthalpy into the contributions from electrostatic, van der Waals, and bonded interactions in the force field shows that the unfavourable binding enthalpy is due to the bonded interactions being more favourable in the dissociated state, suggesting the presence of structural strain in the bound complex. Thus, the simulations provide microscopic explanations for the experimental findings and provide a possible route towards the targeted design of switchable nanocontainers with modified binding properties.

Inclusion Complexation of Organic Micropollutants with β-Cyclodextrin

Recently, β-cyclodextrin (βCD)-based polymers with enhanced adsorption kinetics and high removal capacity of organic micropollutants (OMPs) and uptake rates have been synthesized and tested experimentally. Although the exact physical–chemical mechanisms via which these polymers capture the various types of OMPs are not yet fully understood, it is suggested that the inclusion complex formation of OMPs with βCD is very important. In this study, the inclusion complex formation of OMPs with βCD in an aqueous solution is investigated by using the well-established attach–pull–release method in force field-based molecular dynamics simulations. A representative set of OMPs is selected based on the measured occurrences in surface and ground waters and the directives published by the European Union. To characterize the formation of the inclusion complex, the binding free energies, enthalpies, and entropies are computed and compared to experimental values. It is shown that computations using the q4md-CD/GAFF/Bind3P force field combination yield binding free energies that are in reasonable agreement with the experimental results for all OMPs studied. The binding enthalpies are decomposed into the main contributing interaction types. It is shown that, for all studied OMPs, the van der Waals interactions are favorable for the inclusion complexion and the hydrogen bond formation of the guest with the solvent and βCD plays a crucial role in the binding mechanism. Our findings show that MD simulations can adequately describe the inclusion complex formation of βCD with OMPs, which is the first step toward understanding the underlying mechanisms via which the βCD-based polymers capture OMPs.

Theoretical and Experimental Studies on Inclusion Complexes of Pinostrobin and β-Cyclodextrins

Pinostrobin (PNS) belongs to the flavanone subclass of flavonoids which shows several biological activities such as anti-inflammatory, anti-cancerogenic, anti-viral and anti-oxidative effects. Similar to other flavonoids, PNS has a quite low water solubility. The purpose of this work is to improve the solubility and the biological activities of PNS by forming inclusion complexes with β-cyclodextrin (βCD) and its derivatives, heptakis-(2,6-di-O-methyl)-β-cyclodextrin (2,6-DMβCD) and (2-hydroxypropyl)-β-cyclodextrin (HPβCD). The A_L-type diagram of the phase solubility studies of PNS exhibited the formed inclusion complexes with the 1:1 molar ratio. Inclusion complexes were prepared by the freeze-drying method and were characterized by differential scanning calorimetry (DSC). Two-dimensional nuclear magnetic resonance (2D-NMR) and steered molecular dynamics (SMD) simulation revealed two different binding modes of PNS, i.e., its phenyl- (P-PNS) and chromone- (C-PNS) rings preferably inserted into the cavity of βCD derivatives whilst only one orientation of PNS, where the C-PNS ring is inside the cavity, was detected in the case of the parental βCD. All PNS/βCDs complexes had a higher dissolution rate than free PNS. Both PNS and its complexes significantly exerted a lowering effect on the IL-6 secretion in LPS-stimulated macrophages and showed a moderate cytotoxic effect against MCF-7 and HeLa cancer cell lines in vitro.

Validation and Comparison of Force Fields for Native Cyclodextrins in Aqueous Solution

Molecular dynamics simulations of native α-, β-, and γ-cyclodextrin in aqueous solution have been conducted with the goal to investigate the performance of the CHARMM36 force field, the AMBER-compatible q4md-CD force field, and five variants of the GROMOS force field. The properties analyzed are structural parameters derived from X-ray diffraction and NMR experiments as well as hydrogen bonds and hydration patterns, including hydration free enthalpies. Recent revisions of the torsional-angle parameters for carbohydrate systems within the GROMOS family of force fields lead to a significant improvement of the agreement between simulated and experimental NMR data. Therefore, we recommend using the variant 53A6_GLYC instead of 53A6 and 56A6_{CARBO_R} or 2016H66 instead of 56A6_CARBO to simulate cyclodextrins in solution. The CHARMM36 and q4md-CD force fields show a similar performance as the three recommended GROMOS parameter sets. A significant difference is the more flexible nature of the cyclodextrins modeled with the CHARMM36 and q4md-CD force fields compared to the three recommended GROMOS parameter sets.

Evaluating Force Field Performance in Thermodynamic Calculations of Cyclodextrin Host-Guest Binding: Water Models, Partial Charges, and Host Force Field Parameters

Computational prediction of noncovalent binding free energies with methods based on molecular mechanical force fields has become increasingly routine in drug discovery projects, where they promise to speed the discovery of small molecule ligands to bind targeted proteins with high affinity. Because the reliability of free energy methods still has significant room for improvement, new force fields, or modifications of existing ones, are regularly introduced with the aim of improving the accuracy of molecular simulations. However, comparatively little work has been done to systematically assess how well force fields perform, particularly in relation to the calculation of binding affinities. Hardware advances have made these calculations feasible, but comprehensive force field assessments for protein–ligand sized systems still remain costly. Here, we turn to cyclodextrin host–guest systems, which feature many hallmarks of protein–ligand binding interactions but are generally much more tractable due to their small size. We present absolute binding free energy and enthalpy calculations, using the attach-pull-release (APR) approach, on a set of 43 cyclodextrin-guest pairs for which experimental ITC data are available. The test set comprises both α- and β-cyclodextrin hosts binding a series of small organic guests, each with one of three functional groups: ammonium, alcohol, or carboxylate. Four water models are considered (TIP3P, TIP4Pew, SPC/E, and OPC), along with two partial charge assignment procedures (RESP and AM1-BCC) and two cyclodextrin host force fields. The results suggest a complex set of considerations when choosing a force field for biomolecular simulations. For example, some force field combinations clearly outperform others at the binding enthalpy calculations but not for the binding free energy. Additionally, a force field combination which we expected to be the worst performer gave the most accurate binding free energies – but the least accurate binding enthalpies. The results have implications for the development of improved force fields, and we propose this test set, and potential future elaborations of it, as a powerful validation suite to evaluate new force fields and help guide future force field development.

Free-energy patterns in inclusion complexes: the relevance of non-included moieties in the stability constants

A MD/PMF-based procedure is designed for quantification of the interaction and respective components, guiding complex formation in water between β-CD and several naphthalene derivatives, highlighting the relevance of substituents.

Binding of Pollutants to Biomolecules: A Simulation Study

A number of cases around the world have been reported where animals were found dead or dying with symptoms resembling a thiamine (vitamin B) deficiency, and for some of these, a link to pollutants has been suggested. Here, we investigate whether biomolecules involved in thiamin binding and transport could be blocked by a range of different pollutants. We used in silico docking of five compound classes (25 compounds in total) to each of five targets (prion protein, ECF-type ABC transporter, thi-box riboswitch receptor, thiamin pyrophosphokinase, and YKoF protein) and subsequently performed molecular dynamics (MD) simulations to assess the stability of the complexes. The compound classes were thiamin analogues (control), pesticides, veterinary medicines, polychlorinated biphenyls, and dioxins, all of which are prevalent in the environment to some extent. A few anthropogenic compounds were found to bind the ECF-type ABC transporter, but none binds stably to prion protein. For the riboswitch, most compounds remained in their binding pockets during 50 ns of MD simulation, indicating that RNA provides a promiscuous binding site. In both YKoF and thiamin pyrophosphokinase (TPK), most compounds remain tightly bound. However, TPK biomolecules undergo pollutant-induced conformational changes. Although most compounds are found to bind to some of these targets, a larger data set is needed along with more quantitative methods like free energy perturbation calculations before firm conclusions can be drawn. This study is in part a test bed for large-scale quantitative computational screening of interactions between biological entities and pollutant molecules.

Insights into the Packing Switching of the EphA2 Transmembrane Domain by Molecular Dynamic Simulations

Receptor tyrosine kinases play an important role in mediating cell migration and adhesion associated with various biology processes. With a single-span transmembrane domain (TMD), the activities of the receptors are regulated by the definite packing configurations of the TMDs. For the EphA2 receptor, increasing studies have been conducted to investigate the packing domains that induce its switching TMD dimerization. However, the inherent transformation mechanisms including the interrelations among the involved packing domains remain unclear. Herein, we applied multiple simulation methods to explore the underlying packing mechanisms within the EphA2 TMD dimer. Our results demonstrated that the G(540)xxxG(544) contributed to the formation of the right-handed configuration while the heptad repeat L(535)xxxG(539)xxA(542)xxxV(546)xxL(549)xxxG(553) motif together with the FFxH(559) region mediated the parallel mode. Furthermore, the FF(557) residues packing mutually as rigid riveting structures were found comparable to the heptad repeat motif in maintaining the parallel configuration. In addition, the H(559) residue associated definitely with the lower bilayer leaflet, which was proved to stabilize the parallel mode significantly. The simulations provide a full range of insights into the essential packing motifs or residues involved in the switching TMD dimer configurations, which can enrich our comprehension toward the EphA2 receptor.

Deconvoluting Protein (Un)folding Structural Ensembles Using X-Ray Scattering, Nuclear Magnetic Resonance Spectroscopy and Molecular Dynamics Simulation

The folding and unfolding of protein domains is an apparently cooperative process, but transient intermediates have been detected in some cases. Such (un)folding intermediates are challenging to investigate structurally as they are typically not long-lived and their role in the (un)folding reaction has often been questioned. One of the most well studied (un)folding pathways is that of Drosophila melanogaster Engrailed homeodomain (EnHD): this 61-residue protein forms a three helix bundle in the native state and folds via a helical intermediate. Here we used molecular dynamics simulations to derive sample conformations of EnHD in the native, intermediate, and unfolded states and selected the relevant structural clusters by comparing to small/wide angle X-ray scattering data at four different temperatures. The results are corroborated using residual dipolar couplings determined by NMR spectroscopy. Our results agree well with the previously proposed (un)folding pathway. However, they also suggest that the fully unfolded state is present at a low fraction throughout the investigated temperature interval, and that the (un)folding intermediate is highly populated at the thermal midpoint in line with the view that this intermediate can be regarded to be the denatured state under physiological conditions. Further, the combination of ensemble structural techniques with MD allows for determination of structures and populations of multiple interconverting structures in solution.

Reconciling structural and thermodynamic predictions using all-atom and coarse-grain force fields: the case of charged oligo-arginine translocation into DMPC bilayers

Using the translocation of short, charged cationic oligo-arginine peptides (mono-, di-, and triarginine) from bulk aqueous solution into model DMPC bilayers, we explore the question of the similarity of thermodynamic and structural predictions obtained from molecular dynamics simulations using all-atom and Martini coarse-grain force fields. Specifically, we estimate potentials of mean force associated with translocation using standard all-atom (CHARMM36 lipid) and polarizable and nonpolarizable Martini force fields, as well as a series of modified Martini-based parameter sets. We find that we are able to reproduce qualitative features of potentials of mean force of single amino acid side chain analogues into model bilayers. In particular, modifications of peptide-water and peptide-membrane interactions allow prediction of free energy minima at the bilayer-water interface as obtained with all-atom force fields. In the case of oligo-arginine peptides, the modified parameter sets predict interfacial free energy minima as well as free energy barriers in almost quantitative agreement with all-atom force field based simulations. Interfacial free energy minima predicted by a modified coarse-grained parameter set are -2.51, -4.28, and -5.42 for mono-, di-, and triarginine; corresponding values from all-atom simulations are -0.83, -3.33, and -3.29, respectively, all in units of kcal/mol. We found that a stronger interaction between oligo-arginine and the membrane components and a weaker interaction between oligo-arginine and water are crucial for producing such minima in PMFs using the polarizable CG model. The difference between bulk aqueous and bilayer center states predicted by the modified coarse-grain force field are 11.71, 14.14, and 16.53 kcal/mol, and those by the all-atom model are 6.94, 8.64, and 12.80 kcal/mol; those are of almost the same order of magnitude. Our simulations also demonstrate a remarkable similarity in the structural aspects of the ensemble of configurations generated using the all-atom and coarse-grain force fields. Both resolutions show that oligo-arginine peptides adopt preferential orientations as they translocate into the bilayer. The guiding theme centers on charged groups maintaining coordination with polar and charged bilayer components as well as local water. We also observe similar behaviors related with membrane deformations.

Practical Aspects of Free-Energy Calculations: A Review

Free-energy calculations in the framework of classical molecular dynamics simulations are nowadays used in a wide range of research areas including solvation thermodynamics, molecular recognition, and protein folding. The basic components of a free-energy calculation, that is, a suitable model Hamiltonian, a sampling protocol, and an estimator for the free energy, are independent of the specific application. However, the attention that one has to pay to these components depends considerably on the specific application. Here, we review six different areas of application and discuss the relative importance of the three main components to provide the reader with an organigram and to make nonexperts aware of the many pitfalls present in free energy calculations.

Cooperative Binding of Cyclodextrin Dimers to Isoflavone Analogues Elucidated by Free Energy Calculations

Dimerization of cyclodextrin (CD) molecules is an elementary step in the construction of CD-based nanostructured materials. Cooperative binding of CD cavities to guest molecules facilitates the dimerization process and, consequently, the overall stability and assembly of CD nanostructures. In the present study, all three dimerization modes (head-to-head, head-to-tail, and tail-to-tail) of β-CD molecules and their binding to three isoflavone drug analogues (puerarin, daidzin, and daidzein) were investigated in explicit water surrounding using molecular dynamics simulations. Total and individual contributions from the binding partners and solvent environment to the thermodynamics of these binding reactions are quantified in detail using free energy calculations. Cooperative drug binding to two CD cavities gives an enhanced binding strength for daidzin and daidzein, whereas for puerarin no obvious enhancement is observed. Head-to-head dimerization yields the most stable complexes for inclusion of the tested isoflavones (templates) and may be a promising building block for construction of template-stabilized CD nanostructures. Compared to the case of CD monomers, the desolvation of CD dimers and entropy changes upon complexation prove to be influential factors of cooperative binding. Our results shed light on key points of the design of CD-based supramolecular assemblies. We also show that structure-based calculation of binding thermodynamics can quantify stabilization caused by cooperative effects in building blocks of nanostructured materials.

Coarse-grained versus atomistic simulations: realistic interaction free energies for real proteins

MOTIVATION

To assess whether two proteins will interact under physiological conditions, information on the interaction free energy is needed. Statistical learning techniques and docking methods for predicting protein-protein interactions cannot quantitatively estimate binding free energies. Full atomistic molecular simulation methods do have this potential, but are completely unfeasible for large-scale applications in terms of computational cost required. Here we investigate whether applying coarse-grained (CG) molecular dynamics simulations is a viable alternative for complexes of known structure.

RESULTS

We calculate the free energy barrier with respect to the bound state based on molecular dynamics simulations using both a full atomistic and a CG force field for the TCR-pMHC complex and the MP1-p14 scaffolding complex. We find that the free energy barriers from the CG simulations are of similar accuracy as those from the full atomistic ones, while achieving a speedup of >500-fold. We also observe that extensive sampling is extremely important to obtain accurate free energy barriers, which is only within reach for the CG models. Finally, we show that the CG model preserves biological relevance of the interactions: (i) we observe a strong correlation between evolutionary likelihood of mutations and the impact on the free energy barrier with respect to the bound state; and (ii) we confirm the dominant role of the interface core in these interactions. Therefore, our results suggest that CG molecular simulations can realistically be used for the accurate prediction of protein-protein interaction strength.

AVAILABILITY AND IMPLEMENTATION

The python analysis framework and data files are available for download at http://www.ibi.vu.nl/downloads/bioinformatics-2013-btt675.tgz.