Molecular dynamics simulations of electroosmosis in perfluorosulfonic acid polymer

An atomistic MD simulation method has been developed to study the electroosmotic drag in the hydrated perfluorosulfonic acid polymer. The transport characteristics of the hydroniums and water molecules are evaluated from their velocity distribution functions with an electric field applied. It is shown that the microstructure of the hydrated perfluorosulfonic acid polymer is not perturbed significantly by the electric field up to 2 V/microm, and the velocity distribution functions obey the peak shifted Maxwell velocity distribution functions. The evaluated peak shifting velocities are only about 1% of the average thermal motion. The hydronium flow and water flow are evaluated from the average transport velocities or the peak shifting velocities. The electroosmotic drag coefficients from the MD simulations are in good correspondence with the experimental values. It is also shown that the electroosmotic drag coefficient has no or weak temperature dependence.     

Anisotropic diffusion of water in perfluorosulfonic acid membrane and hydrocarbon membranes

Polymer electrolyte membranes that are applied for polymer electrolyte fuel cell (PEFC) retain water in their three-dimensional network structure. Diffusion behavior of water in the membranes was analyzed by pulsed field gradient (PFG)-NMR method to estimate diffusion coefficient of proton species as water or hydronium ion. The membrane samples were put in a sample tube vertically or horizontally toward to the field gradient axis under determined temperature and humidity conditions. As the results, anisotropic diffusion behavior of water in the membranes was indicated. Anisotropic properties depended on the sample type, preparation conditions of the wet membranes, and measurement conditions. A perfluorosulfonic acid membrane tended to have smaller anisotropy while hydrocarbon membranes showed greater anisotropy.     

Comparison of kinetic Monte Carlo and molecular dynamics simulations of diffusion in a model glass former

We present results from kinetic Monte Carlo (KMC) simulations of diffusion in a model glass former. We find that the diffusion constants obtained from KMC simulations have Arrhenius temperature dependence, while the correct behavior, obtained from molecular dynamics simulations, can be super-Arrhenius. We conclude that the discrepancy is due to undersampling of higher-lying local minima in the KMC runs. We suggest that the relevant connectivity of minima on the potential energy surface is proportional to the energy density of the local minima, which determines the "inherent structure entropy." The changing connectivity with potential energy may produce a correlation between dynamics and thermodynamics.     

On the consequences of side chain flexibility and backbone conformation on hydration and proton dissociation in perfluorosulfonic acid membranes

The flexibility of the side chain and effects of conformational changes in the backbone on hydration and proton transfer in the short-side-chain (SSC) perfluorosulfonic acid fuel cell membrane have been investigated through first principles based molecular modelling studies. Potential energy profiles determined at the B3LYP/6-31G(d,p) level in the two pendant side chain fragments: CF(3)CF(-O(CF(2))(2)SO(3)H)-(CF(2))(7)-CF(-O(CF(2))(2)SO(3)H)CF(3) indicate that the largest CF(2)-CF(2) rotational barrier along the backbone is nearly 28.9 kJ mol(-1) higher than the minimum energy staggered trans conformation. Furthermore, the calculations reveal that the stiffest portion of the side chain is near to its attachment site on the backbone, with CF-O and O-CF(2) barriers of 38.1 and 28.0 kJ mol(-1), respectively. The most flexible portion of the side chain is the carbon-sulfur bond, with a barrier of only 8.8 kJ mol(-1). Extensive searches for minimum energy structures (at the B3LYP/6-311G(d,p) level) of the same polymeric fragment with 4-7 explicit water molecules reveal that the perfluorocarbon backbone may adopt either an elongated geometry, with all carbons in a trans configuration, or a folded conformation as a result of the hydrogen bonding of the terminal sulfonic acids with the water. These electronic structure calculations show that the fragments displaying the latter 'kinked' backbone possessed stronger binding of the water to the sulfonic acid groups, and also undergo proton dissociation with fewer water molecules. The calculations point to the importance of the flexibility in both the backbone and side chains of PFSA membranes to effectively transport protons under low humidity conditions.     

Application of molecular dynamics simulations in molecular property prediction II: diffusion coefficient

In this work, we have evaluated how well the general assisted model building with energy refinement (AMBER) force field performs in studying the dynamic properties of liquids. Diffusion coefficients (D) have been predicted for 17 solvents, five organic compounds in aqueous solutions, four proteins in aqueous solutions, and nine organic compounds in nonaqueous solutions. An efficient sampling strategy has been proposed and tested in the calculation of the diffusion coefficients of solutes in solutions. There are two major findings of this study. First of all, the diffusion coefficients of organic solutes in aqueous solution can be well predicted: the average unsigned errors and the root mean square errors are 0.137 and 0.171 × 10(-5) cm(-2) s(-1), respectively. Second, although the absolute values of D cannot be predicted, good correlations have been achieved for eight organic solvents with experimental data (R(2) = 0.784), four proteins in aqueous solutions (R(2) = 0.996), and nine organic compounds in nonaqueous solutions (R(2) = 0.834). The temperature dependent behaviors of three solvents, namely, TIP3P water, dimethyl sulfoxide, and cyclohexane have been studied. The major molecular dynamics (MD) settings, such as the sizes of simulation boxes and with/without wrapping the coordinates of MD snapshots into the primary simulation boxes have been explored. We have concluded that our sampling strategy that averaging the mean square displacement collected in multiple short-MD simulations is efficient in predicting diffusion coefficients of solutes at infinite dilution.     

Dynamic mechanism of fatty acid transport across cellular membranes through FadL: molecular dynamics simulations

FadL is an important member of the family of fatty acid transport proteins within membranes. In this study, 11 conventional molecular dynamics (CMD) and 25 steered molecular dynamics (SMD) simulations were performed to investigate the dynamic mechanism of transport of long-chain fatty acids (LCFAs) across FadL. The CMD simulations addressed the intrinsically dynamic behavior of FadL. Both the CMD and SMD simulations revealed that a fatty acid molecule can move diffusively to a high-affinity site (HAS) from a low-affinity site (LAS). During this process, the swing motion of the L3 segment and the hydrophobic interaction between the fatty acid and FadL could play important roles. Furthermore, 22 of the SMD simulations revealed that fatty acids can pass through the gap between the hatch domain and the transmembrane domain (TMD) by different pathways. SMD simulations identified nine possible pathways for dodecanoic acid (DA) threading the barrel of FadL. The binding free energy profiles between DA and FadL along the MD trajectories indicate that all of the possible pathways are energetically favorable for the transport of fatty acids; however, one pathway (path VI) might be the most probable pathway for DA transport. The reasonability and reliability of this study were further demonstrated by correlating the MD simulation results with the available mutagenesis results. On the basis of the simulations, a mechanism for the full-length transport process of DA from the extracellular side to the periplasmic space mediated by FadL is proposed.     

Mesoscale simulation of morphology in hydrated perfluorosulfonic acid membranes

Current fuel cell proton exchange membranes rely on a random network of conducting hydrophilic domains to transport protons across the membrane. Despite extensive investigation, details of the structure of the hydrophilic domains in these membranes remain unresolved. In this study a dynamic self-consistent mean field theory has been applied to obtain the morphologies of hydrated perfluorosulfonic acid membranes (equivalent weight of 1100) as a model system for Nafion at several water contents. A coarse-grained mesoscale model was developed by dividing the system into three components: backbone, side chain, and water. The interaction parameters for this model were generated using classical molecular dynamics. The simulated morphology shows phase separated micelles filled with water, surrounded by side chains containing sulfonic groups, and embedded in the fluorocarbon matrix. The size distribution and connectivity of the hydrophilic domains were analyzed and the small angle neutron scattering (SANS) pattern was calculated. At low water content (lambda<6, where lambda is the number of water molecules per sulfonic group) the isolated domains obtained from simulation are nearly spherical with a domain size smaller than that fitted to experimental SANS data. At higher water content (lambda>8), the domains deform into elliptical and barbell shapes as they merge. The simulated morphology, hydrophilic domain size and shape are generally consistent with some experimental observations.     

Predicting the proton conductivity of perfluorosulfonic acid membrane via combining statistical thermodynamics and molecular dynamics simulation

The electrochemical properties of a perfluorosulfonic acid (PFSA) membrane are estimated using a combination of molecular dynamics simulation and statistical thermodynamic model. We obtain all parameters in an ionic conductivity model from an atomistic simulation and remove all adjusted model parameters. From a microscopic point of view, the hydrated PFSA membrane shows micro‐phase segregation which separated into hydrophilic and hydrophobic phases. Our present work originates with this phenomenon and we treat this phase segregation as if it is a continuous phase for each of which the proton (H⁺) is transported inside the PFSA membrane/solvent (water and alcohols) mixture. The chemical potential for a given system is estimated using a molecular simulation technique to predict the van der Waals interaction energy between the polymer and solvent. In addition, the self diffusion coefficients are calculated from the molecular dynamics simulation. We study various polymer/solvent compositions to understand the concentration dependence of self diffusion coefficient. Our self diffusion coefficients and also the predicted final ionic conductivity agree well with previously reported experimental data. © 2011 Wiley Periodicals, Inc. J Polym Sci Part B: Polym Phys 49: 1455–1463, 2011     

Langevin dynamics simulations of the diffusion of molecular knots in tensioned polymer chains

Motivated by recent experiments, in which knots have been tied in individual biopolymer molecules, we use Langevin dynamics simulations to study the diffusion of a knot along a tensioned polymer chain. We find that the dependence of the knot diffusion coefficient on the tension can be non-monotonic. This behavior can be explained by the model, in which the motion of the knot involves cooperative displacement of a local knot region. At low tension, the overall viscous drag force that acts on the knot region is proportional to the number N of monomers that participate in the knot, which decreases as the tension is increased, leading to faster diffusion. At high tension the knot becomes tight and its dynamics are dominated by the chain's internal friction, which increases with the increasing tension, thereby slowing down the knot diffusion. This model is further supported by the observation that the knot diffusion coefficient measured across a set of different knot types is inversely proportional to N. We propose that the lack of tension dependence of the knot diffusion coefficients measured in recent experiments is due to the fact that the experimental values of the tension are close to the turnover between the high- and low-force regimes.     

Predicting hydration free energies using all-atom molecular dynamics simulations and multiple starting conformations

Molecular dynamics simulations in explicit solvent were applied to predict the hydration free energies for 23 small organic molecules in blind SAMPL2 test. We found good agreement with experimental results, with an RMS error of 2.82 kcal/mol over the whole set and 1.86 kcal/mol over all the molecules except several hydroxyl-rich compounds where we find evidence for a systematic error in the force field. We tested two different solvent models, TIP3P and TIP4P-Ew, and obtained very similar hydration free energies for these two models; the RMS difference was 0.64 kcal/mol. We found that preferred conformation of the carboxylic acids in water differs from that in vacuum. Surprisingly, this conformational change is not adequately sampled on simulation timescales, so we apply an umbrella sampling technique to include free energies associated with the conformational change. Overall, the results of this test reveal that the force field parameters for some groups of molecules (such as hydroxyl-rich compounds) still need to be improved, but for most compounds, accuracy was consistent with that seen in our previous tests.     

Refinement of the primary hydration shell model for molecular dynamics simulations of large proteins

A realistic representation of water molecules is important in molecular dynamics simulation of proteins. However, the standard method of solvating biomolecules, that is, immersing them in a box of water with periodic boundary conditions, is computationally expensive. The primary hydration shell (PHS) method, developed more than a decade ago and implemented in CHARMM, uses only a thin shell of water around the system of interest, and so greatly reduces the computational cost of simulations. Applying the PHS method, especially to larger proteins, revealed that further optimization and a partial reworking was required and here we present several improvements to its performance. The model is applied to systems with different sizes, and both water and protein behaviors are compared with those observed in standard simulations with periodic boundary conditions and, in some cases, with experimental data. The advantages of the modified PHS method over its original implementation are clearly apparent when it is applied to simulating the 82 kDa protein Malate Synthase G. © 2009 Wiley Periodicals, Inc. J Comput Chem 2009     

The SAMPL4 hydration challenge: evaluation of partial charge sets with explicit-water molecular dynamics simulations

We used blind predictions of the 47 hydration free energies in the SAMPL4 challenge to test multiple partial charge models in the context of explicit solvent free energy simulations with the general AMBER force field. One of the partial charge models, IPolQ-Mod, is a fast continuum solvent-based implementation of the IPolQ approach. The AM1-BCC, restrained electrostatic potential (RESP) and IpolQ-Mod approaches all perform reasonably well (R² > 0.8), while VCharge, though faster, gives less accurate results (R² of 0.5). The AM1-BCC results are more accurate than those of RESP for tertiary amines and nitrates, but the overall difference in accuracy between these methods is not statistically significant. Interestingly, the IPolQ-Mod method is found to yield partial charges in very close agreement with RESP. This observation suggests that the success of RESP may be attributed to its fortuitously approximating the arguably more rigorous IPolQ approach.     

Effect of acetone accumulation on structure and dynamics of lipid membranes studied by molecular dynamics simulations

The modulation of the properties and function of cell membranes by small volatile substances is important for many biomedical applications. Despite available experimental results, molecular mechanisms of action of inhalants and organic solvents, such as acetone, on lipid membranes remain not well understood. To gain a better understanding of how acetone interacts with membranes, we have performed a series of molecular dynamics (MD) simulations of a POPC bilayer in aqueous solution in the presence of acetone, whose concentration was varied from 2.8 to 11.2 mol%. The MD simulations of passive distribution of acetone between a bulk water phase and a lipid bilayer show that acetone favors partitioning into the water-free region of the bilayer, located near the carbonyl groups of the phospholipids and at the beginning of the hydrocarbon core of the lipid membrane. Using MD umbrella sampling, we found that the permeability barrier of ∼0.5 kcal/mol exists for acetone partitioning into the membrane. In addition, a Gibbs free energy profile of the acetone penetration across a bilayer demonstrates a favorable potential energy well of −3.6 kcal/mol, located at 15–16 Å from the bilayer center. The analysis of the structural and dynamics properties of the model membrane revealed that the POPC bilayer can tolerate the presence of acetone in the concentration range of 2.8–5.6 mol%. The accumulation of the higher acetone concentration of 11.2 mol% results, however, in drastic disordering of phospholipid packing and the increase in the membrane fluidity. The acetone molecules push the lipid heads apart and, hence, act as spacers in the headgroup region. This effect leads to the increase in the average headgroup area per molecule. In addition, the acyl tail region of the membrane also becomes less dense. We suggest, therefore, that the molecular mechanism of acetone action on the phospholipid bilayer has many common features with the effects of short chain alcohols, DMSO, and chloroform.     

Mechanical properties of surfactant bilayer membranes from atomistic and coarse-grained molecular dynamics simulations

We use simulations to predict the stability and mechanical properties of two amphiphilic bilayer membranes. We carry out atomistic MD simulations and investigate whether it is possible to use an existing coarse-grained (CG) surfactant model to map the membrane properties. We find that certain membranes can be represented well by the CG model, whereas others cannot. Atomistic MD simulations of the erucate membrane yield a headgroup area per surfactant a(0) of 0.26 nm(2), an elastic modulus K(A) of 1.7 N/m, and a bending rigidity kappa of 5 k(B)T. We find that the CG model, with the right choice for the size and potential well depth of the head, correctly reproduces a(0), kappa, as well as the fluctuation spectrum over the whole range of q values. Atomistic MD simulations of EHAC, on the other hand, suggest that this membrane is unstable. This is indicated by the fact that kappa is of the order of k(B)T, which means that the interface is extremely flexible and diffuse, and K(A) is close to zero, which means that the surface tension is zero. We argue that the CG model can be used if the headgroups are uncharged, dipolar, or effectively dipolar due to headgroup charge screening induced by counterion condensation.     

FMO-MD simulations on the hydration of formaldehyde in water solution with constraint dynamics

Full-quantum mechanical fragment molecular orbital-based molecular dynamics (FMO-MD) simulations were applied to the hydration reaction of formaldehyde in water solution under neutral conditions. Two mechanisms, a concerted and a stepwise one, were considered with respect to the nucleophilic addition and the proton transfer. Preliminary molecular orbital calculations by means of polarized continuum model reaction field predicted that the hydration prefers a concerted mechanism. Because the calculated activation barriers were too high for free FMO-MD simulations to give reactive trajectories spontaneously, a More O'Ferrall-Jencks-type diagram was constructed from the statistical analysis of the FMO-MD simulations with constraint dynamics. The diagram showed that the hydration proceeds through a zwitterionic-like (ZW-like) structure. The free energy changes along the reaction coordinate calculated by means of the blue moon ensemble for the hydration and the amination of formaldehyde indicated that the hydration proceeds through a concerted process through the ZW-like structure, whereas the amination goes through a stepwise mechanism with a ZW intermediate. In inspection of the FMO-MD trajectories, water-mediated cyclic proton transfers were observed in both reactions on the way from the ZW-like structure to the product. These proton transfers also have an asynchronous character, in which deprotonation from the nucleophilic oxygen atom (or nitrogen atom for amination) precedes the protonation of the carbonyl oxygen atom. The results showed the strong advantage of the FMO-MD simulations to obtain detailed information at a molecular level for solution reactions.     

Potential of hydrogen bond in water. Comparison of the theory with vibrational spectra and results of molecular dynamics simulations

Potential of hydrogen bond is the function which relates its energy to geometrical parameters of hydrogen bridge: its length R(O…O) and angles between direction O…O and OH group [φ (H-O…O)] and/or lone pair of proton accepting oxygen atom [χ(-O…O)]. Previously we have suggested an approach to design such potentials based on the approximate numerical solution of a reverse problem of the spectrum band shape in the frames of the fluctuation theory of hydrogen bonding. In the given work this method is applied to construction of the two-parameter potentials that depend on parameters {R(O…O), φ (H-O…O} or {φ (H-O…O), χ (-O…O)}. Using them, the spectra of OH vibrations of HOD molecules in a liquid phase are calculated theoretically in good agreement with experiment in the temperature range up to 200 °C. Distributions of angles P(φ, T), P(χ, T), and energies P(E) are calculated also. The same distributions and spectra are independently calculated on the basis of the geometrical parameters of the hydrogen bridges obtained from molecular dynamics models of water. The shapes of the spectral contours and their temperature evolution calculated for computer models turned out to be similar to experimental ones only for the potential that includes the length of H-bond R(O…O).     

Dual function of the hydration layer around an antifreeze protein revealed by atomistic molecular dynamics simulations

Atomistic molecular dynamics simulations are used to investigate the mechanism by which the antifreeze protein from the spruce budworm, Choristoneura fumiferana, binds to ice. Comparison of structural and dynamic properties of the water around the three faces of the triangular prism-shaped protein in aqueous solution reveals that at low temperature the water structure is ordered and the dynamics slowed down around the ice-binding face of the protein, with a disordering effect observed around the other two faces. These results suggest a dual role for the solvation water around the protein. The preconfigured solvation shell around the ice-binding face is involved in the initial recognition and binding of the antifreeze protein to ice by lowering the barrier for binding and consolidation of the protein:ice interaction surface. Thus, the antifreeze protein can bind to the molecularly rough ice surface by becoming actively involved in the formation of its own binding site. Also, the disruption of water structure around the rest of the protein helps prevent the adsorbed protein becoming covered by further ice growth.     

Prediction of hydration free energies for aliphatic and aromatic chloro derivatives using molecular dynamics simulations with the OPLS-AA force field

All-atom molecular dynamics computer simulations were used to blindly predict the hydration free energies of a range of chloro-organic compounds as part of the SAMPL3 challenge. All compounds were parameterized within the framework of the OPLS-AA force field, using an established protocol to compute the absolute hydration free energy via a windowed free energy perturbation approach and thermodynamic integration. Three different approaches to deriving partial charge parameters were pursued: (1) using existing OPLS-AA atom types and charges with minor adjustments of partial charges on equivalent connecting atoms; (2) calculation of quantum mechanical charges via geometry optimization, followed by electrostatic potential (ESP) fitting, using Jaguar at the LMP2/cc-pVTZ(-F) level; and (3) via geometry optimization and CHelpG charges (Gaussian03 at the HF/6-31G* level), followed by two-stage RESP fitting. Protocol 3 generated the most accurate predictions with a root mean square (RMS) error of 1.2 kcal mol(-1) for the entire data set. It was found that the deficiency of the standard OPLS-AA parameters, protocol 1 (RMS error 2.4 kcal mol(-1) overall), was mostly due to compounds with more than three chlorine substituents on an aromatic ring. For this latter subset, the RMS errors were 1.4 kcal mol(-1) (protocol 3) and 4.3 kcal mol(-1) (protocol 1), respectively. We propose new OPLS-AA atom types for aromatic carbon and chlorine atoms in rings with ≥4 Cl-substituents that perform better than the best QM-based approach, resulting in an RMS error of 1.2 kcal mol(-1) for these difficult compounds.     

Long-timescale simulations of diffusion in molecular solids

Kinetic processes play a crucial role in the formation and evolution of molecular layers. In this perspective we argue that adaptive kinetic Monte Carlo is a powerful simulation technique for determining key kinetic processes in molecular solids. The applicability of the method is demonstrated by simulating the diffusion of a CO admolecule on a water ice surface, which is an important process for the formation of organic compounds on interstellar dust grains. CO diffusion is found to follow Arrhenius behavior and the corresponding effective activation energy for diffusion is determined to be 50 ± 1 meV. A coarse graining algorithm is applied which greatly enhances the efficiency of the simulations at low temperatures, down to 10 K, without altering the underlying physical processes. Eventually, we argue that a combination of both on- and off-lattice kinetic Monte Carlo techniques is a good way for simulating large-scale processes in molecular solids over long time spans.