Relative Free Energy Differences Research Articles

Ionomer and Carbon Aggregate Structure in Catalyst Ink Using Coarse-Grained Molecular Dynamics Simulations

The cathode catalyst layer (CL) is one of the principal components for proton exchange membrane fuel cells. CL is typically prepared by drying a catalyst ink consisting of Pt-C (platinum on a carbon power support), ionomer, and solvents (water/alcohol mixtures). Under good solvent conditions for the polymer backbone, where solvents dissolve the backbone chains, ionomers exhibit polyelectrolyte-type behavior, forming loose ionic aggregates. The interaction between ionomer and solvents controls the ionomer conformations in the solutions, which determines the final microstructure of ionomer in CLs. The characterization of ionomer structures has been extensively studied using experimental measurements such as small-angle neutron scattering (SANS) and small-angle X-ray scattering (SAXS).1-2 They observed that hydrophobic backbones of Nafion form the cylindrical-like core of self-assembled ionomer structures while ionic groups are located at the periphery of this aggregates. Theoretically, while most of studies have primarily focused on the solid membrane state with water volume fraction significantly below 50%, ionomer and carbon aggregations in solutions are also an important phenomenon to understand the self-assembly behaviors of ionomers as key structure-forming elements in CL. In this study, ionomer and carbon aggregations in a mixture of 1-propanol (NPA) and water have been investigated using coarse-grained molecular dynamics simulations. To perform simulations of larger systems and longer time spans compared with atomistic simulations, the reduced spatial resolution of the bead representation based on the relative free energy difference was used as a coarse-grained model. The dependence of NPA content on the ionomer structures was studied by systematically changing the NPA content in the system. The self-assembly behavior of ionomers into cylindrical bundle-like aggregates was observed for all NPA content solutions (Figure 1). Ionomer aggregates showed different behavior of formation at different NPA contents. The radius of an ionomer bundle was found to become smaller as NPA content increases, which is in good agreement with the trend seen in experimental measurements.3 The effects of carbon particles on the ionomer aggregation for various ionomer concentrations will also be discussed in the presentation. References Gebel, G.; Loppinet, B.; Hara, H.; Hirasawa, E., Small-Angle Neutron Scattering Study of Polyethylene-Co-Methacrylate Ionomer Aqueous Solutions. J. Phys. Chem. B 1997, 101, 3980-3987.Jiang, S.; Xia, K.-Q.; Xu, G., Effect of Additives on Self-Assembling Behavior of Nafion in Aqueous Media. Macromolecules 2001, 34, 7783-7788.Yamaguchi, M.; Matsunaga, T.; Amemiya, K.; Ohira, A.; Hasegawa, N.; Shinohara, K.; Ando, M.; Yoshida, T., Dispersion of Rod-Like Particles of Nafion in Salt-Free Water/1-Propanol and Water/Ethanol Solutions. J. Phys. Chem. B 2014, 118, 14922-14928. Figure 1

Molecular simulation of capillary phase transitions in flexible porous materials.

We used flat-histogram sampling Monte Carlo to study capillary phase transitions in deformable adsorbent materials. Specifically, we considered a pure adsorbate fluid below its bulk critical temperature within a slit pore of variable pore width. The instantaneous pore width is dictated by a number of factors, such as adsorbate loading, reservoir pressure, fluid-wall interaction, and bare adsorbent properties. In the slit pores studied here, the bare adsorbent free energy was assumed to be biparabolic, consisting of two preferential pore configurations, namely, the narrow pore and the large pore configurations. Four distinct phases could be found in the adsorption isotherms. We found a low-pressure phase transition, driven primarily by capillary condensation/evaporation and accompanied by adsorbent deformation in response. The deformation can be a relatively small contraction/expansion as seen in elastic materials, or a large-scale structural transformation of the adsorbent. We also found a high-pressure transition driven by excluded volume effects, which tends to expand the material and thus results in a large-scale structural transformation of the adsorbent. The adsorption isotherms and osmotic free energies can be rationalized by considering the relative free energy differences between the basins of the bare adsorbent free energy.

Theoretical approach to the conformational analyses of dithiazinane, thiadiazinane and triazinane, their N-borane adducts and N–H cations

Computational methods were applied to the study of the different nitrogen and chair conformers of [1,3,5]-dithiazinanes, [1,3,5]-thiadiazinanes and [1,3,5]-triazinanes. Optimisation of heterocycles bearing N–H and N–Me were performed by B3LYP density functional and the 6-311++G** basis set methods, and their free energies were correlated. N-Borane adducts and the protonated cations of these heterocycles were also optimized. The relative free energy differences (ΔG) in gas phase of the chair conformers of eighteen heterocycles were calculated. Calculations agree with experimental data and allow to explain the steric and electronic factors that determine the preferred conformations. The distribution of the electronic density and the NBO atomic charges of the optimized structures support the experimental conformational analyses. The bonding and intramolecular interactions were analysed by using critical point search of the electronic density. They were also evaluated by exchange reduced density gradients associated to overlapping, using the PBE functional with the DZUP basis set. Proton–hydride, proton–proton, hydride–hydride, hydride–sulphur, proton–sulphur and hydride-nitrogen stabilising weak interactions were analysed.

Computational insights into the mechanism of iron carbonyl-catalyzed ethylene hydrosilylation or dehydrogenative silylation

The hydrosilylation of alkenes catalyzed by Fe(CO)5 is an intricate process involving formation of alkylsilane, vinylsilane, and alkane. Herein, we represent a computational study of this reaction to determine the underlying reaction mechanism using density functional theory (DFT) techniques. On the basis of an extensive exploration of the potential energy surfaces, the modified Chalk–Harrod mechanism was found to be competitive with the Chalk–Harrod mechanism in the stoichiometric reaction of ethylene and trimethylsilane. The source of product selectivity was predicted to be determined by the relative stability of transition state TS(4b–5b), which is for the ethylene-insertion into FeSi bond from (H)(Me3Si)Fe(CO)3(C2H4) to (H)(C2H4SiMe)Fe(CO)3, and TS(5a–7a1) for the conversion of ethylene hydrometallation product (Me3Si)Fe(CO)3(η2-H)(C2H4) to (C2H5)(Me3Si)Fe(CO)4. The relative free energy difference of 3.27kcal/mol between these two transition states gives a percentage ratio of 79:21 for C2H5SiMe3 to C2H4SiMe3, which is qualitatively in agreement with the experimental observations. The largely excessive ethylene favors not only the formation of Fe(CO)3(C2H4)(H)(SiMe3) 4b but also the release of vinylsilane from Fe(CO)3(H2)(C2H3SiMe3) or Fe(CO)4(C2H3SiMe3) via the addition of another ethylene molecular to the metal center. Alternatively, in the presence of excessive R3SiH, species Fe(CO)4(H)(SiR3) will be dominate over Fe(CO)4(C2H4). In the catalytic process, Fe(CO)4(H)(SiR3) plays the main role of active species. The experimental findings were rationalized in terms of two reaction pathways, the accessibility of which depended on the ratio of silane/alkene.

Density Functional Theory Based Study on Cis–Trans Isomerism of the Amide Bond in Homodimers of β2,3- and β3-Substituted Homoproline

Preference for a cis/trans peptide bond between residues of dipeptides formed by substituted β(2,3) (I) and β(3) (II) homoproline is investigated using density functional theory (DFT). Potential energy surfaces for monomer and linear dimers are explored at the B3LYP/6-31G(d,p) level of theory. Minimum energy conformations of the dipeptides are optimized using B3LYP, PBE1PBE, B97D, and M06-2X functionals at the 6-31G(d,p) level of basis set in both the gas phase and solvent phase. The relative free energy difference between the selected conformations is marginal. Results obtained using the functionals M06-2X and B97D on dimers of I and II, respectively, agree with experimental results. The lowest energy conformations predicted by B97D/6-31G(d,p) and M06-2X/6-31G(d,p) levels of theory show greater relative MP2 correlation energy. Dipeptides of I with hydrophilic substituents show preference for a trans peptide bond. Support for cis/trans isomerism in dimers of I with hydrophobic substituents comes from potential energy surfaces and free energy data. Although dipeptides of II with hydrophilic substituents show preference for cis peptide bond, the dipeptides with hydrophobic substituent prefer trans bond.

Analyzing machupo virus-receptor binding by molecular dynamics simulations.

In many biological applications, we would like to be able to computationally predict mutational effects on affinity in protein–protein interactions. However, many commonly used methods to predict these effects perform poorly in important test cases. In particular, the effects of multiple mutations, non alanine substitutions, and flexible loops are difficult to predict with available tools and protocols. We present here an existing method applied in a novel way to a new test case; we interrogate affinity differences resulting from mutations in a host–virus protein–protein interface. We use steered molecular dynamics (SMD) to computationally pull the machupo virus (MACV) spike glycoprotein (GP1) away from the human transferrin receptor (hTfR1). We then approximate affinity using the maximum applied force of separation and the area under the force-versus-distance curve. We find, even without the rigor and planning required for free energy calculations, that these quantities can provide novel biophysical insight into the GP1/hTfR1 interaction. First, with no prior knowledge of the system we can differentiate among wild type and mutant complexes. Moreover, we show that this simple SMD scheme correlates well with relative free energy differences computed via free energy perturbation. Second, although the static co-crystal structure shows two large hydrogen-bonding networks in the GP1/hTfR1 interface, our simulations indicate that one of them may not be important for tight binding. Third, one viral site known to be critical for infection may mark an important evolutionary suppressor site for infection-resistant hTfR1 mutants. Finally, our approach provides a framework to compare the effects of multiple mutations, individually and jointly, on protein–protein interactions.

Reconstructing the Most Probable Folding Transition Path from Replica Exchange Molecular Dynamics Simulations.

The characterization of transition pathways between long-lived states, and the identification of the corresponding transition state ensembles are useful tools in the study of rare events such as protein folding. In this work we demonstrate how the most probable transition path between metastable states can be recovered from replica exchange molecular dynamic simulation data by using the dynamic string method. The local drift vector in collective variables is determined via short continuous trajectories between replica exchanges at a given temperature, and points along the string are updated based on this drift vector to produce reaction pathways between the folded and unfolded state. The method is applied to a designed beta hairpin-forming peptide to obtain information on the folding mechanism and transition state using different sets of collective variables at various temperatures. Two main folding pathways differing in the order of events are found and discussed, and the relative free energy differences for each path estimated. Finally, the structures near the transition state are found and described.

Reconstruction of the Most Probable Folding Transition Path from All-Atom Replica Exchange Simulations, using the Dynamic String Method

In this work we demonstrate how the most probable transition path between metastable states can be recovered from replica exchange molecular dynamic simulation data by using the dynamic string method. The local drift vector in collective variables is determined via the short continuous trajectories between replica exchanges at a given temperature. A string is updated based on this drift vector to produce reaction pathways between the folded and unfolded state. The method is applied to a beta hairpin-forming peptide to obtain information on the folding mechanism and transition state using different sets of collective variables. Two main folding pathways differing in the order of events are found and discussed, and the relative free energy differences for each path estimated. The string that shows turn formation before native arrangement of the hydrophobic core has lower free energy difference from the unfolded state and thus is expected to be more probable. Finally, in both cases, the conformations on the peak of the free energy barrier have the tryptophans on the same side of the beta structure suggesting this arrangement to be the rate-limiting step.This work was supported by the NSF MCB-1050966 and CCNI.

Retaining the Self Interactions in Alchemical Free Energy Calculations

Alchemical transformations of chemical species into others provide a common numerical method to compute free energy differences. This method is particularly well suited to compute relative free energy differences formulated as thermodynamic cycles. Moreover, it allows us to compute experimentally measurable quantities, such as the relative binding free energy of small molecules to proteins.The drawback of the alchemical methods is that the annihilation/creation processes of mutated particles present sampling issues and singularities that reduce the accuracy of the calculation. It is therefore important to understand whether the process of annihilation can be simplified and perhaps even avoided.We report a rigorous proof that it is possible to find an alchemical pathway from the native to the mutant system, in which “self interactions” (bonded and non-bonded) are retained without affecting the relative free energy difference. We also provide a molecular dynamics illustrating of this proof in which a cycles of solvated side chains is considered.

Significance of thermodynamic and physical characteristics on permeation of ions during membrane separation: Hydrated radius, hydration free energy and viscous effects

Selected thermodynamic and physical characteristics of ions were compared in relation to their permeation ability during membrane filtration. The characteristics analyzed included ionic radius, hydrated radius, hydration free energy, Jones–Dole viscosity coefficient (B). Correlations between ionic parameters were analyzed to explain the ability of some ions to permeate by rearranging the water molecules (i.e., soft ions) within their hydration shells during transport as well as near membrane inter-ionic interactions. Hydrated size and strength of the hydration shell significantly affect ability of ions to permeate through the thin film membranes. Selectivity of membrane processes depend on the hydration free energy of ions. Ion selectivity of a membrane (or partitioning of ions on the membrane) is governed by differences in relative free energies between the solution and membrane phases. The Jones–Dole coefficient (B) for both cations and anions show strong correlations with selected thermodynamic characteristics including entropy of hydration and hydration free energy. The order of ions with respect to selected parameters show similarities with the Hofmeister series which is often used to explain the ability of ions to interact with hydrogen in proteins. The same ability (i.e., interaction with hydrogen) affects the behavior of ions both near and within membrane, allowing some ions to permeate either by size exclusion or ability to rearrange the hydration waters under shear forces during transport.

Towards Finding a Better Reaction Coordinate: NtrC and DIMS

All atom molecular dynamics is appealing due to the detailed insights towards protein function. But, all-atom molecular dynamics is limited due to the small amount of time spent in the transitions between stable states. Dynamic importance sampling (DIMS) enhances this sampling problem by biasing, with correction, to encourage sampling away from the stable states and to thus improve estimates of barrier crossing times. Nevertheless, DIMS has limitations based on the choice of order parameters that define intermediate states. To improve the representation, in DIMS, of order parameters, we have performed simulations of NtrC and have analyzed the transition populations. Our analysis, using techniques of nonlinear time series analysis and nonlinear dimensionality reduction, has enabled us to determine a new numerical method for finding reduced parameter descriptions of transitions in complex biomolecular systems. The results should enable an iterative approach to finding the best order parameter and at the same time improve convergence for estimates of reaction pathways, kinetics and relative free energy differences.

Efficiency of alchemical free energy simulations. I. A practical comparison of the exponential formula, thermodynamic integration, and Bennett's acceptance ratio method

We investigate the relative efficiency of thermodynamic integration, three variants of the exponential formula, also referred to as thermodynamic perturbation, and Bennett's acceptance ratio method to compute relative and absolute solvation free energy differences. Our primary goal is the development of efficient protocols that are robust in practice. We focus on minimizing the number of unphysical intermediate states (λ-states) required for the computation of accurate and precise free energy differences. Several indicators are presented which help decide when additional λ-states are necessary. In all tests Bennett's acceptance ratio method required the least number of λ-states, closely followed by the "double-wide" variant of the exponential formula. Use of the exponential formula in only strict "forward" or "backward" mode was not found to be competitive. Similarly, the performance of thermodynamic integration in terms of efficiency was rather poor. We show that this is caused by the use of the trapezoidal rule as method of numerical quadrature. A systematic study focusing on the optimization of thermodynamic integration is presented in a companion paper.

Prediction of Relative Stability between TACE/Gelastatin and TACE/Gelastatin Hydroxamate

A gelastatins (1), natural MMP inhibitors, and their hydroxamate analogues (2) in TACE enzyme evaluated for discovery of potent TACE inhibitors. We have employed molecular dynamics simulations to compute the relative free energy of hydration and binding to TACE for gelastatin (1) and its hydroxamate analogue (2). The relative free energy difference is directly described in this article using the free energy perturbation approach as a means to accurately predict the TACE inhibitor of gelastatin analogues. The results show that the good agreement between the experimental and theoretical relative free energies of binding, gelastatin hydroxamate (2) binds stronger to TACE by -3.37 kcal/mol. The desolvation energy costs significantly reduced binding affinity, hydroxamate group associated with high desolvation energy formed strong favorable interactions with TACE with more than compensated for the solvation costs and therefore led to an improvement in relative binding affinity.

The severity of osteogenesis imperfecta: A comparison to the relative free energy differences of collagen model peptides

Molecular dynamics simulations were carried out to calculate free energy differences between the folded and unfolded states of wild type and mutant collagen model peptides. The calculated stability of the collagen models was compared with the severity of osteogenesis imperfecta. Free energy differences of Gly → Xaa (Xaa: Ser, Cys, Glu, and Asp) mutations between the wild type and the mutants at position 15 of the model peptide were 3.8, 4.2, 5.6, and 8.8 kcal/mol, respectively. The corresponding free energy differences of a second Gly mutation at the same position in different chains were, on average, 1.3, 1.5, 2.9, and 5.4 kcal/mol, respectively. Free energy simulations were also performed to estimate the relative stability between an oxidized form and a reduced form of the mutants containing two Cys residues, which indicated that the mutant of the collagen-like peptide containing an intramolecular disulfide bond was more stable than the mutant containing one Cys residue but less stable than the wild type. The calculated free energy differences between an oxidized and a reduced form of the mutants containing two Cys residues are 0.8 and 2.6 kcal/mol for the disulfide bonds between Chains A and B and between Chains A and C, respectively.

Effect of Single Central Base Pair Mismatch on the Conformation and Stability of Surface Immobilized DNA: Molecular Dynamics Studies

We combine molecular dynamics simulations and relative free energy calculation to investigate how single central base mismatch influence the conformation and stability of surface immobilized DNA. For comparison purpose, molecular dynamics simulations have been performed on a central mismatched duplex and the corresponding perfectly matched duplex, respectively. Simulations have been conducted with the generalized Born model. The relative free energy differences of the two systems suggest the single base-pair mismatch may lower the stability of DNA as expectation. The result obtained from the structural analysis shows the average conformations of the two DNA duplexes resemble the B-form, which suggests the central base pair mismatch has little effect on the overall conformation of DNA. Furthermore, our result shows the fluctuation of the perfectly matched duplex is obviously larger than that of the central base pair mismatched duplex. This study gives us a new atomic-level insight into DNA probe-target interaction on the surface. The results presented here will also aid in the design of DNA microarrays for single nucleotide polymorphisms screening.

Elevated polyamines induce c-MYC overexpression by perturbing quadruplex–WC duplex equilibrium

The biological role of quadruplexes and polyamines has been independently associated with cancer. However, quadruplex-polyamine mediated transcriptional regulation remain unaddressed. Herein, using c-MYC quadruplex model, we have attempted to understand quadruplex–polyamine interaction and its role in transcriptional regulation. We initially employed biophysical approach involving CD, UV and FRET to understand the role of polyamines (spermidine and spermine) on conformation, stability, molecular recognition of quadruplex and to investigate the effect of polyamines on quadruplex–Watson Crick duplex transition. Our study demonstrates that polyamines affect the c-MYC quadruplex conformation, perturb its recognition properties and delays duplex formation. The relative free energy difference (ΔΔG°) between the duplex and quadruplex structure indicate that polyamines stabilize and favor c-MYC quadruplex over duplex. Further, we investigated the influence of polyamine mediated perturbation of this equilibrium on c-MYC expression. Our results suggest that polyamines induce structural transition of c-MYC quadruplex to a transcriptionally active motif with distinctive molecular recognition property, which drives c-MYC expression. These findings may allow exploiting quadruplex–polyamines interaction for developing antiproliferative strategies to combat aberrant gene expression.

Graph-theoretical identification of dissociation pathways on free energy landscapes of biomolecular interaction

Biomolecular association and dissociation reactions take place on complicated interaction free energy landscapes that are still very hard to characterize computationally. For large enough distances, though, it often suffices to consider the six relative translational and rotational degrees of freedom of the two particles treated as rigid bodies. Here, we computed the six-dimensional free energy surface of a dimer of water-soluble alpha-helices by scanning these six degrees of freedom in about one million grid points. In each point, the relative free energy difference was computed as the sum of the polar and nonpolar solvation free energies of the helix dimer and of the intermolecular coulombic interaction energy. The Dijkstra graph algorithm was then applied to search for the lowest cost dissociation pathways based on a weighted, directed graph, where the vertices represent the grid points, the edges connect the grid points and their neighbors, and the weights are the reaction costs between adjacent pairs of grid points. As an example, the configuration of the bound state was chosen as the source node, and the eight corners of the translational cube were chosen as the destination nodes. With the strong electrostatic interaction of the two helices giving rise to a clearly funnel-shaped energy landscape, the eight lowest-energy cost pathways coming from different orientations converge into a well-defined pathway for association. We believe that the methodology presented here will prove useful for identifying low-energy association and dissociation pathways in future studies of complicated free energy landscapes for biomolecular interaction.

Computational mutagenesis reveals the role of active-site tyrosine in stabilising a boat conformation for the substrate: QM/MM molecular dynamics studies of wild-type and mutant xylanases

Molecular dynamics simulations have been performed for non-covalent complexes of phenyl beta-xylobioside with the retaining endo-beta-1,4-xylanase from B. circulans (BCX) and its Tyr69Phe mutant using a hybrid QM/MM methodology. A trajectory initiated for the wild-type enzyme-substrate complex with the proximal xylose ring bound at the -1 subsite (adjacent to the scissile glycosidic bond) in the (4)C(1) chair conformation shows spontaneous transformation to the (2,5)B boat conformation, and potential of mean force calculations indicate that the boat is approximately 30 kJ mol(-1) lower in free energy than the chair. Analogous simulations for the mutant lacking one oxygen atom confirm the key role of Tyr69 in stabilizing the boat in preference to the (4)C(1) chair conformation, with a relative free energy difference of about 20 kJ mol(-1), by donating a hydrogen bond to the endocyclic oxygen of the proximal xylose ring. QM/MM MD simulations for phenyl beta-xyloside in water, with and without a propionate/propionic acid pair to mimic the catalytic glutamate/glutamic acid pair of the enzyme, show the (4)C(1) chair to be stable, although a hydrogen bond between the OH group at C2 of xylose and the propionate moiety seems to provide some stabilization for the (2,5)B conformation.

Effect of loop length variation on quadruplex-Watson Crick duplex competition

The effect of loop length on quadruplex stability has been studied when the G-rich strand is present along with its complementary C-rich strand, thereby resulting in competition between quadruplex and duplex structures. Using model sequences with loop lengths varying from T to T5, we carried out extensive FRET to discover the influence of loop length on the quadruplex-Watson Crick duplex competition. The binding data show an increase in the binding affinity of quadruplexes towards their complementary strands upon increasing the loop length. Our kinetic data reveal that unfolding of the quadruplex in presence of a complementary strand involves a contribution from a predominant slow and a small population of fast opening conformer. The contribution from the fast opening conformer increases upon increasing the loop length leading to faster duplex formation. FCS data show an increase in the interconversion between the quadruplex conformers in presence of the complementary strand, which shifts the equilibrium towards the fast opening conformer with an increase in loop length. The relative free-energy difference (ΔΔG°) between the duplex and quadruplex indicates that an increase in loop length favors duplex formation and out competes the quadruplex.

Role of Locked Nucleic Acid Modified Complementary Strand in Quadruplex/Watson−Crick Duplex Equilibrium

In the human genome, the G-rich sequences that form quadruplexes are present along with their C-rich complementary strands; this suggests the existence of equilibrium between a quadruplex and a Watson-Crick duplex which allows the execution of their respective biological functions. We have investigated the sensitivity of this equilibrium to pharmacological agents by employing locked nucleic acid (LNA) modified complementary strands, and demonstrated successful invasion of the stable telomeric quadruplex d[(G(3)TTA)(3)G(3)]. Fluorescence, UV, ITC, and SPR studies were performed to understand the binding process involving the preformed quadruplex and LNA-modified complementary strands compared with that involving the unmodified complementary strand. Our data indicate that LNA modifications in the complementary strand shift the equilibrium toward the duplex state. These modifications confer increased thermodynamic stability to the duplex and increase the magnitude of relative free energy (DeltaDeltaG degrees) difference between duplex and quadruplex, thus favoring the predominance of duplex population over quadruplex. This superior ability of LNA-modified complementary strand can be exploited to pave an exploratory approach in which it hybridizes to a telomeric quadruplex and drives duplex formation, and inhibits the recognition of 3' G-rich overhang by RNA template of telomerase which guides telomere extension.