Calcium waves in differentiated N1E-115 Neuroblastoma cells, experimentally initiated by applying bradykinin (BK), are the result of calcium release from the endoplasmic reticulum (ER) through channels activated by Inositol-1.4.5- triphosphate (InsP₃ ). InsP₃, in turn, is released from the plasma membrane under the stimulus effect and diffuses throughout the cell. The mathematical formulation of the calcium and InsP₃ dynamics incorporates the models of the major components participating in the phenomenon: ER, InsP₃ receptors and SERCA pumps embedded in the ER membrane, plasma membrane extrusion mechanisms, and calcium buffering in the cytosol. Below we describe all the components in detail. To combine them into a complex model, we use a new, redesigned and rigorously tested, version of the Virtual Cell, a general framework for modeling cell biological processes (Schaff et al., 1997). The Virtual Cell allows simulation of the complex dynamics of spatially distributed molecular species while mapping them onto an experimental geometry.

InsP₃ dynamics. InsP₃ dynamics accounts for its production at the plasma membrane, diffusion into the cytosol, and gradual degradation (Allbritton et al., 1992). InsP₃ production is modeled by imposing an appropriate boundary condition at the plasma membrane for the flux density:

(1)

Eq. (1) corresponds to a very rapid (practically instantaneous on the time scale of interest) increase in InsP₃ production immediately after the BK application and the following exponential decay with the time constant k₀. The flux density amplitude  is generally a function of spatial variables because of a non-uniform distribution of BK receptors. Since this is a 2D simulation, this value was also corrected for the difference of surface-to-volume ratios between the soma and neurite regions in the actual 3D cell geometry. See section Reduction of the Model to 2 Dimensions.  The average value of and the time constant k₀ value (see Table 1) were picked to correctly reproduce the experimentally observed time course of the average InsP₃ concentration.

Diffusion and degradation of InsP₃ are described by the equation

, (2)

which, together with the flux jump condition (1), provides the full description of the InsP₃ dynamics in our model. The values of the diffusion constant , the degradation rate , and the basal InsP₃ level at t = 0, [IP₃]₀ , are also given in Table 1.

Note that different membrane conditions and initial conditions are used to model the InsP₃ uncaging experiments. In this case,  = 0 and the InsP₃ concentration at time t = 0 (immediately after UV flash) is assumed to be [IP₃]₀ + [IP₃]_max , where [IP₃]_max is the uncaged amount of InsP₃.

Endoplasmic Reticulum. We use a continuous approximation in our modeling of ER. This assumption should provide an accurate description within the resolution limits of the light microscope used to record Ca²⁺ dynamics. As an internal calcium store, ER is considered to have an infinite capacity and, therefore, is characterized by a constant free calcium concentration .

The ER distribution in a cell is not uniform. The experimental data indicate that the InsP₃ receptor density and the pump density are considerably higher in the soma and experience a sharp decrease at the junction with the neurite. We use these experimental data to modulate the average rate of change of calcium concentration due to fluxes through the calcium channels and pumps. The actual function of spatial variables that was used in the simulations is defined in Section Reduction of the Model to 2 Dimensions.

Calcium Dynamics. The calcium dynamics is governed by the partial differential equation (Sneyd et al., 1995),

, (3)

that takes into account effects of calcium release from ER through the InsP3 sensitive channels, leak, and uptake by the SERCA pumps combined with calcium diffusion and buffering.  are the rates of change of calcium concentration due to the channel release, pump uptake and leak respectively. The factor a is a function of spatial coordinates that modulates calcium fluxes due to the non-uniform ER distribution (see Reduction of the Model to 2 Dimensions).  is the rate of change of calcium concentration due to buffering. Our treatment of calcium buffering is based on the rapid buffer approximation.(Wagner and Keizer, 1994) and is described below in Section Calcium Buffering.

Eq.(3) is accompanied by a jump membrane condition (see Eq. (9) below) that accounts for the plasma membrane extrusion mechanisms. The initial calcium distribution is assumed uniform.

Channel kinetics. We use the simplified version (Li and Rinzel, 1994) of the model proposed in (De Young and Keizer, 1992) to describe the calcium release through an InsP₃ sensitive channel:

(4)

where  is the maximal possible rate;  is the dissociation constant for InsP₃ binding to a channel;  is the dissociation constant for calcium binding to an activation site;  is the calcium concentration in ER; h is the probability of the inhibition site not being occupied, and the governing equation for h is

(5)

where  is the on-rate of calcium binding to the inhibition site and  is the corresponding dissociation constant.

SERCA pump kinetics. The SERCA pump uptake is modeled by a Hill type equation (Gill and Chueh, 1985; Lytton et al., 1992),

(6)

where  is the maximal rate and  is the corresponding dissociation constant. The constant values are given in Table 2.

Leak. The rate of change of the calcium concentration due to calcium leak from ER to the cytosol is modeled as

(7)

where the leak constant L is to be determined from the conditions of a flux balance at a steady state. The initial steady state values of variables are

Calcium Buffering. Calcium buffering has a significant impact on the overall calcium dynamics. The endogenous buffers are usually assumed to be practically stationary with the binding ratio /» 40.0 ( is the total buffer concentration, and is its dissociation constant) and the on rate of the order of 5.0x10⁸ M^-1s^-1 (Xu et al., 1997). Taking into account typical values of , we arrive at very short characteristic times of buffering. The presence of a small parameter makes it possible to apply the pseudosteady approximation (Wagner and Keizer, 1994) when the bound buffer is considered to be in a rapid equilibrium with the local calcium. The fluorescent dyes used in the simulated experiments (Fura-2 in BK experiments and Calcium Green in uncaging experiments), also act as rapid buffers, but are mobile in contrast with the endogenous buffers. Correspondingly, we introduced two types of buffers in our simulations: the stationary buffer B₁, representing endogenous buffers, and the mobile buffer B₂, modeling a fluorescent indicator. The term  from Eq. (3) then takes a form

,

where the variables [B₁], [B₂], [CaB₁], and [CaB₂] are governed by the equations

, (8.1)

, (8.2)

, (8.3)

. (8.3)

Our numerical approach to fast buffering (and, generally, to any fast reactions in biochemical systems) is based on the combination of time splitting and the pseudosteady approximation applied to fast reactions. Time splitting allows us to separate slow and fast processes and use different methods for their treatment. No preliminary analytical manipulation is necessary. Since fast buffering is treated within a pseudosteady approximation, the algorithm does not require information about the kinetic constants of fast buffering. All we need to describe fast reactions are equilibrium constants

, ,

and total buffer concentrations ,  (see Table 3).

Plasma membrane extrusion mechanisms. This type of calcium extrusion (through plasma membrane pumps and - Ca²⁺ exchangers) is important for the correct description of calcium distribution in a cell during the refractory phase of the calcium signal. According to the experimental data (Herrington et al., 1996 ), the corresponding flux, as a function of the cytosolic calcium, exhibits a threshold behavior. There is no flux below the critical concentration [Ca]_c while above this concentration the flux grows linearly as a function of [Ca]. Thus, we introduce the following flux membrane condition for calcium:

(9)

where the rate of flux density g is generally a function of spatial variables depending on the distributions of the plasma membrane pumps and exchangers. These distributions were assumed uniform in our simulations. Similar to the description of InsP₃ fluxes (see Reduction of the Model to 2 Dimensions), the correction of g was made due to the difference of actual 3D surface-to-volume ratios of the soma and the neurite.

Experimental ER, receptor, and pump distributions. To determine the average ER, InsP3-R, BKR, and SERCA2 distributions in N1E-115 neuroblastomas, analysis of immunostained fixed cells was performed. Differentiated N1E-115 neuroblastoma cells were fixed by a 10 minute incubation at room temperature in a 4% paraformaldehyde/0.1% glutaraldehyde solution. After washing with phosphate-buffered saline (PBS), fixed cells were then permeabilized in 0.01-0.1 % Triton-X/PBS solution for 2-5 minutes at room temperature, followed by 2 hours of non-specific blocking in a 5% milk/PBS solution. Following wash, the cells were incubated with the primary antibody (in a 3% bovine serum albumin (BSA)/PBS solution) for up to 19 hours at 4^oC. Primary antibodies used were rabbit anti-IP₃I-receptor (1/100 dilution; generously donated by Barbara Ehrlich), rabbit anti-BK2-receptor (1/50 dilution; generously donated by Werner Esterl-Muller), rabbit anti-BIP (1/100; StressGen), and rabbit anti-SERCA2 (1/10 dilution; generously donated by Sergei Syrbu). For the negative controls, either the primary antibody was omitted or the primary was substituted with the IgG of the appropriate animal. After the wash, the cells were incubated with the secondary antibody (goat anti-mouse or anti-rabbit rhodamine, 1/100 dilution, Molecular Probes) in 3% BSA/PBS for 1 hour at room temperature. For BK imaging, 1/50 goat anti-rabbit Liss-Rhodamine (Jackson Pharmaceuticals) was used as the secondary antibody. Following a final wash, the cover slips were mounted on glass slides in a glycerol/PBS anti-fade solution (SlowFade, Molecular Probes, Inc.) and sealed with nail polish. Z-series were collected on the confocal microscope using 568 nm excitation. These series were imaged with VoxelView (Vital Images, Inc. Fairfield, IA) and analyzed with RatioView (developed by James Schaff) to calculate the fluorescence intensities in 6 regions of each cell: distal soma, proximal soma, proximal neurite, middle neurite, distal neurite, and growth cone. All image analysis was performed on SGI workstations (Silicon Graphics, Inc. Mountain View, CA). These values were exported to an Excel spreadsheet, where all intensities were converted to percentage intensity of each cell?s distal soma. This data is then calibrated according to the calculated convolution occurring in the image, making estimates using 3-D models chosen to best approximate the cellular geometry (Fink et al., 1998). Example images can be found at Biophysical Journal Online.

 The experimental values were then used in the model of the cell.

Reduction of the model to 2 Dimensions. The analysis of the cell geometry indicates that, though the average surface-to-volume ratio of the 2D model s = 0.274 m m^-1 is close to the actual 3D one, the corresponding values for the soma and the neurite are different. Particularly,

while  ( here and in the following the subscripts n and s denote values associated with the neurite and the soma correspondingly). Since the surface-to-volume ratios are important for the overall dynamics, we need to make corresponding corrections of the membrane flux density distributions.

We first consider the InsP₃ flux density. Its amplitude  in Eq. (1) can be described as

(10)

The function b accounts for the experimentally determined non-uniform distribution of the BK receptors. To determine and  in terms of the average flux density J₀, we note that

(11.1)

and , (11.2)

where and are the volume fractions of the neurite and soma in the 2D model correspondingly. Solving Eqs. (11) simultaneously, we obtain

, (12.1)

. (12.2)

Similar corrections should be made for the rate of calcium flux density g in Eq. (9):

(13)

where

, (14.1)

. (14.2)

In Eqs. (12) and (14),  = 0.431,  = 0.569,  0.379 m m^-1, 
0.141m m^-1. The average value of the rate g₀ is given in Table 2.

Finally, we define the function a of Eq. (3) as a= a₀ f where a ₀ = 1.41, and the normalized ER density distribution f is taken from experimental data.

Numerical methods. We use the finite volume method (Patankar, 1980) for solving partial differential equations (PDEs), the forward Euler method for solving ordinary differential equations (ODEs), and the Newton method for solving the system of nonlinear equations for fast buffering (Press at al., 1992). The diffusion terms of PDEs are treated implicitly while the reaction terms - explicitly. We solve the system of equations within the Virtual Cell computational framework (Schaff et al. 1997). Based on the object oriented approach, the framework provides convenient tools for creating all necessary compartments within an arbitrary geometry and all necessary variables within each compartment. It then builds up the discretized version of equations to be solved and employs PDE, ODE or algebraic solvers where they are necessary. The solvers have been rigorously tested against exact solutions. The system of linear algebraic equations, obtained after the discretization of PDEs, is iteratively solved by the line-by-line method (Patankar, 1980) using 6 iterations for each dimension. We decomposed the 2D computational domain into 7.7x10³ control volumes of the 1.2 m m mesh size and ran simulations using a 0.01 sec time step. The computations were performed on an SGI workstation and ordinarily took about 25 min per simulation. Comparing results, obtained with various mesh size and time increments, we estimate that the computational error does not exceed 3%.

TABLE 1. InsP₃ dynamics.
 

Parameter	Symbol	Value	Units	Reference
Average flux density amplitude		9.14	m M m m s-1	fit to experiment
Flux time constant	k0	0.411	s-1	fit to experiment
Degradation rate		0.16	s-1	(Wang et al., 1995), present work
Initial InsP3 concentration	[IP3]0	0.16	m M	experimentally determined
Uncaged amount of InsP3	[IP3]max	1.5	m M	experimentally determined
InsP3 diffusion coefficient		283	m m2 s-1	(Allbritton et al., 1992)

TABLE 2. ER and plasma membrane Ca²⁺ fluxes and diffusion.
 

Parameter	Symbol	Value	Units	Reference
Average amplitude of channel release	Jmax	3500	m M s-1	(Bezprozvanny et al., 1991; Kupferman et al., 1998), present work
Average amplitude of pump uptake	Vmax	1.59	m M s-1	fit to experiment
Leak constant	L	0.00858	m M s-1	determined by the steady state balance
Dissociation constant for Ca2+ binding to an activation site		0.2	m M	constrained by stability conditions
Dissociation constant for Ca2+ binding to an inhibition site		0.11	m M	constrained by stability conditions
On-rate for Ca2+ binding to an inhibition site		2.1	m M-1 s-1	fit to experiment
Dissociation constant for InsP3 binding to a channel		0.8	m M	fit to experiment
Dissociation constant for Ca2+ binding to a pump		0.27	m M	(Gill and Chueh, 1985; Lytton et al., 1992)
Initial Ca2+ concentration	[Ca]0	0.05	m M	from indicator fluorescence
Ca2+ concentration in ER		400	m M	(Miyawaki et al., 1997; Meldolesi and Pozzan, 1998)
Average rate of Ca2+ flux density at the plasma membrane		11.0	m m s-1	(Herrington et al., 1996), fit to experiment
Threshold concentration for calcium extrusion at the plasma membrane	[Ca]c	0.2	m M	(Herrington et al., 1996)
Ca2+ diffusion coefficient		220	m m2 s-1	(Klingauf and Neher, 1997)

 

TABLE 3. Calcium buffering.
 

Parameter	Symbol	Value	Units	Reference
Total concentration of stationary endogenous buffers (bound and unbound)	B1	450	m M	(Xu et al., 1997; Klingauf and Neher, 1997)
Dissociation constant for Ca2+ binding to an endogenous buffer	K1	10	m M	(Klingauf and Neher, 1997)
Total concentration of the exogenous buffer (BK experiments - Fura-2)	B2	75.0	m M	experimentally determined
Dissociation constant for Ca2+ binding to the exogenous buffer (Fura-2)	K2	0.24	m M	(Klingauf and Neher, 1997)
Diffusion coefficient of the exogenous buffer (Fura-2)		50	m m2 s-1	(Klingauf and Neher, 1997)
Total concentration of the exogenous buffer (uncaging experiments - Calcium Green)	B2	15.0	m M	measured [CG-1]cyt
Dissociation constant for Ca2+ binding to the exogenous buffer (Calcium Green)	K2	0.26	m M	(Eberhard and Erne, 1991)
Diffusion coefficient of the exogenous buffer (Calcium Green)		18.4	m m2 s-1	calculated based on diffusion of a sphere

References.

Allbritton, N., T. Meyer, and L. Stryer. 1992. Range of messenger action of calcium ion and inositol 1,4,5-triphosphate. Science. 258: 1812 - 1815.

Bezprozvanny, I., J. Watras., and B. E. Ehrlich. 1991. Bell-shaped calcium-response curves of Ins(1.4.5)P3- and calcium gated channels from endoplasmic reticulum of cerebellum. Nature 351: 751 - 754.

De Young, G. W., and J. Keizer. 1992. A single -pool inositol 1,4,5-triphosphate-receptor-based model for agonist-stimulated oscillations in Ca2+ concentration. Proc. Natl. Acad. Sci. U. S. A. 89:9895 - 9899.

Eberhard, M., and P. Erne. Calcium Binding to fluorescent Calcium Indicators: Calcium Green, Calcium Orange, and Calcium Crimson. Biochem. Biophys. Res. Comm. 180(1): 209 - 215.

Fink, C., F. Morgan, L. Loew. 1998. Intracellular Fluorescent Probe Concentrations by Confocal Microscopy. Biophys. J. 75 (4):1648-1658.

Gill, D., and S.-H. Chueh. 1985. An intracellular (ATP + Mg2+)-dependent calcium pump within the N1E-115neuronal cell line. J. Biol. Chem. 260 (16): 9289 - 9297.

Herrington, J., Y. B. Park, D. F. Babcock, and B. Hille. 1996. Dominant role of mitochondria in clearance of large Ca2+ loads from rat adrenal chromaffin cells. Neuron. 16: 219 - 228.

Klingauf, J., and E. Neher. 1997. Modeling buffered ca2+ diffusion near the membrane: implications for secretion in neuroendocrine cells. Biophys. J. 72:674 - 690.

Kupferman, R., P. P. Mitra, P. C. Hohenberg, and S. S.-H. Wang. 1997. Analytical calculation of calcium wave characteristics. Biophys. J. 72:2430 - 2444.

Lytton, J., M. Westlin, S. E. Burk, G. E. Shull, and D. H. MacLennan. 1992. Functional comparisons between isoforms of the sarcoplasmic or endoplasmic reticulum family of calcium pumps. J. Biol. Chem. 237:14483 - 14489.

Li, Y. X., and J. Rinzel. 1994. Equations for InsP3 receptor mediated [Ca2+]i oscillations derived from a detailed kinetic model: a Hodgkin-Huxley like formalism. J. Theor. Biol. 166:461 - 473.

Meldolesi, J., and T. Pozzan. 1998. The endoplasmic reticulum Ca2+ store: a view from the lumen. TIBS 23: 10 - 14.

Miyawaki, A., J. Llopis, R. Heim, J. M. McCaffery, J. A. Adams, M. Ikura, and R. Y. Tsien. 1997. Fluorescent indicators for Ca2+ based on green fluorescent proteins and calmodulin. Nature 388: 882 - 887.

Patankar, S. V. 1980. Numerical Heat Transfer and Fluid Flow. Taylor&Francis.

Press, W. H., S. A. Tuekolsky, W. T. Vetterling, and B. P. Flannery. 1992. Numerical Recipes in C (The Art of Scientific Computing). Cambridge University Press.

Schaff, J., C. C. Fink, B. Slepchenko, J. H. Carson, L. M. Loew. 1997. A general computational framework for modeling cellular structure and function. Biophys. J. 73:1135 - 1146.

Sneyd, J., J. Keizer, and M. J. Sanderson. 1995. Mechanisms of calcium oscillations and waves: a quantitative analysis. FASEB J. 9:1463 - 1472.

Wagner, J., and J. Keizer. 1994. Effects of rapid buffers on Ca2+ diffusion and Ca2+ oscillations. Biophys. J. 67:447 - 456.

Wang, S. S.-H., A. A. Alousi, and S. H. Thompson. 1995. The lifetime of inositol 1,4,5-triphosphate in single cells. J. Gen. Physiol. 105:149 - 171.

Xu, T., M. Naraghi, H. Kang, and E. Neher. 1997. Kinetic studies of Ca2+ binding and Ca2+ clearance in the cytosol of adrenal chromaphin cells. Biophys. J. 73:532 - 545.