Electrodiffusion: a continuum modeling framework for biomolecular systems with realistic spatiotemporal resolution

Differential behaviors of calcium-induced calcium release in one dimensional dendrite by Nernst-Planck equation, cable model and pure diffusion model

The source and dynamics of calcium is the key factor that regulates dendritic integration. Apart from the voltage-gated and ligand-gated calcium influx, an important source of calcium is from inner store of endoplasmic reticulum with a regenerative process of calcium-induced calcium release (CICR). To trigger this process, inositol 1,4,5-trisphosphate (IP3) and calcium are needed to satisfy certain requirements. The aim of our paper is to investigate how the CICR depends on the dynamics of membrane potential. We utilize one dimensional dendritic model to calculate membrane potential by Nernst-Planck Equation (NPE) and cable model and Pure Diffusion (PD) model, computational simulations are carried out to inject the calcium influx by synaptic stimulation and to predict subsequent CICR and calcium wave propagation. Our results demonstrate that CICR initiation and calcium wave propagation have much difference between electro-diffusion process of NPE and cable model. We find that cable model has lower threshold of IP3 stimulation to trigger CICR but is more difficult for calcium propagation than NPE, PD model requires even higher threshold of IP3 to initiate CICR process and calcium duration is shorter than NPE; the regenerative calcium wave propagates with faster speed in NPE than that in cable model and in PD model. Our work addresses the important role of electro-diffusion dynamics of charged ions in regulating CICR process in dendritic structure; and provides theoretical predictions for neurological process which requires sustaining calcium for downstream signaling processes.

Electrokinetic translocation of a deformable nanoparticle controlled by field effect in nanopores

Nanopores have become a popular single-molecule manipulation and detection technology. In this paper, we have constructed a continuum model of the nanopore; the arbitrary Lagrangian-Eulerian (ALE) method is used to describe the motion of particles and fluid. The mathematical model couples the stress-strain equation for the dynamics of a deformable particle, the Poisson equation for the electric field, the Navier-Stokes equations for the flow field, and the Nernst-Planck equations for ionic transport. Based on the model, the mechanism of field-effect regulation of particles passing through a nanopore is investigated. The results show that the transport of particles which is controlled by the field effect depends on the electroosmotic flow (EOF) generated by the gate electrode in the nanopore and the electrostatic interaction between the nanopore and particles. That also explains the asymmetry of particle transport velocity in the nanopore with a gate electrode. When the gate potential is negative, or the gate electrode length is small, the maximum deformation of the particles is increased. The field-effect regulation in the nanopore provides an active and compatible method for nanopore detection, and provides a convenient method for the active control of the particle deformation in the nanopore.

An electrodiffusive neuron-extracellular-glia model for exploring the genesis of slow potentials in the brain

Within the computational neuroscience community, there has been a focus on simulating the electrical activity of neurons, while other components of brain tissue, such as glia cells and the extracellular space, are often neglected. Standard models of extracellular potentials are based on a combination of multicompartmental models describing neural electrodynamics and volume conductor theory. Such models cannot be used to simulate the slow components of extracellular potentials, which depend on ion concentration dynamics, and the effect that this has on extracellular diffusion potentials and glial buffering currents. We here present the electrodiffusive neuron-extracellular-glia (edNEG) model, which we believe is the first model to combine compartmental neuron modeling with an electrodiffusive framework for intra- and extracellular ion concentration dynamics in a local piece of neuro-glial brain tissue. The edNEG model (i) keeps track of all intraneuronal, intraglial, and extracellular ion concentrations and electrical potentials, (ii) accounts for action potentials and dendritic calcium spikes in neurons, (iii) contains a neuronal and glial homeostatic machinery that gives physiologically realistic ion concentration dynamics, (iv) accounts for electrodiffusive transmembrane, intracellular, and extracellular ionic movements, and (v) accounts for glial and neuronal swelling caused by osmotic transmembrane pressure gradients. The edNEG model accounts for the concentration-dependent effects on ECS potentials that the standard models neglect. Using the edNEG model, we analyze these effects by splitting the extracellular potential into three components: one due to neural sink/source configurations, one due to glial sink/source configurations, and one due to extracellular diffusive currents. Through a series of simulations, we analyze the roles played by the various components and how they interact in generating the total slow potential. We conclude that the three components are of comparable magnitude and that the stimulus conditions determine which of the components that dominate.

A Kirchhoff-Nernst-Planck framework for modeling large scale extracellular electrodiffusion surrounding morphologically detailed neurons

Many pathological conditions, such as seizures, stroke, and spreading depression, are associated with substantial changes in ion concentrations in the extracellular space (ECS) of the brain. An understanding of the mechanisms that govern ECS concentration dynamics may be a prerequisite for understanding such pathologies. To estimate the transport of ions due to electrodiffusive effects, one must keep track of both the ion concentrations and the electric potential simultaneously in the relevant regions of the brain. Although this is currently unfeasible experimentally, it is in principle achievable with computational models based on biophysical principles and constraints. Previous computational models of extracellular ion-concentration dynamics have required extensive computing power, and therefore have been limited to either phenomena on very small spatiotemporal scales (micrometers and milliseconds), or simplified and idealized 1-dimensional (1-D) transport processes on a larger scale. Here, we present the 3-D Kirchhoff-Nernst-Planck (KNP) framework, tailored to explore electrodiffusive effects on large spatiotemporal scales. By assuming electroneutrality, the KNP-framework circumvents charge-relaxation processes on the spatiotemporal scales of nanometers and nanoseconds, and makes it feasible to run simulations on the spatiotemporal scales of millimeters and seconds on a standard desktop computer. In the present work, we use the 3-D KNP framework to simulate the dynamics of ion concentrations and the electrical potential surrounding a morphologically detailed pyramidal cell. In addition to elucidating the single neuron contribution to electrodiffusive effects in the ECS, the simulation demonstrates the efficiency of the 3-D KNP framework. We envision that future applications of the framework to more complex and biologically realistic systems will be useful in exploring pathological conditions associated with large concentration variations in the ECS.

Conic shapes have higher sensitivity than cylindrical ones in nanopore DNA sequencing

Nanopores have emerged as helpful research tools for single molecule detection. Through continuum modeling, we investigated the effects of membrane thickness, nanopore size, and pore shape on current signal characteristics of DNA. The simulation results showed that, when reducing the pore diameter, the amplitudes of current signals of DNA increase. Moreover, we found that, compared to cylindrically shaped nanopores, conical-shaped nanopores produce greater signal amplitudes from biomolecules translocation. Finally, we demonstrated that continuum model simulations for the discrimination of DNA and RNA yield current characteristics approximately consistent with experimental measurements and that A-T and G-C base pairs can be distinguished using thin conical solid-state nanopores. Our study not only suggests that computational approaches in this work can be used to guide the designs of nanopore for single molecule detection, but it also provides several possible ways to improve the current amplitudes of nanopores for better resolution.

Incorporating Born solvation energy into the three-dimensional Poisson-Nernst-Planck model to study ion selectivity in KcsA K^{+} channels

Potassium channels are much more permeable to potassium than sodium ions, although potassium ions are larger and both carry the same positive charge. This puzzle cannot be solved based on the traditional Poisson-Nernst-Planck (PNP) theory of electrodiffusion because the PNP model treats all ions as point charges, does not incorporate ion size information, and therefore cannot discriminate potassium from sodium ions. The PNP model can qualitatively capture some macroscopic properties of certain channel systems such as current-voltage characteristics, conductance rectification, and inverse membrane potential. However, the traditional PNP model is a continuum mean-field model and has no or underestimates the discrete ion effects, in particular the ion solvation or self-energy (which can be described by Born model). It is known that the dehydration effect (closely related to ion size) is crucial to selective permeation in potassium channels. Therefore, we incorporated Born solvation energy into the PNP model to account for ion hydration and dehydration effects when passing through inhomogeneous dielectric channel environments. A variational approach was adopted to derive a Born-energy-modified PNP (BPNP) model. The model was applied to study a cylindrical nanopore and a realistic KcsA channel, and three-dimensional finite element simulations were performed. The BPNP model can distinguish different ion species by ion radius and predict selectivity for K^{+} over Na^{+} in KcsA channels. Furthermore, ion current rectification in the KcsA channel was observed by both the PNP and BPNP models. The I-V curve of the BPNP model for the KcsA channel indicated an inward rectifier effect for K^{+} (rectification ratio of ∼3/2) but indicated an outward rectifier effect for Na^{+} (rectification ratio of ∼1/6).

CHARMM-GUI 10 years for biomolecular modeling and simulation

CHARMM-GUI, http://www.charmm-gui.org, is a web-based graphical user interface that prepares complex biomolecular systems for molecular simulations. CHARMM-GUI creates input files for a number of programs including CHARMM, NAMD, GROMACS, AMBER, GENESIS, LAMMPS, Desmond, OpenMM, and CHARMM/OpenMM. Since its original development in 2006, CHARMM-GUI has been widely adopted for various purposes and now contains a number of different modules designed to set up a broad range of simulations: (1) PDB Reader & Manipulator, Glycan Reader, and Ligand Reader & Modeler for reading and modifying molecules; (2) Quick MD Simulator, Membrane Builder, Nanodisc Builder, HMMM Builder, Monolayer Builder, Micelle Builder, and Hex Phase Builder for building all-atom simulation systems in various environments; (3) PACE CG Builder and Martini Maker for building coarse-grained simulation systems; (4) DEER Facilitator and MDFF/xMDFF Utilizer for experimentally guided simulations; (5) Implicit Solvent Modeler, PBEQ-Solver, and GCMC/BD Ion Simulator for implicit solvent related calculations; (6) Ligand Binder for ligand solvation and binding free energy simulations; and (7) Drude Prepper for preparation of simulations with the CHARMM Drude polarizable force field. Recently, new modules have been integrated into CHARMM-GUI, such as Glycolipid Modeler for generation of various glycolipid structures, and LPS Modeler for generation of lipopolysaccharide structures from various Gram-negative bacteria. These new features together with existing modules are expected to facilitate advanced molecular modeling and simulation thereby leading to an improved understanding of the structure and dynamics of complex biomolecular systems. Here, we briefly review these capabilities and discuss potential future directions in the CHARMM-GUI development project. © 2016 Wiley Periodicals, Inc.

Introducing biomimetic shear and ion gradients to microfluidic spinning improves silk fiber strength

Silkworm silk is an attractive biopolymer for biomedical applications due to its high mechanical strength and biocompatibility; as a result, there is increasing interest in scalable devices to spin silk and recombinant silk so as to improve and customize their properties for diverse biomedical purposes (Vepari and Kaplan 2007 Prog. Polym. Sci. 32 ). While artificial spinning of regenerated silk fibroins adds tunability to properties such as degradation rate and surface functionalization, the resulting fibers do not yet approach the mechanical strength of native silkworm silk. These drawbacks reduce the applicability and attractiveness of artificial silk (Kinahan et al 2011 Biomacromolecules 12 ). Here, we used computational fluid dynamic simulations to incorporate shear in tandem with biomimetic ion gradients by coupling a modular novel glass microfluidic device to our previous co-axial flow device. Fibers spun with this combined apparatus demonstrated a significant increase in mechanical strength compared to fibers spun with the basic apparatus alone, with a three-fold increase in Young's modulus and extensibility and a twelve-fold increase in toughness. These results thus demonstrate the critical importance of ionic milieu and shear stress in spinning strong fibers from solubilized silk fibroin.

Effect of Ionic Diffusion on Extracellular Potentials in Neural Tissue

Recorded potentials in the extracellular space (ECS) of the brain is a standard measure of population activity in neural tissue. Computational models that simulate the relationship between the ECS potential and its underlying neurophysiological processes are commonly used in the interpretation of such measurements. Standard methods, such as volume-conductor theory and current-source density theory, assume that diffusion has a negligible effect on the ECS potential, at least in the range of frequencies picked up by most recording systems. This assumption remains to be verified. We here present a hybrid simulation framework that accounts for diffusive effects on the ECS potential. The framework uses (1) the NEURON simulator to compute the activity and ionic output currents from multicompartmental neuron models, and (2) the electrodiffusive Kirchhoff-Nernst-Planck framework to simulate the resulting dynamics of the potential and ion concentrations in the ECS, accounting for the effect of electrical migration as well as diffusion. Using this framework, we explore the effect that ECS diffusion has on the electrical potential surrounding a small population of 10 pyramidal neurons. The neural model was tuned so that simulations over ∼100 seconds of biological time led to shifts in ECS concentrations by a few millimolars, similar to what has been seen in experiments. By comparing simulations where ECS diffusion was absent with simulations where ECS diffusion was included, we made the following key findings: (i) ECS diffusion shifted the local potential by up to ∼0.2 mV. (ii) The power spectral density (PSD) of the diffusion-evoked potential shifts followed a 1/f2 power law. (iii) Diffusion effects dominated the PSD of the ECS potential for frequencies up to several hertz. In scenarios with large, but physiologically realistic ECS concentration gradients, diffusion was thus found to affect the ECS potential well within the frequency range picked up in experimental recordings.

Charged Substrate and Product Together Contribute Like a Nonreactive Species to the Overall Electrostatic Steering in Diffusion-Reaction Processes

The Debye-Hückel limiting law is used to study the binding kinetics of substrate-enzyme system as well as to estimate the reaction rate of a electrostatically steered diffusion-controlled reaction process. It is based on a linearized Poisson-Boltzmann model and known for its accurate predictions in dilute solutions. However, the substrate and product particles are in nonequilibrium states and are possibly charged, and their contributions to the total electrostatic field cannot be explicitly studied in the Poisson-Boltzmann model. Hence the influences of substrate and product on reaction rate coefficient were not known. In this work, we consider all the charged species, including the charged substrate, product, and mobile salt ions in a Poisson-Nernst-Planck model, and then compare the results with previous work. The results indicate that both the charged substrate and product can significantly influence the reaction rate coefficient with different behaviors under different setups of computational conditions. It is interesting to find that when substrate and product are both considered, under an overall neutral boundary condition for all the bulk charged species, the computed reaction rate kinetics recovers a similar Debye-Hückel limiting law again. This phenomenon implies that the charged product counteracts the influence of charged substrate on reaction rate coefficient. Our analysis discloses the fact that the total charge concentration of substrate and product, though in a nonequilibrium state individually, obeys an equilibrium Boltzmann distribution, and therefore contributes as a normal charged ion species to ionic strength. This explains why the Debye-Hückel limiting law still works in a considerable range of conditions even though the effects of charged substrate and product particles are not specifically and explicitly considered in the theory.

Hybrid finite element and Brownian dynamics method for charged particles

Diffusion is often the rate-determining step in many biological processes. Currently, the two main computational methods for studying diffusion are stochastic methods, such as Brownian dynamics, and continuum methods, such as the finite element method. A previous study introduced a new hybrid diffusion method that couples the strengths of each of these two methods, but was limited by the lack of interactions among the particles; the force on each particle had to be from an external field. This study further develops the method to allow charged particles. The method is derived for a general multidimensional system and is presented using a basic test case for a one-dimensional linear system with one charged species and a radially symmetric system with three charged species.

A Stabilized Finite Element Method for Modified Poisson-Nernst-Planck Equations to Determine Ion Flow Through a Nanopore

The conventional Poisson-Nernst-Planck equations do not account for the finite size of ions explicitly. This leads to solutions featuring unrealistically high ionic concentrations in the regions subject to external potentials, in particular, near highly charged surfaces. A modified form of the Poisson-Nernst-Planck equations accounts for steric effects and results in solutions with finite ion concentrations. Here, we evaluate numerical methods for solving the modified Poisson-Nernst-Planck equations by modeling electric field-driven transport of ions through a nanopore. We describe a novel, robust finite element solver that combines the applications of the Newton's method to the nonlinear Galerkin form of the equations, augmented with stabilization terms to appropriately handle the drift-diffusion processes. To make direct comparison with particle-based simulations possible, our method is specifically designed to produce solutions under periodic boundary conditions and to conserve the number of ions in the solution domain. We test our finite element solver on a set of challenging numerical experiments that include calculations of the ion distribution in a volume confined between two charged plates, calculations of the ionic current though a nanopore subject to an external electric field, and modeling the effect of a DNA molecule on the ion concentration and nanopore current.

Electrodiffusive model for astrocytic and neuronal ion concentration dynamics

The cable equation is a proper framework for modeling electrical neural signalling that takes place at a timescale at which the ionic concentrations vary little. However, in neural tissue there are also key dynamic processes that occur at longer timescales. For example, endured periods of intense neural signaling may cause the local extracellular K⁺-concentration to increase by several millimolars. The clearance of this excess K⁺ depends partly on diffusion in the extracellular space, partly on local uptake by astrocytes, and partly on intracellular transport (spatial buffering) within astrocytes. These processes, that take place at the time scale of seconds, demand a mathematical description able to account for the spatiotemporal variations in ion concentrations as well as the subsequent effects of these variations on the membrane potential. Here, we present a general electrodiffusive formalism for modeling of ion concentration dynamics in a one-dimensional geometry, including both the intra- and extracellular domains. Based on the Nernst-Planck equations, this formalism ensures that the membrane potential and ion concentrations are in consistency, it ensures global particle/charge conservation and it accounts for diffusion and concentration dependent variations in resistivity. We apply the formalism to a model of astrocytes exchanging ions with the extracellular space. The simulations show that K⁺-removal from high-concentration regions is driven by a local depolarization of the astrocyte membrane, which concertedly (i) increases the local astrocytic uptake of K⁺, (ii) suppresses extracellular transport of K⁺, (iii) increases axial transport of K⁺ within astrocytes, and (iv) facilitates astrocytic relase of K⁺ in regions where the extracellular concentration is low. Together, these mechanisms seem to provide a robust regulatory scheme for shielding the extracellular space from excess K⁺.

When neurons generate electrical signals they release potassium ions (K⁺) into the extracellular space. During periods of intense neural activity, the local extracellular K⁺ may increase drastically. If it becomes too high, it can lead to neural dysfunction. Astrocytes (a kind of glial cells) are involved in preventing this from happening. Astrocytes can take up excess K⁺, transport it intracellularly, and release it in regions where the concentration is lower. This process is called spatial buffering, and a full mechanistic understanding of it is currently lacking. The aim of this work is twofold: First, we develop a formalism for modeling ion concentration dynamics in the intra- and extracellular space. The formalism is general, and could be used to simulate many cellular processes. It accounts for ion transports due to diffusion (along concentration gradients) as well as electrical migration (along voltage gradients). It extends previous, related formalisms, which have focused only on intracellular dynamics. Secondly, we apply the formalism to model how astrocytes exchange ions with the extracellular space. We conclude that the membrane mechanisms possessed by astrocytes seem optimal for shielding the extracellular space from excess K⁺, and provide a full mechanistic description of the spatial (K⁺) buffering process.

A molecular level prototype for mechanoelectrical transducer in mammalian hair cells

The mechanoelectrical transducer (MET) is a crucial component of mammalian auditory system. The gating mechanism of the MET channel remains a puzzling issue, though there are many speculations, due to the lack of essential molecular building blocks. To understand the working principle of mammalian MET, we propose a molecular level prototype which constitutes a charged blocker, a realistic ion channel and its surrounding membrane. To validate the proposed prototype, we make use of a well-established ion channel theory, the Poisson-Nernst-Planck equations, for three-dimensional (3D) numerical simulations. A wide variety of model parameters, including bulk ion concentration, applied external voltage, blocker charge and blocker displacement, are explored to understand the basic function of the proposed MET prototype. We show that our prototype prediction of channel open probability in response to blocker relative displacement is in remarkable accordance with experimental observation of rat cochlea outer hair cells. Our results appear to suggest that tip links which connect hair bundles gate MET channels.

Progress in developing Poisson-Boltzmann equation solvers

This review outlines the recent progress made in developing more accurate and efficient solutions to model electrostatics in systems comprised of bio-macromolecules and nano-objects, the last one referring to objects that do not have biological function themselves but nowadays are frequently used in biophysical and medical approaches in conjunction with bio-macromolecules. The problem of modeling macromolecular electrostatics is reviewed from two different angles: as a mathematical task provided the specific definition of the system to be modeled and as a physical problem aiming to better capture the phenomena occurring in the real experiments. In addition, specific attention is paid to methods to extend the capabilities of the existing solvers to model large systems toward applications of calculations of the electrostatic potential and energies in molecular motors, mitochondria complex, photosynthetic machinery and systems involving large nano-objects.

Hybrid finite element and Brownian dynamics method for diffusion-controlled reactions

Diffusion is often the rate determining step in many biological processes. Currently, the two main computational methods for studying diffusion are stochastic methods, such as Brownian dynamics, and continuum methods, such as the finite element method. This paper proposes a new hybrid diffusion method that couples the strengths of each of these two methods. The method is derived for a general multidimensional system, and is presented using a basic test case for 1D linear and radially symmetric diffusion systems.

Poisson-Nernst-Planck equations for simulating biomolecular diffusion-reaction processes II: size effects on ionic distributions and diffusion-reaction rates

The effects of finite particle size on electrostatics, density profiles, and diffusion have been a long existing topic in the study of ionic solution. The previous size-modified Poisson-Boltzmann and Poisson-Nernst-Planck models are revisited in this article. In contrast to many previous works that can only treat particle species with a single uniform size or two sizes, we generalize the Borukhov model to obtain a size-modified Poisson-Nernst-Planck (SMPNP) model that is able to treat nonuniform particle sizes. The numerical tractability of the model is demonstrated as well. The main contributions of this study are as follows. 1), We show that an (arbitrarily) size-modified PB model is indeed implied by the SMPNP equations under certain boundary/interface conditions, and can be reproduced through numerical solutions of the SMPNP. 2), The size effects in the SMPNP effectively reduce the densities of highly concentrated counterions around the biomolecule. 3), The SMPNP is applied to the diffusion-reaction process for the first time, to our knowledge. In the case of low substrate density near the enzyme reactive site, it is observed that the rate coefficients predicted by SMPNP model are considerably larger than those by the PNP model, suggesting both ions and substrates are subject to finite size effects. 4), An accurate finite element method and a convergent Gummel iteration are developed for the numerical solution of the completely coupled nonlinear system of SMPNP equations.

Adaptive Skin Meshes Coarsening for Biomolecular Simulation

In this paper, we present efficient algorithms for generating hierarchical molecular skin meshes with decreasing size and guaranteed quality. Our algorithms generate a sequence of coarse meshes for both the surfaces and the bounded volumes. Each coarser surface mesh is adaptive to the surface curvature and maintains the topology of the skin surface with guaranteed mesh quality. The corresponding tetrahedral mesh is conforming to the interface surface mesh and contains high quality tetrahedral that decompose both the interior of the molecule and the surrounding region (enclosed in a sphere). Our hierarchical tetrahedral meshes have a number of advantages that will facilitate fast and accurate multigrid PDE solvers. Firstly, the quality of both the surface triangulations and tetrahedral meshes is guaranteed. Secondly, the interface in the tetrahedral mesh is an accurate approximation of the molecular boundary. In particular, all the boundary points lie on the skin surface. Thirdly, our meshes are Delaunay meshes. Finally, the meshes are adaptive to the geometry.

Poisson-Nernst-Planck Equations for Simulating Biomolecular Diffusion-Reaction Processes I: Finite Element Solutions

In this paper we developed accurate finite element methods for solving 3-D Poisson-Nernst-Planck (PNP) equations with singular permanent charges for electrodiffusion in solvated biomolecular systems. The electrostatic Poisson equation was defined in the biomolecules and in the solvent, while the Nernst-Planck equation was defined only in the solvent. We applied a stable regularization scheme to remove the singular component of the electrostatic potential induced by the permanent charges inside biomolecules, and formulated regular, well-posed PNP equations. An inexact-Newton method was used to solve the coupled nonlinear elliptic equations for the steady problems; while an Adams-Bashforth-Crank-Nicolson method was devised for time integration for the unsteady electrodiffusion. We numerically investigated the conditioning of the stiffness matrices for the finite element approximations of the two formulations of the Nernst-Planck equation, and theoretically proved that the transformed formulation is always associated with an ill-conditioned stiffness matrix. We also studied the electroneutrality of the solution and its relation with the boundary conditions on the molecular surface, and concluded that a large net charge concentration is always present near the molecular surface due to the presence of multiple species of charged particles in the solution. The numerical methods are shown to be accurate and stable by various test problems, and are applicable to real large-scale biophysical electrodiffusion problems.

Solutions to a reduced Poisson-Nernst-Planck system and determination of reaction rates

We study a reduced Poisson-Nernst-Planck (PNP) system for a charged spherical solute immersed in a solvent with multiple ionic or molecular species that are electrostatically neutralized in the far field. Some of these species are assumed to be in equilibrium. The concentrations of such species are described by the Boltzmann distributions that are further linearized. Others are assumed to be reactive, meaning that their concentrations vanish when in contact with the charged solute. We present both semi-analytical solutions and numerical iterative solutions to the underlying reduced PNP system, and calculate the reaction rate for the reactive species. We give a rigorous analysis on the convergence of our simple iteration algorithm. Our numerical results show the strong dependence of the reaction rates of the reactive species on the magnitude of its far field concentration as well as on the ionic strength of all the chemical species. We also find non-monotonicity of electrostatic potential in certain parameter regimes. The results for the reactive system and those for the non-reactive system are compared to show the significant differences between the two cases. Our approach provides a means of solving a PNP system which in general does not have a closed-form solution even with a special geometrical symmetry. Our findings can also be used to test other numerical methods in large-scale computational modeling of electro-diffusion in biological systems.

Kinetics of diffusion-controlled enzymatic reactions with charged substrates

The Debye-Hückel limiting law (DHL) has often been used to estimate rate constants of diffusion-controlled reactions under different ionic strengths. Two main approximations are adopted in DHL: one is that the solution of the linearized Poisson-Boltzmann equation for a spherical cavity is used to estimate the excess electrostatic free energy of a solution; the other is that details of electrostatic interactions of the solutes are neglected. This makes DHL applicable only at low ionic strengths and dilute solutions (very low substrate/solute concentrations). We show in this work that through numerical solution of the Poisson-Nernst-Planck equations, diffusion-reaction processes can be studied at a variety of conditions including realistically concentrated solutions, high ionic strength, and certainly with non-equilibrium charge distributions. Reaction rate coefficients for the acetylcholine-acetylcholinesterase system are predicted to strongly depend on both ionic strength and substrate concentration. In particular, they increase considerably with increase of substrate concentrations at a fixed ionic strength, which is open to experimental testing. This phenomenon is also verified on a simple model, and is expected to be general for electrostatically attracting enzyme-substrate systems.PACS Codes: 82.45.Tv, 87.15.VvMSC Codes: 92C30.

A list-based method for fast generation of molecular surfaces

In this paper we present a fast and reliable method for generation of molecular surfaces. While the method is readily applicable to van der Waals surface generation, we shall be focusing on solvent-accessible surfaces (SAS) and solvent-excluded surfaces (SES) of a molecule. A list-based method is utilized to represent and generate the union of multiple spheres with arbitrary radii, with which we are able to produce the SAS and SES of a molecule in a very efficient way. The surface generated is represented by a quadrilateral mesh, which can be easily converted into a triangular mesh if needed. Some results will be provided to demonstrate the speed and quality of our mesh generation algorithms.

Enzymatic activity versus structural dynamics: the case of acetylcholinesterase tetramer

The function of many proteins, such as enzymes, is modulated by structural fluctuations. This is especially the case in gated diffusion-controlled reactions (where the rates of the initial diffusional encounter and of structural fluctuations determine the overall rate of the reaction) and in oligomeric proteins (where function often requires a coordinated movement of individual subunits). A classic example of a diffusion-controlled biological reaction catalyzed by an oligomeric enzyme is the hydrolysis of synaptic acetylcholine (ACh) by tetrameric acetylcholinesterase (AChEt). Despite decades of efforts, the extent to which enzymatic efficiency of AChEt (or any other enzyme) is modulated by flexibility is not fully determined. This article attempts to determine the correlation between the dynamics of AChEt and the rate of reaction between AChEt and ACh. We employed equilibrium and nonequilibrium electro-diffusion models to compute rate coefficients for an ensemble of structures generated by molecular dynamics simulation. We found that, for the static initial model, the average reaction rate per active site is approximately 22-30% slower in the tetramer than in the monomer. However, this effect of tetramerization is modulated by the intersubunit motions in the tetramer such that a complex interplay of steric and electrostatic effects either guides or blocks the substrate into or from each of the four active sites. As a result, the rate per active site calculated for some of the tetramer structures is only approximately 15% smaller than the rate in the monomer. We conclude that structural dynamics minimizes the adverse effect of tetramerization, allowing the enzyme to maintain similar enzymatic efficiency in different oligomerization states.

Molecular surface-free continuum model for electrodiffusion processes

Incorporation of van der Waals interactions enables the continuum model of electrodiffusion in biomolecular system to avoid the artifacts of introducing a molecular surface and the painful task of the surface mesh generation. Calculation examples show that the electrostatics, diffusion-reaction kinetics, and molecular surface defined as an isosurface of a certain density distribution can be extracted from the solution of the Poisson-Nernst-Planck equations using this model. The molecular surface-free model enables a wider usage of some modern numerical methodologies such as finite element methods for biomolecular modeling, and sheds light on a new paradigm of continuum modeling for biomolecular systems.

Continuum simulations of acetylcholine consumption by acetylcholinesterase: a Poisson-Nernst-Planck approach

The Poisson-Nernst-Planck (PNP) equation provides a continuum description of electrostatic-driven diffusion and is used here to model the diffusion and reaction of acetylcholine (ACh) with acetylcholinesterase (AChE) enzymes. This study focuses on the effects of ion and substrate concentrations on the reaction rate and rate coefficient. To this end, the PNP equations are numerically solved with a hybrid finite element and boundary element method at a wide range of ion and substrate concentrations, and the results are compared with the partially coupled Smoluchowski-Poisson-Boltzmann model. The reaction rate is found to depend strongly on the concentrations of both the substrate and ions; this is explained by the competition between the intersubstrate repulsion and the ionic screening effects. The reaction rate coefficient is independent of the substrate concentration only at very high ion concentrations, whereas at low ion concentrations the behavior of the rate depends strongly on the substrate concentration. Moreover, at physiological ion concentrations, variations in substrate concentration significantly affect the transient behavior of the reaction. Our results offer a reliable estimate of reaction rates at various conditions and imply that the concentrations of charged substrates must be coupled with the electrostatic computation to provide a more realistic description of neurotransmission and other electrodiffusion and reaction processes.