Spatial Buffering Mechanism : Mathematical Model and Computer Simulations

It is generally accepted that the spatial buffering mechanism is important to buffer extracellular-space potassium in the brain-cell microenvironment. In the past, this phenomenon, generally associated with glial cells, has been treated analytically and numerically using a simplified one-dimensional description. The present study extends the previous research by using a novel numerical scheme for the analysis of potassium buffering mechanisms in the extracellular brain-cell microenvironment. In particular, a lattice-cellular automaton was employed to simulate a detailed two-compartment model of a two-dimensional brain-cell system. With this numerical approach, the present study elaborates upon previous theoretical work on spatial buffering (SB) by incorporating a more realistic structure of the brain cell microenvironment, which was not feasible earlier. We use the experimental paradigm consisting of iontophoretic injection of KCl to study the SB mechanism. Our simulations confirmed the results reported in the literature obtained by an averaged model. The results also show that the additional effects captured by a simplified twodimensional geometry do not alter significantly the conclusions obtained from the averaged model. The details of applying such a numerical method to the study of ion movements in cellular environments, as well as its potential for future study, are discussed. 2000 Mathematics Subject Classification. 92C05.

1. Introduction.The brain consists of cellular components (neurons and glial cells) and chemical ions (sodium, Na + , potassium, K + , and chloride, Cl − ).Regulation of the extracellular space (ECS) and intracellular space (ICS) concentrations of these ions is necessary for proper functioning of the brain through the electrical activity of neurons.For example, increases in extracellular K + concentration ([K + ] o ) from a resting level of 3 mM to 5 mM have been shown to lead to hyperexcitability of the surrounding neurons [17].Further increases in [K + ] o also may affect the efficiency of synaptic transmitter release [2].Accordingly, [K + ] o is tightly regulated, normally never exceeding 11 mM.
It was suggested that temporary K + accumulation can be caused by prolonged neuronal activity [29,30], iontophoretic K + injection [19,22], application of drugs, or pathological conditions such as hypoxia and ischemia [23,26].To restore resting K + levels, mechanisms by which K + concentration gradients are dissipated, superseding simple passive diffusion, have been demonstrated.In particular, three primary mechanisms have been emphasized for extracellular K + removal around glial cells [6]: enhanced K + transport by the Na + / K + -ATPase, passive KCl cotransport, and current-mediated K + entry by means of K + channels.The latter mechanism, termed spatial buffering (SB) by K + channels, will be the primary focus of this paper.
First introduced by Kuffler and colleagues in their studies on Necturus glial cells [25], SB is the direct consequence of basic biophysical properties.In a glial syncytium connected via gap junctions, the membrane of connected cells tend to remain isopotential.An increased [K + ] o causes a local depolarization across the glial cell membrane, which by the cable properties of the membrane, will spread electrotonically down the cells to more remote regions of the syncytium.Local currents are formed by the unbalanced spatial distribution of membrane potential differences and completed through the ECS.Since glial cell membranes exhibit a high K + conductance, the asymmetric potential differences result in an influx of extracellular K + at the site of local depolarization and an efflux of K + into the ECS at a more distant region where the [K + ] o is still near resting levels.The K + is passively transported from an ECS location of high [K + ] o and dispersed to the ECS of regions with low [K + ] o .The circuit current loop is closed by intracellular and extracellular ionic current flows, primarily those of Na + and Cl − .It has been shown that this passive buffering mechanism is energy-independent and generally more efficient than diffusion for transporting K + through the interstitial space [14,16].
At rest, the glial membrane is selectively permeable to K + ions and passively obeys the Nernst equation [6].The K + conductance can be attributed to at least four different voltage-dependent K + channels: the inward rectifier, Kir; delayed rectifier, Kd; transient A-type, K A ; and Ca 2+ -activated channel, K Ca [28].Of particular interest to the study of SB is the Kir channel, as it maintains an open configuration at rest and can be activated at more negative potentials.In particular, the inward rectifier channel exhibits a different response to a hyperpolarization in contrast to a depolarization.In the case of the former, the channel retains a high open probability and thus has a high conductance compared to that during depolarization.This result for the inward rectification properties of the channel allows a preferential influx of K + into the glial cell.The efficiency of SB via Kir channels is enhanced by increases in the influx of K + in regions of high [K + ] o [21] and by the spreading of the membrane depolarization along the glial cell and syncytium [1].
It should be noted that the distribution of the Kir channels in certain glial cells, such as the Müller cells in the retina of amphibians, has been shown to be nonuniform with a preponderance at the endfoot of the vitreous humor [5,20,21,27].This nonuniformity may manifest itself in a higher K + conductance at the endfoot processes.It has been suggested that the asymmetric distribution of K + conductances most efficiently directs a K + efflux through the high-conductance areas in situations involving only a few space constants [4,13].The effects of these asymmetries are addressed in the present study.
Specifically, we study K + dynamics and homeostasis in the brain-cell microenvironment.We build a theoretical model to examine the migration of increased [K + ] o through the system in response to the injection of KCl, dealing principally with the SB mechanism.Generally, the ECS and ICS of the brain-cell microenvironment have complicated geometrical structures.This intricate array of neurons and glial cells constrain the migration of ions in the ECS to a tortuous and hindered pathway.In the study of K + dynamics, the influence of volume fraction and the detailed structure of the ECS, which is described macroscopically by tortuosity, have been analyzed by treating the brain as a porous medium [24].In previous studies, the glial syncytium geometry was modelled as a one-dimensional equivalent cable [6,15].While this simplification is conducive to an analytical examination of glial SB using continuous descriptions, it does not directly address the complicated geometry of a glial syncytium.Presumably, the geometrical component of SB may influence the buffer mechanism, as it dictates the form of the possible trajectories of the migrating particles.The various buffering mechanisms will dissipate the elevated K + .Any such dissipating mechanisms will require the movement of the ions, whether through the ECS, the ICS, or both.As such, the structure of both compartments dictates, in part, the evolution of the migrating particles.Accordingly, a model that effectively allows for a more accurate study of the geometrical components of the brain-cell microenvironment while incorporating the appropriate ICS, ECS, and membrane properties-such as those described in previous theoretical work [6,15]-would allow for a more comprehensive and accurate study of K + buffering mechanisms.In the present work, a two-dimensional model of the glial syncytium is generated, taking into account volume fractions of the ECS and ICS.Through our computer simulations using a lattice cellular automata, we are able to study more realistic geometries of the brain-cell microenvironment as well as incorporate an increased number of components of the K + homeostatic mechanisms.
At a macroscopic level, the migration of ions in the ECS and ICS can be viewed as a diffusion process coupled with nonlinear sources and sinks governed by a nonlinear diffusion equation.In contrast, at a microscopic level, the migration of ions corresponds to a random walk of the ions as they collide with surrounding molecules and cell membranes.This microscopic property allows us to make use of a different approach, namely, a lattice cellular automata (LCA), in the study of ion diffusion [7,9,11,12].The concept behind the LCA is to design a self-reproducing, spatially discrete system that evolves in discrete time by following simple rules of evolution.This generates an artificial microscopic structure that reproduces the system's behavior on a macroscopic scale.
There are several advantages to using an LCA to study the migration of ions in the brain.An important advantage is that it can deal with complicated geometries and interface conditions.Due to the microscopic interpretation of the ion dynamics in the LCA, the membrane fluxes can be taken into account more naturally than when a continuous description is used.The continuous description loses the basic intuition of the microscopic phenomena, particularly in mimicking ionic currents through geometrically complex membrane boundaries.Another advantage is that the evolution of several different ions can be simulated concurrently, with each ion being updated according to its own set of rules in the ECS and ICS.The ion concentrations and fluxes are coupled by the dependence of the membrane currents on the membrane potential and concentration profiles.The LCA is a self-reproducing system, and it will evolve to its homeostatic state similar to a realistic brain-cell system after various external stimuli and/or changes in the environment have been introduced.The LCA model can readily be improved by adding more realistic rules of evolution.A major advantage of the LCA is that it is a stable computational algorithm, whereas other computational schemes for continuum models, e.g., finite difference methods, have temporal and spatial step restrictions for stability.
In section 2, the mathematical modelling of the diffusion of ions in the ICS and the ECS and movements of ions across membranes are addressed.First, a discussion of the lattice structure representative of the brain-cell microenvironment is presented.Then, the evolution rules of the lattice dynamics, as well as the membrane conditions, are explained.Finally, we define the governing equations of the membrane potential and currents.Each membrane current is modelled by Ohm's law in addition to a passive pumping flux.In section 3, the simulation results are described.We begin with the migration of ions for a single cell environment.This simulation is used to explicitly demonstrate the properties of our virtual glial cell, how the cell maintains its resting state and how the system returns to its stable resting state after imposing a transient perturbation from its resting state.Next, we address the issue of the spatial buffer.This study proceeds through the analysis of two different geometries: a multicellular glial system constructed with a realistic volume fraction and elliptical glial cells arranged in an array.In the case studies, the K + buffering mechanisms, including the cable and endfoot effects, are explored using the iontophoretic paradigm.The paper concludes with a discussion in section 4.

2.1.
The brain-cell microenvironment.Initially, we consider a general twocompartment model, which could represent a homogeneous ECS and glial ICS, for example.By neglecting the detailed components within these two compartments, ions are assumed to undergoes free diffusion in the ECS and ICS.The concentrations of each distinct ion are governed by the diffusion equation where t is time, C is the ionic concentration, and D is the diffusion coefficient of the ion.We assign diffusion coefficients to each kind of ion in aqueous solution and assume that they are the same in the ICS and ECS.The diffusion of ions in the ECS of brain tissue is analogous to diffusion in porous media composed of one permeant phase and one impermeant phase.
In addition, we assume that all electrophysiological features for each cell are uniformly distributed over the cell membrane for most of the simulations.In particular, the densities of the relevant ion channels are assumed to be evenly distributed over each cell membrane.However, a nonuniform distribution of Kir channels can be accounted for easily by modifying the description of the LCA.
The model system then has three components: the ICS, the ECS, and the membrane that separates the two compartments.Although the brain-cell microenvironment has a complex three-dimensional structure, we will consider only the twodimensional case, analogous to a thin section through the brain tissue.For each cell, the membrane is a closed curve that encircles a set of ICS lattice nodes.The ions diffuse freely in the ECS and in the ICS and can move across the membrane in the form of membrane currents.Two factors that influence the migration of ions in the brain are the volume fraction of the ECS and the geometric shapes of the cells.In the current study, we consider the glial cells as the only cellular components of the brain-cell microenvironment, ignoring the neuronal membrane contributions.Experiments indicate ECS and glial volume fractions are approximately 0.2 and 0.4, respectively [6,10,18].A detailed description of the cellular geometries in this study are given in section 3.
2.2.Migration of ions in the ICS and ECS.We design a two-dimensional LCA by defining a square lattice representation of our model system.The lattice system consists of M + 1 sublattices, with M being the number of cells.One lattice is assigned to the uninterrupted ECS, while each of the other independent lattices is assigned to the ICS of each individual glial cell.The sublattices are disjoint, but the ECS lattice is coupled to all of the other sublattices through the membranes.The membranes themselves do not occupy any lattice nodes, but they are assumed to lie between adjacent ICS and ECS nodes.If there are no membrane currents, then the ionic particles placed on the nodes of a given sublattice will migrate independently of particles on other sublattices.However, if there are membrane currents, then ionic particles may move between the ICS and ECS lattices.The complicated geometrical structure of the brain-cell microenvironment is embodied in the shape of the membrane, that is, the shapes of the ECS and ICS lattice structures.This shape helps define the boundary conditions that influence the migration of ions across the interface between the ICS and ECS.In this study, the distance between neighboring nodes, λ, in each square lattice is taken to be 0.1 µm.
As described in Dai and Miura [9], there are three iteration steps involved in the evolution of the LCA: injection, collision, and propagation.During the initial injection step, ions are placed at each node of the lattice with a certain density distribution to mimic the realistic microenvironment of the brain.Particles that may move out of the system boundary are compensated for by the constant injection of particles at a certain rate through the four boundaries such that the density of ions in the ICS and the ECS are held at the resting values when the whole system is at rest.Furthermore, in the study of the spatial buffering phenomenon, additional ionic particles are injected at certain nodes in the ECS lattice for a given amount of time, mimicking the iontophoresis paradigm.In the second step, at each node, to replicate the random walk of the diffusion process, the migrating ions undergo collisions subject to specified collision rules.The third step is the migration or propagation of the particles within a given lattice and between adjacent lattices via the membrane currents.In this step, a particle may move from its present node to one of its nearest neighboring nodes.If a particle moves towards the boundary between the ECS and ICS, it is possible for it to migrate across the membrane and enter the adjacent lattice.The probability that this will occur is determined by the membrane current.In each time step, the above three steps operate at each node, and the distribution of particles is updated.
There are four directions, v j (j = 1,2,3,4), in which an ion travelling on the lattice can move, as shown schematically in Figure 1.Also, for a given time step, a particle could remain at a given node.As a particle can enter or leave a node in one of five ways, including not moving at all, there are 25 combinations of incoming-outgoing motions that exist at each node for a particle at every two time steps.Figure 1.A particle migrating on the lattice can travel in four directions, given by the directional vectors v j for j = 0,1,2,3,4.The direction v 0 corresponds to a particle that remains still at a given position.
Each node is labeled by the discrete spatial vector R = (x, y).We define the number of the particles of a given ion at each lattice node by where N j is the number of particles of a given ion that are moving in the direction of the unit vector v j with v 0 = 0.With this notation, the dynamics of our LCA model are given by two steps: collision Here τ is the time step of evolution.The p l,j is the probability that a particle which entered the node at R with direction v j is permuted to a particle with direction v l .The sum of these probabilities must satisfy For the sake of simplicity, we restrict ourselves to cyclic permutations.Namely, the probability of an incoming particle being permuted to direction v l is p l , which is independent of its incoming direction, that is, p l,j = p l .Furthermore, it is reasonable to assume that the lattice spacing is uniform, and hence, for a given ion species, p l , l = 1, 2, 3, 4, is the same for all nodes and directions in the same lattice The value of p 0 is determined from the diffusion coefficient of each ion species in a compartment.When both τ and λ are small and τ = O(λ 2 ) [7,9], it can be proved by employing the Chapman-Enskog expansion that the macroscopic law governing the evolution of the LCA is the diffusion equation.The diffusion coefficient of a given ion species in a specific compartment is given by which allows for a mapping of parameters between the LCA and the macroscopic diffusion process.We are able to fit realistic values of the diffusion coefficients for each ion species by defining the time scale, τ , and value of p 0 .In this study, where we treat only K + and Cl − movements, we select p 0,Cl = 0 since chloride ions have a larger diffusion coefficient than potassium ions.It then follows that the diffusive time step is determined by where D Cl is given in Table 1.The probability, p 0,K + , then can be determined from (7).
In determining the size of the time step used in the numerical experiments, it is important to account for the differences between the speeds of the particles as a result of diffusion and the speed of propagation of the membrane potential.We can compare these speeds by comparing the associated diffusion coefficients.For the membrane potential, the effective diffusion coefficient is characterized by u/C m in (16).Note that u/C m is equal to Λ 2 /Γ where Λ and Γ are the electrical space and time constants, respectively, of the glial cell membrane under resting conditions.It follows from the parameters presented in Table 1 that the membrane potential diffusion coefficient u/C m is approximately 4000 times the particle diffusion coefficient.Consequently, the above-determined time step must be modified to account for this difference in speeds by dividing by this factor of approximately 4000.Thus, each numerical iteration should correspond to a time step That is, in terms of the LCA, a diffusion step would occur once in every 4000 iterations of the membrane dynamics.If this differential between the membrane and diffusive time scales were implemented, it would increase the computational time dramatically.To run the system in a more realistic computational scenario, we allow for diffusion every 40 iterations of membrane dynamics.While the above LCA is valid away from the cell membrane, that is, in the interiors of the ECS and ICS, we need to describe boundary conditions at two types of boundaries: the boundary of the whole system and the boundary between the ICS and the ECS.Though the four border sides of the whole system may end up in either the ICS or ECS, in the simulations the simplified model cells were chosen such that the surrounding border nodes were completely in the ECS.This 73 mM [6] Resting 146 mM Donnan equil Static K + permutation probability p 0,K 0.082 eqn 7 Static Cl − permutation probability p 0,Cl 0

mV Nernst eqn
Apparent ECS resistivity r 0 1000 Ω cm [6] Apparent ICS resistivity r i 4000 Ω cm [6] Glial membrane surface density a 4884 cm −1 [6] KCl uptake time constant Note: The effective diffusion coefficient for the membrane potential is given by u/C m in equation ( 16) and is equal to Λ 2 /Γ where Λ and Γ are space and time constants, respectively.Since Λ is given by [ Compared with the diffusion coefficient for K + , we see that the membrane propagation velocity is approximately 4000 times that of diffusion.implies that there were no broken cells in the system since all ICS lattices were fully encircled by membrane.Accordingly, the density of each kind of ion along the bordering ECS was set to the resting ECS levels.By employing this requirement, we stipulate that there should be a constant injection of particles to compensate for the leakage of ions from the four sides.
The boundary condition between the ICS and ECS lattices is determined by the shape of the cell membrane.There are six groups of boundaries between the ICS and ECS due to the geometric shape of the ICS lattice node boundary.One ECS lattice boundary node may border one, two, or three ICS lattice nodes, labeled I, II, and III, respectively, in Figure 2. Similarly, one ICS lattice boundary node may have one, two, or three ECS neighbors.In group I, there are four possible configurations in which an ECS node may be bordered by one ICS node.These configurations correspond to one of the four neighboring nodes being an ICS node.Similarly, there are six and four arrangements for groups II and III, respectively.Therefore, to account for the detailed structure of the ICS and ECS on the membrane current, a total of 28 configurations must be taken into account.
To deal with the complicated boundary conditions, we mimic the movement of ionic particles in the brain.To do this, we assume that a particle is sitting on an ECS boundary node.At the next time step, it can remain at this node or move in one of four directions.If the direction points to an ECS node, meaning a node within its own lattice, it will move freely to its neighbor.Conversely, if the direction points to the boundary, then there are two possible outcomes.If there is no net membrane current for this ion species, then the particle cannot propagate through the membrane and is bounced back.Alternatively, if there is a nonzero membrane current for this ion species, then the particle may bounce back or propagate through the membrane and enter the ICS.The probability of being bounced back or moving into the ICS is determined by the biophysical properties of the membrane, the local membrane potential, and the local ion concentration in the ECS and ICS.Whatever the shape of the boundary, if one direction from the ECS boundary node points to the cell boundary, then there must exist a corresponding opposite direction at the matching ICS boundary node.Once a particle has moved from the ECS to the ICS, it becomes a particle that may freely move within the ICS lattice.
The LCA rule at the boundary can be explicitly determined from the above description.We first consider case I in Figure 2, where the boundary is in the direction of v 4 .Note that due to the assumption of cyclic permutation, we have N j (r, t) = p j (r)C(r, t).If there is no membrane current, the LCA rule for the node reads as where R is the ECS node of interest (see Figure 2).The first term reflects the sum of the particles that migrate toward node R from adjacent nodes.We note that since there is no membrane current, any particles migrating from the ICS node R are not considered.The second term corresponds to the density of the particles that would have otherwise migrated in the direction of R but have bounced back.Similarly, the equation for the corresponding ICS node, R , is If there is a net migration of ions across the membrane, then the above LCA needs to be modified to incorporate the passage of these ions through the membrane.As elaborated on below, we consider two types of ionic currents: current through the potassium inward rectifier Kir channel, I Kir ; and a pump current, I p .Combining these additional currents with those described above, the density at the boundary ECS node evolves according to (10) Explicitly including the j = 2 term in the summation, equation ( 10) can be rewritten as We now make the simplifying assumption that the discrete membrane lies midway between the extracellular and intracellular nodes, that is, the membrane node lies at a position given by R− v 2 2 .To ensure that this does give the correct macroscopic dynamics, namely, that C y = ∂C/∂y = I Kir + I p , we first expand C(R, t + τ ) and p j (R − v j )C(R − v j , t) for j = 1, 3, 4 at node R and at time t in Taylor series for small = |v j | and τ : Finally, we expand , and time t: Hence, from the normalization condition (6) for the probabilities, equation ( 10) is equivalent to Assuming that the time unit τ is proportional to 2 , that is, τ = k 2 for some positive constant k, and solving for ∂C/∂y(R, t), we obtain In the limit where → 0, equation ( 13) becomes the desired membrane boundary condition.
Similarly, for the ICS node R , the boundary LCA becomes ( 14) Following an analogous procedure as above for the ICS node, we obtain which also approaches the membrane boundary condition in the appropriate limit.
The corresponding results for the other three configurations of case I follow similarly.In the case of a type II boundary, we can separate it into two pairs of directions.Each pair of directions has the same rule as that of a type I boundary given above.Similarly, the type III boundary can be separated into three pairs.As a result, the complicated boundary can be determined by direct, simple rules; however, this requires that the exact pairs of 28 different configurations be known.
2.3.Membrane potential and currents.The number of particles that move across the boundary between the ICS and the ECS is governed by the membrane properties, the membrane potential, and the ionic concentrations.In terms of modelling, the membrane potential and the ionic currents are assumed to be continuous both in time and space.As treated above, the boundary conditions are determined by the concentrations of ions on its two sides, that is, by their density values at the boundary nodes in the ICS and ECS.Thus, the membrane serves to couple the ICS and ECS.
The glial membrane potential, v m , is defined by convention as the difference between the intracellular and extracellular potentials (v i −v o ).Membrane potential dynamics can be described by cable theory: where C m is the membrane capacitance, and I i is the membrane current through the ion channels of a given species i, determined by the respective channel variables and membrane potential.
To study the effects of the membrane potential dynamics, we use a finite difference method to discretize the continuous membrane.The membrane now is described by nodes with finite spacing.Each node of the membrane corresponds to one boundary node in the ICS and one in the ECS.We use the same time step τ as that of the LCA.We note that the diffusive time scale for the membrane potential is much faster (4000 times faster) than that of the ion diffusion in the ICS and ECS.We could use the pseudosteady solution of the membrane potential equation (16).However, since ( 16) is nonlinear (via the membrane current), iterative methods must be used.An alternative way is to solve the time-dependent problem, which is used in our simulation.This method may not be as fast as iterative methods, such as Newton's method; however, it is more robust.In addition, our method can be used readily to study other phenomena where dynamics of the membrane potential are important.
Moreover, because of the nearly exclusive permeability of glial cells to K + , we assume that the transmembrane current I is carried solely by K + .Although passive Cl − membrane flux is considered, we presume that it offers a negligible contribution to changes in the membrane current [6].For the constitutive relationship between the membrane current and membrane potential, we assume the Hodgkin-Huxley form of an ohmic relationship, namely, where E K is the glial equilibrium potential described by the Nernst equation for K + .The specific glial membrane K + conductance is given by g K .As outlined above, we consider only the Kir channel and thus make use of the empirically derived K + conductivity of the Kir channels given by Newman [21]: The resting membrane conductance is given by g o K and corresponds to the conductance when and ∆v = 0.The first term in braces describes the influence of ∆v on the rectification, and the second term reflects the impact of the membrane potential on the channel's opening probabilities.
In addition to current through the inward rectifier channels, passive flux of KCl is considered.According to the model of Boyle and Conway [3], variation in the concentration of both K + and Cl − in either the ICS and ECS will result in passive redistribution of KCl across the membrane.In their model, the final concentrations in both compartments obey a Donnan equilibrium [6].As such, we take the passive KCl movement to be proportional to the concentration difference.Note that the values of the resting densities have been chosen to satisfy a Donnan equilibrium.The change in the concentration, P , of a given ion species due to passive flux is given by where τ is the first-order uptake constant; C K,i and C Cl,i represent [K + ] i and [Cl − ] i concentrations, respectively, and C ∞ K,i is the [K + ] i concentration at the Donnan equilibrium.
Note that our treatment of diffusion through the LCA assumes electroneutrality.It then follows that in our synthetic brain-cell microenvironment, the Cl − ions act purely as a spectator species relative to the K + ions.Specifically, the transmembrane flux of Cl − follows the passive flux of K + and has no impact on the ultimate evolution of the system.In particular, it is solely the K + dynamics that govern the Cl − membrane current, and, therefore, the actual state of the Cl − ions with regard to their distribution across the membrane will not drive the chloride ions to reestablish their initial resting conditions, satisfying the Donnan equilibrium.Consider, for example, the situation where the intracellular and extracellular K + concentrations return to their equilibrium values ahead of the Cl − .In this case, the system has effectively returned to its resting state, and further iteration would show no change.However, the intracellular and extracellular Cl − densities are not necessarily at their equilibrium levels.They are held at these values because their movement depends on a K + concentration gradient, of which there is none.Considering the passive nature of the anion, it does not seem necessary to account for the Cl − ions; however, by monitoring the Cl − movement, we easily obtain the component of the total K + membrane current that consists of passive movement.Additionally, coupling the passive flux of K + with an anion gives a more reasonable biological treatment.As we are primarily concerned with potassium spatial buffering, whether the Cl − returns fully to the resting value does not have a significant impact on the system and the mechanisms of this study.However, to correct for the possibility of having the system evolve to a state where the Cl − particles are not at the initial resting values, we can add an additional Cl − pumping current.For simplicity, we can represent the pump current I Cl,p using an ohmic relationship, with constant conductance g Cl,p and a driving force determined by the difference between the present Cl − Nernst potential e Cl and the predetermined resting Cl − Nernst potential E o Cl : I Cl,p = g Cl,p (e Cl − E o Cl ).

Note that by the Donnan equilibrium condition, E o
Cl is equal to E o K .Hence, by incorporating this additional pumping current, the system will move toward its initial resting state.
In summary, we consider the membrane current of two ionic species, K + and Cl − .The K + current is the result of flux through Kir channels as well as passive redistribution governed by local fluctuations in concentration.The Cl − current, on the other hand, is assumed to result from the passive flux due to coupling with the K + current.An additional Cl − pumping current is incorporated to ensure that the Cl − concentration returns to its resting value.Moreover, the concentrations also can be changed by iontophoretic injection, which can be represented by a source or sink term.Combining the above currents, we can describe the changes in intracellular and extracellular concentrations of these two species.
Finally, returning to our treatment of the membrane boundary, if we know the membrane potential and total ionic current, then at each time step, the number of particles exchanged across the membrane is given by where ∆S is the area corresponding to a node of membrane, and F is the Faraday constant.Because the brain is modelled as a two-dimensional system, this parameter cannot be determined from realistic neurophysiological parameters.In the present study, we do not know exactly how many particles can be transported through the membrane.As an approximation, we assume the maximum number of ions passing through the membrane to be constant.

3.
Results.Table 1 lists the parameter values used for all our simulations.
3.1.Single circular cell.We begin our investigation of the SB properties of our two-dimensional LCA by analyzing how a single glial cell maintains a stable resting membrane potential and responds to perturbations.We construct a 100 x 30 lattice system with one circular cell placed at the center of the ECS environment (see Figure 3; note the different scales on the two axes).The ECS and ICS densities of both K + and Cl + are set to their respective resting levels.Accordingly, the system is initially in its resting state.The evolution of the system in response to particle injection can be monitored by recording the particle densities at nodes of the system as the system is iterated.The ICS and ECS recording sites, as well as the injection sites, are shown in Figure 3.In the single cell simulations, the K + ions were injected at the nodes indicated by diamonds in Figure 3 at a rate of 64 mM/msec for 3.125 msec.
3.1.1.Spatial buffering mechanism.The simulations were run with the above injection protocol.The injected K + ions can diffuse around the single glial cell or can enter into the ICS through the membrane current.As described above, the mechanisms by which K + ions may enter the cell are passive flux, driven by a transmembrane density gradient, and ohmic membrane current through inward rectifier channels.While the excess K + ions must still migrate to the borders for the system to return to its resting state, there is now the potential for other mechanisms, in addition to simple diffusion, to accelerate the diffusion of the particles.
As K + ions are injected, the local [K + ] o increases.At a segment of membrane near this site of increased [K + ] o , there is an increased (more positive) resting equilibrium potential.The local depolarization creates an electrical driving force on the K + ions, and because of the permeability of the membrane to potassium, the ions move across the membrane into the intracellular lattice.Concurrently, the deviation of the intracellular K + density results in the passive cotransport of KCl across the membrane.The cotransport acts in the direction opposite to the direction of the influx of K + through the inward rectifier channels, damping the impact of the inward channel current.Furthermore, the local depolarization moves electrotonically along the membrane away from the site of local depolarization.In so doing, segments of the membrane not immediately adjacent to the injection site are depolarized, although not to the same degree as at the injection site since there is some spatial dissipation of the potential as it propagates along the membrane (Figure 4).Because a glial cell has the tendency to remain isopotential, at segments of membrane further removed from the injection site, the resting membrane potential is now more negative than the actual membrane potential.As a result, there is an outwardly directed driving force on the K + .The K + permeability through the Kir channels allows for the cation to move through the membrane into the ECS lattice.Effectively, the initial iontophoretic injection of K + ions formed a current loop: in regions of high extracellular K + , there is positive current into the glial cell and a subsequent outward current away from the location of increased K + density.In our system, the current is closed by an ECS return current carried by Cl − .Throughout the cell, the movement of K + through the Kir channels driven by membrane depolarization is countered by the passive flux of K + .This description of the spatial buffering mechanism can be determined from analyzing the evolution of the system.First, we consider extracellular nodes located along the membrane.As expected, the K + ion concentrations at these nodes increase as a result of the injected ions being diffused throughout the system (Fig- ures 5a and 5c).The nodes in closer proximity to the stimulus are the first ones subjected to higher potassium levels.In our system, two components cause the increase in [K + ] o .First and predominantly, for the nodes relatively close to the injection point, direct diffusion of the particles results in the observed increase.As we move further away from the injection site, the aforementioned current loop plays a more predominant role.Because the propagation of the membrane potential depolarization along the membrane surface is faster than that of ion diffusion, the [K + ] o at further-removed sites is increased by the movement of K + out of the cell through Kir channels.Effectively, the cell has taken in K + ions from the region of increased ion concentration and then expelled K + ions at a distant region.Note that in Figure 5c, there appears to be a slight undershoot of [K] o at the sites distant from the location of injection.This may be caused by the nonlinear interaction between the membrane potential and cross-membrane transport and diffusion in ICS and ECS.Due to depolarization of the membrane potential, the potassium ions are driven out of the cell (flow from ICS to ECS) and carried away by diffusion in the ECS.When the membrane potential returns to its resting level, the ECS potassium concentration becomes lower than its equilibrium value, caused by the flux of potassium back into the cell as well as by the amount diffused away in the Figure 4.The local depolarization resulting from the K + ion injection spreads electrotonically along the membrane surface.Segments not immediately adjacent to the injection site thus are depolarized.As you move further from the source, the depolarization decreases.The indicated nodes refer to the geometry shown in Figure 3.
ECS during the depolarized stage.This appears as an undershoot when its value eventually comes back to the rest value at a slower diffusive time scale.While the potassium concentration increases at all extracellular nodes to varying degrees in SB, depending on location, both the magnitude and direction of the K + ion concentration changes at the intracellular nodes depend on their location relative to the injection site (Figures 5b and 5d).As described above, the increased K + ion levels result in the flux of the cation into the cell, raising the intracellular concentration accordingly.Further away, the membrane is depolarized while the concentrations remain transiently at the resting levels.This forces the intracellular K + ions to move out of the cell across the membrane, lowering the intracellular concentration, such that this depolarized potential becomes the new resting membrane potential.For this reason, we observe that at intracellular boundary nodes away from the injection site, there is a decrease in K + .The magnitude of this decrease is larger at the segments of membrane where diffusion has not already raised the extracellular concentration level.
3.1.2.Spatial buffering mechanism with increased endfoot conductance.So far, we have assumed that the cell membrane has uniform properties.In reality, however, the electrical properties of a glial cell membrane are not uniformly distributed.As described above, a specific example is the increased conductance at the endfeet of the Müller cells in the retina of amphibians [5,20,21,27].Obviously, it will be interesting to find out whether the increased endfoot conductance enhances or retards the diffusion of ions via the spatial buffering mechanism.To simulate the situation of increased endfoot conductance, we increased the Kir channel conductance by a factor of 25 from the value given in Table 1 along the region indicated by the solid curve in Figure 3.The simulation was run with the same initial conditions and injection protocol as in the previous section.In Figure 6, we have compared the ion concentration for the system with and without In descending order, the frames of the left column are the same as Figure 5.(b) The right column shows recordings at the same lattice sites, but in a simulation run with an increase in endfoot conductance.With increased endfeet conductance, the amount of K + moved through the glial cell is increased.Also, [K + ] ions are directed preferentially from one side of the cell to the other as shown by the third panels from the top.
The larger Kir current causes a greater local depolarization of the membrane which spreads along the cell.At the other end of the cell, there is an observable difference between the two conditions.The increased conductance allows more intracellular K + to be moved from the cell to the extracellular space.This is shown by the enhanced increase in [K + ] o and decrease in [K + ] i at lattice sites bordering the membrane (see third and fourth rows of Figure 6).Effectively, the larger endfeet conductance allows for a more directed and increased movement of K + away from the injection site.
3.2.Syncytium of multicell environment.To look at the impact of geometry on the evolution of the system, a cellular environment was artificially generated and is shown in Figure 7.When each cell is in its resting state, K + ions are injected at a rate of 64 mM/msec for 3.125 msec at an ECS node located approximately at the center of the system (Figure 7).At the site of injection, however, the [K + ] o increases only by a factor of about two, reflecting the balance between the injection process and the concurrent dilution by diffusion away from the source (Figure 8).The migration processes of K + and Cl − ions then are recorded and analyzed.It should be noted that during intense neuronal stimulation or during seizure, the ECS [K + ] o may reach a ceiling level of 12 mM; while in the case of injury-such as hypoxia, ischemia, trauma, or hypoglycemia-the ceiling level can be exceeded and the ECS K + concentration can reach levels of 25 mM [31].Accordingly, the amount of K + ions injected at the point sources in our simulations remain well within the normal ceiling level as well as the pathophysiological range.The number of injected particles in these simulations represents a balance between using physiologically accurate injection rates and being able to capture the qualitative component of the results.As we are primarily concerned with the implementation of a novel method to study SB mechanisms and the resulting qualitative effects, we use the above injection conditions.In this multicellular system, the SB mechanism is observed as follows.As K + ions are injected into the center of the system, the cations begin to diffuse away from the injection point toward the cell membranes.A narrow ECS surrounds the glial cells, giving little space for the diffusion to move the excess K + ions through the system.When the K + ions arrive at a cell, the membrane depolarizes.This depolarization spreads electronically around the cell, and K + ions are extruded from the ICS to the ECS at a location distant from the initial site of depolarization.In this way, K + ions are moved radially outward from the site of injection by the cells surrounding the injection point.In turn, there is an accumulation of K + ions at the distal sides of the cells that surround the injection point, which causes a local depolarization in the surrounding cells.This mechanism persists until the initial accumulation of K + has been dissipated.
The SB in the multicellular system is depicted in Figure 9, which shows recordings of [K + ] in the ICS (Figure 9a) and ECS (Figure 9b) and the membrane potential (Figure 9c) at two locations each for the three cells indicated in Figure 7.At nodes facing the injection site, [K + ] i increases as K + ions move into the cell, whereas at a node distal from the injection site, [K + ] i decreases as K + ions are extruded.The ECS [K + ] increases at each recording site as expected.The response is greatest in Cell 1 and decreases as we move away from the injection site to Cell 2 and then to Cell 3.This is because, as the K + front moves away from the center of the system, the accumulation is distributed over a larger area which degrades the response.
With the above injection protocol, the multicellular system was simulated with and without membrane flux.By removing the membrane currents, we can determine how effective diffusion alone is in dissipating the accumulated K + ions.This then can be compared with the dispersion of K + ions when the membrane currents (i.e., the spatial buffer) are included.Monitoring the cumulative flux through the four borders of the system gives a measure of how quickly K + ions can move from the center of the system to the periphery.As shown in Figure 10, at each border, the increase in flux is initially both more rapid and larger when the spatial buffer mechanism is included, as compared to simple diffusion.In the latter case, the injected particles must navigate a tortuous path towards the periphery, while in the former, some of the K + ions are moved through the glial cells as described above.As a result, a larger amount of K + ions are moved more quickly to the system's borders.In each case, the flux reaches a maximum by a subsequent decrease.The maximum flux occurs later in the system when only diffusion is included.During the decrease to the resting level where SD is not included, the flux through the border for the system becomes larger than when SB is included.This occurs because SB has sequestered the excess K + intracellularly and can return the ECS K + concentration more quickly to its resting value.These observations can be seen in all four cases, with the slight variations in magnitude due to the differences in the cell membrane geometry.Again, it is important to note that in these simulations, the discrepancy between membrane potential propagation velocity and ion diffusion  In the system with SB, there is a more rapid and greater increase in border flux as compared to the system without SB.This reflects a more efficient movement of excess K + ions to the system's periphery.
has been artificially decreased.The actual time taken to disperse the accumulated ions is expected to be longer.Observations analogous to those above can be seen when comparing the time course of K + ion migration through nodes that are equidistant from the injection source (Figures 11a and b).In the system with SB, the movement of K + ions occurs more quickly than without SB.The differences in the magnitudes of the K + ion number profiles between nodes equidistant from the injection source reflect the differences in geometry.For example, consider Figure 11a.The path from the injection source (node A in Figure 7) to node B is relatively straight except for one small cell that is blocking the direct path.Injected K + ions can move from node A unimpeded in the direction of node B. When some of the particles arrive at the cell blocking the direct path, this causes a depolarization of the cell which propagates electrotonically in both directions around the cell, causing K + to be extruded on the far side in the direction of node B. Other particles can navigate around the small cell directly to node B. In contrast, the most direct path connecting the injection source A with node C is hindered by parts of cells.While some K + ions will move by diffusion towards node C, the spatial buffer mechanism will not move K + directly in the direction of node C due to the orientation of the impeding cells relative to the injection point source.As a result the increase in K + at node C will be less than at node B, even though they are both equidistant from the point of injection.Also, note that in both Figures 11a and 11b, the maxima of the ECS K + number profile is greater in the system without SB as compared to that with SB.By having K + ions enter the glial cells, the cells themselves can hold onto some of the K + , thereby keeping the ECS K + lower.

Endfoot effect.
As in the case for the one-cell system, we examine the effect of an increased endfoot conductance on spatial buffering.Elliptical cells were arranged in a 100 x 180 regular array.The injection protocol used for this geometry was a constant injection of K + ions along the left boundary at a rate of 32 mM/sec for 3.125 msec.Simulations were run for the case with a uniform distribution of Kir conductance and with increased conductance at the ends of the cells.
Figure 12 shows the flux of K + through the right-most boundary for simulations with and without the increased endfoot conductances.It can be seen that the increased endfoot conductance enhances the spatial buffering mechanism.Simulations indicate that by having an increased conductance at opposing ends of the cells, the accumulated K + is dissipated more quickly.After diffusing toward the cellular array, a greater amount of K + is moved by the Kir channels into the cell.This results in a larger depolarization, which then spreads along the cellular membrane.At the other end of the cell, the depolarization causes K + to be moved out of the cell as described previously.However, with the increased conductance, more K + is moved.Effectively, the movement of K + is directed away from the injection site more efficiently.The larger conductances bias the dispersion of ions to occur linearly away from the injection site.

Conclusion.
In this study, a lattice cellular automaton approach was used to investigate the brain-cell microenvironment comprised of ECS and glial cells.This numerical method can readily deal with the complicated geometry and boundary conditions of neuronal tissue because of the microscopic interpretation of the LCA.To demonstrate the usefulness of the LCA in modelling neuronal processes, we studied the spatial buffering mechanism of glial cells.The model considered K + and Cl − dynamics.The ECS and ICS of the system were assumed to be uniform with constant diffusion coefficients for both ionic species.The interactions between the ECS and ICS were through the complicated membrane boundary conditions.K + ions were able to move through the membrane via Kir channels and via a passive flux that is coupled to Cl − .The Cl − also could move through the membrane by an ohmic pump.
Simulations of a single glial cell showed that our system does exhibit the spatial buffering mechanism.To study the buffering mechanism, K + ions were injected near the cell and the evolution of the ion numbers was traced.The ions move by diffusion toward the membrane, where they cause a local depolarization, causing K + to move into the intracellular lattice via the Kir current.The local depolarization would concurrently spread along the membrane away from the site of local injection.As a result, at locations away from the injection, there was an outwardly directed driving force on the K + ions, which would move K + ions out of the cell.In this way, the glial cell could disperse the accumulated K + from the site of injection.This process describes the spatial buffering mechanism.
To further reveal how readily the LCA can be adapted to suit more specific membrane boundary conditions, the Kir conductance along the region of membrane facing away from the injection site was increased.The purpose was to give a simple representation of the increased Kir conductance seen in the endfoot region of amphibian retina glial cells, where this increased conductance has been suggested as a mechanism to help siphon K + ions into the vitreous humor.By comparing simulations with and without the increased conductance, it was shown that in a single cell, the increased conductance effectively directs the K + away from the injection site toward the endfoot membrane.
An objective of this study was to demonstrate that this numerical method is able to treat complicated geometries associated with cellular environments.The representative two-dimensional geometry used in this preliminary investigation was not intended to mimic exactly a particular cellular geometry seen in slices of neuronal tissue with the accompanying membrane properties.The geometry was generated to contain approximately circular cells, resulting in a volume fraction of approximately 0.23, which is near the measured value of 0.20 [18].The specific aim of this system was to account for the geometrical parameters that influence the movement of particles in the brain-cell microenvironment, including the volume fraction and tortuosity.As such, implementing a more detailed representation of the individual glia morphologies was not considered, even though it can be implemented easily.A circular cellular morphology was used only to approximate the irregular and elongated shapes of glial cells.Again, the increased efficiency of the spatial buffer mechanism could be traced through the evolution of the system.We also were able to look more closely at the influence of geometry on K + ion buffering.The simulations indicate that while differing geometries do impact the movement of K + ions on the scale of a single cell, this effect does not appear to be significant on a larger scale.The results of the average models used in previous studies in the literature [6,15] appear to capture the essential elements of the spatial buffer mechanism.
We note that the time scales for diffusion in the ICS and ECS are approximately 4000 times longer than the time scale for propagation of membrane potential along the cell via electrotonic spread.In our simulations, the diffusion time scale was set to be approximately 100 times shorter than it actually is.Therefore, in reality the buffering mechanism should be slightly more effective than our results have shown.The reason for choosing this diffusion time scale was to increase the computational efficiency in running numerous simulations.Also, we only considered an isotropic tissue, although anisotropic tissue can be handled easily.
Our simplified model can be extended to incorporate other mechanisms that may play a role in glial buffering of K + ions.For example, incorporating Na + dynamics into the model would allow consideration of the Na + /K + pump on the glial buffering mechanism.Also, the LCA can account for the different types of K + channels located at different anatomical locations on the cell membranes.
A more detailed investigation as to the importance of geometry also can be carried out.In particular, glial cells are arranged in syncytia of glial cells connected by gap junctions.The electrical coupling of the cells may further improve the efficiency of the buffering by allowing for electrical activity to more effectively spread a considerable distance away from the site of K + ion accumulation.Without coupling, glial cells have both a short physical length and length constant which can prohibit the spatial buffer mechanism from working efficiently on larger space scales.The extensive coupling through gap junctions may enhance buffering, and this can be implemented in a straightforward fashion using the current approach.
What also would be interesting is to add a neuronal compartment to the model.In the brain-cell microenvironment, glial cells are interposed between the neurons.An accurate representation of brain-cell tissue would include both cell types.Such a model could be used to study K + ion buffering, but also could be applied to other problems such as mathematical treatments of spreading depression, a slowly propagating chemical and electrical wave in cortical neuronal tissue.
In summary, we have presented a versatile LCA numerical scheme that can be used effectively to model ion dynamics in realistic neuronal tissue.There are numerous potential applications of this numerical method in the study of physiological systems, e.g., a three-compartment model can be used to study spreading depression, which will be the subject of a future project.

Figure 2 .
Figure 2.Each boundary ECS lattice node may border one, two, or three ICS lattice nodes, labeled I, II, and III, respectively.Similarly, one boundary ICS lattice node may have one, two, or three ECS neighboring nodes (not shown).This results in a total of 28 membrane configurations that need to be taken into account when defining the boundary conditions for the LCA.In the above, • and • represent ECS and ICS nodes, respectively.

Figure 3 .
Figure 3. Schematic diagram of the single cell microenvironment with injection and recording sites.The ECS and ICS nodes at which the particle concentrations were recorded are depicted by the triangles and squares, respectively.The injection sites are shown by diamonds.The solid line indicates the region of membrane at which the Kir conductance was increased.The increased conductance was used in simulations only when specified.

Figure 5 .
Figure 5.Time courses of [K + ] concentrations at (a) extracellular and (b) intracellular nodes facing toward the injection site and at (c) extracellular and (d) intracellular nodes facing away from the injection site.The indicated nodes refer to the geometry shown in Figure 3.

Figure 6 .
Figure 6.Comparison of the time courses of [K + ] at various lattice sites with and without increases in endfoot conductances.(a) In descending order, the frames of the left column are the same as Figure 5.(b) The right column shows recordings at the same lattice sites, but in a simulation run with an increase in endfoot conductance.With increased endfeet conductance, the amount of K + moved through the glial cell is increased.Also, [K + ] ions are directed preferentially from one side of the cell to the other as shown by the third panels from the top.

Figure 7 .
Figure 7. Schematic diagram of the multicellular environment.The injection site-labeled as A-is indicated by a diamond; ECS recording sites-labeled as B, C, D, and E-are indicated by triangles.

Figure 10 .
Figure 10.Cumulative K + flux through the system's four borders with (solid line) and without (dashed line) spatial buffering.(a) Right.(b) Top.(c) Left.(d) Bottom.In the system with SB, there is a more rapid and greater increase in border flux as compared to the system without SB.This reflects a more efficient movement of excess K + ions to the system's periphery.

Figure 11 .
Figure 11.Time courses of [K + ] at nodes equidistant from the injection site.The labels correspond to the nodes indicated in Figure7.

Figure 12 .
Figure 12.Time courses of the changes in flux through the rightmost boundary of the system for simulations run with and without increased endfoot conductances.With increased endfoot conductance, K + ions are moved more quickly to the far boundary, away from the injection site.

Table 1 .
Typical parameter values used in the LCA simulations