One-dimensional Potential Of Mean Force Research Articles

Counter-Intuitive Features of Particle Dynamics in Nanopores.

Using the framework of a continuous diffusion model based on the Smoluchowski equation, we analyze particle dynamics in the confinement of a transmembrane nanopore. We briefly review existing analytical results to highlight consequences of interactions between the channel nanopore and the translocating particles. These interactions are described within a minimalistic approach by lumping together multiple physical forces acting on the particle in the pore into a one-dimensional potential of mean force. Such radical simplification allows us to obtain transparent analytical results, often in a simple algebraic form. While most of our findings are quite intuitive, some of them may seem unexpected and even surprising at first glance. The focus is on five examples: (i) attractive interactions between the particles and the nanopore create a potential well and thus cause the particles to spend more time in the pore but, nevertheless, increase their net flux; (ii) if the potential well-describing particle-pore interaction occupies only a part of the pore length, the mean translocation time is a non-monotonic function of the well length, first increasing and then decreasing with the length; (iii) when a rectangular potential well occupies the entire nanopore, the mean particle residence time in the pore is independent of the particle diffusivity inside the pore and depends only on its diffusivity in the bulk; (iv) although in the presence of a potential bias applied to the nanopore the "downhill" particle flux is higher than the "uphill" one, the mean translocation times and their distributions are identical, i.e., independent of the translocation direction; and (v) fast spontaneous gating affects nanopore selectivity when its characteristic time is comparable to that of the particle transport through the pore.

On the Nonaggregation of Normal Adult Hemoglobin and the Aggregation of Sickle Cell Hemoglobin.

Sickle cell disease is a genetic disorder associated with a single mutation (Glu-β6 → Val-β6) in the β chains of hemoglobin, causing the polymerization of deoxygenated sickle cell hemoglobin (deoxy-HbS). The deoxy-HbS binding free energy was recently studied through molecular simulations, and a value of -14 ± 1 kcal mol-1 was found. Here, we studied the binding free energy of normal adult hemoglobin (deoxy-HbA), which does not polymerize at normal physiological conditions, with the aim of elucidating the importance of the presence of Val-β6 and of the absence of Glu-β6 on the aggregation of deoxy-HbS. A binding free energy of -4.4 ± 0.5 kcal mol-1 was found from a one-dimensional potential of mean force. Hydrophobic interactions are shown to represent less than 20% of the interactions in the contact interface, and despite similarly strong hydrogen-bonded ion pairs (i.e., salt bridges) and water bridged electrostatic interactions are found for deoxy-HbA and deoxy-HbS, a large repulsive potential energy is associated with Glu-β6, whereas a mild attractive potential energy is connected with Val-β6. Interestingly, Asp-β73 switches from forming a major electrostatic repulsive pair with Glu-β6 in deoxy-HbA, to forming a major attractive residue pair with Val-β6 in deoxy-HbS, consistent with the view that damping of electrostatic repulsions involving Glu-β6, namely, those associated with Asp-β73, could be responsible for the polymerization of deoxy-HbA at high potassium phosphate concentrations. Solvation analysis shows that functional groups forming salt bridges and water bridged interactions preserve a nearly intact first hydration sphere, avoiding a complete dewetting free energy penalty. These results support the view that the absence of Glu-β6 is more important than the presence of Val-β6, and that although hydrophobic effects, associated with the Val-β6 dehydration and interaction with the hydrophobic pocket in the neighbor tetramer, are important, electrostatic interactions are dominant, opposite to a picture where HbS association is driven by hydrophobic interactions.

Molecular Dynamics Simulations of Ion Selectivity in a Claudin-15 Paracellular Channel.

Claudins are tissue-specific transmembrane proteins able to form junctions between two cells and regulate the flow of physiological solutes parallel to the cell walls, that is, the paracellular transport. Claudin-15 is highly expressed in the intestine where it forms efficient Na+ channels and Cl- barriers. However, the molecular details of these biological complexes are still unclear. Here, the permeation process of Na+, K+, and Cl- ions inside a refined structural model of a claudin-15 paracellular channel is investigated using all-atom molecular dynamics simulations in a double-bilayer and explicit solvent. One-dimensional potential of mean force (PMF) profiles, calculated using umbrella sampling (US) simulations, show that the channel allows the passage of the two physiological cations while excluding chloride. These features are generated by the action of several acidic residues, in particular the ring of D55 residues which is located at the narrowest region of the pore, in correspondence with the energy minimum for cations and the peak for chloride. We also used the Voronoi-tessellated milestoning method to obtain additional PMF profiles and the permeation timescale of the three ions. The milestoning PMFs agree well with those obtained by US, and the rate calculation reveals that the passage of chloride is almost 30 times slower than that of sodium. Our results are consistent with the known ability of claudin-15 to regulate tight junction selectivity and with the experimentally determined role of the acidic residues. This further validates our structural model and provides insights into the atomistic details of ion transport in paracellular channels that could be shared by other claudin-based architectures.

On the Binding Free Energy and Molecular Origin of Sickle Cell Hemoglobin Aggregation.

Protein aggregation is associated with various diseases, including Alzheimer and Parkinson as well as sickle cell disease (SCD). From a molecular point of view, protein aggregation depends on a complex balance of electrostatic and hydrophobic interactions mediated by water. An impressive manifestation of the importance of this balance concerns the human hemoglobin (HbA) mutant, HbS (sickle cell Hb), where a single substitution at the 6th position of HbA β-chains, from glutamic acid to valine, causes the polymerization of deoxygenated HbS (deoxy-HbS), responsible for SCD. HbS polymerization is believed to occur via a double nucleation mechanism initiated by the formation of HbS fibers (homogeneous nucleation), followed by fiber growth. Furthermore, it was proposed that homogeneous nucleation proceeds through a two-step mechanism, where metastable dense clusters play the role of nucleation precursors. Thus, hindering or delaying the formation of such precursors could represent a potential SCD therapeutic route. Here, we study, through molecular dynamics, the binding free energy and protein-protein contacts involved in the deoxy-HbS dimer aggregation and stabilization process. A binding free energy of ∼-14.0 ± 1 kcal/mol is estimated from a one-dimensional potential of mean force. Analysis of protein-protein interactions shows that both electrostatic and van der Waals interactions play an important role on the aggregation of HbS. With respect to the former, our results indicate that aggregation is largely favored by the formation of salt bridges (SB), mostly, Lys-Glu, Lys-Asp, and Heme-Lys SB, which outweigh electrostatic repulsions involving similar residues. Thus, our results suggest that a potential antisickling drug could be one with the ability to weaken or hinder the formation of a few SB between carboxylate and ammonium groups.

Dimerization of the EphA1 receptor tyrosine kinase transmembrane domain: Insights into the mechanism of receptor activation.

EphA1 is a receptor tyrosine kinase (RTK) that plays a key role in developmental processes, including guidance of the migration of axons and cells in the nervous system. EphA1, in common with other RTKs, contains an N-terminal extracellular domain, a single transmembrane (TM) α-helix, and a C-terminal intracellular kinase domain. The TM helix forms a dimer, as seen in recent NMR studies. We have modeled the EphA1 TM dimer using a multiscale approach combining coarse-grain (CG) and atomistic molecular dynamics (MD) simulations. The one-dimensional potential of mean force (PMF) for this system, based on interhelix separation, has been calculated using CG MD simulations. This provides a view of the free energy landscape for helix–helix interactions of the TM dimer in a lipid bilayer. The resulting PMF profiles suggest two states, consistent with a rotation-coupled activation mechanism. The more stable state corresponds to a right-handed helix dimer interacting via an N-terminal glycine zipper motif, consistent with a recent NMR structure (2K1K). A second metastable state corresponds to a structure in which the glycine zipper motif is not involved. Analysis of unrestrained CG MD simulations based on representative models from the PMF calculations or on the NMR structure reveals possible pathways of interconversion between these two states, involving helix rotations about their long axes. This suggests that the interaction of TM helices in EphA1 dimers may be intrinsically dynamic. This provides a potential mechanism for signaling whereby extracellular events drive a shift in the repopulation of the underlying TM helix dimer energy landscape.

Unrestrained Computation of Free Energy along a Path

We apply the adaptive biasing potential (ABP) method to optimize the principal curve defining a conformational transition between two known end states and to subsequently compute the one-dimensional potential of mean force as a function of arc length along the principal curve. This approach allows the use of the ABP method in a collective variable space of arbitrary dimension and offers several advantages over line-search methods. First, configurations are neither generated along an initial path for the transition nor equilibrated during evolution of the path. Second, and most importantly, the powerful sampling provided by the ABP serves to accelerate the dynamics during the optimization and computation of the free energy. Finally, the free energy is formulated as a potential of mean force that captures changes in the reaction channel along the principal curve, in contrast to the free energy profile evaluated from the local free-energy gradient in restrained path optimization methods. We first demonstrate the ABP formulation of path optimization using a two-dimensional potential surface and then with a more complex system of Src protein tyrosine kinase. The method is shown to be efficient and robust in the case of rugged, free-energy landscapes.

Theory and Simulation of Ion Conduction in the Pentameric GLIC Channel

GLIC is a bacterial member of the large family of pentameric ligand-gated ion channels. To study ion conduction through GLIC and other membrane channels, we combine the one-dimensional potential of mean force for ion passage with a Smoluchowski diffusion model, making it possible to calculate single-channel conductance in the regime of low ion concentrations from all-atom molecular dynamics (MD) simulations. We then perform MD simulations to examine sodium ion conduction through the GLIC transmembrane pore in two systems with different bulk ion concentrations. The ion potentials of mean force, calculated from umbrella sampling simulations with Hamiltonian replica exchange, reveal a major barrier at the hydrophobic constriction of the pore. The relevance of this barrier for ion transport is confirmed by a committor function that rises sharply in the barrier region. From the free evolution of Na(+) ions starting at the barrier top, we estimate the effective diffusion coefficient in the barrier region, and subsequently calculate the conductance of the pore. The resulting diffusivity compares well with the position-dependent ion diffusion coefficient obtained from restrained simulations. The ion conductance obtained from the diffusion model agrees with the value determined via a reactive-flux rate calculation. Our results show that the conformation in the GLIC crystal structure, with an estimated conductance of ~1 picosiemens at 140 mM ion concentration, is consistent with a physiologically open state of the channel.

Curcumin Prevents Aggregation in α-Synuclein by Increasing Reconfiguration Rate

α-Synuclein is a protein that is intrinsically disordered in vitro and prone to aggregation, particularly at high temperatures. In this work, we examined the ability of curcumin, a compound found in turmeric, to prevent aggregation of the protein. We found strong binding of curcumin to α-synuclein in the hydrophobic non-amyloid-β component region and complete inhibition of oligomers or fibrils. We also found that the reconfiguration rate within the unfolded protein was significantly increased at high temperatures. We conclude that α-synuclein is prone to aggregation because its reconfiguration rate is slow enough to expose hydrophobic residues on the same time scale that bimolecular association occurs. Curcumin rescues the protein from aggregation by increasing the reconfiguration rate into a faster regime.

Failure of one-dimensional Smoluchowski diffusion models to describe the duration of conformational rearrangements in floppy, diffusive molecular systems: A case study of polymer cyclization

Motivated by recent experimental efforts to measure the duration of individual folding∕unfolding transitions in proteins and RNA, here we use simulations to study the duration of a simple transition mimicking an elementary step in biopolymer folding: the closure of a loop in a long polymer chain. While the rate of such a transition is well approximated by a one-dimensional Smoluchowski model that views the end-to-end distance dynamics of a polymer chain as diffusion governed by the one-dimensional potential of mean force, the same model fails rather dramatically to describe the duration of such transitions. Instead, the latter timescale is well described by a model where the chain ends diffuse freely, uninfluenced by the average entropic force imposed by the polymer chain. The effective diffusion coefficient then depends on the length scale of the loop closure transition. Our findings suggest that simple one-dimensional models, when applied to estimate the duration of reactive events in complex molecular systems, should be used with caution.

Downhill binding energy surface of the barnase-barstar complex

We have employed biased molecular dynamics simulations in explicit solvent to characterize the one-dimensional potential of mean force for the dissociation process of the barnase-barstar protein-protein complex. Unbinding of barstar from wild-type barnase was compared with dissociation from four charge-deletion mutants of barnase. Interestingly, we find in all cases that unbinding of barnase and barstar is an uphill process on a smooth, tilted energy landscape. The total free energy difference between the dissociated and bound state was similar for wild-type barnase-barstar and for the R87A mutant of barnase. The values for the three other mutant barnase mutants K27A, R59A, and R83Q were only about half as much. Besides, we have analyzed the conformational dynamics of important residues at the barnase-barstar interface. In the bound state, their conformational fluctuations are reduced relatively to the free state because of the formation of intermolecular contacts. Interestingly, we find that some residues also show decreased mobility at intermediate stages of the unbinding process suggesting that these residues may be involved in the first contacts being formed on binding. © 2010 Wiley Periodicals, Inc. Biopolymers 93: 977-985, 2010.

Standard Free Energy of Binding from a One-Dimensional Potential of Mean Force.

A practical approach that enables one to calculate the standard free energy of binding from a one-dimensional potential of mean force (PMF) is proposed. Umbrella sampling and the weighted histogram analysis method are used to generate a PMF along the reaction coordinate of binding. At each point, a restraint is applied orthogonal to the reaction coordinate to make possible the determination of the volume sampled by the ligand. The free energy of binding from an arbitrary unbound volume to the restrained bound form is calculated from the ratio of the PMF integrated over the bound region to that of the unbound. Adding the free energy changes from the standard-state volume to the unbound volume and from the restrained to the unrestrained bound state gives the standard free energy of binding. Exploration of the best choice of binding paths is also made. This approach is first demonstrated on a model binding system and then tested on the benzamidine-trypsin system for which reasonable agreement with experiment is found. A comparison is made with other methods to obtain the standard free energy of binding from the PMF.

Lipid Bilayer Deformation and the Free Energy of Interaction of a Kv Channel Gating-Modifier Toxin

A number of membrane proteins act via binding at the water/lipid bilayer interface. An important example of such proteins is provided by the gating-modifier toxins that act on voltage-gated potassium (Kv) channels. They are thought to partition to the headgroup region of lipid bilayers, and so provide a good system for probing the nature of interactions of a protein with the water/bilayer interface. We used coarse-grained molecular dynamics simulations to compute the one-dimensional potential of mean force (i.e., free energy) profile that governs the interaction between a Kv channel gating-modifier toxin (VSTx1) and model phospholipid bilayers. The reaction coordinate sampled corresponds to the position of the toxin along the bilayer normal. The course-grained representation of the protein and lipids enabled us to explore extended time periods, revealing aspects of toxin/bilayer dynamics and energetics that would be difficult to observe on the timescales currently afforded by atomistic molecular dynamics simulations. In particular, we show for this model system that the bilayer deforms as it interacts with the toxin, and that such deformations perturb the free energy profile. Bilayer deformation therefore adds an additional layer of complexity to be addressed in investigations of membrane/protein systems. In particular, one should allow for local deformations that may arise due to the spatial array of charged and hydrophobic elements of an interfacially located membrane protein.

Diffusion model of solute dynamics in a membrane channel: Mapping onto the two-site model and optimizing the flux

The steady-state flux through a singly occupied membrane channel is found for both discrete and continuum models of the solute dynamics in the channel. The former describes the dynamics as nearest-neighbor jumps between N sites, while the latter assumes that the molecule diffuses in a one-dimensional potential of mean force. For both models it is shown that the flux is the same as that for a simple two-site model with appropriately chosen rate constants, which contain all the relevant information about the more detailed dynamics. An interesting consequence of single occupancy is that the flux has a maximum as a function of the channel-solute interaction. If this interaction is too attractive, the molecule will never leave the channel, thus blocking it for the passage of other molecules. If it is too repulsive, the solute molecule will never enter the channel. Thus the flux vanishes in the two limits and, hence, has a maximum somewhere in-between. In the framework of the diffusion model, we find the optimal intrachannel potential of mean force that maximizes the flux using the calculus of variations. For a symmetric channel this potential is flat and occupies the entire channel. In the general case of an asymmetric channel, the optimal potential is obtained by tilting the optimal flat potential for the corresponding symmetric channel around the channel center, so that the solute is driven towards the reservoir with the lower solute concentration by a constant force. This implies that the flux is higher when the solute binding near the channel exit is stronger than that near the entrance.

Transient violations of the second law of thermodynamics in protein unfolding examined using synthetic atomic force microscopy and the fluctuation theorem

The synthetic atomic force microscopy (AFM) method is developed to simulate a periodically replicated atomistic system subject to force and length fluctuations characteristic of an AFM experiment. This new method is used to examine the forced-extension and subsequent rupture of the alpha-helical linker connecting periodic images of a spectrin protein repeat unit. A two-dimensional potential of mean force (PMF) along the length and a reaction coordinate describing the state of the linker was calculated. This PMF reveals that the basic material properties of the spectrin repeat unit are sensitive to the state of linker, an important feature that cannot be accounted for in a one-dimensional PMF. Furthermore, nonequilibrium simulations were generated to examine the rupture event in the context of the fluctuation theorem. These atomistic simulations demonstrate that trajectories which are in apparent violation of the second law can overcome unfolding barriers at significantly reduced rupture forces.

A theoretical study of intra-molecular vibrational effects on fractionation factors for molecules containing intra-molecular low-barrier hydrogen bonds

The fractionation factor is defined as the equilibrium constant for the reaction; R–H+DOH ⇌ R–D+HOH. Of interest are values of fractionation factors that obtain for reactions where reactants and/or products form intra-molecular low-barrier hydrogen bonds. Experimentally measured isotopic fractionation factors are usually interpreted via a one-dimensional potential energy surface along the intrinsic proton hydrogen bond coordinate. Here we show that coupling the intrinsic proton coordinate to an intra-molecular vibrational mode has large effects on the values of isotopic fractionation factors and offer a picture for the fractionation factor for low-barrier hydrogen bonds based on the one-dimensional potential of mean force on the intrinsic proton coordinate.

Simulation of the kinetics of ligand binding to a protein by molecular dynamics: geminate rebinding of nitric oxide to myoglobin.

We have begun to use molecular dynamics to simulate the kinetics of nitric oxide rebinding to myoglobin after photodissociation. Rebinding was simulated using a potential function that switches smoothly between a nonbinding potential and a binding potential as a function of the position and orientation of the ligand, with no barrier arising from the crossing of potential surfaces of different electron spin. In 96 of 100 trajectories, the ligand rebound in < 15 ps. The kinetic progress curve was obtained by determining the time in each trajectory at which the ligand rebound and then calculating the fraction of unbound ligands as a function of time. The curve can be well reproduced by a simple model based on the dynamics of a Langevin particle moving on a one-dimensional potential of mean force calculated from nonreactive protein trajectories. The rate of escape from the energy well adjacent to the heme is in good agreement with the value calculated from experimental data, suggesting that a multiple-well model provides a plausible explanation for the nonexponential rebinding kinetics. A transition-state analysis suggests that protein conformational relaxation coupled to the displacement of the iron from the heme plane is an unlikely cause for the nonexponential rebinding of nitric oxide.