Polymer-grafted surfaces and channels are increasingly used for the design of responsive materials and sensors due to robust performance and ease of use. Various strategies for the control of the nanoscale morphologies of such materials and devices are being tested. Entropic repulsion between the polymer chains in a grafted brush of sufficient density causes the chains to extend in the direction perpendicular to the grafting surface in comparison to the position of unattached polymers. When nanoparticles having attractive interactions with the polymers are introduced into the solvent, these nanoparticles tend to infiltrate into the brush and reduce its extension. Under certain conditions, a sharp reduction in brush height extension can occur over a narrow range of nanoparticle concentrations in solution. We describe a way of controlling transport through polymer-functionalized nanochannels with nanoparticle additives, relying on the physics of nanoparticles and polymer brushes under confinement, and we suggest a blueprint for the creation of a tunable nanovalve. The nanovalve is modeled as a cylinder with a polymer brush grafted on its inside surface. Brush properties such as the chain length and the grafting density are chosen so that the brush chains extend into the center of the cylinder in the absence of nanoparticles, occluding the flux of analyte molecules through the pore. When nanoparticles that are attracted to the polymers are introduced into solution, they infiltrate into the brush and partially collapse it against the cylindrical grafting surface, opening space in the center of the cylinder through which analyte molecules can flow. The operation of such a nanovalve is analyzed via self-consistent field theory calculations in the strong-stretching approximation. Self-consistent field analysis is supported by Langevin dynamics simulations of the underlying coarse-grained model of the polymer–nanoparticle system.

Influx of ferrous ions from the cytoplasm through 3-fold pores in the shell of ferritin protein is computed using a 3-dimensional Poisson–Nernst–Planck electrodiffusion model, with inputs such as the pore structure and the diffusivity profile of permeant Fe2+ ions extracted from all-atom molecular dynamics (MD) simulations. These calculations successfully reproduce experimental estimates of the transit time of Fe2+ through the ferritin coat, which is on the millisecond time scale and hence much too long to be directly simulated via all-atom MD. This is also much longer than the typical time scale for ion transit in standard membrane spanning ion channels whose pores bear structural similarity to that of the 3-fold ferritin pore. The slow time scale for Fe2+ transport through ferritin pores is traced to two features that distinguish the ferritin pore system from standard ion channels, namely, (i) very low concentration of cytoplasmic Fe2+ under physiological conditions and (ii) very small internal diffusion coefficients for ions inside the ferritin pore resulting from factors that include the divalent nature of Fe2+ and two rings of negatively charged amino acids surrounding a narrow geometric obstruction within the ferritin pore interior.

We study the mechanism of vacancy migration and phase transitions of 3D crystalline colloidal arrays (CCA) using Langevin dynamics simulations. We calculate the self-diffusion coefficient of the colloid particles and the diffusion constant for vacancies as a function of temperature and DLVO potential parameters. We investigate the phase behavior of several systems with different interaction potential parameters using Voronoi analysis. Voronoi polyhedra tessellation, which is a useful method for characterizing the nearest neighbor environment around each atom, provides an efficient and effective way to identify phase transitions as well as geometrical changes in crystals. Using Voronoi analysis, we show that several neighboring particles are involved in the vacancy migration process that causes the vacancy to diffuse.

A coarse-graining method based on the partitioning of atoms into compact flexible clusters is used to formulate the dynamics of the nonequilibrium response of a protein to ligand dissociation. The α-carbon positions are used as the degrees of freedom. The net stiffness between each pair of neighboring α-carbons is calculated for the quasi-static, overdamped regime within the harmonic (quadratic potential energy surface) using the equivalent stiffness matrix of the network of atoms occupying the intervening space within the locally interacting region. This localized approach realizes a divide and conquer strategy that results in a substantial reduction in computational complexity while accurately predicting relaxations under general loading conditions. A close correlation between the shapes and time scales of the relaxation curves of the coarse-grained and all-atom instances of two medium-sized proteins, T4 lysozyme and ferric binding protein (each of which having known apo and holo structures), was observed for the holo to the apo transitions. Furthermore, for both proteins the dominant modes of motion and the decay rates of the temporal relaxation profiles monitoring the separation distance between select amino acid pairs were found to be nearly identical when calculated on the coarse-grained and all-atom scales.

An invertebrate glutamate-gated chloride channel (GluCl) has recently been crystallized in an open-pore state. This channel is homologous to the human Cys-loop receptor family of pentameric ligand-gated ion channels, including anion-selective GlyR and GABAR and cation-selective nAChR and 5HT3. We implemented molecular dynamics (MD) in conjunction with an elastic network model to perturb the X-ray structure of GluCl and investigated the open channel stability and its ion permeation characteristics. Our study suggests that TM2 helical tilting may close GluCl near the hydrophobic constriction L254 (L9′), similar to its cation-selective homologues. Ion permeation characteristics were determined by Brownian dynamics simulations using a hybrid MD/continuum electrostatics approach to evaluate the free energy profiles for ion transport. Near the selectivity filter region (P243 or P-2′), the free energy barrier for Na+ transport is over 4 kBT higher than that for Cl–, indicating anion selectivity of the channel. Furthermore, three layers of positivity charged rings in the extracellular domain also contribute to charge selectivity and facilitate Cl– permeability over Na+. Collectively, the charge selectivity of GluCl may be determined by overall electrostatic and ion dehydration effects, perhaps not deriving from a single region of the channel (the selectivity filter region near the intracellular entrance).

Langevin dynamics is used to compute the time evolution of the nonequilibrium motion of the atomic coordinates of a protein in response to ligand dissociation. The protein potential energy surface (PES) is approximated by a harmonic basin about the minimum of the unliganded state. Upon ligand dissociation, the protein undergoes relaxation from the bound to the unbound state. A coarse graining scheme based on rotation translation blocks (RTB) is applied to the relaxation of the two domain iron transport protein, ferric binding protein. This scheme provides a natural and efficient way to freeze out the small amplitude, high frequency motions within each rigid fragment, thereby allowing for the number of dynamical degrees of freedom to be reduced. The results obtained from all flexible atom (constraint free) dynamics are compared to those obtained using RTB-Langevin dynamics. To assess the impact of the assumed rigid fragment clustering on the temporal relaxation dynamics of the protein molecule, three distinct rigid block decompositions were generated and their responses compared. Each of the decompositions was a variant of the one-block-per-residue grouping, with their force and friction matrices being derived from their fully flexible counterpart. Monitoring the time evolution of the distance separating a selected pair of amino acids, the response curves of the blocked decompositions were similar in shape to each other and to the control system in which all atomic degrees of freedom are fully independent. The similar shape of the blocked responses showed that the variations in grouping had only a minor impact on the kinematics. Compared with the all atom responses, however, the blocked responses were faster as a result of the instantaneous transmission of force throughout each rigid block. This occurred because rigid blocking does not permit any intrablock deformation that could store or divert energy. It was found, however, that this accelerated response could be successfully corrected by scaling each eigenvalue in the appropriate propagation matrix by the least-squares fitted slope of the blocked vs nonblocked eigenvalue spectra. The RTB responses for each test system were dominated by small eigenvalue overdamped Langevin modes. The large eigenvalue members of each response dissipated within the first 5 ps, after which the long time response was dominated by a modest set of low energy, overdamped normal modes, that were characterized by highly cooperative, functionally relevant displacements. The response assuming that the system is in the overdamped limit was compared to the full phase space Langevin dynamics results. The responses after the first 5 ps were nearly identical, confirming that the inertial components were significant only in the initial stages of the relaxation. Since the propagator matrix in the overdamped formulation is real-symmetric and does not require the inertial component in the propagator, the computation time and memory footprint was reduced by 1 order of magnitude.

Analytical estimation of state-to-state rate constants is carried out for a recently developed discrete state model of chloride ion motion in a CLC chloride channel (Coalson and Cheng, J. Phys. Chem. B2010, 114, 1424). In the original presentation of this model, the same rate constants were evaluated via three-dimensional Brownian dynamics simulations. The underlying dynamical theory is an appropriate single- or multiparticle three-dimensional Smoluchowski equation. Taking advantage of approximate geometric symmetries (based on the details of the model channel geometry), well-known formulas for state-to-state transition rates are appealed to herein and adapted as necessary to the problem at hand. Rates of ionic influx from a bulk electrolyte reservoir to the nearest binding site within the channel pore are particularly challenging to compute analytically because they reflect multi-ion interactions (as opposed to single-ion dynamics). A simple empirical correction factor is added to the single-ion rate constant formula in this case to account for the saturation of influx rate constants with increasing bulk Cl– concentration. Overall, the agreement between all analytically estimated rate constants is within a factor of 2 of those computed via three-dimensional Brownian dynamics simulations, and often better than this. Current–concentration curves obtained using rate constants derived from these two different computational approaches agree to within 25%.

Bacterial Gloeobacter violaceus pentameric ligand-gated ion channel (GLIC) is activated to cation permeation upon lowering the solution pH. Its function can be modulated by anesthetic halothane. In the present work, we integrate molecular dynamics (MD) and Brownian dynamics (BD) simulations to elucidate the ion conduction, charge selectivity, and halothane modulation mechanisms in GLIC, based on recently resolved X-ray crystal structures of the open-channel GLIC. MD calculations of the potential of mean force (PMF) for a Na+ revealed two energy barriers in the extracellular domain (R109 and K38) and at the hydrophobic gate of transmembrane domain (I233), respectively. An energy well for Na+ was near the intracellular entrance: the depth of this energy well was modulated strongly by the protonation state of E222. The energy barrier for Cl− was found to be 3−4 times higher than that for Na+. Ion permeation characteristics were determined through BD simulations using a hybrid MD/continuum electrostatics approach to evaluate the energy profiles governing the ion movement. The resultant channel conductance and a near-zero permeability ratio (PCl/PNa) were comparable to experimental data. On the basis of these calculations, we suggest that a ring of five E222 residues may act as an electrostatic gate. In addition, the hydrophobic gate region may play a role in charge selectivity due to a higher dehydration energy barrier for Cl− ions. The effect of halothane on the Na+ PMF was also evaluated. Halothane was found to perturb salt bridges in GLIC that may be crucial for channel gating and open-channel stability, but had no significant impact on the single ion PMF profiles.

A discrete-state model of chloride ion motion in a ClC chloride channel is constructed, following a previously developed multi-ion continuous space model of the same system (Cheng, M. H.; Mamonov, A. B.; Dukes, J. W.; Coalson, R. D. J. Phys. Chem. B 2007, 111, 5956) that included a simplistic representation of the fast gate in this channel. The reducibility of the many-body continuous space to the eight discrete-state model considered in the present work is examined in detail by performing three-dimensional Brownian dynamics simulations of each allowed state-to-state transition in order to extract the appropriate rate constant for this process, and then inserting the pairwise rate constants thereby obtained into an appropriate set of kinetic master equations. Experimental properties of interest, including the rate of Cl− ion permeation through the open channel and the average rate of closing of the fast gate as a function of bulk Cl− ion concentrations in the intracellular and extracellular electrolyte reservoirs are computed. Good agreement is found between the results obtained via the eight discrete-state model versus the multi-ion continuous space model, thereby encouraging continued development of the discrete-state model to include more complex behaviors observed experimentally in these channels.

Dynamic linear response theory is adapted to the problem of computing the time evolution of the atomic coordinates of a protein in response to the unbinding of a ligand molecule from a binding pocket within the protein. When the ligand dissociates from the molecule, the protein molecule finds itself out of equilibrium and its configuration begins to change, ultimately coming to a new stable configuration corresponding to equilibrium in a force field that lacks the ligand-protein interaction terms. Dynamic linear response theory (LRT) relates the nonequilibrium motion of the protein atoms that ensues after the ligand molecule dissociates to equilibrium dynamics in the force field, or equivalently, on the potential energy surface (PES) relevant to the unliganded protein. In general, the connection implied by linear response theory holds only when the ligand−protein force field is small. However, in the case where the PES of the unliganded protein system is a quadratic (harmonic oscillator) function of the coordinates, and the force of the ligand upon the protein molecule in the ligand-bound conformation is constant (the force on each atom in the protein is independent of the location of the atom), dynamic LRT is exact for any ligand−protein force field strength. An analogous statement can be made for the case where the atoms in the protein are subjected to frictional and random noise forces in accord with the Langevin equation (to account for interaction of the protein with solvent, for example). We numerically illustrate the application of dynamic LRT for a simple harmonic oscillator model of the ferric binding protein, and for an analogous model of T4 lysozyme. Using a physically appropriate value of the viscosity of water to guide the choice of friction parameters, we find relaxation time scales of residue−residue distances on the order of several hundred ps. Comparison is made to relevant experimental measurements.

We conduct computational analyses of ion permeation characteristics in a model glycine receptor (GlyR) modified by photo-sensitive compounds. In particular, we consider hypothetical attachment to the channel of charge-neutral chemical groups which can be photo-activated by shining light of an appropriate wavelength on the system. After illumination, the attached molecules become charged via a photodissociation process or excited into a charge-separated state (thus generating a significant electric dipole). We carry out Brownian Dynamics simulations of ion flow through the channel in the presence of the additional charges generated in this fashion. Based on these calculations, we predict that photo-activation of appropriately positioned photo-sensitive compounds near the channel mouth can significantly modify the rate of ion permeation and the current rectification ratio. Possible implications for GlyR-based device designs are briefly discussed.

An improved approach to updating the electric field in simulations of Coulomb gases using the local lattice technique introduced by Maggs and Rossetto [A.C. Maggs, V. Rossetto, Phys. Rev. Lett. 88 (2002) 196402] is described and tested. Using the Fast Fourier Transform (FFT) an independent configuration of electric fields subject to Gauss' law constraint can be generated in a single update step. This FFT based method is shown to outperform previous approaches to updating the electric field in the simulation of a basic test problem in electrostatics of strongly correlated systems.

The local diffusion constant of K+ inside the Gramicidin A (GA) channel has been calculated using four computational methods based on molecular dynamics (MD) simulations, specifically: Mean Square Displacement (MSD), Velocity Autocorrelation Function (VACF), Second Fluctuation Dissipation Theorem (SFDT) and analysis of the Generalized Langevin Equation for a Harmonic Oscillator (GLE-HO). All methods were first tested and compared for K+ in bulk water—all predicted the correct diffusion constant. Inside GA, MSD and VACF methods were found to be unreliable because they are biased by the systematic force exerted by the membrane-channel system on the ion. SFDT and GLE-HO techniques properly unbias the influence of the systematic force on the diffusion properties and predicted a similar diffusion constant of K+ inside GA, namely, ca. 10 times smaller than in the bulk. It was found that both SFDT and GLE-HO methods require extensive MD sampling on the order of tens of nanoseconds to predict a reliable diffusion constant of K+ inside GA.

Early crystal structures of prokaryotic CLC proteins identified three Cl– binding sites: internal (Sint), central (Scen), and external (Sext). A conserved external GLU (GLUex) residue acts as a gate competing for Sext. Recently, the first crystal structure of a eukaryotic transporter, CmCLC, revealed that in this transporter GLUex competes instead for Scen. Here, we use molecular dynamics simulations to investigate Cl– transport through CmCLC. The gating and Cl–/H+ transport cycle are inferred through comparative molecular dynamics simulations with protonated and deprotonated GLUex in the presence/absence of external potentials. Adaptive biasing force calculations are employed to estimate the potential of mean force profiles associated with transport of a Cl– ion from Sext to Sint, depending on the Cl– occupancy of other sites. Our simulations demonstrate that protonation of GLUex is essential for Cl– transport from Sext to Scen. The Scen site may be occupied by two Cl– ions simultaneously due to a high energy barrier (∼8 Kcal/mol) for a single Cl– ion to translocate from Scen to Sint. Binding two Cl– ions to Scen induces a continuous water wire from Scen to the extracellular solution through the side chain of the GLUex gate. This may initiate deprotonation of GLUex, which then drives the two Cl– ions out of Scen toward the intracellular side via two putative Cl– transport paths. Finally, a conformational cycle is proposed that would account for the exchange stoichiometry.