A theoretical study of polymorphism in VQIVYK fibrils.

The VQIVYK fragment from the Tau protein, also known as PHF6, is essential for aggregation of Tau into neurofibrillary lesions associated with neurodegenerative diseases. VQIVYK itself forms amyloid fibrils composed of paired ß-sheets. Therefore, the full Tau protein and VQIVYK fibrils have been intensively investigated. A central issue in these studies is polymorphism, the ability of a protein to fold into more than one structure. Using all-atom molecular simulations, we generate five stable polymorphs of VQIVYK fibrils, establish their relative free energy with umbrella sampling methods, and identify the side chain interactions that provide stability. The two most stable polymorphs, which have nearly equal free energy, are formed by interdigitation of the mostly hydrophobic VIY "face" sides of the ß-sheets. Another stable polymorph is formed by interdigitation of the QVK "back" sides. When we turn to examine structures from cryo-electron microscopy experiments on Tau filaments taken from diseased patients or generated in vitro, we find that the pattern of side chain interactions found in the two most stable face-to-face as well as the back-to-back polymorphs are recapitulated in amyloid structures of the full protein. Thus, our studies suggest that the interactions stabilizing PHF6 fibrils explain the amyloidogenicity of the VQIVYK motif within the full Tau protein and provide justification for the use of VQIVYK fibrils as a test bed for the design of molecules that identify or inhibit amyloid structures.

Mechanism of surface freezing of alkanes.

Using molecular dynamics simulation of octane (C8) and nonadecane (C19), we probe the mechanism of n-alkane surface freezing, the appearance of a crystalline monolayer above the liquid at a temperature Tsf above the bulk freezing point Tf. Formation of a crystalline monolayer occurs robustly in these systems. When Tf > Tsf, the surface frozen phase is metastable with respect to the solid but persists for long periods for study in simulations. Surface freezing of both C8 and C19 is driven by significant energy-lowering when alkane chains become ordered along the surface normal, and we elucidate the origins of this phenomenon. The degree of configurational disorder in the surface frozen layer relative to the solid is much larger for C8 compared to C19. From the Gibbsian viewpoint, we extract the excess energy and entropy of the liquid and surface frozen phases. We also consider the surface frozen layer as an intervening third phase, the viewpoint taken in previous theoretical analyses. Here, we find significantly increased entropy of the surface frozen phase of C8 associated with configurational disorder, while the energy and entropy of the surface frozen phase of C19 are marginally different from the bulk solid. Finally, by combining our previously determined solid-vapor surface free energies of C8 and C19 with liquid-vapor surface tensions from this work, we eliminate wetting as a possible mechanism for C8 surface freezing, but it remains a possibility for C19. We analyze the molecular structure of the liquid, surface frozen, and solid surfaces and discuss its relevance to thermodynamic properties.

On the determination of the crystal-vapor surface free energy, and why a Gaussian expression can be accurate for a system far from Gaussian.

The crystal-vapor surface free energy γ is an important physical parameter governing physical processes, such as wetting and adhesion. We explore exact and approximate routes to calculate γ based on cleaving an intact crystal into non-interacting sub-systems with crystal-vapor interfaces. We do this by turning off the interactions, ΔV, between the sub-systems. Using the soft-core scheme for turning off ΔV, we find that the free energy varies smoothly with the coupling parameter λ, and a single thermodynamic integration yields the exact γ. We generate another exact method, and a cumulant expansion for γ by expressing the surface free energy in terms of an average of e(-ßΔV) in the intact crystal. The second cumulant, or Gaussian approximation for γ is surprisingly accurate in most situations, even though we find that the underlying probability distribution for ΔV is clearly not Gaussian. We account for this fact by developing a non-Gaussian theory for γ and find that the difference between the non-Gaussian and Gaussian expressions for γ consist of terms that are negligible in many situations. Exact and approximate methods are applied to the (111) surface of a Lennard-Jones crystal and are also tested for more complex molecular solids, the surface of octane and nonadecane. Alkane surfaces were chosen for study because their crystal-vapor surface free energy has been of particular interest for understanding surface freezing in these systems.

Polarization charge: Theory and applications to aqueous interfaces.

When an electric field is applied across an interface, a dielectric will acquire a polarization charge layer, assumed infinitely thin in the theory of macroscopic dielectrics and also in most treatments of electrokinetic phenomena in nanoscale structures. In this work we explore the polarization charge layer in molecular detail. Various formal relations and a linear response theory for the polarization charge are presented. Properties of the polarization charge layer are studied for three aqueous interfaces: air-water, a crystalline silica surface with water, and an amorphous silica surface with water. The polarization charge is calculated from equilibrium simulations via linear response theory and from non-equilibrium simulations, and the results are within statistical error. The polarization charge is found to be distributed within a region whose width is on the order of a nanometer.

Experimental evidence for surface freezing in supercooled n-alkane nanodroplets.

Intermediate chain length (16 ≤i≤ 50) n-alkanes are known to surface freeze at temperatures that are up to three degrees higher than the equilibrium melting point [B. M. Ocko et al., Phys. Rev. E, 1997, 55, 3164-3182]. Our recent experimental results suggest that highly supercooled nanodroplets of n-octane and n-nonane also surface freeze, and subsequently bulk crystallization occurs. The data yield surface and bulk nucleation rates on the order of ~10(15) cm(-2) s(-1) and ~10(22) cm(-3) s(-1), respectively, at temperatures between 180 K and 200 K. Molecular dynamics simulations at the united atom level were used to follow the freezing of a supercooled n-octane drop and show that an ordered monolayer develops on the surface of the droplet almost immediately, and the rest of the droplet then freezes in a layer-by-layer manner.

A New Approach for Estimating the Free Energy Differences among Multiple Thermodynamic States in Statistical Simulations.

In this letter, a new approach to compute free energy differences (FEDs) between multiple thermodynamics states is introduced. The method directly uses energy probability densities, which can be extracted with high accuracy from equilibrium simulations to obtain FEDs. Methods in current use, such as Bennett acceptance ratio (BAR), its multistate generalization (MBAR), or the weighted histogram analysis method (WHAM), require iterative solution of nonlinear equations which are known to be slowly convergent. The equations providing MBAR FEDs are identical to those derived earlier by Souaille and Roux in a method that has become known informally as "binless WHAM". In contrast, we obtain FEDs by solution of linear equations. For the classic two-state problem, the statistical error of our method, solving linear equations, is shown analytically to match that of BAR under common conditions.

Analysis of the subcritical carbon dioxide-water interface.

We follow the evolution of the H(2)O/CO(2) interface at 300 K from the low pressure limit to near-critical pressures in molecular dynamics simulations using the SPC water and EPM2 carbon dioxide models. The intrinsic structure of the interface is elucidated by accumulating density profiles relative to the fluctuating capillary wave surface. Our main finding is that a carbon dioxide film of increasing density and thickness grows in two stages at the interface while the structure of the water surface barely changes. At low density, the entire film density profile grows linearly with the bulk CO(2) density. This regime continues up to a bulk CO(2) density of roughly 0.00095 Å(-3). At pressures above this point, we observe a distinct second peak in the CO(2) density, along with a tail of excess density that decays exponentially with distance from the interface. The decay length of the exponential tail diverges with increasing CO(2) pressure according to an inverse power law decay. Over the entire range of pressures, the CO(2) film had no detectable effect on the orientational order of the water surface. As expected, when the film of excess CO(2) at the interface grows, we find that the surface tension drops with increasing pressure. This is in qualitative accord with existing measurements, although the rate at which the surface tension falls with increasing pressure according to the SPC and EPM2 models is too small, indicating that the surface excess of CO(2) is underestimated by these models.

The water-amorphous silica interface: analysis of the Stern layer and surface conduction.

To explain why dynamical properties of an aqueous electrolyte near a charged surface seem to be governed by a surface charge less than the actual one, the canonical Stern model supposes an interfacial layer of ions and immobile fluid. However, large ion mobilities within the Stern layer are needed to reconcile the Stern model with surface conduction measurements. Modeling the aqueous electrolyte-amorphous silica interface at typical charge densities, a prototypical double layer system, the flow velocity does not vanish until right at the surface. The Stern model is a good effective model away from the surface, but cannot be taken literally near the surface. Indeed, simulations show no ion mobility where water is immobile, nor is such mobility necessary since the surface conductivity in the simulations is comparable to experimental values.

Dielectric Behavior near a Spherical Ion.

The dielectric response of a polar solvent to an ion is analyzed in terms of the bound charge, the net charge that accumulates near the ion as a consequence of the inhomogeneous polarization of the surrounding solvent. We demonstrate that the total bound charge arising in a full molecular treatment is identical to the total bound charge from standard continuum theory. In continuum theory, the bound charge resides in an infinitely thin layer, while in a molecular description the bound charge is spread over a region of finite width. Near simple atomic ions, the width of the bound charge distribution is roughly 1.3 nm. By simulating a sequence of ion charges from 0.1 to 2 e, where e is the magnitude of the electron charge, we analyze the applicability of linear response theory, which has been used by several authors. With increasing charge, the nonlinear response extends to an increasing distance from the ion. However, outside the region containing bound charge, the response is linear and in accord with continuum theory. Previous attempts to assign a dielectric constant to a solvent in the interfacial region are analyzed.

Site disorder in ice VII arising from hydrogen bond fluctuations.

Recently, multisite models used for the refinement of neutron diffraction data have suggested that the structure of ice VII is quite unlike that of its ordered counterpart, ice VIII. We investigate the oxygen site disorder by modeling the site displacement, obtained from periodic DFT calculations, as a function of the local hydrogen bond network. Then, using graph invariants to describe hydrogen bond fluctuations in the thermodynamic limit, we perform statistical mechanical calculations using the oxygen site displacement model developed here. We find that the probability distribution of the oxygen atom about its perfect lattice site more closely resembles the 100 model rather than the recently suggested 111 model, although both models represent a simplification of the actual site distribution. We also find a unimodel distribution for the hydrogen bonded oxygen-oxygen distance and a trimodel distribution for the nearest nonbonded oxygen-oxygen distance with a peak separation of approximately 0.1 A.

Molecular Dynamics Study of the Electric Double Layer and Nonlinear Spectroscopy at the Amorphous Silica-Water Interface.

The electrical double layer (EDL) at the amorphous silica-aqueous electrolyte interface is of long-standing scientific interest and current technological relevance. Using extensive molecular dynamics simulations, we have studied this EDL as a function of salt concentration for a silica surface charge density of -0.82e/nm2 (e = electron charge). The simulation results can be captured with a simple model by breaking the double-layer region into three zones: an inner region in which the Na+ counterion population is independent of [NaCl] and there are no Cl- co-ions, an intermediate region which hosts a population of nonexchangeable Na+ plus another group of Na+ and Cl- ions whose population is described by a Langmuir adsorption model, and an outer region where the ion distribution is well-described using the Poisson-Boltzmann theory. When the asymptotic [NaCl] >0.17 M, the adsorption of Na+ in the intermediate zone leads to an overcompensation of the negatively charged silica surface. Nonlinear spectroscopic experiments on the water-amorphous silica interface have been interpreted by others using the Gouy-Chapman model at low salt concentration and the constant capacitance model at high salt concentration. We discuss the applicability of these and other models and the implications for interpretation of the results of second harmonic and sum frequency generation experiments.

Origin of slow relaxation following photoexcitation of W7 in myoglobin and the dynamics of its hydration layer.

Molecular dynamics simulations are used to calculate the time-dependent Stokes shift following photoexcitation of Trp-7 (W7) in myoglobin. In agreement with experiment, a long time (approximately 60 ps) component is observed. Since the long time Stokes shift component is absent when we repeat the calculation with protein frozen at the instant of photoexcitation, we firmly establish that protein flexibility is required to observe slow Stokes shift dynamics in this case. A transition between sub-states near the middle of a 30 ns ground-state trajectory gave us an opportunity to compare solvation dynamics in two different environments. While some of the superficial features are different, we find that the underlying dynamics are shared by the two isomers. It is necessary to look beyond a decomposition of the Stokes shift into protein and water contributions and probe the underlying dynamics of protein side groups, backbone, and water dynamics to obtain a full picture of the relaxation process. We analyze water residence times, diffusion, and reorientation dynamics in the hydration layer. We find slow components in each of these quantities and critically examine their origin and how they affect the observed Stokes shift.

Hydrogen bond ordering in ice V and the transition to ice XIII.

The proton ordered version of ice V, ice XIII, was recently identified using Raman spectroscopy and neutron diffraction techniques. The transformation, between 108 and 117 K, only occurred in the presence of a small amount of dopant, similar to the proton ordering transition of ice Ih/XI. In this work, we investigate the hydrogen bond fluctuations in ice V and XIII with statistical mechanical techniques that use results from periodic electronic density functional theory calculations as input. We find a number of low-lying hydrogen bond configurations, approximately 20 within 10 K/water above the ground state state configuration, the structure of which agrees with fully ordered ice XIII. Using an analytic theory, graph invariants, we developed effective spin-lattice Hamiltonians governing hydrogen bond fluctuations to perform statistical mechanical calculations for a large simulation cell containing 6048 water molecules. Two models were constructed, one more elaborate than the first, to gauge the robustness of our methods when the unit cell is very complex and a large number of configurations lie close in energy to the ground state. The predicted proton ordering transitions, 62 and 72 K for the two models, are in qualitative agreement with experiment. Occupation probabilities, obtained from our simulations, compare well with values from recent neutron diffraction experiments and help verify our effective Hamiltonians. In both models, we find that a second order phase transition intervenes about 10 K above the transition to ice XIII, but its effect is negligible on the behavior of thermodynamic functions near the transition to ice XIII.

Model for the water-amorphous silica interface: the undissociated surface.

The physical and chemical properties of the amorphous silica-water interface are of crucial importance for a fundamental understanding of electrochemical and electrokinetic phenomena, and for various applications including chromatography, sensors, metal ion extraction, and the construction of micro- and nanoscale devices. A model for the undissociated amorphous silica-water interface reported here is a step toward a practical microscopic model of this important system. We have extended the popular BKS and SPC/E models for bulk silica and water to describe the hydrated, hydroxylated amorphous silica surface. The parameters of our model were determined using ab initio quantum chemical studies on small fragments. Our model will be useful in empirical potential studies, and as a starting point for ab initio molecular dynamics calculations. At this stage, we present a model for the undissociated surface. Our calculated value for the heat of immersion, 0.3 J x m(-2), falls within the range of reported experimental values of 0.2-0.8 J x m(-2). We also study the perturbation of water properties near the silica-water interface. The disordered surface is characterized by regions that are hydrophilic and hydrophobic, depending on the statistical variations in silanol group density.

Biomolecules at the amorphous silica/water interface: Binding and fluorescence anisotropy of peptides.

We investigate binding of the tripeptides Lys-Trp-Lys (KWK) and Glu-Trp-Glu (EWE) to the amorphous silica surface using atomistic simulations. These peptides were chosen because they were previously utilized in experiments measuring binding affinity and steady-state fluorescence anisotropy from the indole chromophore of the tryptophan residue. Our simulations were performed using silica with surface change density of -0.8 elementary charges per square nanometer, which is expected at neutral pH. Even though positive charged KWK binds more strongly to the negatively charged silica surface, EWE also forms bound complexes with the surface that are stable for at least 15ns of simulation, in agreement with the experiments which revealed evidence for binding of both KWK and EWE to silica. Binding mechanisms include a wide variety of electrostatic interactions, as well as hydrophobic interactions between the indole group and hydrophobic areas of the heterogeneous silica surface. The long-time limit of the fluorescence anisotropy of tryptophan is calculated from the simulations in order to help interpret the recent experiments. We identify several factors which control the magnitude of the fluorescence anisotropy for each binding configuration.

A molecular dynamics study of Lys-Trp-Lys: structure and dynamics in solution following photoexcitation.

We report studies of the structure and dynamics of a tripeptide Lys-Trp-Lys (KWK) in aqueous solution following photoexcitation by molecular dynamics simulations. For ground-state KWK, we observe three stable conformations with free energy differences of less than 5.2 kJ/mol. Each conformer is stabilized by a pi-cation interaction between one of three protonated amino groups and the indole moiety. For the excited state of tryptophan in KWK, the simulated molecular dynamics of the three isomers are similar, all in good agreement with recent femtosecond experiments (J. Phys. Chem. B 2005, 109, 16901). Specifically, we observe: (1) the fluorescence anisotropy is dominated by a single-exponential component and decays in approximately 130 ps, (2) the total dynamic Stokes shift reaches approximately 2700 cm(-1), and (3) the excited state relaxation dynamics occurs on several time scales ranging from femtoseconds to tens of picoseconds. The relaxation dynamics involve rapid initial response of neighboring water, followed by local motions of flexible peptide chains. These processes drive global restructuring of the tripeptide on a rather flat energy surface, inducing slower dynamics evident in both the water and protein contributions to the stabilization energy of the photoexcited chromophore. The water and protein dynamics are strongly correlated. On a still longer time scale, we observe isomerization of two excited state conformers to the other most stable one, an analogue for evolution of trajectories along the funnel on the rugged free energy landscape to the final "native" state. Our studies suggest new experiments to detect this unique dynamics.

Hydrogen bond topology and the ice VII/VIII and Ih/XI proton ordering phase transitions.

Ice Ih, ordinary ice at atmospheric pressure, is a proton-disordered crystal that when cooled under special conditions is believed to transform to ferroelectric proton-ordered ice XI, but this transformation is still subject to controversy. Ice VII, also proton disordered throughout its region of stability, transforms to proton-ordered ice VIII upon cooling. In contrast to the ice Ih/XI transition, the VII/VIII transition and the crystal structure of ice VIII are well characterized. In order to shed some light on the ice Ih proton ordering transition, we present the results of periodic electronic density functional theory calculations and statistical simulations. We are able to describe the small energy differences among the innumerable H-bond configurations possible in a large simulation cell by using an analytic theory to extrapolate from electronic DFT calculations on small unit cells to cells large enough to approximate the thermodynamic limit. We first validate our methods by comparing our predictions to the well-characterized ice VII/VIII proton ordering transition, finding agreement with respect to both the transition temperature and structure of the low-temperature phase. For ice Ih, our results indicate that a proton-ordered phase is attainable at low temperatures, the structure of which is in agreement with the experimentally proposed ferroelectric structure. The predicted transition temperature of is in qualitative agreement with the observed transition at on KOH-doped ice samples.

Prediction of a phase transition to a hydrogen bond ordered form of ice VI.

Ice VI is a hydrogen bond disordered crystal over its known region of stability. In this work, we predict that ice VI will transform into a hydrogen bond ordered phase near 108 K, and have identified the likely low-temperature phase as ferroelectric (space group Cc) with an antiferroelectric structure (space group P2(1)2(1)2(1)) close by in energy. Electronic density functional theory calculations provide input to our calculations, which are extended to cells large enough for statistical simulations by using graph invariants. A significant decrease in the configurational entropy is predicted as hydrogen bonds exhibit partial order above the transition, provided that the hydrogen bonds can equilibrate on an experimental time scale. Conversely, partial disorder is predicted at temperatures below the transition. Although some evidence for ordering of ice VI has been observed in experiments, a low-temperature proton ordered phase has not been identified experimentally.

Electro-osmotic flow of a model electrolyte.

Electro-osmotic flow is studied by nonequilibrium molecular dynamics simulations in a model system chosen to elucidate various factors affecting the velocity profile and facilitate comparison with existing continuum theories. The model system consists of spherical ions and solvent, with stationary, uniformly charged walls that make a channel with a height of 20 particle diameters. We find that hydrodynamic theory adequately describes simple pressure-driven (Poiseuille) flow in this model. However, Poisson-Boltzmann theory fails to describe the ion distribution in important situations, and therefore continuum fluid dynamics based on the Poisson-Boltzmann ion distribution disagrees with simulation results in those situations. The failure of Poisson-Boltzmann theory is traced to the exclusion of ions near the channel walls resulting from reduced solvation of the ions in that region. When a corrected ion distribution is used as input for hydrodynamic theory, agreement with numerical simulations is restored. An analytic theory is presented that demonstrates that repulsion of the ions from the channel walls increases the flow rate, and attraction to the walls has the opposite effect. A recent numerical study of electro-osmotic flow is reanalyzed in the light of our findings, and the results conform well to our conclusions for the model system.

DNA Binding to the Silica Surface.

We investigate the DNA-silica binding mechanism using molecular dynamics simulations. This system is of technological importance, and also of interest to explore how negatively charged DNA can bind to a silica surface, which is also negatively charged at pH values above its isoelectric point near pH 3. We find that the two major binding mechanisms are attractive interactions between DNA phosphate and surface silanol groups and hydrophobic bonding between DNA base and silica hydrophobic region. Umbrella sampling and the weighted histogram analysis method (WHAM) are used to calculate the free energy surface for detachment of DNA from a binding configuration to a location far from the silica surface. Several factors explain why single-stranded DNA (ssDNA) has been observed to be more strongly attracted to silica than double-stranded (dsDNA): (1) ssDNA is more flexible and therefore able to maximize the number of binding interactions. (2) ssDNA has free unpaired bases to form hydrophobic attachment to silica while dsDNA has to break hydrogen bonds with base partners to get free bases. (3) The linear charge density of dsDNA is twice that of ssDNA. We devise a procedure to approximate the atomic forces between biomolecules and amorphous silica to enable large-scale biomolecule-silica simulations as reported here.