INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUDING REMARKS APPENDIX A APPENDIX B REFERENCES

Ca²⁺-triggered secretion by neuroendocrine cells is known to be a relatively slow process and with longer latencies when compared with the secretion of neurotransmitters in synapses (Augustine et al., 1985; Llinás et al., 1981). In chromaffin cells, for instance, it is known that secretion continues during tens of milliseconds after a short pulse (Chow et al., 1992). In pancreatic -cells, latency between elevation of [Ca²⁺] and exocytosis has recently been reported (Eliasson et al., 1997).

A plausible explanation for such slow mechanisms of secretion may be found in the existence of cytosolic buffers that delay the response by slowing down the free calcium transient. The way in which Ca²⁺ exogenous chelators interfere with secretion, in both chromaffin (Chow et al., 1996) and -cells (Bokvist et al., 1995; Pertusa et al., submitted for publication), strongly supports such a possibility.

Speaking in general terms, one can summarize the basic ingredients of calcium signaling (Clapham, 1995) in neuroendocrine cells as consisting of transient of calcium ions through the cell membrane, uptake and release of Ca²⁺ by internal stores, diffusion of calcium, and binding/unbinding by endogenous (or exogenously added) buffers. When secretion is studied, the dynamics of the release granules responsible for secretion as well as their spatial distribution also have to be taken into account.

Several studies have previously considered the problem of buffered calcium diffusion in excitable cells, modeling the system by means of diffusion-reaction differential equations, solved numerically using finite differences schemes (Sala and Hernández-Cruz, 1990; Nowycky and Pinter, 1993; Klingauf and Neher, 1997). Other approaches, to mathematically simplify the models, considered different kinds of approximations valid either for rapid buffers (Wagner and Keizer, 1994; Smith et al., 1996) or in the linear regime (Naraghi and Neher, 1997). Such previous studies addressed the importance of the study of Ca²⁺-buffered diffusion both to gain insight into the secretory response (Klingauf and Neher, 1997) as well as to correctly interpret data from fluorescence experiments (Nowycky and Pinter, 1993; Smith et al., 1996) that use exogenous buffers as indicators of Ca²⁺ intracellular activities.

In this study we have adopted a Monte Carlo method to model the buffered diffusion of calcium. Monte Carlo simulations have been previously used successfully to study reaction and diffusion processes in biological systems (Saxton, 1994, 1996; Riley et al., 1995; Kruk et al., 1997). Our motivations for adopting such a scheme are several. First, no a priori assumptions about the symmetry of the problem are needed to reduce the number of spatial dimensions; in this way, we can perform a three-dimensional simulation of diffusion. Second, the boundary conditions of the problem are more easily taken into account, and one can vary such conditions with ease. For instance, with a Monte Carlo simulation it is a simple task to model the entrance of ions through channel pores (not necessarily regularly distributed) and to consider the possibility of clustered channels. Third, Monte Carlo techniques are robust methods for solving complex problems in which different agents come into play together, interacting among themselves. The robustness and versatility of the method enable us to safely study tendencies when the parameters affected by large uncertainties are varied, which is indeed the case of the diffusional and kinetic parameters of our system. Last, but not least: as we shall see, given the relatively small number of Ca²⁺ ions and buffer molecules involved, a stochastic approach to the problem may be more appropriate and safer than the solution of the reaction-diffusion differential equations (in terms of concentrations).

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUDING REMARKS APPENDIX A APPENDIX B REFERENCES

We have built a Monte Carlo simulation to treat the problem of 3-D diffusion and kinetic reactions of calcium ions and buffers. The main ingredients of the simulation are

1.	Selection of the geometrical configuration and distribution of channel pores
2.	Selection of initial concentrations
3.	Diffusion of Ca2+ and mobile buffers
4.	Binding and unbinding of Ca2+ ions by buffers
5.	Entry of Ca2+ through channel pores

One can also consider extrusion mechanisms through the membrane (Sala and Hernández-Cruz, 1990; Klingauf and Neher, 1997); however, the contribution to Ca²⁺ concentration profiles for time intervals smaller than 100 ms seems to be irrelevant (Sala and Hernández-Cruz, 1990). Other phenomena not considered in the present simulation are the effect of obstacles (release granules, mitochondria) in diffusion (Saxton, 1994; Ölveczky and Verckman, 1998), sequestration, and release by internal stores (like mitochondriasee Xu et al., 1997; or granules, sarcoplasmic reticulum, etc.). In the present study we limited our field of investigation to the problem of buffered calcium diffusion in a subdomain with uniformly (but not necessarily regularly) distributed calcium channels (or channel clusters), with a Ca²⁺ influx discretized in space (the entry is restricted to the pores).

Let us now discuss the different ingredients of the simulation of Ca²⁺-buffered diffusion and how the sequence of the simulation works.

Geometrical configuration

The only a priori geometrical simplification considered in our study, to save computational time, is the restriction to a conical subdomain of the whole cell, assuming zero flux on the lateral boundaries of the cone. This is a reasonable approximation if one assumes that the conical subdomain (which, typically, will have from 15 to 45 channel pores) is embedded in a larger submembrane region where the channels are more or less uniformly distributed. In such a case, one can always find, approximately, a conical boundary with zero flux. When a spherical cell with uniformly distributed channels is considered, the conical section should reach the center of the cell.

The selection of such a geometrical configuration considerably relaxes the assumption of symmetry around a channel pore assumed by Klingauf and Neher (1997) and thus allows us to consider randomly but uniformly distributed channels or channel clusters, not necessarily regularly distributed.

We take a conical section with R = 5 µm of depth and a base with radius r = 1 µm (Fig. 1), and we assume that the net flux of ions and mobile buffer molecules is zero on the lateral side of the cone. The conical section would reach the center of a typical pancreatic -cell. We approximate the spherical cell membrane by a plane surface (upper side of the cone), discarding negligible curvature effects (its surface increases by only 1% when curvature is considered).

View larger version (57K): [in this window] [in a new window]   FIGURE 1   Conical domain considered for the simulation. The dots represent the calcium channels.

A 3-D and orthogonal grid mapping the conical section is considered, with length of the sides l = 0.105 µm (see Fig. 1). The total number of cubic compartments is 4450, 253 of which are located at the upper slice of the cone (submembrane domain); the number decreases until the greatest depth is reached, at the vertex of the cone, where there is only one compartment.

On the upper side of the cone we distribute the channel pores uniformly and randomly. Because of the assumption of symmetry needed to balance the number of ions on the fictitious walls of the conical subdomain, the distribution of channels (or channel clusters) should preferably not break the (approximate) uniformity in the surface distribution.

The typical number of calcium channel pores ranges from 15 to 45, corresponding to a density of channels ranging from 5 to 15 pores per µm², which seems to be in the range of physiological values (Klingauf and Neher, 1997).

Notice that when a cell of radius R = 5 µm with uniformly distributed calcium channels is considered, the relation between the whole-cell current (I_wc) and the incoming current for our conical domain (I_d) is given by the ratio of the surface of the submembrane domain over the surface of the whole cell, that is, I_d/I_wc = r²/(2R)² = 10². In our simulations, we always refer to the domain current (I_d); the whole-cell current would be a factor 100 larger in the situation described.

Diffusion

The 3-D diffusion of ions and mobile buffers is modeled as a random walk process. During every time step t and for each dimension, a particle has a 50% chance of remaining in its initial position and a 25% chance of moving in either the positive or the negative direction.

The time step corresponding to such probabilities is given by the relation (Kruk et al., 1997)

(1)

where D = 220 µm²/s the diffusion coefficient for Ca²⁺ in cytoplasm (Allbritton et al., 1992). To model 3-D diffusion we have to take into account the probability of staying in the same cubic compartment of the grid and the probabilities of moving to the surrounding 26 cubic compartments (Fig. 2). After each diffusional time interval t we decide the fate of each ion with the probabilities of staying or moving in each compartment of the grid as described in Fig. 2.

View larger version (53K): [in this window] [in a new window]   FIGURE 2   Probabilities for a particle in the center of the cube (corresponding to a given compartment of the grid) of staying in this same compartment or moving to a neighbor compartment in a 3-D random walk step. Notice that each vertex, middle point of an edge, center of a side, or center of the cube in the diagram corresponds to a given cubic compartment of the grid.

Depending on the mobile buffer considered (EGTA, Fura-2, etc.), the corresponding diffusion coefficient can be typically about two to five times smaller than that for Ca²⁺, so that the characteristic diffusion times are about two to five times bigger. Thus the diffusion of buffer particles can be simulated by moving them in exactly the same way that Ca²⁺ ions are moved, but only after two to five simulation steps (Kruk et al., 1997). By simply tuning the probabilities of staying or moving for the buffer particles we could also account for noninteger ratios of diffusion coefficients; however, given the considerable uncertainties in the actual values of diffusion coefficients in cytoplasm (Klingauf and Neher, 1997), considering integer ratios seems to be enough.

Typically, the total number of buffer molecules per compartment in our simulations is as "large" as 300. Thus to save time it is convenient not to select the next position of each of them, but to redistribute the bulk of buffer particles according to the probabilities of going to the different 27 positions (see Fig. 2). When we have N > 64 buffer particles to diffuse, we take N = 64c + R, with R < 64 (c and R are integers). Then we consider the "bulk" diffusion of the 64c buffer particles according to the following distribution:

c	buffer particles for each corner of the cube
2c	buffer particles for each middle point of an edge
4c	buffer particles for each center of a side
8c	buffer particles in the center of the cube,

and we decide probabilistically the position of the remaining R particles according to the probabilities in Fig. 2.

For the free calcium ions, given the number of incoming ions and the moderate expected rise of free calcium (for instance, for 10 µM there are seven ions per compartment in average), the bulk redistribution is absent on most occasions and only the probabilistic selection comes into play. This fact suggests that, for the size of the grid considered, the calcium diffusion is better described probabilistically.

Of course, we must also consider the reflection in the boundaries of the subdomain. This effect is taken into account by rejecting those events that result in the positioning of a particle outside the subdomain; in other words, when a particle encounters a wall, we do not move it until the next diffusion time step (in case this further step takes the ion inside the conical subdomain).

The reflection takes place both in the physical upper boundary (cell membrane) as well as in the fictitious walls (lateral side of the cone); reflection in this last case is equivalent to considering that the number of outgoing particles equals that of incoming ones.

Incoming Ca²⁺ current

As previously commented, the value of the whole-cell current I_wc is on the order of 100 times the domain current, considering a spherical cell of radius R = 5 µm with uniformly distributed calcium channels.

Typically, for a current of I_d = 10 pA the number of ions entering the subdomain would be ~400 ions per unit of diffusion time (12.5 µs). Thus, for a density of channels _cha = 15 channels/µm², corresponding to 48 channels for the conical domain, the unitary current will be ~200 fA, corresponding to an approximate number of two ions per channel every 3 µs.

We adopt a simple scheme for simulating the entry of ions. We consider a pulse with a given initial current, typically I_d = 2-10 pA in the conical domain, and decreasing exponentially. We take

(2)

and I = 0 if t > T. We will also consider constant pulses.

Among other selections, in our simulation we use the values I_d⁰ = 10 pA,  = 20 ms, T = 50 ms, which corresponds to a time-averaged current I_d = 3.7 pA and a unitary current of 76 fA, within the range of physiological values (Klingauf and Neher, 1997).

At each simulation time step, the number of incoming ions is evaluated using the exponential (or any another selection), and then the entries are randomly selected from among all available channels. In this way, the number of ions entering through each pore is obtained and the number of free calcium ions in the corresponding compartment beneath the pore is incremented.

This description of the incoming calcium current implies a spatial discretization of the entry points, which, as we have already pointed out, considerably relaxes the assumptions on the symmetry of the distribution of calcium channels. However, it is well known that calcium channels open and close stochastically, so that, for a more realistic description of calcium currents, a discretization in time, considering open/close probabilities for the channels, would be desirable. The time dependence of the calcium current is determined by the variation of the open/close/inactivation probabilities for the calcium channels. Nevertheless, as commented by Klingauf and Neher (1997), channel gating can be expected not to be of importance for the understanding of the (relatively slow) secretory response of neuroendocrine cells. On the other hand, it is likely that fluctuating channels could be important in the study of neurotransmitter release or, in general, in the investigation of rapid responses to calcium signals.

Our Monte Carlo simulation naturally allows a stochastic time description of the calcium current, and as an illustration we will present an example showing the effect of channel gating (see section Channel gating: a simple model).

The opening and closing of calcium channels is modeled as a Markovian process: given the mean open (_o) and closed (_c) times for the channels, we first select probabilistically the initial state of the calcium channels with probabilities P_o = _o/(_o + _c) of being initially open and P_c = 1  P_o of being initially closed. Then, for each closed channel we choose the time it will remain closed, using an exponential distribution with expected value _c (we proceed similarly with the open channels). After the time for a channel to remain closed has expired, the channel opens and we decide the time it will remain open (and similarly for the transitions open  closed).

We are not considering inactivated states; nor are we taking into account the possible variation of the open/closed (and inactivation) times. The study of the dynamics of calcium channels and its dependence on voltage, inactivation by [Ca²⁺], and so on lies beyond the scope of this work and deserves a separate analysis.

Kinetics

During each diffusional time step t and in each compartment of the grid the kinetic equations that describe the buffering processes are applied and the calcium and buffer concentrations evolve. In terms of concentrations, one can describe the first-order kinetics of the process
by the system of equations

(3)

where k₊ⁱ is the forward binding rate for the buffer i and kⁱ is the unbinding rate.

We take the values of the kinetic constants from Klingauf and Neher (1997), except for the endogenous buffer. We only consider two buffering systems (as in Sala and Hernández-Cruz, 1990; Nowycky and Pinter, 1993): an endogenous fixed buffer (500 µM), with the kinetic constants given by Klingauf and Neher (1997), plus an exogenous (and mobile) buffer (Fura-2 or EGTA). Of course, we could have included more buffering systems, but we preferred to keep the model within the minimum number of parameters, given the considerable uncertainties in the kinetic constants for the buffers (in particular for the endogenous buffer(s)).

Because the simulation is performed in terms of the number of ions and buffer molecules, we adopt the first-order kinetic equations in terms of the number of particles of each species in each small cubic compartment of the grid. In terms of the number of free calcium ions N_Ca, free buffer molecules N_Bi, and bound buffer molecules N_CaBi, Eq. 3 reads

(4)

where N_Bi^T represents the total (free + bound) number of buffer particles of type i, N_Ca^T is the total number of calcium ions, and N_Ca represents the variation of free calcium ions in a given cubic compartment during a time interval dt  t/n. ₊ⁱ, ⁱ are defined below.

To integrate the kinetic equations we subdivide the interval t into n subintervals of length dt, evaluating probabilistically the changes in the free calcium and the free and bound buffer concentrations for each dt. During the integration of the kinetic equations inside each time interval t, the total (free + bound) number of buffer particles of type i and the total number of calcium ions N_Ca^T remain constant. These numbers can only change because of diffusion, which is taking place for each t (> dt).

The effective constants ⁱ and ₊ⁱ for binding and unbinding in each time interval dt are given by

(5)

where t is given by Eq. 1, n is the number of subdivisions of the interval t used to integrate the kinetic equations, and V is the volume of each cubic compartment.

Starting from the initial conditions for N_Ca, N_Bi, and N_CaBi and choosing n large enough (dt small enough), one can obtain the final N_Ca after each interval dt. Not only the number of calcium ions is required; the changes in the number of free and bound buffer molecules must also be evaluated in each iteration of Eq. 4 and actualized for the next step dt. Because the effects of binding and unbinding for the two buffer species are coupled, the solution to the safe integration of Eq. 4 seems to consist of taking a very small dt and splitting Eq. 4 into different equations for each of the effects (binding-unbinding) and the buffers. Indeed, for small enough dt the small variations due to each of the effects (binding/unbinding) and buffers should have a small impact on the rest. However, given the small number of ions and buffer molecules in each compartment (in particular, the almost nonexistent Ca²⁺ before influx), the variations in each dt could be too small to record any single change when N_Ca is considered to be an integer number and finite differences are used to integrate Eq. 4. In fact, using finite differences implies that N_Ca (and N_Ca, N_Bi, etc.) is a real and not necessarily integer number. In other words, N_Ca, N_Ca, etc., should be interpreted as average values (or expected values); in this way, a macroscopic description in terms of concentrations (Eq. 3) would be equivalent to Eq. 4. However, given the small number of ions and buffer molecules per compartment of the grid, the more straightforward interpretation of N_Ca, N_Bi, etc., as actual numbers of ions and buffer molecules in a given compartment of the grid is a feasible one.

To implement such a microscopic interpretation, one should describe the system probabilistically, evaluating the probabilities associated with each of the processes and generating events according to such probabilities.

Taking a large enough n one can then split Eq. 4 as follows.

Binding: For each buffer species i one can write the variation in Ca²⁺ ions (ⁱN_Ca) due to binding with such buffer species i as

(6)

which in terms of the binding probability means that each Ca²⁺ ion has a probability

(7)

of being bound by any of the N_Bi free buffer molecules during the time dt. Such binding probabilities can change after each step dt because the number of free buffer molecules N_Bi can vary.

Unbinding: On the other hand, bound buffer molecules can suffer unbinding and increase the number of calcium ions by the amount ⁱN_CaBi; but this variation equals N_CaBi, so that

(8)

in other words, each bound buffer particle of type i has a probability

(9)

of becoming unbound in the interval dt. Notice that the unbinding probability is constant in time, in contrast to the binding probability (arising from the nonlinear terms of the system of equations).

Restrictions: During each diffusional time step, the total number of calcium ions (bound + free) or buffer mole- cules (bound + free) of each kind is held constant in each compartment.

After evaluating the probabilities of binding or unbinding in a time dt, one can decide the number of ions that undergo any of the processes described. Given a probability p for each ion to be bound (or for a buffer molecule to become unbound), the probability for the binding of k ions of a total of N ions (or unbinding of k buffer molecules of a total of N bound buffer molecules) in a time interval dt is given by the binomial distribution

(10)

Then we can select the number of ions and buffer molecules that become bound or unbound in each time interval dt according to the binomial distribution (for algorithms to generate random deviates that follow binomial distributions, see, for instance, Press et al., 1992, pp. 285-286).

To safely integrate the kinetic equations, a large enough n is selected for every diffusion time t. Such a value n is evaluated every time and for each compartment before the kinetic equations are integrated. Therefore, the effective constants ⁱ and ₊ⁱ can be different at different time intervals of the simulation, though keeping the same ratio. The selection of the appropriate number of subdivisions of the interval is made by requiring that all binding (and unbinding) probabilities are smaller enough than 1; that is, we set n by bisectioning the interval as many times as required to fulfill the inequalities:

(11)

Let us stress that this subdivision is performed locally, both in time (at each diffusion time step) and in space (at each compartment of the grid).

Numerical experiments suggest that P_max = 0.1 is a safe selection that guarantees proper convergence of the method. Indeed, we do not find any changes in our results by limiting the probabilities by smaller values of P_max.

In summary, we have considered that, in each time interval t, each bound buffer molecule has a constant probability of becoming unbound while each free calcium ion has a probability of becoming bound, by a given buffer species, proportional to the number of free buffer molecules in a region surrounding the calcium ions (in our simulation, a compartment of the grid). The differential equation (Eqs. 3 and 4) can thus be interpreted as giving statistical averages for the variation of the number calcium ions. In the Monte Carlo simulation we stay at the probabilistic level. The interval t is divided into n subintervals dt = t/n in such a way that all probabilities are small enough. The number of calcium ions that become bound in each interval dt or become unbound is evaluated using the corresponding probabilities and following the distribution in Eq. 10.

The algorithm gives priority to the fastest effects during each t (a waiting time can be considered as done when the buffers are diffused). For instance, unbinding is normally so slow that the probabilistic selection of the number of buffer molecules that undergo unbinding can be made at the end of each interval t and not necessarily each dt  t. In case two or more effects have the same priority in a given time dt, the order in which they act is generated considering equal probabilities for those effects.

Initial conditions

Initially, we consider the number of Ca²⁺ ions and buffer molecules needed to achieve equilibrium for each of the buffer species; in each compartment the condition of equilibrium is given by

(12)

where N_Ca should be understood as the average number of ions in each compartment of the grid. Assuming that the initial concentration of free Ca²⁺ is the basal concentration, [Ca²⁺] = 0.1 µM, which is equivalent to only 310 ions in the whole conical subdomain and 0.07 ions per compartment, Eq. 12 can only be interpreted as giving the average number of free and bound buffer molecules in each compartment. As the expected average number of Ca²⁺ ions in each compartment is by far smaller than 1, locally (in each compartment) there is no possible equilibrium; however, in the whole conical subdomain, equilibrium exists. Of course, because we are following the fate of each ion, there can be one or more ions in each compartment or none at all, but never fractional values.

Regarding the buffers: the typical concentration considered for the fixed endogenous buffer will be 0.5 mM, which, in terms of the number of molecules, means ~350 buffer molecules (free + bound) in each compartment.

These estimations clearly suggest that, at the beginning of a Ca²⁺ pulse, starting from the basal concentration for free Ca²⁺, there is no way local equilibrium of Ca²⁺ with the different buffer species can be reached in each compartment of the grid. Equilibrium has to be interpreted as an average over the whole conical domain, regardless of how fast the binding/unbinding of each buffer species takes place. Only when the number of Ca²⁺ ions is large enough could local equilibrium possibly be found. Therefore, the assumptions of the rapid buffer approximation (RBA) are not expected to be met on the submicron scale at short times after the beginning of a pulse; however, RBA succeeds in estimating local average concentrations for fast in-mobile buffers, as we will later discuss.

As input data, we consider the basal concentration of calcium and the concentration of the different buffer species as given in Table 1. We distribute the corresponding number of Ca²⁺ ions and buffer molecules (free and bound) randomly and uniformly in the conical subdomain. This is the starting point of our Monte Carlo simulation.

Sequence of the simulation

The sequence of the simulation is carried out as follows.

Presimulation: The number of ions and buffer molecules given by the equilibrium conditions are distributed randomly and uniformly along the conical domain. To ensure that global equilibrium is reached, the distributions of ions and buffers develop during 1 ms according to the kinetic equations; during this presimulation lapse all channels remain closed. The system follows its course by means of the same procedures that apply to the actual simulation (except the entrance of ions).

Simulation: For each simulation step, the following actions are taken and repeated in this same order until the end of the simulation.

	Entry of ions: The incoming calcium flux is modeled as previously discussed. On most occasions, except when explicitly stated otherwise, the time dependence of the Ca2+ is varied continuously according to an exponential decay.
	Kinetics: The binding and unbinding are calculated in each compartment, following the previously discussed scheme.
	Diffusion of Ca2+: Ca2+ ions and mobile buffers are diffused by using the previously discussed random walk scheme. As we mentioned before, diffusion of buffers is delayed with respect to ions, that is, we diffuse the mobile buffers (assuming the ratio of diffusion coefficients is k) after k diffusion cycles have been completed.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUDING REMARKS APPENDIX A APPENDIX B REFERENCES

In the following, we discuss the effect of the different parameters involved in the simulation. Unfortunately, there is considerable uncertainty as to the actual values of the physical constants describing the endogenous buffers (Klingauf and Neher, 1997; Xu et al., 1997; Naraghi and Neher, 1997). However, our simulation, being valid quite independently of the possible values, can be used safely to study the tendencies when the parameters are varied.

To illustrate the kind of results that can be obtained from our Monte Carlo method, we adopt as standard values those given by Klingauf and Neher (1997) for chromaffin cells. These values, together with the rest of the parameters of the simulation, are summarized in Table 1. We will vary these values to study the effect of the parameters defining the model.

We will show the time courses for free calcium, endogenous (fixed) buffer, and exogenous (mobile) buffer corresponding to short pulses with a duration of 50 ms and different shapes and total currents. Although the simulation is carried out in terms of the numbers of calcium ions and buffer molecules, the results are presented in terms of more easily readable concentrations. However, the fact that we are dealing with a discrete stochastic simulation that moves particles becomes apparent in the dispersion observed when we plot the data generated by the simulation (every 12.5 µs). Each such point corresponds to the concentration at a given time step of the simulation; the number of each type of particle in any given compartment of the grid suffers stochastic fluctuations but gives overall smooth variations when average values are considered. The Monte Carlo simulation not only gives average values for concentrations but also shows the uncertainty of the simulations. The time fluctuations can be interpreted as dispersions in the concentrations when results averaged over time intervals are considered, possibly limited by the temporal resolution of a given experiment (which is expected to be larger than the time step of the simulation).

To display our results, we consider two kinds of presentation. The most immediate one consists of directly showing the output of the simulation, with the appearance of scatter plots; the second presentation consists of averaging the time courses (with binning of 10 points) and giving the average value and, eventually, error bars (2   confidence level). We choose the first kind of presentation when considering local studies and the second one for showing data averaged over given regions.

Dependence on the intensity of the Ca²⁺ current: saturation levels

Fig. 3 A shows the free Ca²⁺ concentrations for the values listed in Table 1 with varying initial domain currents I_d⁰ and considering Fura-2 as an exogenous buffer. We consider a current given by Eq. 2, with  = 20 ms, T = 50 ms, and I_d⁰ = 10, 5, 2.5 pA. Fig. 3, B and C, shows, respectively, the concentration of free endogenous buffer and exogenous buffer (Fura-2). All of the figures show the averages over the submembrane domain; that is, the results are averages down to 0.105 µm  100 nm of depth.

View larger version (25K): [in this window] [in a new window]   FIGURE 3   Effects of the amplitude of the incoming current on the average concentrations of free calcium and free buffer molecules at depths between 0 and 100 nm. Average values with binning of 10 points, corresponding to time intervals of 10 × 12.5 µs, are shown. Error bars (2- confidence level) for some selected points are displayed for the [Ca2+] time courses; the uncertainties for the buffer time courses are small (lower than 1%), and they are not shown. (A) [Ca2+] time course (with 100 µM Fura-2 added) and for three different exponential currents (see Table 1): Id0 = 10, 5, 2.5 pA. (B) Time course for the free endogenous buffer for the same conditions of A. (C) As in B, for the time course of Fura-2. (D) As in A, without Fura-2.

The selection of an exponentially decaying pulse for the simulation is only one possibility among many others. It has the theoretical benefit that, as our results will show, the Ca²⁺ concentration near the membrane peaks before the pulse ends. This fact is interesting when different time courses for different kinetic parameters are compared. Constant intensities (as in Sala and Hernández-Cruz, 1990; Nowycky and Pinter, 1993; Klingauf and Neher, 1997) will also be considered.

In Fig. 3 A one observes that the free calcium concentration is far from responding linearly when currents in the range of 5-10 pA are considered; in fact, a reduction by a factor of 2 in the incoming current involves a reduction by a factor of 4 in the increase of free Ca²⁺. This fact can be understood by looking at Fig. 3 C: for 10 pA the exogenous buffer saturates rapidly and then the increase in [Ca²⁺] becomes faster. Under such saturation conditions, the [Ca²⁺] can be expected to depart from the linearized buffered approximation (Naraghi and Neher, 1997), which assumes low variations of buffer concentrations. For lower currents the response of Ca²⁺ gradually becomes more linear, as can be seen by comparing the rise of [Ca²⁺] for incoming currents of 5 pA and 2.5 pA; observe also that Fura-2 does not saturate for such currents.

The effect of the exogenous buffer saturation is also shown in Fig. 3 D, where the [Ca²⁺] time course (without Fura-2) is plotted. As can be seen in the figure, the nonlinear effects shown in Fig. 3 A have now become attenuated.

Notice also the steep decrease in free Fura-2 concentration (Fig. 3 C), in contrast to the smoother relative variation of endogenous buffer concentration (Fig. 3 B). The endogenous buffering system rapidly becomes accommodated to the concentration of free Ca²⁺, as can be inferred by noticing that [Ca²⁺] and [B_endo^free] follow a similar pattern of temporal variation (compare Fig. 3, A and B); [B_endo^bound] also demonstrates this pattern (not shown), reaching its maximum when [Ca²⁺] does. In contrast, the exogenous buffer is not able to keep track of the rapid variation in [Ca²⁺] (Fig. 3 C).

Shape of the Ca²⁺ current

Let us now briefly discuss the importance of the shape of the calcium pulses. Until now, we have chosen an exponentially decaying current as described in Eq. 2, with  = 20 ms and T = 50 ms, which for I_d⁰ = 10 pA gives a time-averaged current I_d = 3.7 pA. Let us compare such a result with the corresponding one for a constant pulse with a same duration T = 50 ms and a same total influx charge; that is, we take I = 3.7 pA for T < 50 ms, and we switch off the pulse at T = 50 ms. The comparison between the result of these two different pulses is shown in Fig. 4, A and B (considering Fura-2), and Fig. 4 C (for EGTA).

View larger version (14K): [in this window] [in a new window]   FIGURE 4   Effects of the shape of the current. Average values with binning of 10 points are shown with error bars (2- confidence level) for some selected points. (A) [Ca2+] time course averaged down to 100 nm of depth, considering an exponentially decaying calcium current and a constant current (Fura-2 added). (B) As in A, but for depths between 300 and 400 nm. (C) As in A, but considering EGTA (parameters given in Table 1) instead of Fura-2.

Fig. 4 A shows the average calcium concentrations, varying in time, for the submembrane domain (0-100 nm). After the pulse has been disconnected, one can observe that the same equilibrium is reached for both the exponential and the constant pulses, which clearly indicates that the same amount of Ca²⁺ ions has entered the cell, and thus the same whole-cell equilibrium has to be reached. However, for the exponential pulse the peak is reached at ~20 ms, while the concentration for the constant pulse (slightly larger in the peak) increases steadily until the end of the pulse (50 ms). The importance of the shape of the macroscopic calcium current becomes manifestly important in determining when the peak is reached and, to a lesser extent, in determining how high the peak is. Of course, as deeper regions are considered (for instance, 300-400 nm, Fig. 4 B), the peak for the exponential current is delayed and approaches the peak for the constant current. In any case, for determining concentrations near the membrane a detailed simulation of the calcium current seems to be necessary; not only it is important to know the magnitude of the current, but a good description of the shape of the current is also needed. This will be the subject of further study.

Ca²⁺ concentrations as a function of depth

In Fig. 5 A, we can observe the Ca²⁺ concentrations at different depths. Notice how, as expected, the increase in [Ca²⁺] becomes smaller as the depth grows. In addition, the [Ca²⁺] peak becomes more delayed as deeper regions are considered. As we go deeper into the cell, the peak of [Ca²⁺] begins to disappear. In such deep regions, Ca²⁺ ions arrive with considerable delay, and the Ca²⁺ concentration increases steadily to reach the values given by the equilibrium in the whole subdomain.

View larger version (21K): [in this window] [in a new window]   FIGURE 5   Time courses (100 µM Fura-2 added) at different depths (0-100 nm, 100-200 nm, 200-300 nm, 500-600 nm, 1000-1100 nm) for an exponential calcium current lasting 50 ms. A time interval of 1 s is simulated. Averages with binning of 10 points are shown in the bottom figures, while in the top figures the direct output generated by the simulation is displayed. (A) (Bottom) [Ca2+] time courses at different depths. The [Ca2+] averaged over the whole conical domain is also shown. As the depth becomes larger the rise in [Ca2+] is smaller and the calcium peak is delayed. (Top) detail for long times at 0-100 nm and 1000-1100 nm. (B) As in A, for the free endogenous buffer. (C) As in A and B, for free Fura-2.

In Fig. 5 the calcium and free buffer time course for an exponential pulse (I_d⁰ = 10 pA,  = 20 ms, T = 50 ms) and a total simulated time of 1 s are shown. As can be observed, the concentrations at different depths tend to the same value as the time elapsed becomes longer. Notice, however, that for such "long times" the effect of extrusion of Ca²⁺ should be noticeable. We have not implemented this effect in our simulation because we are mainly concerned with the study of calcium time courses in time intervals on the order of 50-100 ms. The results shown for times of ~1 s, without extrusion, will be used in Appendix A to show that our simulation converges properly at long times. We will show that the final concentrations are in full agreement with the solution of the reaction diffusion problem at long times.

The average values for different depths (Figs. 4 and 5) display the behavior described by Sala and Hernández-Cruz (1990) and Nowycky and Pinter (1993). However, such values are not sufficient to explain the secretory response of cells to calcium influx: as discussed by Klingauf and Neher (1997), the study of submembrane gradients is essential. For such a purpose, the versatility of the Monte Carlo simulation is very useful. We will consider such gradients later. First, let us study the validity of the RBA, focusing on the submembrane domain.

Study of local equilibrium and the RBA hypothesis

Let us consider the following four different kinds of submembrane compartments:

1.	Type A: Compartments beneath a channel pore with an extra neighbor pore
2.	Type B: Compartments beneath a pore with no extra neighbor pores
3.	Type C: Compartments not corresponding to a pore but with a neighbor compartment having a pore
4.	Type D: Compartments that have no pore and no neighboring pore

Such characterization of the submembrane compartments enables us to study the local equilibria of the buffering systems with the calcium ions. The percentage deviation from equilibrium of a given species of buffer (i) at given time (and spatial location) reads

(13)

where [B_i⁰] is the concentration of free buffer i that would be in equilibrium with the calcium concentration [Ca]. By considering the equilibrium equation k₊[Ca][B_i]  k([B_i^T]  [B_i]) = 0, Eq. 13 can be rewritten in terms of the concentration of free buffer, total buffer, and free calcium. To study the local equilibrium in the different types of compartments (k = A, B, C, D), we use Eq. 13 in terms of the number of free calcium ions and buffer molecules:

(14)

where NBik, NBiTk, and NCak represent, respectively, the number of free endogenous buffer molecules, total (free + bound) endogenous buffer molecules, and free Ca2+ ions, summed over all of the compartments of the kind k (k = A, B, C, D). Ncomk is the number of compartments of the type k. We apply Eq. 14 to study local equilibrium for both the endogenous buffer and the exogenous buffers; a negative deviation will mean that there is an excess of free buffer molecules.

Fig. 6, A and B, shows that, regardless of the type of compartment considered, the endogenous buffer reaches equilibrium rapidly and "locally" (but averaging over homologous compartments); this fact suggests that RBA is valid for such a rapid buffer, at least in the average over compartments of the same kind. Notice, however, the dispersion in the plots, which becomes larger as the number of compartments that enter into the average is smaller (A- and B-type compartments). Such stochastic fluctuations will have a larger impact on [Ca²⁺] concentrations as lower influx currents are considered. Let us stress again that for typical free calcium concentrations of 10 µM an average value of around seven free calcium ions is to be expected in each compartment of the grid; thus, for smaller concentrations there are not enough free calcium ions to meet a stable local equilibrium at scales on the order of 100 nm or smaller. The hypothesis of rapid and local equilibrium for the endogenous buffer is thus seen to be valid as an average over homologous compartments, with typical deviations larger as fewer compartments enter into the average and as lower calcium concentrations and incoming currents are considered. In Appendix B we will discuss how the finding of average local equilibrium of the endogenous buffer is explained by estimating the time the endogenous buffer requires to reach equilibrium, starting from a situation away from equilibrium.

View larger version (5521K): [in this window] [in a new window]   FIGURE 6   Tests of equilibrium for the endogenous buffer as well as for Fura-2 (100 µM). For every 10 steps of simulation the output given by the simulation is plotted. (A) Deviation from equilibrium (in %) for the endogenous buffer for I0 = 500 pA and Fura-2 (100 µM) added; the equilibrium is tested both for A-type and B-type compartments (see text). (B) As in A, for compartments of type C and D. (C) Percentage deviation time course of the exogenous buffer (Fura-2) with the conditions of A and B. The deviation is tested for compartments of type A, B, C, and D; the largest deviations occur for the compartments nearest the channels (A and B compartments, which are almost indistinguishable in the figure); the deviation for C-type is larger than for D-type. (D) Percentage deviation from equilibrium of the endogenous buffer as an average at depths between 0 and 100 nm. Two incoming currents are considered: Id0 = 10 pA (largest deviation at small ts) and Id0 = 5 pA. (E) Deviation from equilibrium for Fura-2 in the first 15 ms of the pulse for compartments of type B and two choices of the endogenous buffer kinetic parameters: k+ = 5 × 108 M1 s1, KD = 10 µM, and k+ = 1 × 108 M1 s1, KD = 50 µM. The largest deviation occurs as k+ becomes larger (k is the same in both cases). (F) As in E, for compartments of type C. (G) Percentage deviation time course of EGTA for compartments of type A, B, C, and D. For A and B compartments the deviation is the largest and the plots overlap; the smallest deviation is for D compartments.

Although equilibrium is reached locally, at least on the average, in general there is no average equilibrium in the whole submembrane domain; this is an effect of the discrete distribution of ion channels that limits the validity of shell models (which assume that Ca²⁺ enters continuously all along the cell membrane). In Fig. 6 D, the equilibrium test for the endogenous buffer was performed by counting calcium ions and buffer molecules in the whole first slice of the conical subdomain and for two different currents (I_d⁰ = 10, 5 pA). The observed displacement is in fact nothing but a manifestation of the [Ca²⁺] gradients along the submembrane domain (0-100 nm of depth). Such gradients are believed to be fundamental to explain the secretory response of neuroendocrine cells (release-ready granules are located beneath the membrane). Notice that the higher the current the more pronounced the deviation becomes at the beginning of the pulse, which indicates that higher submembrane gradients appear. Later, we will discuss such submembrane gradients in more detail.

In contrast to the endogenous buffer, Fura-2 (with the parameters from Table 1) is not "fast enough" to reach "locally" instantaneous equilibrium. Initially, the exogenous buffer deviates considerably from equilibrium, even in the average over homologous compartments. Such local deviation reflects the fact that Fura-2 (100 µM) has relatively slow kinetics, and it cannot account for all of the incoming calcium current (see Appendix B for more details on equilibration times and the RBA hypothesis). In addition, because Fura-2 is a mobile buffer, bound and unbound Fura-2 can approach the membrane coming from deeper regions, changing the ratio of bound/unbound and hence displacing equilibrium, especially when high [Ca²⁺] gradients are present. In this way, the mobility of Fura-2 helps to maintain and enhance the deviation of Fura-2 from equilibrium.

It is also worth noting that the kinetics of the endogenous buffer affects the way Fura-2 deviates from equilibrium, as shown in Fig. 6, E and F. The equilibrium condition is shown for Fura-2; one can observe that displacement from equilibrium changes considerably in the first 10 ms, changing K_D (taking a fixed k) for the endogenous buffer. The free Ca²⁺ concentration profiles, of course, become smaller as K_D is taken smaller (not shown). As expected, displacement from equilibrium for Fura-2 persists even when the endogenous buffer is not considered in the simulation (not shown), because, as discussed in Appendix B (see also Naraghi and Neher, 1997), the origin for such a displacement is to be found in the relatively slow equilibration time for Fura-2.

Equilibrium for Fura-2, particularly near the membrane, was shown to depend strongly on the diffusion coefficient. When we performed a simulation considering a buffering system with the same properties as Fura-2, but taking a null diffusion coefficient, we again observed an initial departure from local equilibrium. However, in this case, equilibrium for such a buffer was reached in less than 2 ms in the presence of endogenous buffer (not shown).

Regarding EGTA, because this buffer is slower than Fura-2 (see Appendix B) and it has a considerably larger mobility, even more deviation from local equilibrium is expected, as indeed happens (Fig. 6 G). Considering a buffer with the same properties as EGTA (Table 1) but with a zero diffusion coefficient, we can see in Fig. 7 D (as happened with Fura-2 with the same change) that the initial displacement from equilibrium is much smaller and has a smaller duration compared to the real case in which EGTA moves. It becomes clear that mobility is acting against local equilibrium when concentration gradients are present; such an effect is, of course, more important, as the buffer considered has slower kinetics.

View larger version (27K): [in this window] [in a new window]   FIGURE 7   Effects of the diffusion coefficient of buffers; averages with binning of 10 points are shown, except in D. (A) Averaged [Ca2+] time course considering two different diffusion coefficients for Fura-2 (DFura-2 = 42 µm2/s and DFura-2 = 200 µm2/s) at different distances from the pores (from top to bottom: 0-100 nm, 100-200 nm, 200-300 nm). The uncertainties are similar to those displayed in Fig. 3 A, and they are not shown here. (B) [Ca2+] time course averaged over the submembrane domain (0-100 nm), considering two diffusion coefficients for EGTA: DEGTA = 200 µm2/s and DEGTA = 0 µm2/s. (C) As in B, for the concentration of free EGTA. (D) Percentage deviation time course for free EGTA, considering the diffusion coefficients of B and C. The test of equilibrium is performed over A-type compartments; the points plotted correspond to the output of the simulation (every 10 steps).

Notice how EGTA approaches saturation levels when a zero diffusion coefficient is assumed (Fig. 7 C). When EGTA is forced to be fixed in the simulation, no unbound EGTA can approach, and, given the small K_D and k for EGTA, it rapidly saturates near membrane.

Given that exogenous buffers are used as indicators of [Ca²⁺], the departure from equilibrium of exogenous buffer concentrations indicates that the measurement of [Ca²⁺] in fluorescence experiments (assuming equilibrium of such buffer with calcium) could grossly underestimate the actual free calcium concentration (as suggested by Smith et al., 1996). Indeed, the equilibrium test is displaced to negative values, which shows that there is less exogenous bound buffer than expected from equilibrium. Such displacement is related to the previously discussed "blindness" of both Fura-2 and EGTA to the calcium time course (Nowycky and Pinter, 1993).

Mobility of buffers

In Fig. 7, A and B, the effect of the mobility of buffers on the [Ca²⁺] time course can be observed.

In Fig. 7 A, we consider the choice of parameters from Table 1 for Fura-2 and the possibility of having a buffer with the same parameters as Fura-2, except that its diffusion coefficient is taken to be approximately equal to the coefficient for calcium; the calcium peak is then observed to be lower as the diffusion coefficient is taken to be larger. This reflects the fact that Fura-2 can now transport calcium to other regions much more rapidly because the diffusion coefficient is larger. The effect of considering higher mobilities is then to produce a more rapid washing out of the regions of high Ca²⁺ concentration. In Fig. 7 B, we repeat the simulation, but now for EGTA (with the values of Table 1). We also consider the possibility of having immobile EGTA (zero diffusion coefficient), to show the effect of mobility. As before, the calcium peak becomes smaller as the mobility is higher.

As previously discussed, the mobility of buffers has a considerable impact on the ratio of unbound to bound buffer molecules, especially when high concentration gradients are present (as happens at the start of the pulse). The relatively slow exogenous buffers (Fura-2, EGTA) show larger deviations from local equilibrium at the beginning of the pulse as their mobility is taken higher (Fig. 7 D). Notice, furthermore, that this fact can affect the saturation levels of the buffer (Fig. 7 C). Thus, if the mobility of exogenous buffers was restricted in cytoplasm, the time course of the free buffer could change considerably in the submembrane region, where high gradients appear.

Submembrane [Ca²⁺] gradients and dependence on the density and distribution of calcium channels

As was previously suggested (Klingauf and Neher, 1997), submembrane gradients depend on the density and distribution of channel pores. One of the advantages of our Monte Carlo simulation lies in the easy and natural implementation of the entry of calcium ions through channel pores and the flexibility with which one can distribute the calcium channels over the membrane. On the contrary, the solution of differential equations can be a hard task when the boundary conditions do not show some regularity or symmetry; besides, there exists the intrinsic difficulty of describing discrete entries for the calcium ions when continuous quantities (as concentrations) are being considered. Such a difficulty becomes even greater if one wishes to discretize the entry of calcium in time, that is, if one wants to consider the probabilities of the opening and closing of the channels. A Monte Carlo simulation is able to deal with such situations. However, although in reality channels open and close stochastically, it is believed that this is not a major effect in neuroendocrine cells (see Klingauf and Neher, 1997). Notwithstanding, examples will be given for a simple model of stochastic opening and closing of calcium channels.

At this stage, we keep to a continuous time description of the calcium current and a discrete spatial treatment of the influx that copes with different distributions of channels, regularly distributed or not.

Let us now study the [Ca²⁺] time courses, both as an average over the submembrane domain (depth 0-100 µm) and at different locations in such subdomain. Let us consider an incoming current given by Eq. 2 with I_d⁰ = 10 pA; thus, the time-averaged current per channel will be 76 fA for a density of channels _cha = 15 µm² (48 channels). Fura-2 (100 µM) will be the exogenous buffer.

Different distributions of channels with different densities, keeping the total current fixed, will be taken into account. Together with the selection of our standard (Table 1) value _cha = 15 µm², we will consider _cha = 5 µm², corresponding to 16 channels. In this last case, the unitary current through each pore will be three times bigger. This last distribution can also be interpreted as a distribution of small groups of three channels, clustered in regions smaller than 0.105 × 0.105 µm², with a unitary average current 76 fA through each of the 48 channels. We also consider the possibility of having the 48 channels clustered in groups of four channels, giving rise to 12 clusters; we place each of the four channels for each cluster in squares of side 2 × 0.105 µm (one channel in each of the four compartments forming such square). In all cases we distribute the location of the channels or channel clusters randomly and uniformly along the membrane of our conical domain, checking that the resulting distributions are statistically representative and approximately uniform. We do not allow overlapping of channels or clusters (notice, however, that the case _cha = 5 µm² corresponds to the overlapping of three channels in a region smaller than 0.105 × 0.105 µm²).

Fig. 8 A shows the [Ca²⁺] time course, averaged over the submembrane domain, corresponding to _cha = 15 µm² and _cha = 5 µm². Notice that the [Ca²⁺] time course averaged over the submembrane domain is very similar in the two cases. When the distribution of 12 channel clusters is considered, the resulting averaged [Ca²⁺] is very similar and lies between the two curves in Fig. 8 A (not shown). As a result, the average submembrane concentrations do not depend crucially on the distribution of calcium channels, which validates the use of radial models (Sala and Hernández-Cruz, 1990; Nowicky and Pinter, 1993) to evaluate such averages.

View larger version (37K): [in this window] [in a new window]   FIGURE 8   Effect of the density of channels. (A) Averaged [Ca2+] time course (0-100 nm) for cha = 5 channels/µm2 and cha = 15 channels/µm2 (average values with binning of 10 points are shown). The [Ca2+] plot for cha = 5 channels/µm2 peaks slightly before that for cha = 15 channels/µm2. (B) [Ca2+] time course at compartments of type A, B, C, D (from higher to lower concentrations) and cha = 15 channels/µm2. In this figure, as well as in C and D, we show the direct output of the simulation each 10 steps of simulation. (C) As in B, for cha = 5 channels/µm2. A and B compartments are indistinguishable in this case. (D) As in B, but for clusters of four channels (see text). In this case there are no B-type compartments.

However, for the different types of submembrane compartments (Fig. 8, B and C), the differences in the [Ca²⁺] time course are evident when _cha = 15 µm² (48 entries) with _cha = 5 µm² (16 entries) are compared, not surpr