Charge Inversion in 1:1 Electrolytes: Analyzing the Energetics.

The effects of bulk concentration, surface charge density, ionic diameter, and bulk dielectric constant on charge inversion in 1:1 electrolyte systems are investigated. The framework of the classical density functional theory is used to describe the mean electrostatic potential and the volume and electrostatic correlations, which combine to define the adsorption of ions at a positively charged surface. Our results show that a decrease in the dielectric constant, in particular, creates charge inversion for 1:1 electrolytes by amplifying both the electrostatic potential and the screening component (which is generally much larger than the excluded-volume component). Local electrical potential inversion can occur even for moderate concentrations and surface charges. These findings are especially significant for ionic liquids and systems with organic molecules as solvents, as these generally have a dielectric constant much smaller than water.

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.

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.

Revisiting the Charged Shell Model: A Density Functional Theory for Electrolytes.

Classical density functional theory (DFT) has proven to be a sophisticated and efficient approach for investigating charge systems. In DFT, the excess free energy functional for inhomogeneous charged hard-sphere fluids consists of hard-core interactions and charge-charge electrostatic interactions. The former component can be precisely described by well-established fundamental measure theory (FMT). The latter component is usually computed using the Poisson equation combined with the mean spherical approximation (MSA). In order to predict accurate density profiles of ions and satisfy some thermodynamics sum rules, Roth and Gillespie [J. Phys.: Condens. Matter 2016, 28, 244006] proposed a DFT combining a functional-based version of MSA and an approximated charged shell model. Here, we rebuild the DFT based on the exact charged shell model, and the analytic expressions for the shell interaction potential and the corresponding thermodynamic quantities are provided. The structural and thermodynamic properties of both bulk and inhomogeneous electrolyte systems are analyzed. Moreover, the software named Atif (an advanced theoretical tool for inhomogeneous fluids) is released to the public via this work.

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.

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.

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.

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.

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.

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.

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.

The binding interactions that maintain excitation-contraction coupling junctions in skeletal muscle.

Calcium for contraction of skeletal muscles is released via tetrameric ryanodine receptor (RYR1) channels of the sarcoplasmic reticulum (SR), which are assembled in ordered arrays called couplons at junctions where the SR abuts T tubules or plasmalemma. Voltage-gated Ca2+ (CaV1.1) channels, found in tubules or plasmalemma, form symmetric complexes called CaV tetrads that associate with and activate underlying RYR tetramers during membrane depolarization by conveying a conformational change. Intriguingly, CaV tetrads regularly skip every other RYR tetramer within the array; therefore, the RYRs underlying tetrads (named V), but not the voltage sensor-lacking (C) RYRs, should be activated by depolarization. Here we hypothesize that the checkerboard association is maintained solely by reversible binary interactions between CaVs and RYRs and test this hypothesis using a quantitative model of the energies that govern CaV1.1-RYR1 binding, which are assumed to depend on number and location of bound CaVs. A Monte Carlo simulation generates large statistical samples and distributions of state variables that can be compared with quantitative features in freeze-fracture images of couplons from various sources. This analysis reveals two necessary model features: (1) the energy of a tetramer must have wells at low and high occupation by CaVs, so that CaVs positively cooperate in binding RYR (an allosteric effect), and (2) a large energy penalty results when two CaVs bind simultaneously to adjacent RYR protomers in adjacent tetramers (a steric clash). Under the hypothesis, V and C channels will eventually reverse roles. Role reversal justifies the presence of sensor-lacking C channels, as a structural and functional reserve for control of muscle contraction.

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.

Assessing the accuracy of three classical density functional theories of the electrical double layer.

Classical density functional theory (DFT) is a useful tool to compute the structure of the electrical double layer because it includes ion-ion correlations due to excluded-volume effects (i.e., steric correlations) and ion screening effects (i.e., electrostatic correlations beyond the electrostatic mean-field potential). This paper systematically analyzes the accuracies of three different electrostatic excess free-energy functionals, as compared to Monte Carlo (MC) simulations of the planar electrical double layer, over a large parameter space. Specifically, we tested the reference fluid density (RFD) [Gillespie et al., J. Phys.: Condens. Matter 14, 12129 (2002)10.1088/0953-8984/14/46/317], functionalized mean spherical approximation (fMSA) [Roth and Gillespie, J. Phys.: Condens. Matter 28, 244006 (2016)10.1088/0953-8984/28/24/244006], and bulk fluid (BF) [Kierlik and Rosinberg, Phys. Rev. A 44, 5025 (1991)10.1103/PhysRevA.44.5025; Y. Rosenfeld, J. Chem. Phys. 98, 8126 (1993)10.1063/1.464569] functionals. Previous work compared these DFT methods to MC simulations only for a small set of parameters. Here, a total of twelve different cations were studied, with valences of +1, +2, and +3 and ion diameters of 0.15, 0.30, 0.60, and 0.90 nm at bulk concentrations between 1 µM and 1 M. The anion always had valence -1 and diameter 0.30 nm. The structure of the double layer of these charged, hard-sphere ions was computed for surface charges ranging from 0 to -0.50C/m^{2}. All the DFTs were compared against each other for all the parameters, as well as to 378 MC simulations. Overall, RFD was the best of the three functionals, while BF was the least accurate. fMSA performed significantly better than BF, making it a reasonable choice that is less computationally expensive than RFD. For monovalent cations, all three functionals worked reasonably well, except BF was qualitatively different from MC at very low surface charges. For multivalent cations, BF underestimated charge inversion while fMSA overestimated it. All DFTs performed poorly for small multivalent ions.

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.

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.

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.