Surface charge regulation using classical density functional theory: the effect of divalent potential determining ions.

The charge regulation approach has been used to describe the charge of surfaces susceptible to the presence of protons and other ions. Conventionally, this model is used with the Poisson-Boltzmann equation, which generally neglects the finite size of the ions and the electrostatic correlations. Recently, progress has been made by coupling charge regulation with classical density functional theory (DFT), which explicitly includes these correlations. However, little is known about charge regulation at surfaces with both acid-base equilibria and complexation with multivalent ions. The main purpose of this work is to investigate the role divalent ions play in charge regulation. Using DFT, we show that the size of the divalent ion has significant consequences on the surface charge density and it should not be neglected. For the surface reactions investigated, the larger the size of the divalent cation, the greater the charge on the surface due to higher divalent concentration there. At low divalent concentration, the ion correlations play a second-order but non-negligible role; using Poisson-Boltzmann theory with point ions cannot recover the DFT surface charge. At high concentrations, ion correlations play a dominant role by creating charge inversion.

Numerical solution of compartment-based reaction/diffusion models with DABOSS algorithm.

The intracellular diffusive movement of molecular substances that are buffered by diffusing chelators is often modeled as movement between compartments with constant concentrations within which the buffering occurs. Here, an algorithm to solve such a system of time-dependent differential equations is presented. This Dynamic and Balanced Operator Splitting Scheme (DABOSS) combines dynamic time stepping and operator splitting techniques. The time stepping minimizes the number of time steps while bounding local errors. The balanced operator splitting separates the diffusion and reaction components (each of which is solved efficiently) in a way that preserves the correct steady-state behavior. Analysis shows that DABOSS scales almost linearly in the number of compartments and remains accurate over very long simulations. Moreover, DABOSS works efficiently for nanometer-sized compartments with sources of flux, showing that it is applicable to situations where more spatial resolution is desired.

Individual ion species chemical potentials in the Mean Spherical Approximation.

The Mean Spherical Approximation (MSA) is a commonly used thermodynamic theory for computing the energetics of ions in the primitive model (i.e., charged hard-sphere ions in a background dielectric). For the excess chemical potential, however, the early MSA formulations (which were widely adopted) only included the terms needed to compute the mean excess chemical potential (or the mean activity coefficient). Other terms for the chemical potential µi of individual species i were not included because they sum to 0 in the mean chemical potential. Here, we derive these terms to give a complete MSA formulation of the chemical potential. The result is a simple additive term for µi that we show is a qualitative improvement over the previous MSA version. In addition, our derivation shows that the MSA's assumption of global charge neutrality is not strictly necessary, so that the MSA is also valid for systems close to neutrality.

Simulating cardiac Ca2+ release units: effects of RyR cluster size and Ca2+ buffers on diastolic Ca2+ leak.

Leak of Ca2+ out of the cardiac sarcoplasmic reticulum (SR) via ryanodine receptors (RyRs) during diastole is vital to regulate SR Ca2+ levels. This leak can become deleterious when large spontaneous RyR-mediated Ca2+ release events evoke proarrhythmic Ca2+ waves that can lead to delayed after-depolarizations. Here, we model diastolic SR Ca2+ leak at individual SR Ca2+ release sites using computer simulations of RyR arrays like those in the dyadic cleft. The results show that RyR arrays size has a significant effect on SR Ca2+ leak, with bigger arrays producing larger and more frequent Ca2+ release events. Moreover, big RyR arrays are more susceptible to small changes in the levels of dyadic Ca2+ buffers. Such changes in buffering shift Ca2+ leak from small Ca2+ release events (involving few open RyRs) to larger events (with many open RyRs). Moreover, by analyzing a large parameter space of possible buffering and SR Ca2+ loads, we find further evidence for the hypothesis that SR Ca2+ leak by RyR arrays can undergo a sudden phase transition.

Electrostatic correlations in electrolytes: Contribution of screening ion interactions to the excess chemical potential.

A new theory for the electrostatic component of the chemical potential for homogeneous electrolytes modeled with the primitive model is developed. This Mean Countershell Approximation (MCSA) is an analytic theory derived by including the interactions between the ions' screening clouds. At molar concentrations, these contribute substantially to the excess chemical potential but are absent in classical Debye-Hückel and Mean Spherical Approximation (MSA) theories. Simulations show that the MCSA is highly accurate, including at the low dielectric constants of ionic liquids. While sharing a mathematical framework with the MSA, the MCSA has simpler formulas and is qualitatively more accurate when there is ion size asymmetry.

Recruiting RyRs to Open in a Ca2+ Release Unit: Single-RyR Gating Properties Make RyR Group Dynamics.

In cardiac myocytes, clusters of type-2 ryanodine receptors (RyR2s) release Ca2+ from the sarcoplasmic reticulum (SR) via a positive feedback mechanism in which fluxed Ca2+ activates nearby RyRs. Although the general principles of this are understood, less is known about how single-RyR gating properties define the RyR group dynamics in an array of many channels. Here, we examine this using simulations with three models of RyR gating that have identical open probabilities: the commonly used two-state Markov gating model, one that utilizes multiple exponentials to fit single-channel open time (OT) and closed time (CT) distributions, and an extension of this multiexponential model that also includes experimentally measured correlations between single-channel OTs and CTs. The simulations of RyR clusters that utilize the multiexponential gating model produce infrequent Ca2+ release events with relatively few open RyRs. Ca2+ release events become even smaller when OT/CT correlations are included. This occurs because the correlations produce a small but consistent bias against recruiting more RyRs to open during the middle of a Ca2+ release event, between the initiation and termination phases (which are unaltered compared to the uncorrelated simulations). In comparison, the two-state model produces frequent, large, and long Ca2+ release events because it had a recruitment bias in favor of opening more RyRs. This difference stems from the two-state model's single-RyR OT and CT distributions being qualitatively different from the experimental ones. Thus, the details of single-RyR gating can profoundly affect SR Ca2+ release even if open probability and mean OTs and CTs are identical. We also show that Ca2+ release events can terminate spontaneously without any reduction in SR [Ca2+], luminal regulation, Ca2+-dependent inactivation, or physical coupling between RyRs when Ca2+ flux is below a threshold value. This supports and extends the pernicious attrition/induction decay hypothesis that SR Ca2+ release events terminate below a threshold Ca2+ flux.

Simulating diffusion from a cluster of point sources using propagation integrals.

A computational methodology to simulate the diffusion of ions from point sources (e.g., ion channels) is described. The outlined approach computes the ion concentration from a cluster of many ion channels at pre-specified locations as a function of time using the theory of propagation integrals. How the channels' open/closed states evolve in time does not need to be known at the start of the simulation, but can be updated on-the-fly as the simulation goes along. The technique uses analytic formulas for the solutions of the diffusion equation for three common geometries: (1) ions diffusing from a membrane (planar symmetry); (2) ions diffusing into a narrow cleft for effective two-dimensional diffusion (cylindrical symmetry); and (3) ions diffusing into open space like the cytosol (spherical symmetry). Because these formulas are exact solutions valid for arbitrarily long timesteps, no spatial or time discretizations are necessary. The only discrete locations are where the ion concentration is computed, and the only discrete timesteps are when the channels' open/closed states are updated. Beyond pure diffusion, the technique is generalized to the Excess Buffer Approximation of ion chelation to give an analytic solution of this approximation of the full reaction/diffusion system. Both the pure diffusion and the diffusion/buffering algorithms scale linearly with the number of channels and the number of ion concentration locations.

Calsequestrin depolymerizes when calcium is depleted in the sarcoplasmic reticulum of working muscle.

Calsequestrin, the only known protein with cyclical storage and supply of calcium as main role, is proposed to have other functions, which remain unproven. Voluntary movement and the heart beat require this calcium flow to be massive and fast. How does calsequestrin do it? To bind large amounts of calcium in vitro, calsequestrin must polymerize and then depolymerize to release it. Does this rule apply inside the sarcoplasmic reticulum (SR) of a working cell? We answered using fluorescently tagged calsequestrin expressed in muscles of mice. By FRAP and imaging we monitored mobility of calsequestrin as [Ca2+] in the SR--measured with a calsequestrin-fused biosensor--was lowered. We found that calsequestrin is polymerized within the SR at rest and that it depolymerized as [Ca2+] went down: fully when calcium depletion was maximal (a condition achieved with an SR calcium channel opening drug) and partially when depletion was limited (a condition imposed by fatiguing stimulation, long-lasting depolarization, or low drug concentrations). With fluorescence and electron microscopic imaging we demonstrated massive movements of calsequestrin accompanied by drastic morphological SR changes in fully depleted cells. When cells were partially depleted no remodeling was found. The present results support the proposed role of calsequestrin in termination of calcium release by conformationally inducing closure of SR channels. A channel closing switch operated by calsequestrin depolymerization will limit depletion, thereby preventing full disassembly of the polymeric calsequestrin network and catastrophic structural changes in the SR.

Modeling the Device Behavior of Biological and Synthetic Nanopores with Reduced Models.

Biological ion channels and synthetic nanopores are responsible for passive transport of ions through a membrane between two compartments. Modeling these ionic currents is especially amenable to reduced models because the device functions of these pores, the relation of input parameters (e.g., applied voltage, bath concentrations) and output parameters (e.g., current, rectification, selectivity), are well defined. Reduced models focus on the physics that produces the device functions (i.e., the physics of how inputs become outputs) rather than the atomic/molecular-scale physics inside the pore. Here, we propose four rules of thumb for constructing good reduced models of ion channels and nanopores. They are about (1) the importance of the axial concentration profiles, (2) the importance of the pore charges, (3) choosing the right explicit degrees of freedom, and (4) creating the proper response functions. We provide examples for how each rule of thumb helps in creating a reduced model of device behavior.

Electric Double Layers with Surface Charge Regulation Using Density Functional Theory.

Surprisingly, the local structure of electrolyte solutions in electric double layers is primarily determined by the solvent. This is initially unexpected as the solvent is usually a neutral species and not a subject to dominant Coulombic interactions. Part of the solvent dominance in determining the local structure is simply due to the much larger number of solvent molecules in a typical electrolyte solution.The dominant local packing of solvent then creates a space left for the charged species. Our classical density functional theory work demonstrates that the solvent structural effect strongly couples to the surface chemistry, which governs the charge and potential. In this article we address some outstanding questions relating double layer modeling. Firstly, we address the role of ion-ion correlations that go beyond mean field correlations. Secondly we consider the effects of a density dependent dielectric constant which is crucial in the description of a electrolyte-vapor interface.

Energetics of counterion adsorption in the electrical double layer.

The energetics of the electrical double layer (EDL) is studied in a systematic way to define how different components of the chemical potential help or hinder cation adsorption at a negatively charged wall. Specifically, the steric (i.e., excluded-volume interactions), mean electrostatic, and screening (i.e., electrostatic correlations beyond the mean field) components were computed using classical density functional theory of the primitive model of ions (i.e., ions as charged, hard spheres in a background dielectric). The reduced physics of the primitive model allows for an extensive analysis over a large parameter space: cation valences +1, +2, and +3, cation diameters 0.15, 0.30, 0.60, and 0.90 nm, bulk concentrations ranging from 1 µM to 1M, and surface charges ranging from 0 to -0.50 C/m2. Our results show that all components are necessary to understand the physics of the EDL. The screening component is always significant; for small monovalent cations such as K+, it is generally much larger than the steric component, and for multivalent ions, charge inversion cannot occur without it. At moderate surface charges, the screening component makes the electrostatic potential less negative than in classical Poisson-Boltzmann theory, sometimes even positive (charge inversion). At high surface charges, this is overcome by the repulsive potential of the steric component as the first ion layer becomes extremely crowded. Large negative electrostatic potentials counteract this to draw even more cations into the first layer. Although we used an approximate model of the EDL, the physics inherent in these trends appears to be quite general.

Multiscale analysis of the effect of surface charge pattern on a nanopore's rectification and selectivity properties: From all-atom model to Poisson-Nernst-Planck.

We report a multiscale modeling study for charged cylindrical nanopores using three modeling levels that include (1) an all-atom explicit-water model studied with molecular dynamics, and reduced models with implicit water containing (2) hard-sphere ions studied with the Local Equilibrium Monte Carlo simulation method (computing ionic correlations accurately), and (3) point ions studied with Poisson-Nernst-Planck theory (mean-field approximation). We show that reduced models are able to reproduce device functions (rectification and selectivity) for a wide variety of charge patterns, that is, reduced models are useful in understanding the mesoscale physics of the device (i.e., how the current is produced). We also analyze the relationship of the reduced implicit-water models with the explicit-water model and show that diffusion coefficients in the reduced models can be used as adjustable parameters with which the results of the explicit- and implicit-water models can be related. We find that the values of the diffusion coefficients are sensitive to the net charge of the pore but are relatively transferable to different voltages and charge patterns with the same total charge.

An Experimental Approach to Systematically Probe Charge Inversion in Nanofluidic Channels.

Charge inversion of the surfaces of nanofluidic channels occurs in systems with high-surface charge and/or highly charged ions and is of particular interest because of applications in biological and energy conversion systems. However, the details of such charge inversion have not been clearly elucidated. Specifically, although we can experimentally and theoretically show charge inversion, understanding at what conditions charge inversion occurs, as well how much the charge-inverting ions change the surface, are not known. Here, we show a novel experimental approach for uniquely finding both the ζ-potential and adsorption time of charge inverting ions in aqueous nanofluidic systems.

Ryanodine Receptor Open Times Are Determined in the Closed State.

The ryanodine receptor (RyR) ion channel releases Ca2+ from intracellular stores by conducting Ca2+ but also by recruiting neighboring RyRs to open, as RyRs are activated by micromolar levels of cytosolic Ca2+. Using long single-RyR recordings of the cardiac isoform (RyR2), we conclude that Ca2+ binding to the cytosolic face of RyR while the channel is closed determines the distribution of open times. This mechanism explains previous findings that RyR is not activated by its own fluxed Ca2+. Our measurements also bolster previous findings that luminal [Ca2+] can affect both the cytosolic activation and inactivation sites and that RyR has different gating modes for the same ionic conditions.

Sarcoplasmic Reticulum Ca2+ Release Uses a Cascading Network of Intra-SR and Channel Countercurrents.

In muscle, Ca2+ release from the sarcoplasmic reticulum (SR) into the cytosol is mediated through the ryanodine receptors (RyRs) and sustained by countercurrents that keep the SR membrane potential near 0 mV. Likewise, Ca2+ reuptake by the sarco/endoplasmic reticulum Ca2+ ATPase pump requires countercurrent. Although evidence has suggested that TRIC K+ channels and/or RyR K+ influx provide these countercurrents, the exact sources have not yet been determined. We used an equivalent circuit compartment model of a cardiac SR, the surrounding cytosol, and the dyadic cleft to probe the sources of countercurrent during a complete cardiac cycle. By removing and relocating TRIC K+ channels, as well as limiting when they are active, we explored the various possible sources of SR countercurrent under many conditions. Our simulations indicate that no single channel type is essential for countercurrent. Rather, a cascading network of countercurrents is present with anion fluxes within the SR redistributing charges throughout the full SR volume. This allows ion channels in the entire SR membrane, far from the Ca2+ fluxes through the RyRs in the junctional SR and sarco/endoplasmic reticulum Ca2+ ATPase pump in the nonjunctional SR, to mediate countercurrents that support Ca2+ release and reuptake. This multifactorial network of countercurrents allows Ca2+ release to be remarkably robust.

Sarcoplasmic reticulum Ca2+, Mg2+, K+, and Cl- concentrations adjust quickly as heart rate changes.

During systole, Ca2+ is released from the sarcoplasmic reticulum (SR) through ryanodine receptors (RyRs) while, simultaneously, other ions (specifically K+, Mg2+, and Cl-) provide counter-ion flux. These ions move back into the SR during diastole through the SERCA pump and SR K+ and Cl- channels. In homeostasis, all ion concentrations in different cellular regions (e.g., junctional and non-junctional SR, dyadic cleft, and cytosol) are the same at the beginning and end of the cardiac cycle. Here, we used an equivalent circuit compartment model of the SR and the surrounding cytoplasm to understand the heart rate dependence of SR ion homeostasis. We found that the Ca2+, Mg2+, K+, and Cl- concentrations in the SR and the cytoplasm self-adjust within just a few heartbeats with only very small changes in Mg2+, K+, and Cl- concentrations and membrane voltages (just a few percent). However, those small changes were enough to compensate for the large heart-rate-dependent changes in SR and cytoplasmic Ca2+ concentrations in the new steady state. The modeling suggests that ion adaptation to increases in heart rate is inherent to the system and that physiological changes that increase contractility and cardiac output are accommodated by the same self-adjusting mechanism of producing small changes in ion driving forces. Our findings also support the long-held hypothesis that SR membrane potentials are small (~1-2mV).

Two-Dimensional Electric Double Layer Structure with Heterogeneous Surface Charge.

In this work we present a systematic study of the lateral (parallel to the wall) and normal (perpendicular to the wall) nanostructure of the electric double layer at a heterogeneous interface between two regions of different surface charges, often found in nanoscale electrochemical devices. Specifically, classical density functional theory (DFT) is used to probe a cation concentration range of 10 mM to 1 M, for valences of +1, + 2, and +3, and a diameter range of 0.15-0.9 nm over widely varying surface charges (between -0.15 and +0.15 C/m2). The DFT results predict significant lateral and normal nanostructure in the form of ion concentration oscillations. These results are directly compared with those from Poisson-Boltzmann theory, showing significant deviation between the two theories, not only in the concentration profiles, but also in the sign of the electrostatic potential.

Multiscale modeling of a rectifying bipolar nanopore: explicit-water versus implicit-water simulations.

In a multiscale modeling approach, we present computer simulation results for a rectifying bipolar nanopore at two modeling levels. In an all-atom model, we use explicit water to simulate ion transport directly with the molecular dynamics technique. In a reduced model, we use implicit water and apply the Local Equilibrium Monte Carlo method together with the Nernst-Planck transport equation. This hybrid method makes the fast calculation of ion transport possible at the price of lost details. We show that the implicit-water model is an appropriate representation of the explicit-water model when we look at the system at the device (i.e., input vs. output) level. The two models produce qualitatively similar behavior of the electrical current for different voltages and model parameters. Looking at the details of concentration and potential profiles, we find profound differences between the two models. These differences, however, do not influence the basic behavior of the model as a device because they do not influence the z-dependence of the concentration profiles which are the main determinants of current. These results then address an old paradox: how do reduced models, whose assumptions should break down in a nanoscale device, predict experimental data? Our simulations show that reduced models can still capture the overall device physics correctly, even though they get some important aspects of the molecular-scale physics quite wrong; reduced models work because they include the physics that is necessary from the point of view of device function. Therefore, reduced models can suffice for general device understanding and device design, but more detailed models might be needed for molecular level understanding.

(Almost) Stationary Isotachophoretic Concentration Boundary in a Nanofluidic Channel Using Charge Inversion.

The present work is an experimental study of a new means to induce a quasi-stationary boundary for concentration or separation in a nanochannel induced by charge inversion. Instead of using pressure-driven counter-flow to keep the front stationary, we exploit charge inversion by a highly charged electrolyte, Ru(bpy)3Cl2, that changes the sign of the zeta potential in part of the channel from negative to positive. Having a non-charge inverting electrolyte (MgCl2) in the other part of the channel and applying an electric field can create a standing front at the interface between them without added dispersion due to an externally applied pressure-driven counterflow. The resulting slow moving front position can be easily imaged optically since Ru(bpy)3Cl2 is fluorescent. A simple analytical model for the velocity field and front axial position that reproduces the experimental location of the front shows that the location can be tuned by changing the concentration of the electrolytes (and thus local zeta potential). Both of these give the charge inversion-mediated boundary significant advantages over current methods of concentration and separation and the method is, therefore, of particular importance to chemical and biochemical analysis systems such as chromatography and separations and for enhancing the stacking performance of field amplified sample injection and isotachophoresis. By choosing a non-charge inverting electrolyte other than MgCl2, either this electrolyte or the Ru(bpy)3Cl2 solution can be made to be the leading or trailing electrolyte.

Selecting ions by size in a calcium channel: the ryanodine receptor case study.

Many calcium channels can distinguish between ions of the same charge but different size. For example, when cations are in direct competition with each other, the ryanodine receptor (RyR) calcium channel preferentially conducts smaller cations such as Li(+) and Na(+) over larger ones such as K(+) and Cs(+). Here, we analyze the physical basis for this preference using a previously established model of RyR permeation and selectivity. Like other calcium channels, RyR has four aspartate residues in its GGGIGDE selectivity filter. These aspartates have their terminal carboxyl group in the pore lumen, which take up much of the available space for permeating ions. We find that small ions are preferred by RyR because they can fit into this crowded environment more easily.