## Abstract

We describe a software package for the simulation of exocytotic events from readily releasable pools of secretory vesicles in neuroendocrine cells and presynaptic terminals. The visual package Ca3D_Exolab simulates the entry of Ca^{2+} through the calcium channels, the kinetic reactions of calcium with buffers, the diffusion of calcium and mobile buffers, and the kinetic reactions of calcium with the secretory vesicles. The location of both channels and secretory vesicles can be set by using a graphical interface. Calcium and buffer concentrations at different depths from the cellular membrane and capacitance time courses are obtained as outputs. The software package also provides a descriptive statistical data analysis of the different output data.

- modeling of exocytosis
- Monte Carlo methods

the quantitative modeling of calcium dynamics and related processes is a useful tool for confronting theory with experimental studies of the secretory machinery in neuroendocrine cells and presynaptic terminals. Most computational models of calcium dynamics in cells are based on the solution of differential equations assuming, for instance, a spherical model of the cell with uniform time-dependent distributions of calcium below the membrane (see, for instance, Refs. 7–9 for related software). Although these models provide an average description of the calcium processes, it is extremely difficult to build reliable simulations based on differential equations for describing aspects related to the geometry of the system. For instance, the influence of the spatial distribution of the calcium channels on the calcium and buffer time courses is difficult to model with a differential approach; a spatial organization of the calcium channels in clusters and a random distribution of channels will lead, in general, to very different calcium time courses.

When the secretory machinery of vesicles enters into play, the influence of the geometrical aspects becomes even more important; in this case, the relative position of the calcium channels (organized or not in clusters) and the secretory sensors of the vesicles are crucial parameters of the system. To incorporate information on the geometry of the system, the use of microscopic models, such as those based on Monte Carlo (MC) methods (see Ref. 2, for example), is particularly appropriate. However, because of the complexity of the algorithms involved, MC methods are probably not as popular as differential methods. A “raw” MC code is not an easy thing to use, whereas many mathematical packages are equipped with friendly differential equation solvers. In trying to overcome this limitation, we have built a user-friendly software environment based on the use of MC methods for simulating calcium-triggered secretory processes (from readily releasable vesicles) in cells.

A first version of the software (Calcium3D) for simulating calcium dynamics was presented previously (1). The program now has been enhanced by incorporating the simulation of secretory events from readily releasable pools of vesicles; in this way, accumulated secretory event time courses (proportional to capacitance) can be obtained as outputs. The resulting software (Ca3D_Exolab) includes the following features: different geometrical configurations for both channels and vesicles can be chosen by using a graphical interface, and calcium and buffer concentrations at different depths from the cellular membrane and accumulated secretory event time courses (proportional to capacitance time courses) are obtained as outputs, as well as a descriptive statistical data analysis of the different output data.

## MATHEMATICAL MODELS

The mathematical models behind Ca3D_Exolab were introduced, analyzed, and explained in detail in previous publications (2, 10). We now summarize the main ingredients of the model.

### Domain of Simulation

A three-dimensional (3-D) orthogonal grid maps the (by now) two possible choices of the domain for the simulations: a conical domain and a cylindrical domain. The size Δ*x* of the lateral sides of the cubic compartments of the grid can be selected by the user of the program. This size represents the spatial resolution for the simulations.

The choice of a conical domain is appropriated for simulating secretory events in spherical cells, assuming that the portion of the membrane represented by the base of the cone is embedded inside a larger portion of the membrane where the calcium channels are, on average, uniformly distributed. In this way, one can consider that the net flux of ions and buffers through the lateral sides of the cone is zero. This approximation is also safe for smoothly nonuniform configurations of channels in which the conical domain is placed inside a region where the distribution of channels is sufficiently uniform. The choice of a cylindrical domain seems appropriated for simulating calcium processes in cilliar cells or presynaptic terminals.

### Calcium Current

We consider a simple model for calcium influx in which the calcium current is constant during a pulse. The selection of different shapes of the current and the stochastic modeling of channel gating will be considered in future versions of the software.

### Diffusion of Mobile Particles

For modeling the 3-D diffusion of calcium and mobile buffers, we use a random walk algorithm (2). In this way, calcium ions and mobile buffers move in intervals of time: Δ*t* = (Δ*x*)^{2}/4*D*, where Δ*x* is the spatial resolution for the simulations and *D* is the diffusion coefficient of the mobile particle.

### Secretory Vesicles and Buffers

For simulating the kinetic reactions of the secretory sensors of vesicles and buffers with calcium, we consider a discrete (stochastic) modeling (2, 10). For the secretory vesicles, a noncooperative scheme with a variable number of binding sites (*n*) is considered: where X represents the calcium sensors of the secretory vesicles, *k*_{on} is the on rate of the reaction, *k*_{off} is the off rate of the reaction, and *p* is the fusion rate. The kinetic reactions of calcium ions to binding sites (of both buffers and secretory vesicles) are described by first-order equations: , where [Ca], [B_{i}], and [CaB_{i}] are, respectively, the concentrations of calcium, unbound buffer, and bound calcium. In our simulation scheme, these equations are interpreted stochastically, evaluating the probabilities associated with each of the processes and generating events according to such probabilities (2).

The spatial distribution of the binding sites of buffers is randomly and uniformly distributed over the 3-D domain, whereas the binding sites of the secretory vesicles are spatially restricted to the location of the vesicles. The mechanism for secretion is simulated by placing the number of binding sites, initially free, corresponding to the kinetic scheme considered, in each of the submembrane compartments where the protein lies. In Fig. 1, a schematic representation of the location of vesicles below the membrane is shown. In this representation, a five-binding sites kinetics scheme per vesicle is considered. The calcium channels are also shown.

In our model, the fate of each individual vesicle is recorded. The finite size of the vesicles is neglected: we simulate the kinetics for vesicle fusion, but we do not consider the fact that calcium ions and buffer molecules are not free to diffuse through the space occupied by the vesicles of the readily releasable pool (RRP). This is justified given the small number of vesicles that typically constitute the RRP.

Replenishment of the RRP is not considered in this model. The simulations are restricted to trains of pulses in the subsecond range. Endocytosis and the mobilization of additional vesicles can then be thought to be too slow to be noticed.

## DESCRIPTION OF THE SOFTWARE

Using the theoretical ingredients described above, we have built a software package (Ca3D_Exolab) for simulating secretory events in neuroendocrine cells and presynaptic terminals. The package can be downloaded at http://personales.unican.es/gila/bbva.html.

The MC code for simulating calcium dynamics that, in an improved version, form the core of Ca3D_Exolab, was extensively tested, and its results were proved to match with theoretical expectations at long times (2). Explicit tests for the routines implementing the random walk algorithm and the kinetic reactions of calcium with buffers can be found in Ref. 3.

The outputs of Ca3D_Exolab are calcium and buffer concentration time courses at different depths from the cellular membrane and the number of secretory events as a function of time. For the different output data, a descriptive statistical data analysis is available.

### Two Illustrative Examples

As a practical example, consider two different case studies of secretion: chromaffin cells (see Ref. 11, for example) and calyx of Held (see Refs. 5, 6).

For chromaffin cells, differential methods are generally able to describe the time course of average calcium concentrations below the membrane. Concentrations are high enough so that the averages make sense, and so does the differential modeling of these average values. When high concentrations occur, the MC method implemented uses average descriptions when suitable (see Refs. 2 and 3 for further details).

The calyx of Held is a good example showing when it is not so clear that average descriptions are valid, as commented in Ref. 6. Indeed, calcium concentrations are so low that, for the spatial resolutions considered, very few or no calcium ions may be found, as the results show. This fact, together with the fact that the binding of calcium to vesicle sensors also affects calcium concentration, poses a serious limitation on differential models, which are overcome by the MC approach.

### Chromaffin Cells

We followed descriptions in Refs. 2 and 10 for modeling of these cells. An appropriate domain for simulating localized submembrane time courses is a conical domain. We chose, for example, a cone with a 1-μm radius of base (representing the radius of a portion of the membrane) and 5 μm in height. The spatial resolution Δ*x* was set to 50 nm.

#### Initial conditions and reaction-diffusion parameters.

We considered the following input data for the initial conditions and reaction-diffusion parameters (2, 10): for calcium, initial concentration ([Ca^{2+}]_{rest}) = 0.1 μM, and diffusion coefficient (*D*_{Ca}) = 220 μm^{2}/s; for the endogenous fixed buffer (EFB), total concentration ([EFB]_{total}) = 500 μM, affinity (*K*_{d}) = 10 μM, and forward Ca^{2+} binding rate (*k*_{+}) = 5 × 10^{8} M^{−1}·s^{−1}; for the exogenous mobile buffer fura-2, total concentration ([fura-2]_{total}) = 100 μM, diffusion coefficient (*D*_{Fura-2}) = 42 μm^{2}/s, *K*_{d} = 0.24 μM, and *k*_{+} = 5 × 10^{8} M^{−1}·s^{−1}; and for the kinetic release model for vesicles, there were five independent binding sites, where *K*_{d} = 10 μM and *k*_{+} = 3 × 10^{8} M^{−1}·s^{−1}. In addition, we considered a fusion rate of *p* = 40,000 s^{−1}. The parameters of the system are summarized in Table 1.

#### Topography of channels and vesicles.

Let us, for example, consider two clusters of channels (Fig. 2*A*) and a RRP composed of three vesicles (Fig. 2*B*). For this topographical configuration of channels and vesicles, typical outputs of Ca3D_Exolab are shown in Figs. 3 and 4. Figure 3*A* shows the average [Ca^{2+}] time course from 0 to 50 nm from the cellular membrane.

Figure 3*B* shows the time course of the number of accumulated secretory events. In Fig. 4, the time courses of the average free endogenous and exogenous buffer (fura-2) concentrations (from 0 to 50 nm from the cellular membrane) are shown.

### Calyx of Held

The parameters for the simulation of the calyx of Held presynaptic terminal were taken from Ref. 5. Because the release topography of the system cannot be inferred from experimental studies, the modeling can be used for testing different topographies.

#### Domain of simulation.

In the case of the calyx of Held, we chose as domain for the simulations a cylindrical domain of 400 nm in height and 130 nm in radius (5). This last parameter corresponds approximately to the mean radius of the active zones in the presynaptic membrane, which are the regions of interest.

#### Initial conditions and reaction-diffusion parameters.

We considered the following input data for the initial conditions and reaction-diffusion parameters (5): for calcium, [Ca^{2+}]_{rest} = 50 nM and *D*_{Ca} = 220 μm^{2}/s; for endogenous fixed buffer, [EFB]_{total} = 80 μM, *K*_{d} = 2 μM, and *k*_{+} = 5 × 10^{8} M^{−1}·s^{−1}; for ATP, total concentration ([ATP]_{total}) = 580 μM, diffusion coefficient (*D*_{ATP}) = 220 μm^{2}/s, *K*_{d} = 200 μM, and *k*_{+} = 5 × 10^{8} M^{−1}·s^{−1}; and for the kinetic release model for vesicles, there were five independent binding sites, where *K*_{d} = 10 μM and *k*_{+} = 3 × 10^{8} M^{−1}·s^{−1}. In addition, we considered a fusion rate of *p* = 40,000 s^{−1}.

#### Topography of channels and vesicles.

As mentioned previously, the release site topography at most synapses, including the calyx of Held, cannot be measured experimentally. According to the theoretical study in Ref. 5, at each active zone, the Ca^{2+} that controls phasic release is supplied by 10 or more calcium channels. These channels are grouped into one or a few clusters, with a channel cluster covering an area of <50 nm across. Let us consider a RRP composed of three vesicles associated to this active zone.

We simulated the secretory response of the RRP under two different topographical configurations of vesicles (Fig. 5, *A* and *B*) and a topographical configuration of channels consisting in a single cluster of 12 channels (Fig. 6). The spatial resolution Δ*x* was set to 10 nm.

In Table 2 we summarize the set of input parameters used in the simulations. For simplicity, we assumed a constant current pulse per channel of 0.02 pA and a duration of 1 ms. For the topographical configuration of vesicles shown in Fig. 5*A*, typical outputs of Ca3D_Exolab are shown in Figs. 7–9. Figure 7 shows the time course of the number of accumulated secretory events; Fig. 8 shows the average [Ca^{2+}] time course from 0 to 10 nm from the cellular membrane. The descriptive statistical data analysis for these data is also shown. In Fig. 9, the time courses of the average free endogenous buffer and ATP concentrations (from 0 to 10 nm from the cellular membrane) are shown.

A different secretory scenario appeared when the topographical configuration of vesicles in Fig. 5*B* was considered; in this case, no secretory events were recorded during the simulation time. This illustrates the necessity of being able to simulate different topographies for obtaining information on the feasibility of a given set of kinetic and diffusive parameters. Average descriptions will fail to provide this type of information.

## GRANTS

We acknowledge financial support from Fundación BBVA under the project “Understanding the secretory machinery of living cells: mathematical models and computational tools.” A. Gil acknowledges financial support from the Ministerio de Educación y Ciencia (Programa Ramón y Cajal).

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2007 the American Physiological Society