- Email: [email protected]
CA3D: A Monte Carlo code to simulate 3D buffered diffusion of ions in sub-membrane domains
Computer Physics Communications 136 (2001) 269–293 www.elsevier.nl/locate/cpc
CA3D: A Monte Carlo code to simulate 3D buffered diffusion of ions in sub-membrane domains Amparo Gil ∗ , Javier Segura 1 Instituto de Bioingeniería, Universidad Miguel Hernández, 03202 Elche, Alicante, Spain Received 26 June 2000; accepted 21 December 2000
Abstract CA3D is a simulation program, based on Monte Carlo techniques, which simulates the buffered diffusion of ions in the submembrane domain of cells. Two geometries of submembrane domains are considered. The first, a conical domain, applies to spherical cells with a locally uniform distribution of channels while the second, a cylindrical domain, is intended to model ciliar cells and presynaptic terminals. Other domains can be implemented in the code with ease. The entry of ions takes place through discrete entry points (channel pores). Ions diffuse inside the cell and react with endogenous or exogenous buffers, mobile or fixed. The diffusion of ions and binding molecules is modeled as a 3D random walk process by moving each individual ion and molecule. The reaction of ions with the binding molecules is also modeled probabilistically. The outputs of the code are concentration time courses and 2D maps of concentration of ions in the sub-membrane domain of cells. 2001 Published by Elsevier Science B.V. PACS: 87.10.+e; 87.22.As; 02.70.Lq Keywords: Ion diffusion; Random walk; Buffering; Reaction-diffusion
PROGRAM SUMMARY
Computers for which the program is designed an others on which it is operable: SUN Enterprise 3000; PC Pentium III 450 MHz with 128 MB RAM; PC Pentium II 350 MHz with 64 MB RAM
Title of program: CA3D Catalogue identifier: ADNW Program Summary URL: http://cpc.cs.qub.ac.uk/summaries/ADNW Program obtainable from: CPC Program Library, Queen’s University of Belfast, N. Ireland
Operating systems under which the program has been tested: Solaris (SunOS 5.6), GNU/Linux, Win32 Programming language used: FORTRAN 77 Compilers: f77 version SC4.2 for SunOS, g77 v0.5.24 for Linux, g77 v0.2.95 for MinGW32, Compaq Visual Fortran for Windows No. of bytes in distributed program, including test data, etc.: 15 357
* Corresponding author.
E-mail address: [email protected] (A. Gil). 1 Present address: Departamento de Matemáticas, Universidad Carlos III de Madrid, Leganés Madrid, Spain.
0010-4655/01/$ – see front matter 2001 Published by Elsevier Science B.V. PII: S 0 0 1 0 - 4 6 5 5 ( 0 1 ) 0 0 1 4 4 - 8
270
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
Distribution format: tar gzip file Supplementary material: Pseudocolored animated images of calcium 2D maps obtained from the simulation Keywords: Ion diffusion, random walk, buffering, reaction-diffusion Nature of the physical problem The buffered diffusion of ions in the sub-membrane domain of cells is simulated by means of Monte Carlo methods. Two geometries of submembrane domains are considered. The first, a conical domain, applies to spherical cells with a locally uniform distribution of channels while the second, a cylindrical domain, is intended to model ciliar cells and presynaptic terminals. Other domains can be implemented in the code with ease. The entry of ions takes place through discrete entry points (channel pores) which are distributed over the compartments of the first slice of the domain (base of the cone or the cylinder). Ions diffuse inside the cell and react with endogenous or exogenous buffers, mobile or fixed.
Method of solution The diffusion of ions and binding molecules is modeled as a 3D random walk process by moving each individual ion and molecule. The reaction of ions with the binding molecules is also modeled probabilistically.
Restrictions on the complexity of the problem The number of radial grid points is restricted to 100 (diameter 200 points). The depth is limited to 100 grid points. In case the grid resolution and the sizes of the domain are chosen in such a way that this limit is surpassed, a warning message is generated and the program gives the option of continuing the simulation by taking the largest possible sizes of the domain, maintaining the spatial resolution. The number of grid points available could be smaller in other (generally older) computer systems, depending on the RAM available. In all the systems checked (with a RAM of 64 MB or higher) the restrictions on the number of grid points described above were always within computable ranges. The dimension of the arrays limit the resolution and size of the domain of simulation that can be taken into account. In this way, for domains with a typical size of L, the best resolution achievable is of the order of L/100. Typical running time Depends on the buffers, grid resolution, unitary current, etc. For example, for a conical domain of 1 µm radius and 5 µm deep, a calcium time course lasting 1 ms takes 34 s of CPU time for a resolution of 100 nm and 378 s for a resolution of 50 nm; in both cases a standard endogenous buffering is considered with the addition of an exogenous buffer (Fura-2 100 µM). See section “CPU times” for more details.
1. Introduction The diffusion of different types of ions in the cytosol of living cells play a key role in their function. For instance, it is well known that the fluxes of Na2+ and K+ across the membranes of excitable cells are responsible for the propagation of action potentials, described by the celebrated Hodgkin–Huxley cable equations. On the other hand Ca2+ is known to act as a second messenger which triggers a number of cellular functions such as: secretion of hormones (like, for instance, insulin), muscle contraction, release of neurotransmitter and others. In particular, much attention has been paid to the diffusion of Ca2+ in intracellular media and there is a vast number of biological publications which rely on differential methods (see for instance [1,2]). Ca2+ diffuses inside the cytosol while at the same time first order kinetic reactions take place by which the Ca2+ ions can bind with a binding site of a given buffer B to give rise to a bound state BCa: Ca2+ + B BCa.
(1)
The diffusion of calcium in the cytosol is also referred to as calcium buffered diffusion. Previous models to describe buffered diffusion were based on differential equations, which need drastic geometrical simplifications in order to be solved with ease. Depending on the geometrical assumptions considered the differential models can lie within two categories: shell models (see for instance [1]), in which it is assumed that ions enter uniformly over the cell membrane and only the radial diffusion is of relevance (for spherical cells) and microdomain models (see [2]), which solve the reaction diffusion equations in the vicinity of a channel pore, assuming azimuthal symmetry around the pore; in the latter case it is necessary to assume that all channel pores are equidistant. However, it is well known that ions do not enter uniformly given that the channel pores are not regularly spaced. These departures from symmetry may have an important impact on the functionality of the cells. In [3] it was described how the clustering of channels gives rise to calcium profiles which can not be obtained as a
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
271
superposition of calcium concentrations developed by each individual channel. This non-linearity amplifies the effect of irregular distributions of calcium channels on the formation of submembrane calcium concentration profiles. As an example of the physiological relevance of such non-uniform configurations we can mention the case of secretion in neuroendocrine cells, which is a process triggered by calcium. There are indications that the vesicles responsible for secretion seem to be non-uniformly distributed [2,4,5]; this fact, together with the complexity of the calcium concentration profiles below the membrane clearly points towards the necessity of a modeling which is free of strong simplifying assumptions regarding the geometry of the system [6]. With the idea in mind to overcome the geometrical simplifications imposed by differential methods, we developed a Monte Carlo simulation for the problem of 3D buffered calcium diffusion [3]. This scheme has proved to be successful in the study of the influence of the geometry in the exocytotic response of neuroendocrine cells [6]. For the spatial resolutions that are relevant (of the order of 100 nm in neuroendocrine cells) and the typical concentrations of calcium (order 10 µM) and calcium buffers (order 100 µM) we can expect few calcium ions inside each cube of 100 nm of side (∼7 ions for 10 µM) and also a moderate number of calcium binding sites. Of course, as we demand a higher spatial resolution, less ions and buffer molecules are expected. Then, it is more appropriate to keep track of each individual ion and molecule instead of using a continuous modeling in terms of concentrations. The code CA3D is a microscopic simulation in which the fundamental variables are the number of ions and buffers; in contrast, the differential methods are macroscopic and deal with concentrations. The average values of the output of our simulation converge to macroscopic simulations when considering symmetrical configurations. The outputs of the code are transformed back into macroscopic variables (concentrations). The program gives as output ion and buffer concentration time courses at different depths as well as 2D maps of concentrations of ions in the submembrane domain. Pseudocolored images of calcium 2D maps obtained from the simulation can be found in http://bessel.umh.es. 2. Description of the method 2.1. Geometrical configuration We consider two basic configurations: a cylindrical volume and a conical volume. As described in Ref. [3], the conical domain is appropriate to describe buffered diffusion in the submembrane domain of spherical cells. The base of the cone represents the membrane of the cell. Assuming that the patch of membrane considered is embedded inside a larger portion of the membrane where calcium channels are, on average, uniformly distributed, we can assume that the flux of ions and buffer molecules through the lateral sides of the cone is zero. This simplification imposes a simple boundary condition in the walls of the membrane: ions bounce back when they encounter a wall. At the base (membrane) this is interpreted as reflection while, on the lateral sides, reflection is equivalent to considering that the outgoing flux equals the incoming flux, due to the symmetry of the problem. An orthogonal 3D regular grid maps the domain of simulation with a distance between grid points chosen by the user. Typical operative ranges for this spatial resolution (l) lie between 10 and 100 nm. Each point of the grid is associated with a cubic compartment of volume (l)3 . Therefore, grid points or cubic compartments are equivalent concepts. Ions and buffer molecules move from one compartment to another compartment of the grid due to diffusion, which is modeled as a random walk process. The time scale for the random walk is given by the spatial resolution (which can be selected by the user) and by the higher diffusion coefficient of the system, which is supposed to be the diffusion coefficient for the ions. The relation between those scales is given by: t = (l)2 /4D
(2)
being l the spatial resolution and D the diffusion coefficient for the ions. As an example, for a spatial resolution of l = 100 nm and for calcium ions (DCa2+ = 220 µm2 /s in cytoplasm [7]), Eq. (2) gives a time resolution
272
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
Fig. 1. Flux diagram of the program.
t = 12.5 µs. The time resolution is the time step of the simulation: ions will be moved every time step t using a random walk algorithm in which each ion will have a 50% chance to remain in its initial position and a 25% chance to move either in the positive or negative direction for each spatial dimension. The cylindrical domain fits to study Ca2+ signals in ciliar cells or in presynaptic terminals. We also impose total reflection on the boundaries of the cylinder. 2.2. Sequence of the simulations The sequence of the simulations follows the flux diagram shown in Fig. 1. Before entering the simulation loop, the main program calls the subroutines makegrid, data, districa and chandis for setting a discretization of the simulation domain, transforming the macroscopic input variables into microscopic quantities, setting the initial conditions and distributing the channel pores over the upper slice of the simulation domain. By upper slice we mean the submembrane domain (base of the cone or the cylinder). After the initial data have been processed and the simulation domain and initial conditions have been built, the program enters in the simulations loop. The main ingredients of this loop are the routines walk and kine, combined with the entry of ions through the channels. This last task is performed by the main program which also handles the outputs. 2.3. Input data processing First, the orthogonal grid is built by the subroutine makegrid which, as previously commented, maps the simulation domain to a 3D orthogonal grid. The resolution demanded is one of the input parameters of the program together with the size of the domain. After building the grid the resolution is slightly adjusted in order to fit the grid with the dimensions of the domain. This small correction (smaller as the chosen resolution is better) reduces the possible underestimation of the membrane surface by the rounding in the process of making the grid. The output files of the program show the actual resolution used. After building the grid, the program calls the subroutine data in order to transform the macroscopic input data (concentrations, kinetic rates, . . .) into microscopic data (number of molecules in each compartment, reaction probabilities per molecule, . . .). The subroutine data also evaluates the initial conditions for buffers and ions; it is assumed that ions and buffers are initially at equilibrium. Among other parameters, data evaluates the total number of free and bound buffers
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
273
from the equilibrium condition taking as input the total concentrations of buffers as well as the basal (initial) concentration of free ions. After data conversion by the routine data, districa distributes uniformly the initial (basal) ion concentrations and the concentration of buffers (free and bound) over the whole simulation domain. Regarding the distributions of the entry points for the ions (ionic channels), when the user selects a randomly generated distribution, the program calls the routine chandis to distribute the number of channels specified over the upper layer of the domain (representing the cell membrane). The user has also the option of determining the position of the channels in the file “channel.dat” by giving the grid coordinates of the channels in the upper layer of the grid. 2.4. Simulation loop In the simulation loop three consecutive actions are repeated for each simulation step t: (1) During the time intervals for which there is an influx of ions (as determined in “input.dat”) , the main program sets, for each time step, a number of extra ions in the submembrane compartments corresponding to the position of a channel. (2) Ions and buffer molecules suffer a random walk step. (3) The kinetic equations for the binding of ions by buffers evolves during the time lapse t. 2.4.1. Influx of ions The influx of ions is controlled by the main program. By default, the program simulates trains of pulses with constant incoming current. The user should determine in the file “input.dat” the number and duration of the pulses, the frequency at which they are delivered as well as the (constant) magnitude of the current through each channel pore. Then, using the microscopic current (number of incoming ions in each t) evaluated in data 2 the program sets the extra number of ions “below” each channel in each t. Given that the number of incoming ions in each t could be smaller than 1, the program chooses probabilistically the number of ions which enter through each channel to avoid rounding errors. The user has also the possibility of choosing different shapes of the unitary current. This has to be specified in the following piece of the main program: cccccccccccccccccccccccccccccccccc c For non constant pulses enter c c here your choice of function c c for the unitary current, c c scaling the currents and the c c time variable with the factors c c scei and scet (see below). c c c c 1) SET fnio= YOUR FUNCTION c c c c Here is an example for a c c unitary current cccccccccccc c I(t)=0.213exp(-t/20) pA (t in ms) c c lasting 50 ms c c c 2 Note: the program is prepared to simulate calcium buffered diffusion. For other ions, the parameter charge in data should be modified to
account for the number of charges per ion.
274
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
c if(i.lt.(50*scet))then c fnio=0.213*scei*exp(-i/(20*scet)) c else c fnio=0 c endif c c c (the time variable is i) cccccccccccc c c c 2) UNCOMMENT THE LINES: c c c c npulse=10000000 c ifnio=int(fnio) c fnio=fnio-ifnio cccccccccccccccccccccccccccccccccc c WARNING: you can not simulate c c trains of pulses unless you c c explicitly specify this fact c c in your function c cccccccccccccccccccccccccccccccccc 2.4.2. Subroutine walk: diffusion The 3D diffusion of ions and mobile buffers is modeled as a random walk process. After each time step t = (l)2 /4D and for each dimension, a particle has a 50% chance of remaining in its initial position and a 25% chance of moving either in the positive or the negative direction. Another possible prescription to implement the random walk is considering a time step t = (l)2 /2D and assigning 50% probability of movement in either direction for each spatial dimension. 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. Depending on the mobile buffer considered (EGTA, Fura-2, . . .), the corresponding diffusion coefficient can be typically about 2 to 5 times smaller than that for Ca2+ so that the characteristic diffusion times are about 2 to 5 times bigger. Thus, the diffusion of buffer particles is simulated by moving them in exactly the same way that ions are moved but only after 2 to 5 simulation steps. This means that the program takes integer ratios between the largest diffusion coefficient (ions) and the smallest (buffers). Given the experimental uncertainties associated with diffusion coefficients in the cytoplasm this approximation is accurate enough and saves a considerable amount of computational time. An additional way of saving CPU time is by considering “bulk” diffusion when the number of particles N to diffuse in a given compartment is greater than 64 [3]. In such a case, we take N = 64c + R, with R < 64 (c and R integer numbers). Then, we consider the “bulk” diffusion of the 64c buffer particles according to the probabilities of moving in each direction of the 3D grid [3] and we decide individually and probabilistically the position of the remaining R particles. It is preferable to consider the 3-point implementation of the random walk (with t = (l)2 /4D) to make the bulk diffusion smoother and without large local calcium fluctuations. Regarding the reflection in the boundaries of the sub-domain we take this effect into account by rejecting those events which result in the positioning of a particle outside the sub-domain; 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 sub-domain). 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 being equivalent to considering that the number of outgoing particles equals the number of incoming particles.
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
275
2.4.3. Subroutine kine: kinetics During each diffusional time step t and in each compartment of the grid, the kinetic equations which 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 d[Ca] i i = k− [CaBi ] − k+ [Ca][Bi ] , dt i (3) d[Bi ] d[Ca] d[CaBi ] d[Bi ] = ; = − , dt dt dt dt i
i k+
i , the unbinding rate. being the forward binding rate for the buffer i and k− Since 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 specie in each small cubic compartment of the grid. In terms of the number of free calcium ions NCa , free buffer molecules NBi and bound buffer molecules NCaBi , Eq. (3) reads: i i NCa = K− NCaBi − K+ NCa NBi , i (4) + N NCaBi + NCa = NCaT , N Bi = NBT ; CaBi i i
where NBT represents the total (free + bound) number of buffer particles of type i, NCaT is the total number of (cali cium) ions and NCa represents the variation of free calcium ions in a given cubic compartment during a time interi and Ki for unbinding and binding in each time interval dt are given by val dt ≡ t/n. The effective constants K− + t i i t ; K+ , (5) = k+ n nV where n is the number of sub-divisions of the interval t used to integrate the kinetic equations and V is the volume of each cubic compartment. In order to implement a microscopic interpretation we describe the system probabilistically, evaluating the probabilities associated with each of the processes and generating events according to such probabilities [3]. During the time dt, the unitary probabilities (probabilities for the binding of a single ion or unbinding of a single bound i N for binding and P i = Ki for unbinding, provided we have taken a sufficiently buffer molecule) are Pbi = K+ Bi − u small dt (n large enough). n is chosen in such a way that all probabilities are small enough. Numerical experiments show that taking 0.1 as the maximal probability (for binding or unbinding) is an accurate selection. The number of calcium ions that become bound in each interval dt or become unbound are evaluated using the corresponding probabilities and following a binomial distribution N (6) P (k; p, N) = pk (1 − p)N−k , k i i = k− K−
where P (k; p, N) is 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 molecules) in a time dt, p being the unitary probability (Pbi or Pui ). The code BTPEC [8], obtainable from NETLIB repository [9], is used to generate binomial random deviates. This program uses the uniform random number generator RAND1 [10], which is also used in the rest of the simulation program. The algorithm for kinetics here described gives priority to the fastest effects during each t. For instance, unbinding is normally so slow that the probabilistic selection of the number of buffer molecules that undergo unbinding can be done 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.
276
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
2.5. Description of the inputs, outputs and parameters Inputs: The input data for the program must be written in the file ‘input.dat’. Their meaning is the following: • icondit=logical flag which controls the timing option (uses the non-standard function secnds). The value by default is .false. and the program lines for timing appear as commentaries. To use this function set icondit to the value .true. and uncomment the corresponding lines in the program (search for the word “TIMING” in the program to locate these lines). This function will estimate the duration of a run of the program by counting the number of iterations of the simulation loop during the first 10 s of CPU time and extrapolating the time for the whole run. • imicro=kind of microdomain. If imicro=1, the selected microdomain is conical while imicro=2 corresponds to a cylindrical domain. • rbas=radius of the base of the microdomain in µm. • rdepth=depth of the microdomain in µm. • ncha=number of ionic channels. • idncha=distribution of ionic channels. If idncha=1, the selected distribution is random and uniform. If idncha=2, the distribution of channels can be specified by the user: the program reads the coordinates of the entry channels from the file ‘channel.dat’. • deltax=spatial resolution in µm. • npulses=number of calcium pulses (or ionic pulses in general). • tpulse=duration of each pulse in ms (assumed equal for all the pulses). • fhz=frequency at which the pulses are delivered. • tsimul=duration of the simulation in ms. • di0=unitary current in pA (current of ions through each channel). • difca=diffusion coefficient of ions in the cytosol in µm2 /s. • difbf=diffusion coefficient of buffer 1 in µm2 /s. • difbm=diffusion coefficient of buffer 2 in µm2 /s. • caini=basal (initial) concentration of ions in µM. • bfixed=concentration of buffer 1 in µM. • bmobil=concentration of buffer 2 in µM. 1 ) in M−1 s−1 . • dkbf=forward binding rate of buffer 1 (k+ 1 = k 1 /k 1 ) in µM. • dkdf=dissociation constant of buffer 1 (KD − + 2 • dkbm=forward binding rate of buffer 2 (k+ ) in M−1 s−1 . 2 = k 2 /k 2 ) in µM. • dkdm=dissociation constant of buffer 2 (KD − + • nslice=number of layers of the grid (at progressive depths) where the averaged [Ca2+ ] is required. Maximum: 9. • ibuf=choice for averaged buffer concentrations. If ibuf=0(1), the averaged buffer concentrations at the nslice layers are not (are) calculated. • nmaps=number of 2D ion concentration maps to be calculated. Maximum: 25. • tmap(1:nmaps)=times at which an ion concentration map is obtained in ms (include as many lines in the input file as concentration times demanded) Outputs: The outputs provided by the program are averaged ion concentrations as a function of time (for each of the first nslice layers of the microdomain) or 2-dimensional ion concentration maps (with a temporal bining given as a parameter). The output files are the following: • Averaged ion concentrations in the first nslice layers: files “caslixx”, where xx= 01, 02, 03, . . . gives the depth of the slice (in grid units, 1 grid unit=resolution).
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
277
• Averaged free and bound buffer concentrations in the first nslice layers (for ibuf=1): files “b1fslxx”, “b1bslxx”,“b2fslxx”,“b2fslxx”,“b2bslxx” with xx= 01, 02, . . .. • 2D ion concentration maps: files “camapxx”, xx= 01, 02, . . . . Parameters: • ibin=number of temporal points for the bining of the 2D calcium maps (default: 10 points). 2.6. Description of the subroutines and dependencies The main program calls the subroutines: makegrid, data, chandis, districa, walk, kine, with a calling sequence as described in Fig. 1. We now describe the purpose of these subroutines as well as their dependencies: makegrid(imicro,ncomps,ncomph,jl,il,numcom,nt) Purpose: Routine to make the grid for a conical or a cylindrical domain. Arguments: Inputs: imicro: selects the geometry of the domain. If imicro=1, the selected domain is conical while imicro=2 stands for a cylindrical domain. ncomps: number of grid points along the radius of the base of the domain. ncomph: number of grid points along the depth of the domain. Outputs: jl(1:ncomph), il(1:ncomph,jl): range of the coordinates of the grid points as a function of depth. numcom(1:ncomph): number of compartments as a function of depth. nt: total number of compartments of the domain. This routine is called only once and at the start of the simulation. After the call to the routine makegrid the domain of simulation becomes defined by the following set of grid points: {{{i |il(k, j )|}, j |j l(k)|}, k = 1, 2, . . . , ncomph}, where k measures the depth and i, j are the “horizontal” coordinates. k = 1 represents the sub-membrane domain. data(deltax,difca,difbf,difbm,di0,fhz,tpulse,tsimul,nmaps,tmap) Purpose: transformation of the input data to microscopic quantities used in the simulation and setting of the initial conditions. Arguments: Inputs: deltax: selected resolution for the simulation. difca: diffusion coefficient of Ca2+ . difbf: diffusion coefficient of buffer 1. difbm: diffusion coefficient of buffer 2. di0: unitary current (per channel). fhz: frequency of delivery of ionic pulses. tpulse: duration of each pulse. tsimul: duration of the simulation. nmaps: number of 2D calcium maps. tmap(1:nmap): times at which the 2D calcium maps are obtained. Parameters: charge: number of charges of one ion (default: charge=2). The outputs of this routine, which is called only once by the main program, are passed to the main using common declarations. This routine gives as output the total number of free ions and free and bound buffers
278
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
of both types, together with the time step for the simulation in ms. The routine data also renormalizes the rest of macroscopic constants (kinetic and geometrical constants) transforming them to microscopic quantities. chandis(iseed,ncha,ie,je) Purpose: generation of a uniform and random distribution of channels. Arguments: Inputs: iseed: seed for the simulation. ncha: number of calcium channels to be distributed. Outputs: ie(1:ncha),je(1:ncha): coordinates of the channels (the depth coordinate being k = 1). This routine is called only once and at the start of the simulation. districa(iseed,numcom,inica,nt,ni) Purpose: distributes a given specie of particle (ions, free buffer molecules, . . .) all over the domain of simulation, being all compartments of the grid equally probable. Arguments: Inputs: iseed: seed for the simulation. numcom(1:ncomph): number of compartments of the domain as a function of the depth. inica: number of particles to be distributed. nt: total number of compartments in the domain. Outputs: ni(-il(ncomph,jl):il(ncomph,jl),-jl(ncomph):jl(ncomph),ncomph): initial number of particles in each compartment of the grid. This routine is called only once and at the start of the simulation. kine(iseed,ni,nb1,nb2,nbm1,nbm2) Purpose: implementation of the first order kinetics reactions in each compartment of the domain. Arguments: Inputs: iseed: seed for the simulation. ni(-il(ncomph,jl):il(ncomph,jl),-jl(ncomph):jl(ncomph),ncomph): number of calcium ions in each compartment of the domain. nb1(-il(ncomph,jl):il(ncomph,jl),-jl(ncomph):jl(ncomph),ncomph): number of free buffer molecules of type 1 in each compartment of the domain. nb2(-il(ncomph,jl):il(ncomph,jl),-jl(ncomph):jl(ncomph),ncomph):number of free buffer molecules of type 2 in each compartment of the domain. nbm1(-il(ncomph,jl):il(ncomph,jl),-jl(ncomph):jl(ncomph),ncomph): number of bound buffer molecules of type 1 in each compartment of the domain. nbm2(-il(ncomph,jl):il(ncomph,jl),-jl(ncomph):jl(ncomph),ncomph):number of bound buffer molecules of type 2 in each compartment of the domain. Outputs: The same arrays as for the inputs but actualized after the changes made by the kinetic reactions. walk(iseed,ni) Purpose: implementation of a 3D random walk algorithm for the movement of particles. Arguments:
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
279
Fig. 2. Distributions for a group of 106 particles, initially at k = 0, after iterating 500, 1000 and 2000 times the routine for random-walk diffusion (walk) (from top to bottom). The comparison with Eq. (8) is also shown.
Inputs: iseed: seed for the simulation. ni(-il(ncomph,jl):il(ncomph,jl),-jl(ncomph):jl(ncomph),ncomph): number of particles in each compartment of the domain. Outputs: The same array as for the input but actualized after the random walk has taken place. 2.7. Tests of the program The Monte Carlo code has been extensively checked and its results proved to match with theoretical expectation [3] at long times. We are showing now two explicit tests of the two main routines of the code: the routine walk, which performs the random-walk diffusion of the particles and the routine kine,which is used to integrate the kinetic equations. The validity of the random walk routine (walk) was established by testing that free diffusion of particles resulted in correct Gaussian distributions. As previously described this subroutine, given an input distribution of particles in the domain of simulation n(i, j, k), gives as output the actualized array n(i, j, k) after a 3D random walk (time) step. The test of random walk diffusion is shown in Fig. 2. The test consisted in placing 106 particles initially at the center of a rectangular grid with 40 × 40 × 1000 grid points and evaluating the total number of particles for the different values of the last coordinate (range −500 : 500) after a number of iterations of the routine walk. It is well known [11,12] that the solution of this 1-dimensional diffusion problem is given by: N 2 n(z, t) = √ e−z /4Dt 4πDt
(7)
with N = 106 and, as previously commented t = (z)2 /4D being z and t the distance between points of the grid and the time step corresponding to each iteration of the random walk. Then, after i iterations, we will have a distribution of particles: N 2 n(k, i) = √ e−k / i πi
(8)
with k = −500, −499, . . ., 500. In Fig. 2 we show the simulated distributions for i = 500, 1000, 2000. We can observe that the similarity with the Gaussian distribution is more accurate as i is taken larger.
280
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
Fig. 3. Evolution of free calcium for a single buffering system (Fura-2) obtained from the simulation and the analytical solution (Eq. (9)). There is no bound buffer in the system at t = 0. (A) BT ≡ [Fura − 2] = 100 µM, [Ca](0) = CT = 81.64 µM; (B) BT ≡ [Fura − 2] = 200 µM, [Ca](0) = CT = 195.5 µM.
The routine which implements the first-order kinetic reactions (kine) was tested by checking that starting from concentrations far from equilibrium, equilibrium was reestablished and maintained. By using the kinetic equations, it can be easily proved that the evolution of the free calcium for a single buffering system and non-diffusible species is given by [Ca](t) = where λ1/2 =
C1 ± AC2 e−k+ λ 1 ∓ Ae−k+ λ
1/2 t
1/2 t
,
(9)
2 + 2[B ]K + 2K [C ] − 2[B ][C ], [BT ]2 + [CT ]2 + [C]KD T D D T T T
λ1/2 ∓ ([BT ] − [CT ] + KD ) , C1,2 = 2 [Ca](0) − C1 . A = [Ca](0) + C2 The upper sign in Eq. (9) corresponds to the case [Ca](0) > C1 while the lower sign applies when [Ca](0) < C1 . [Ca](0), [BT ], [CT ], KD and k+ are the free calcium concentration at t = 0, the total buffer concentration (free+bound), the total calcium concentration (free+bound), the dissociation constant (KD = k− /k+ ) and the forward binding rate, respectively. C1 = limt →∞ [Ca](t) is the [Ca] concentration at equilibrium. We checked that the routine kine reproduces the dependence of Eq. (9) for two different buffering systems: Fura-2 (k+ = 5 × 108 M−1 s−1 , KD = 0.24 µM) and EGTA (k+ = 1 × 107 M−1 s−1 , KD = 0.15 µM). In this test, the buffers are considered to be immobile. Figs. 3 and 4 show the solution (9) for the case in which [Ca](0) > C1 ; in particular, we considered that there is no bound buffer at t = 0. In Fig. 3 we show the plot for Eq. (9) together with the simulated result for two different concentrations of Fura-2. In Fig. 4 we repeat the test but for EGTA as buffer. The simulated results show the concentration time course in the first layer of a conical domain of 1 µm radius and 5 µm deep (with initial uniform concentrations).
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
281
Fig. 4. Evolution of the free calcium concentration for a single immobile buffering system (EGTA). The simulated results are compared with the analytical solution (Eq. (9), upper signs). There is no bound buffer at t = 0. (A) BT ≡ [EGTA] = 1000 µM, [Ca](0) = CT = 870.56 µM; (B) BT ≡ [EGTA] = 500 µM, [Ca](0) ≡ CT = 490.44 µM.
Fig. 5. Temporal evolution of free calcium for a single buffering system obtained from the simulation and comparison with the analytical solution (Eq. (9), lower signs). There is no free buffer at t = 0 ([Ca](0) = 0). (A) CT = BT = [Fura − 2] = 100 µM; (B) CT = BT = [EGTA] = 500 µM.
Fig. 5 compares Eq. (9) with the simulated results in the case [Ca](0) < C1 . We consider that initially the buffer is completely bound and that [Ca](0) = 0. As can be seen in the figures, the agreement between the simulated time course and the analytical solution is very good in all the cases considered.
282
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293 Table 1 CPU times spent in a Pentium III 450 MHz for the simulation of a calcium time course lasting 0.5 ms. The grid resolution is 100 nm Buffer
k+ (M−1 s−1 )
KD (µM)
DBuffer (µm2 /s)
[Buffer] (µM)
CPU (s)
Fura-2
5 × 108
0.24
42
100
17.11
Fura-2
5 × 108
0.24
42
200
18.75
EGTA
1 × 107
0.15
200
100
34.38
EGTA
1 × 107
0.15
200
200
36.92
BAPTA
5 × 108
0.2
200
100
35.26
BAPTA
5 × 108
0.2
200
200
39.17
Table 2 CPU times spent in a Pentium III 450 MHz for the simulation of a calcium time course lasting 0.5 ms. The grid resolution is 50 nm Buffer
k+ (M−1 s−1 )
KD (µM)
DBuffer (µm2 /s)
[Buffer] (µM)
CPU (s)
Fura-2
5 × 108
0.24
42
100
189.27
Fura-2
5 × 108
0.24
42
200
215.85
EGTA
1 × 107
0.15
200
100
303.83
EGTA
1 × 107
0.15
200
200
426.43
BAPTA
5 × 108
0.2
200
100
291.90
BAPTA
5 × 108
0.2
200
200
413.68
3. CPU times The CPU time spent by the code depends considerably on the parameters of the simulation. In Tables 1, 2 we show the CPU times spent in a Pentium III of 450 MHz (256 MB of RAM) for simulating a calcium time course 0.5 ms long. Results are shown for two different spatial resolutions. We consider 3 exogenous buffers commonly used in fluorescence experiments (Fura-2, EGTA and BAPTA) together with a “standard” immobile endogenous buffer (500 µM, k+ = 5 × 108 M−1 s−1 , KD = 10 µM [2]). The simulation domain is conical, 5 µm deep and 1 µm base radius. The values of the kinetics constants and diffusion coefficients for each buffer are taken from reference [2]. The unitary current is 0.1 pA and the number of channels considered is 30.
Acknowledgements The authors acknowledge the use of the routine BTPEC (algorithm 678 collected from ACM) and function RAND1 [ACM Transactions on Mathematical Software 5 (1979) 132–138], both obtained from NETLIB repository. The authors also acknowledge financial support from Generalitat Valenciana (GV99-146-1-04) and Fundació La Marató de TV3.
References [1] F. Sala, A. Hernández-Cruz, Biophys. J. 57 (2) (1990) 313–324. [2] J. Klingauf, E. Neher, Biophys. J. 72 (2) (1997) 674–690.
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293 [3] [4] [5] [6] [7] [8] [9] [10] [11] [12]
A. Gil, J. Segura, J.A.G. Pertusa, B. Soria, Biophys. J. 78 (1) (2000) 13–33. M. Oheim, D. Loerke, W. Stühmer, R.H. Chow, Eur. Biophys. J. 28 (2) (1999) 91–101. T. Voets, E. Neher, T. Moser, Neuron 23 (3) (1999) 607–615. J. Segura, A. Gil, B. Soria, Biophys. J. 79 (4) (2000) 1771–1786. N.L. Allbritton, T. Meyer, L. Stryer, Science 258 (1992) 1812–1815. V. Kachitvichyanukul, B. Schmeiser, Comm. ACM 31 (2) (1988) 216–222. http://www.netlib.org. L. Schrage, ACM Trans. Math. Software 5 (1979) 132–138. R.D. Present, Kinetic Theory of Gases, McGraw-Hill, 1958, pp. 65-69. J. Crank, The Mathematics of Diffusion, Oxford at Clarendon Press, 1956, pp. 9-11.
283
284
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
TEST INPUT DATA .false. 1 1. 5. 30 1 0.1 1 0.3 0. 0.3 0.2 220. 0. 42 0.1 500. 100. 5.e+8 10. 5.e+8 0.24 3 0 1 0.1
!logical flag for timing (default .false.). See details in ca3d.f !Type of domain. 1 = Conical. 2 = Cylindrical. !Radius of the base of the domain (in microns). !Height of the domain (in microns). !Number of entry calcium channels. !Distribution of channels. 1 = Random. 2 = Specified by the user. !Spatial resolution in microns !Number of calcium pulses. !Duration of each pulse (in ms). !Frequency at which the pulses are delivered (in Hz). !Duration of the simulation (in ms). !Unitary current per channel (in pA). !Diffusion coefficient of ions (in micron**2/s). !Diffusion coefficient of buffer 1 (in micron**2/s). !Diffusion coefficient of buffer 2 (in micron**2/s). !Basal concentration of calcium (in microMol). !Concentration of buffer 1 (in microMol). !Concentration of buffer 2 (in microMol). !Forward binding rate of the endogenous buffer (in Mol**-1s**-1). !Dissociation constant of the endogenous buffer (in microMol). !Forward binding rate of the exogenous buffer (in Mol**-1s**-1). !Dissociation constant of the exogenous buffer (in microMol) !Number of slices at which the averaged [Ca2+] is required. !Are the buffer concentrations also required?. 0-->no, 1-->yes. !Number of 2D calcium maps to be calculated. !time at which the first calcium map is obtained (in ms). !time at which the second calcium map is obtained (in ms). !time at which the third calcium map is obtained (in ms). !....... !time at which the 25th calcium map is obtained (in ms).
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
TEST OUTPUT FILE: ‘casli01’ #********************************* # Averaged calcium concentration # at distances from 0.nm # to 105.nm to the # membrane #********************************* # # time (ms) Calcium (microMol) # 0.013 0.253 0.025 0.281 0.038 0.349 0.050 0.355 0.063 0.332 0.076 0.338 0.088 0.304 0.101 0.332 0.113 0.422 0.126 0.360 0.139 0.400 0.151 0.343 0.164 0.422 0.176 0.388 0.189 0.439 0.201 0.428 0.214 0.445 0.227 0.467 0.239 0.495 0.252 0.450 0.264 0.658 0.277 0.580 0.290 0.456
285
286
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
FILE: ‘casli02’ #********************************* # Averaged calcium concentration # at distances from 105.nm # to 210.nm to the # membrane #********************************* # # time (ms) Calcium (microMol) # 0.013 0.158 0.025 0.191 0.038 0.146 0.050 0.180 0.063 0.174 0.076 0.253 0.088 0.259 0.101 0.281 0.113 0.135 0.126 0.248 0.139 0.186 0.151 0.214 0.164 0.180 0.176 0.264 0.189 0.225 0.201 0.248 0.214 0.276 0.227 0.191 0.239 0.264 0.252 0.197 0.264 0.253 0.277 0.281 0.290 0.281
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
FILE: ‘casli03’ #********************************* # Averaged calcium concentration # at distances from 210.nm # to 315.nm to the # membrane #********************************* # # time (ms) Calcium (microMol) # 0.013 0.186 0.025 0.118 0.038 0.118 0.050 0.141 0.063 0.118 0.076 0.129 0.088 0.107 0.101 0.135 0.113 0.113 0.126 0.096 0.139 0.152 0.151 0.118 0.164 0.152 0.176 0.118 0.189 0.090 0.201 0.073 0.214 0.124 0.227 0.186 0.239 0.203 0.252 0.135 0.264 0.141 0.277 0.113 0.290 0.101
287
288
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
FILE: ‘camap01’ #***************** # 2D calcium map #***************** # # at time 0.100 ms # # # x y Calcium (microMol) # 0 -9 0.000 -4 -8 2.278 -3 -8 0.854 -2 -8 0.142 -1 -8 0.285 0 -8 0.142 1 -8 0.285 2 -8 0.000 3 -8 0.000 4 -8 0.285 -5 -7 0.285 -4 -7 0.285 -3 -7 0.142 -2 -7 0.427 -1 -7 0.712 0 -7 0.285 1 -7 0.285 2 -7 0.569 3 -7 0.854 4 -7 0.142 5 -7 0.142 -6 -6 0.285 -5 -6 0.000 -4 -6 0.142 -3 -6 0.142 -2 -6 0.285 -1 -6 0.142 0 -6 0.285 1 -6 0.142 2 -6 0.285 3 -6 0.427 4 -6 0.142 5 -6 0.142 6 -6 0.000 -7 -5 2.136 -6 -5 0.712 -5 -5 0.000 -4 -5 0.142 -3 -5 0.000 -2 -5 0.427 -1 -5 0.000
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
0 1 2 3 4 5 6 7 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 -8 -7 -6 -5 -4 -3 -2 -1 0 1
-5 -5 -5 -5 -5 -5 -5 -5 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -3 -3 -3 -3 -3 -3 -3 -3 -3 -3 -3 -3 -3 -3 -3 -3 -3 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2
0.000 0.285 0.142 0.142 1.139 0.712 0.000 0.000 0.285 0.285 0.142 0.142 0.142 0.000 0.285 0.142 0.000 0.427 0.142 0.569 0.142 1.851 0.285 0.000 0.142 0.142 0.142 0.142 0.142 0.142 0.000 0.000 0.142 0.427 0.569 1.139 0.997 2.136 0.427 0.569 0.142 0.000 0.285 0.000 0.142 0.285 0.142 0.142 0.427 0.285 1.566 1.566
289
290
2 3 4 5 6 7 8 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 -8 -7 -6 -5 -4 -3 -2 -1 0
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
-2 -2 -2 -2 -2 -2 -2 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1
0.997 0.712 0.569 0.285 0.000 0.000 0.142 0.000 0.000 0.142 0.142 0.285 0.142 1.424 1.281 1.139 0.000 0.142 0.000 0.142 0.142 0.285 0.000 0.142 0.142 0.000 0.000 0.142 0.142 1.139 0.712 0.569 0.142 0.142 0.427 0.142 0.000 0.000 0.142 0.285 0.427 0.427 0.427 0.142 0.000 0.000 0.000 0.285 0.000 0.997 0.285 0.285
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
1 2 3 4 5 6 7 8 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 -8 -7 -6 -5 -4 -3 -2 -1 0 1
1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4
0.000 0.000 0.427 0.000 0.142 0.142 0.142 1.139 0.000 0.000 0.142 0.142 0.142 0.000 0.285 0.285 0.000 0.142 0.427 0.142 0.285 0.427 0.142 0.427 0.427 0.569 0.854 0.000 0.142 0.000 0.142 0.000 0.000 0.142 0.000 0.000 0.427 0.427 1.708 0.569 1.993 0.854 2.420 0.142 0.285 0.000 0.285 0.000 0.569 0.427 0.000 0.000
291
292
2 3 4 5 6 7 8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 -5 -4 -3 -2 -1 0 1 2 3 4 5 -4 -3 -2 -1 0 1
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 7 7 7 7 7 7 7 7 7 7 7 8 8 8 8 8 8
0.142 1.424 0.142 0.427 0.142 1.424 0.427 0.000 0.142 0.142 0.997 2.278 1.424 0.854 0.569 0.142 0.569 0.142 0.285 0.142 0.712 0.142 0.000 0.285 1.993 1.566 0.569 0.285 0.569 1.708 0.712 0.285 0.427 0.142 0.285 0.142 0.427 0.712 0.569 0.142 0.285 0.569 0.427 0.997 0.285 0.285 0.142 0.569 1.708 0.569 0.000 0.142
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
2 3 4 0
8 8 8 9
0.427 0.569 0.000 0.000
293
CA3D: A Monte Carlo code to simulate 3D buffered diffusion of ions in sub-membrane domains Amparo Gil ∗ , Javier Segura 1 Instituto de Bioingeniería, Universidad Miguel Hernández, 03202 Elche, Alicante, Spain Received 26 June 2000; accepted 21 December 2000
Abstract CA3D is a simulation program, based on Monte Carlo techniques, which simulates the buffered diffusion of ions in the submembrane domain of cells. Two geometries of submembrane domains are considered. The first, a conical domain, applies to spherical cells with a locally uniform distribution of channels while the second, a cylindrical domain, is intended to model ciliar cells and presynaptic terminals. Other domains can be implemented in the code with ease. The entry of ions takes place through discrete entry points (channel pores). Ions diffuse inside the cell and react with endogenous or exogenous buffers, mobile or fixed. The diffusion of ions and binding molecules is modeled as a 3D random walk process by moving each individual ion and molecule. The reaction of ions with the binding molecules is also modeled probabilistically. The outputs of the code are concentration time courses and 2D maps of concentration of ions in the sub-membrane domain of cells. 2001 Published by Elsevier Science B.V. PACS: 87.10.+e; 87.22.As; 02.70.Lq Keywords: Ion diffusion; Random walk; Buffering; Reaction-diffusion
PROGRAM SUMMARY
Computers for which the program is designed an others on which it is operable: SUN Enterprise 3000; PC Pentium III 450 MHz with 128 MB RAM; PC Pentium II 350 MHz with 64 MB RAM
Title of program: CA3D Catalogue identifier: ADNW Program Summary URL: http://cpc.cs.qub.ac.uk/summaries/ADNW Program obtainable from: CPC Program Library, Queen’s University of Belfast, N. Ireland
Operating systems under which the program has been tested: Solaris (SunOS 5.6), GNU/Linux, Win32 Programming language used: FORTRAN 77 Compilers: f77 version SC4.2 for SunOS, g77 v0.5.24 for Linux, g77 v0.2.95 for MinGW32, Compaq Visual Fortran for Windows No. of bytes in distributed program, including test data, etc.: 15 357
* Corresponding author.
E-mail address: [email protected] (A. Gil). 1 Present address: Departamento de Matemáticas, Universidad Carlos III de Madrid, Leganés Madrid, Spain.
0010-4655/01/$ – see front matter 2001 Published by Elsevier Science B.V. PII: S 0 0 1 0 - 4 6 5 5 ( 0 1 ) 0 0 1 4 4 - 8
270
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
Distribution format: tar gzip file Supplementary material: Pseudocolored animated images of calcium 2D maps obtained from the simulation Keywords: Ion diffusion, random walk, buffering, reaction-diffusion Nature of the physical problem The buffered diffusion of ions in the sub-membrane domain of cells is simulated by means of Monte Carlo methods. Two geometries of submembrane domains are considered. The first, a conical domain, applies to spherical cells with a locally uniform distribution of channels while the second, a cylindrical domain, is intended to model ciliar cells and presynaptic terminals. Other domains can be implemented in the code with ease. The entry of ions takes place through discrete entry points (channel pores) which are distributed over the compartments of the first slice of the domain (base of the cone or the cylinder). Ions diffuse inside the cell and react with endogenous or exogenous buffers, mobile or fixed.
Method of solution The diffusion of ions and binding molecules is modeled as a 3D random walk process by moving each individual ion and molecule. The reaction of ions with the binding molecules is also modeled probabilistically.
Restrictions on the complexity of the problem The number of radial grid points is restricted to 100 (diameter 200 points). The depth is limited to 100 grid points. In case the grid resolution and the sizes of the domain are chosen in such a way that this limit is surpassed, a warning message is generated and the program gives the option of continuing the simulation by taking the largest possible sizes of the domain, maintaining the spatial resolution. The number of grid points available could be smaller in other (generally older) computer systems, depending on the RAM available. In all the systems checked (with a RAM of 64 MB or higher) the restrictions on the number of grid points described above were always within computable ranges. The dimension of the arrays limit the resolution and size of the domain of simulation that can be taken into account. In this way, for domains with a typical size of L, the best resolution achievable is of the order of L/100. Typical running time Depends on the buffers, grid resolution, unitary current, etc. For example, for a conical domain of 1 µm radius and 5 µm deep, a calcium time course lasting 1 ms takes 34 s of CPU time for a resolution of 100 nm and 378 s for a resolution of 50 nm; in both cases a standard endogenous buffering is considered with the addition of an exogenous buffer (Fura-2 100 µM). See section “CPU times” for more details.
1. Introduction The diffusion of different types of ions in the cytosol of living cells play a key role in their function. For instance, it is well known that the fluxes of Na2+ and K+ across the membranes of excitable cells are responsible for the propagation of action potentials, described by the celebrated Hodgkin–Huxley cable equations. On the other hand Ca2+ is known to act as a second messenger which triggers a number of cellular functions such as: secretion of hormones (like, for instance, insulin), muscle contraction, release of neurotransmitter and others. In particular, much attention has been paid to the diffusion of Ca2+ in intracellular media and there is a vast number of biological publications which rely on differential methods (see for instance [1,2]). Ca2+ diffuses inside the cytosol while at the same time first order kinetic reactions take place by which the Ca2+ ions can bind with a binding site of a given buffer B to give rise to a bound state BCa: Ca2+ + B BCa.
(1)
The diffusion of calcium in the cytosol is also referred to as calcium buffered diffusion. Previous models to describe buffered diffusion were based on differential equations, which need drastic geometrical simplifications in order to be solved with ease. Depending on the geometrical assumptions considered the differential models can lie within two categories: shell models (see for instance [1]), in which it is assumed that ions enter uniformly over the cell membrane and only the radial diffusion is of relevance (for spherical cells) and microdomain models (see [2]), which solve the reaction diffusion equations in the vicinity of a channel pore, assuming azimuthal symmetry around the pore; in the latter case it is necessary to assume that all channel pores are equidistant. However, it is well known that ions do not enter uniformly given that the channel pores are not regularly spaced. These departures from symmetry may have an important impact on the functionality of the cells. In [3] it was described how the clustering of channels gives rise to calcium profiles which can not be obtained as a
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
271
superposition of calcium concentrations developed by each individual channel. This non-linearity amplifies the effect of irregular distributions of calcium channels on the formation of submembrane calcium concentration profiles. As an example of the physiological relevance of such non-uniform configurations we can mention the case of secretion in neuroendocrine cells, which is a process triggered by calcium. There are indications that the vesicles responsible for secretion seem to be non-uniformly distributed [2,4,5]; this fact, together with the complexity of the calcium concentration profiles below the membrane clearly points towards the necessity of a modeling which is free of strong simplifying assumptions regarding the geometry of the system [6]. With the idea in mind to overcome the geometrical simplifications imposed by differential methods, we developed a Monte Carlo simulation for the problem of 3D buffered calcium diffusion [3]. This scheme has proved to be successful in the study of the influence of the geometry in the exocytotic response of neuroendocrine cells [6]. For the spatial resolutions that are relevant (of the order of 100 nm in neuroendocrine cells) and the typical concentrations of calcium (order 10 µM) and calcium buffers (order 100 µM) we can expect few calcium ions inside each cube of 100 nm of side (∼7 ions for 10 µM) and also a moderate number of calcium binding sites. Of course, as we demand a higher spatial resolution, less ions and buffer molecules are expected. Then, it is more appropriate to keep track of each individual ion and molecule instead of using a continuous modeling in terms of concentrations. The code CA3D is a microscopic simulation in which the fundamental variables are the number of ions and buffers; in contrast, the differential methods are macroscopic and deal with concentrations. The average values of the output of our simulation converge to macroscopic simulations when considering symmetrical configurations. The outputs of the code are transformed back into macroscopic variables (concentrations). The program gives as output ion and buffer concentration time courses at different depths as well as 2D maps of concentrations of ions in the submembrane domain. Pseudocolored images of calcium 2D maps obtained from the simulation can be found in http://bessel.umh.es. 2. Description of the method 2.1. Geometrical configuration We consider two basic configurations: a cylindrical volume and a conical volume. As described in Ref. [3], the conical domain is appropriate to describe buffered diffusion in the submembrane domain of spherical cells. The base of the cone represents the membrane of the cell. Assuming that the patch of membrane considered is embedded inside a larger portion of the membrane where calcium channels are, on average, uniformly distributed, we can assume that the flux of ions and buffer molecules through the lateral sides of the cone is zero. This simplification imposes a simple boundary condition in the walls of the membrane: ions bounce back when they encounter a wall. At the base (membrane) this is interpreted as reflection while, on the lateral sides, reflection is equivalent to considering that the outgoing flux equals the incoming flux, due to the symmetry of the problem. An orthogonal 3D regular grid maps the domain of simulation with a distance between grid points chosen by the user. Typical operative ranges for this spatial resolution (l) lie between 10 and 100 nm. Each point of the grid is associated with a cubic compartment of volume (l)3 . Therefore, grid points or cubic compartments are equivalent concepts. Ions and buffer molecules move from one compartment to another compartment of the grid due to diffusion, which is modeled as a random walk process. The time scale for the random walk is given by the spatial resolution (which can be selected by the user) and by the higher diffusion coefficient of the system, which is supposed to be the diffusion coefficient for the ions. The relation between those scales is given by: t = (l)2 /4D
(2)
being l the spatial resolution and D the diffusion coefficient for the ions. As an example, for a spatial resolution of l = 100 nm and for calcium ions (DCa2+ = 220 µm2 /s in cytoplasm [7]), Eq. (2) gives a time resolution
272
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
Fig. 1. Flux diagram of the program.
t = 12.5 µs. The time resolution is the time step of the simulation: ions will be moved every time step t using a random walk algorithm in which each ion will have a 50% chance to remain in its initial position and a 25% chance to move either in the positive or negative direction for each spatial dimension. The cylindrical domain fits to study Ca2+ signals in ciliar cells or in presynaptic terminals. We also impose total reflection on the boundaries of the cylinder. 2.2. Sequence of the simulations The sequence of the simulations follows the flux diagram shown in Fig. 1. Before entering the simulation loop, the main program calls the subroutines makegrid, data, districa and chandis for setting a discretization of the simulation domain, transforming the macroscopic input variables into microscopic quantities, setting the initial conditions and distributing the channel pores over the upper slice of the simulation domain. By upper slice we mean the submembrane domain (base of the cone or the cylinder). After the initial data have been processed and the simulation domain and initial conditions have been built, the program enters in the simulations loop. The main ingredients of this loop are the routines walk and kine, combined with the entry of ions through the channels. This last task is performed by the main program which also handles the outputs. 2.3. Input data processing First, the orthogonal grid is built by the subroutine makegrid which, as previously commented, maps the simulation domain to a 3D orthogonal grid. The resolution demanded is one of the input parameters of the program together with the size of the domain. After building the grid the resolution is slightly adjusted in order to fit the grid with the dimensions of the domain. This small correction (smaller as the chosen resolution is better) reduces the possible underestimation of the membrane surface by the rounding in the process of making the grid. The output files of the program show the actual resolution used. After building the grid, the program calls the subroutine data in order to transform the macroscopic input data (concentrations, kinetic rates, . . .) into microscopic data (number of molecules in each compartment, reaction probabilities per molecule, . . .). The subroutine data also evaluates the initial conditions for buffers and ions; it is assumed that ions and buffers are initially at equilibrium. Among other parameters, data evaluates the total number of free and bound buffers
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
273
from the equilibrium condition taking as input the total concentrations of buffers as well as the basal (initial) concentration of free ions. After data conversion by the routine data, districa distributes uniformly the initial (basal) ion concentrations and the concentration of buffers (free and bound) over the whole simulation domain. Regarding the distributions of the entry points for the ions (ionic channels), when the user selects a randomly generated distribution, the program calls the routine chandis to distribute the number of channels specified over the upper layer of the domain (representing the cell membrane). The user has also the option of determining the position of the channels in the file “channel.dat” by giving the grid coordinates of the channels in the upper layer of the grid. 2.4. Simulation loop In the simulation loop three consecutive actions are repeated for each simulation step t: (1) During the time intervals for which there is an influx of ions (as determined in “input.dat”) , the main program sets, for each time step, a number of extra ions in the submembrane compartments corresponding to the position of a channel. (2) Ions and buffer molecules suffer a random walk step. (3) The kinetic equations for the binding of ions by buffers evolves during the time lapse t. 2.4.1. Influx of ions The influx of ions is controlled by the main program. By default, the program simulates trains of pulses with constant incoming current. The user should determine in the file “input.dat” the number and duration of the pulses, the frequency at which they are delivered as well as the (constant) magnitude of the current through each channel pore. Then, using the microscopic current (number of incoming ions in each t) evaluated in data 2 the program sets the extra number of ions “below” each channel in each t. Given that the number of incoming ions in each t could be smaller than 1, the program chooses probabilistically the number of ions which enter through each channel to avoid rounding errors. The user has also the possibility of choosing different shapes of the unitary current. This has to be specified in the following piece of the main program: cccccccccccccccccccccccccccccccccc c For non constant pulses enter c c here your choice of function c c for the unitary current, c c scaling the currents and the c c time variable with the factors c c scei and scet (see below). c c c c 1) SET fnio= YOUR FUNCTION c c c c Here is an example for a c c unitary current cccccccccccc c I(t)=0.213exp(-t/20) pA (t in ms) c c lasting 50 ms c c c 2 Note: the program is prepared to simulate calcium buffered diffusion. For other ions, the parameter charge in data should be modified to
account for the number of charges per ion.
274
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
c if(i.lt.(50*scet))then c fnio=0.213*scei*exp(-i/(20*scet)) c else c fnio=0 c endif c c c (the time variable is i) cccccccccccc c c c 2) UNCOMMENT THE LINES: c c c c npulse=10000000 c ifnio=int(fnio) c fnio=fnio-ifnio cccccccccccccccccccccccccccccccccc c WARNING: you can not simulate c c trains of pulses unless you c c explicitly specify this fact c c in your function c cccccccccccccccccccccccccccccccccc 2.4.2. Subroutine walk: diffusion The 3D diffusion of ions and mobile buffers is modeled as a random walk process. After each time step t = (l)2 /4D and for each dimension, a particle has a 50% chance of remaining in its initial position and a 25% chance of moving either in the positive or the negative direction. Another possible prescription to implement the random walk is considering a time step t = (l)2 /2D and assigning 50% probability of movement in either direction for each spatial dimension. 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. Depending on the mobile buffer considered (EGTA, Fura-2, . . .), the corresponding diffusion coefficient can be typically about 2 to 5 times smaller than that for Ca2+ so that the characteristic diffusion times are about 2 to 5 times bigger. Thus, the diffusion of buffer particles is simulated by moving them in exactly the same way that ions are moved but only after 2 to 5 simulation steps. This means that the program takes integer ratios between the largest diffusion coefficient (ions) and the smallest (buffers). Given the experimental uncertainties associated with diffusion coefficients in the cytoplasm this approximation is accurate enough and saves a considerable amount of computational time. An additional way of saving CPU time is by considering “bulk” diffusion when the number of particles N to diffuse in a given compartment is greater than 64 [3]. In such a case, we take N = 64c + R, with R < 64 (c and R integer numbers). Then, we consider the “bulk” diffusion of the 64c buffer particles according to the probabilities of moving in each direction of the 3D grid [3] and we decide individually and probabilistically the position of the remaining R particles. It is preferable to consider the 3-point implementation of the random walk (with t = (l)2 /4D) to make the bulk diffusion smoother and without large local calcium fluctuations. Regarding the reflection in the boundaries of the sub-domain we take this effect into account by rejecting those events which result in the positioning of a particle outside the sub-domain; 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 sub-domain). 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 being equivalent to considering that the number of outgoing particles equals the number of incoming particles.
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
275
2.4.3. Subroutine kine: kinetics During each diffusional time step t and in each compartment of the grid, the kinetic equations which 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 d[Ca] i i = k− [CaBi ] − k+ [Ca][Bi ] , dt i (3) d[Bi ] d[Ca] d[CaBi ] d[Bi ] = ; = − , dt dt dt dt i
i k+
i , the unbinding rate. being the forward binding rate for the buffer i and k− Since 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 specie in each small cubic compartment of the grid. In terms of the number of free calcium ions NCa , free buffer molecules NBi and bound buffer molecules NCaBi , Eq. (3) reads: i i NCa = K− NCaBi − K+ NCa NBi , i (4) + N NCaBi + NCa = NCaT , N Bi = NBT ; CaBi i i
where NBT represents the total (free + bound) number of buffer particles of type i, NCaT is the total number of (cali cium) ions and NCa represents the variation of free calcium ions in a given cubic compartment during a time interi and Ki for unbinding and binding in each time interval dt are given by val dt ≡ t/n. The effective constants K− + t i i t ; K+ , (5) = k+ n nV where n is the number of sub-divisions of the interval t used to integrate the kinetic equations and V is the volume of each cubic compartment. In order to implement a microscopic interpretation we describe the system probabilistically, evaluating the probabilities associated with each of the processes and generating events according to such probabilities [3]. During the time dt, the unitary probabilities (probabilities for the binding of a single ion or unbinding of a single bound i N for binding and P i = Ki for unbinding, provided we have taken a sufficiently buffer molecule) are Pbi = K+ Bi − u small dt (n large enough). n is chosen in such a way that all probabilities are small enough. Numerical experiments show that taking 0.1 as the maximal probability (for binding or unbinding) is an accurate selection. The number of calcium ions that become bound in each interval dt or become unbound are evaluated using the corresponding probabilities and following a binomial distribution N (6) P (k; p, N) = pk (1 − p)N−k , k i i = k− K−
where P (k; p, N) is 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 molecules) in a time dt, p being the unitary probability (Pbi or Pui ). The code BTPEC [8], obtainable from NETLIB repository [9], is used to generate binomial random deviates. This program uses the uniform random number generator RAND1 [10], which is also used in the rest of the simulation program. The algorithm for kinetics here described gives priority to the fastest effects during each t. For instance, unbinding is normally so slow that the probabilistic selection of the number of buffer molecules that undergo unbinding can be done 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.
276
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
2.5. Description of the inputs, outputs and parameters Inputs: The input data for the program must be written in the file ‘input.dat’. Their meaning is the following: • icondit=logical flag which controls the timing option (uses the non-standard function secnds). The value by default is .false. and the program lines for timing appear as commentaries. To use this function set icondit to the value .true. and uncomment the corresponding lines in the program (search for the word “TIMING” in the program to locate these lines). This function will estimate the duration of a run of the program by counting the number of iterations of the simulation loop during the first 10 s of CPU time and extrapolating the time for the whole run. • imicro=kind of microdomain. If imicro=1, the selected microdomain is conical while imicro=2 corresponds to a cylindrical domain. • rbas=radius of the base of the microdomain in µm. • rdepth=depth of the microdomain in µm. • ncha=number of ionic channels. • idncha=distribution of ionic channels. If idncha=1, the selected distribution is random and uniform. If idncha=2, the distribution of channels can be specified by the user: the program reads the coordinates of the entry channels from the file ‘channel.dat’. • deltax=spatial resolution in µm. • npulses=number of calcium pulses (or ionic pulses in general). • tpulse=duration of each pulse in ms (assumed equal for all the pulses). • fhz=frequency at which the pulses are delivered. • tsimul=duration of the simulation in ms. • di0=unitary current in pA (current of ions through each channel). • difca=diffusion coefficient of ions in the cytosol in µm2 /s. • difbf=diffusion coefficient of buffer 1 in µm2 /s. • difbm=diffusion coefficient of buffer 2 in µm2 /s. • caini=basal (initial) concentration of ions in µM. • bfixed=concentration of buffer 1 in µM. • bmobil=concentration of buffer 2 in µM. 1 ) in M−1 s−1 . • dkbf=forward binding rate of buffer 1 (k+ 1 = k 1 /k 1 ) in µM. • dkdf=dissociation constant of buffer 1 (KD − + 2 • dkbm=forward binding rate of buffer 2 (k+ ) in M−1 s−1 . 2 = k 2 /k 2 ) in µM. • dkdm=dissociation constant of buffer 2 (KD − + • nslice=number of layers of the grid (at progressive depths) where the averaged [Ca2+ ] is required. Maximum: 9. • ibuf=choice for averaged buffer concentrations. If ibuf=0(1), the averaged buffer concentrations at the nslice layers are not (are) calculated. • nmaps=number of 2D ion concentration maps to be calculated. Maximum: 25. • tmap(1:nmaps)=times at which an ion concentration map is obtained in ms (include as many lines in the input file as concentration times demanded) Outputs: The outputs provided by the program are averaged ion concentrations as a function of time (for each of the first nslice layers of the microdomain) or 2-dimensional ion concentration maps (with a temporal bining given as a parameter). The output files are the following: • Averaged ion concentrations in the first nslice layers: files “caslixx”, where xx= 01, 02, 03, . . . gives the depth of the slice (in grid units, 1 grid unit=resolution).
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
277
• Averaged free and bound buffer concentrations in the first nslice layers (for ibuf=1): files “b1fslxx”, “b1bslxx”,“b2fslxx”,“b2fslxx”,“b2bslxx” with xx= 01, 02, . . .. • 2D ion concentration maps: files “camapxx”, xx= 01, 02, . . . . Parameters: • ibin=number of temporal points for the bining of the 2D calcium maps (default: 10 points). 2.6. Description of the subroutines and dependencies The main program calls the subroutines: makegrid, data, chandis, districa, walk, kine, with a calling sequence as described in Fig. 1. We now describe the purpose of these subroutines as well as their dependencies: makegrid(imicro,ncomps,ncomph,jl,il,numcom,nt) Purpose: Routine to make the grid for a conical or a cylindrical domain. Arguments: Inputs: imicro: selects the geometry of the domain. If imicro=1, the selected domain is conical while imicro=2 stands for a cylindrical domain. ncomps: number of grid points along the radius of the base of the domain. ncomph: number of grid points along the depth of the domain. Outputs: jl(1:ncomph), il(1:ncomph,jl): range of the coordinates of the grid points as a function of depth. numcom(1:ncomph): number of compartments as a function of depth. nt: total number of compartments of the domain. This routine is called only once and at the start of the simulation. After the call to the routine makegrid the domain of simulation becomes defined by the following set of grid points: {{{i |il(k, j )|}, j |j l(k)|}, k = 1, 2, . . . , ncomph}, where k measures the depth and i, j are the “horizontal” coordinates. k = 1 represents the sub-membrane domain. data(deltax,difca,difbf,difbm,di0,fhz,tpulse,tsimul,nmaps,tmap) Purpose: transformation of the input data to microscopic quantities used in the simulation and setting of the initial conditions. Arguments: Inputs: deltax: selected resolution for the simulation. difca: diffusion coefficient of Ca2+ . difbf: diffusion coefficient of buffer 1. difbm: diffusion coefficient of buffer 2. di0: unitary current (per channel). fhz: frequency of delivery of ionic pulses. tpulse: duration of each pulse. tsimul: duration of the simulation. nmaps: number of 2D calcium maps. tmap(1:nmap): times at which the 2D calcium maps are obtained. Parameters: charge: number of charges of one ion (default: charge=2). The outputs of this routine, which is called only once by the main program, are passed to the main using common declarations. This routine gives as output the total number of free ions and free and bound buffers
278
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
of both types, together with the time step for the simulation in ms. The routine data also renormalizes the rest of macroscopic constants (kinetic and geometrical constants) transforming them to microscopic quantities. chandis(iseed,ncha,ie,je) Purpose: generation of a uniform and random distribution of channels. Arguments: Inputs: iseed: seed for the simulation. ncha: number of calcium channels to be distributed. Outputs: ie(1:ncha),je(1:ncha): coordinates of the channels (the depth coordinate being k = 1). This routine is called only once and at the start of the simulation. districa(iseed,numcom,inica,nt,ni) Purpose: distributes a given specie of particle (ions, free buffer molecules, . . .) all over the domain of simulation, being all compartments of the grid equally probable. Arguments: Inputs: iseed: seed for the simulation. numcom(1:ncomph): number of compartments of the domain as a function of the depth. inica: number of particles to be distributed. nt: total number of compartments in the domain. Outputs: ni(-il(ncomph,jl):il(ncomph,jl),-jl(ncomph):jl(ncomph),ncomph): initial number of particles in each compartment of the grid. This routine is called only once and at the start of the simulation. kine(iseed,ni,nb1,nb2,nbm1,nbm2) Purpose: implementation of the first order kinetics reactions in each compartment of the domain. Arguments: Inputs: iseed: seed for the simulation. ni(-il(ncomph,jl):il(ncomph,jl),-jl(ncomph):jl(ncomph),ncomph): number of calcium ions in each compartment of the domain. nb1(-il(ncomph,jl):il(ncomph,jl),-jl(ncomph):jl(ncomph),ncomph): number of free buffer molecules of type 1 in each compartment of the domain. nb2(-il(ncomph,jl):il(ncomph,jl),-jl(ncomph):jl(ncomph),ncomph):number of free buffer molecules of type 2 in each compartment of the domain. nbm1(-il(ncomph,jl):il(ncomph,jl),-jl(ncomph):jl(ncomph),ncomph): number of bound buffer molecules of type 1 in each compartment of the domain. nbm2(-il(ncomph,jl):il(ncomph,jl),-jl(ncomph):jl(ncomph),ncomph):number of bound buffer molecules of type 2 in each compartment of the domain. Outputs: The same arrays as for the inputs but actualized after the changes made by the kinetic reactions. walk(iseed,ni) Purpose: implementation of a 3D random walk algorithm for the movement of particles. Arguments:
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
279
Fig. 2. Distributions for a group of 106 particles, initially at k = 0, after iterating 500, 1000 and 2000 times the routine for random-walk diffusion (walk) (from top to bottom). The comparison with Eq. (8) is also shown.
Inputs: iseed: seed for the simulation. ni(-il(ncomph,jl):il(ncomph,jl),-jl(ncomph):jl(ncomph),ncomph): number of particles in each compartment of the domain. Outputs: The same array as for the input but actualized after the random walk has taken place. 2.7. Tests of the program The Monte Carlo code has been extensively checked and its results proved to match with theoretical expectation [3] at long times. We are showing now two explicit tests of the two main routines of the code: the routine walk, which performs the random-walk diffusion of the particles and the routine kine,which is used to integrate the kinetic equations. The validity of the random walk routine (walk) was established by testing that free diffusion of particles resulted in correct Gaussian distributions. As previously described this subroutine, given an input distribution of particles in the domain of simulation n(i, j, k), gives as output the actualized array n(i, j, k) after a 3D random walk (time) step. The test of random walk diffusion is shown in Fig. 2. The test consisted in placing 106 particles initially at the center of a rectangular grid with 40 × 40 × 1000 grid points and evaluating the total number of particles for the different values of the last coordinate (range −500 : 500) after a number of iterations of the routine walk. It is well known [11,12] that the solution of this 1-dimensional diffusion problem is given by: N 2 n(z, t) = √ e−z /4Dt 4πDt
(7)
with N = 106 and, as previously commented t = (z)2 /4D being z and t the distance between points of the grid and the time step corresponding to each iteration of the random walk. Then, after i iterations, we will have a distribution of particles: N 2 n(k, i) = √ e−k / i πi
(8)
with k = −500, −499, . . ., 500. In Fig. 2 we show the simulated distributions for i = 500, 1000, 2000. We can observe that the similarity with the Gaussian distribution is more accurate as i is taken larger.
280
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
Fig. 3. Evolution of free calcium for a single buffering system (Fura-2) obtained from the simulation and the analytical solution (Eq. (9)). There is no bound buffer in the system at t = 0. (A) BT ≡ [Fura − 2] = 100 µM, [Ca](0) = CT = 81.64 µM; (B) BT ≡ [Fura − 2] = 200 µM, [Ca](0) = CT = 195.5 µM.
The routine which implements the first-order kinetic reactions (kine) was tested by checking that starting from concentrations far from equilibrium, equilibrium was reestablished and maintained. By using the kinetic equations, it can be easily proved that the evolution of the free calcium for a single buffering system and non-diffusible species is given by [Ca](t) = where λ1/2 =
C1 ± AC2 e−k+ λ 1 ∓ Ae−k+ λ
1/2 t
1/2 t
,
(9)
2 + 2[B ]K + 2K [C ] − 2[B ][C ], [BT ]2 + [CT ]2 + [C]KD T D D T T T
λ1/2 ∓ ([BT ] − [CT ] + KD ) , C1,2 = 2 [Ca](0) − C1 . A = [Ca](0) + C2 The upper sign in Eq. (9) corresponds to the case [Ca](0) > C1 while the lower sign applies when [Ca](0) < C1 . [Ca](0), [BT ], [CT ], KD and k+ are the free calcium concentration at t = 0, the total buffer concentration (free+bound), the total calcium concentration (free+bound), the dissociation constant (KD = k− /k+ ) and the forward binding rate, respectively. C1 = limt →∞ [Ca](t) is the [Ca] concentration at equilibrium. We checked that the routine kine reproduces the dependence of Eq. (9) for two different buffering systems: Fura-2 (k+ = 5 × 108 M−1 s−1 , KD = 0.24 µM) and EGTA (k+ = 1 × 107 M−1 s−1 , KD = 0.15 µM). In this test, the buffers are considered to be immobile. Figs. 3 and 4 show the solution (9) for the case in which [Ca](0) > C1 ; in particular, we considered that there is no bound buffer at t = 0. In Fig. 3 we show the plot for Eq. (9) together with the simulated result for two different concentrations of Fura-2. In Fig. 4 we repeat the test but for EGTA as buffer. The simulated results show the concentration time course in the first layer of a conical domain of 1 µm radius and 5 µm deep (with initial uniform concentrations).
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
281
Fig. 4. Evolution of the free calcium concentration for a single immobile buffering system (EGTA). The simulated results are compared with the analytical solution (Eq. (9), upper signs). There is no bound buffer at t = 0. (A) BT ≡ [EGTA] = 1000 µM, [Ca](0) = CT = 870.56 µM; (B) BT ≡ [EGTA] = 500 µM, [Ca](0) ≡ CT = 490.44 µM.
Fig. 5. Temporal evolution of free calcium for a single buffering system obtained from the simulation and comparison with the analytical solution (Eq. (9), lower signs). There is no free buffer at t = 0 ([Ca](0) = 0). (A) CT = BT = [Fura − 2] = 100 µM; (B) CT = BT = [EGTA] = 500 µM.
Fig. 5 compares Eq. (9) with the simulated results in the case [Ca](0) < C1 . We consider that initially the buffer is completely bound and that [Ca](0) = 0. As can be seen in the figures, the agreement between the simulated time course and the analytical solution is very good in all the cases considered.
282
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293 Table 1 CPU times spent in a Pentium III 450 MHz for the simulation of a calcium time course lasting 0.5 ms. The grid resolution is 100 nm Buffer
k+ (M−1 s−1 )
KD (µM)
DBuffer (µm2 /s)
[Buffer] (µM)
CPU (s)
Fura-2
5 × 108
0.24
42
100
17.11
Fura-2
5 × 108
0.24
42
200
18.75
EGTA
1 × 107
0.15
200
100
34.38
EGTA
1 × 107
0.15
200
200
36.92
BAPTA
5 × 108
0.2
200
100
35.26
BAPTA
5 × 108
0.2
200
200
39.17
Table 2 CPU times spent in a Pentium III 450 MHz for the simulation of a calcium time course lasting 0.5 ms. The grid resolution is 50 nm Buffer
k+ (M−1 s−1 )
KD (µM)
DBuffer (µm2 /s)
[Buffer] (µM)
CPU (s)
Fura-2
5 × 108
0.24
42
100
189.27
Fura-2
5 × 108
0.24
42
200
215.85
EGTA
1 × 107
0.15
200
100
303.83
EGTA
1 × 107
0.15
200
200
426.43
BAPTA
5 × 108
0.2
200
100
291.90
BAPTA
5 × 108
0.2
200
200
413.68
3. CPU times The CPU time spent by the code depends considerably on the parameters of the simulation. In Tables 1, 2 we show the CPU times spent in a Pentium III of 450 MHz (256 MB of RAM) for simulating a calcium time course 0.5 ms long. Results are shown for two different spatial resolutions. We consider 3 exogenous buffers commonly used in fluorescence experiments (Fura-2, EGTA and BAPTA) together with a “standard” immobile endogenous buffer (500 µM, k+ = 5 × 108 M−1 s−1 , KD = 10 µM [2]). The simulation domain is conical, 5 µm deep and 1 µm base radius. The values of the kinetics constants and diffusion coefficients for each buffer are taken from reference [2]. The unitary current is 0.1 pA and the number of channels considered is 30.
Acknowledgements The authors acknowledge the use of the routine BTPEC (algorithm 678 collected from ACM) and function RAND1 [ACM Transactions on Mathematical Software 5 (1979) 132–138], both obtained from NETLIB repository. The authors also acknowledge financial support from Generalitat Valenciana (GV99-146-1-04) and Fundació La Marató de TV3.
References [1] F. Sala, A. Hernández-Cruz, Biophys. J. 57 (2) (1990) 313–324. [2] J. Klingauf, E. Neher, Biophys. J. 72 (2) (1997) 674–690.
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293 [3] [4] [5] [6] [7] [8] [9] [10] [11] [12]
A. Gil, J. Segura, J.A.G. Pertusa, B. Soria, Biophys. J. 78 (1) (2000) 13–33. M. Oheim, D. Loerke, W. Stühmer, R.H. Chow, Eur. Biophys. J. 28 (2) (1999) 91–101. T. Voets, E. Neher, T. Moser, Neuron 23 (3) (1999) 607–615. J. Segura, A. Gil, B. Soria, Biophys. J. 79 (4) (2000) 1771–1786. N.L. Allbritton, T. Meyer, L. Stryer, Science 258 (1992) 1812–1815. V. Kachitvichyanukul, B. Schmeiser, Comm. ACM 31 (2) (1988) 216–222. http://www.netlib.org. L. Schrage, ACM Trans. Math. Software 5 (1979) 132–138. R.D. Present, Kinetic Theory of Gases, McGraw-Hill, 1958, pp. 65-69. J. Crank, The Mathematics of Diffusion, Oxford at Clarendon Press, 1956, pp. 9-11.
283
284
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
TEST INPUT DATA .false. 1 1. 5. 30 1 0.1 1 0.3 0. 0.3 0.2 220. 0. 42 0.1 500. 100. 5.e+8 10. 5.e+8 0.24 3 0 1 0.1
!logical flag for timing (default .false.). See details in ca3d.f !Type of domain. 1 = Conical. 2 = Cylindrical. !Radius of the base of the domain (in microns). !Height of the domain (in microns). !Number of entry calcium channels. !Distribution of channels. 1 = Random. 2 = Specified by the user. !Spatial resolution in microns !Number of calcium pulses. !Duration of each pulse (in ms). !Frequency at which the pulses are delivered (in Hz). !Duration of the simulation (in ms). !Unitary current per channel (in pA). !Diffusion coefficient of ions (in micron**2/s). !Diffusion coefficient of buffer 1 (in micron**2/s). !Diffusion coefficient of buffer 2 (in micron**2/s). !Basal concentration of calcium (in microMol). !Concentration of buffer 1 (in microMol). !Concentration of buffer 2 (in microMol). !Forward binding rate of the endogenous buffer (in Mol**-1s**-1). !Dissociation constant of the endogenous buffer (in microMol). !Forward binding rate of the exogenous buffer (in Mol**-1s**-1). !Dissociation constant of the exogenous buffer (in microMol) !Number of slices at which the averaged [Ca2+] is required. !Are the buffer concentrations also required?. 0-->no, 1-->yes. !Number of 2D calcium maps to be calculated. !time at which the first calcium map is obtained (in ms). !time at which the second calcium map is obtained (in ms). !time at which the third calcium map is obtained (in ms). !....... !time at which the 25th calcium map is obtained (in ms).
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
TEST OUTPUT FILE: ‘casli01’ #********************************* # Averaged calcium concentration # at distances from 0.nm # to 105.nm to the # membrane #********************************* # # time (ms) Calcium (microMol) # 0.013 0.253 0.025 0.281 0.038 0.349 0.050 0.355 0.063 0.332 0.076 0.338 0.088 0.304 0.101 0.332 0.113 0.422 0.126 0.360 0.139 0.400 0.151 0.343 0.164 0.422 0.176 0.388 0.189 0.439 0.201 0.428 0.214 0.445 0.227 0.467 0.239 0.495 0.252 0.450 0.264 0.658 0.277 0.580 0.290 0.456
285
286
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
FILE: ‘casli02’ #********************************* # Averaged calcium concentration # at distances from 105.nm # to 210.nm to the # membrane #********************************* # # time (ms) Calcium (microMol) # 0.013 0.158 0.025 0.191 0.038 0.146 0.050 0.180 0.063 0.174 0.076 0.253 0.088 0.259 0.101 0.281 0.113 0.135 0.126 0.248 0.139 0.186 0.151 0.214 0.164 0.180 0.176 0.264 0.189 0.225 0.201 0.248 0.214 0.276 0.227 0.191 0.239 0.264 0.252 0.197 0.264 0.253 0.277 0.281 0.290 0.281
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
FILE: ‘casli03’ #********************************* # Averaged calcium concentration # at distances from 210.nm # to 315.nm to the # membrane #********************************* # # time (ms) Calcium (microMol) # 0.013 0.186 0.025 0.118 0.038 0.118 0.050 0.141 0.063 0.118 0.076 0.129 0.088 0.107 0.101 0.135 0.113 0.113 0.126 0.096 0.139 0.152 0.151 0.118 0.164 0.152 0.176 0.118 0.189 0.090 0.201 0.073 0.214 0.124 0.227 0.186 0.239 0.203 0.252 0.135 0.264 0.141 0.277 0.113 0.290 0.101
287
288
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
FILE: ‘camap01’ #***************** # 2D calcium map #***************** # # at time 0.100 ms # # # x y Calcium (microMol) # 0 -9 0.000 -4 -8 2.278 -3 -8 0.854 -2 -8 0.142 -1 -8 0.285 0 -8 0.142 1 -8 0.285 2 -8 0.000 3 -8 0.000 4 -8 0.285 -5 -7 0.285 -4 -7 0.285 -3 -7 0.142 -2 -7 0.427 -1 -7 0.712 0 -7 0.285 1 -7 0.285 2 -7 0.569 3 -7 0.854 4 -7 0.142 5 -7 0.142 -6 -6 0.285 -5 -6 0.000 -4 -6 0.142 -3 -6 0.142 -2 -6 0.285 -1 -6 0.142 0 -6 0.285 1 -6 0.142 2 -6 0.285 3 -6 0.427 4 -6 0.142 5 -6 0.142 6 -6 0.000 -7 -5 2.136 -6 -5 0.712 -5 -5 0.000 -4 -5 0.142 -3 -5 0.000 -2 -5 0.427 -1 -5 0.000
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
0 1 2 3 4 5 6 7 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 -8 -7 -6 -5 -4 -3 -2 -1 0 1
-5 -5 -5 -5 -5 -5 -5 -5 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -3 -3 -3 -3 -3 -3 -3 -3 -3 -3 -3 -3 -3 -3 -3 -3 -3 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2
0.000 0.285 0.142 0.142 1.139 0.712 0.000 0.000 0.285 0.285 0.142 0.142 0.142 0.000 0.285 0.142 0.000 0.427 0.142 0.569 0.142 1.851 0.285 0.000 0.142 0.142 0.142 0.142 0.142 0.142 0.000 0.000 0.142 0.427 0.569 1.139 0.997 2.136 0.427 0.569 0.142 0.000 0.285 0.000 0.142 0.285 0.142 0.142 0.427 0.285 1.566 1.566
289
290
2 3 4 5 6 7 8 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 -8 -7 -6 -5 -4 -3 -2 -1 0
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
-2 -2 -2 -2 -2 -2 -2 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1
0.997 0.712 0.569 0.285 0.000 0.000 0.142 0.000 0.000 0.142 0.142 0.285 0.142 1.424 1.281 1.139 0.000 0.142 0.000 0.142 0.142 0.285 0.000 0.142 0.142 0.000 0.000 0.142 0.142 1.139 0.712 0.569 0.142 0.142 0.427 0.142 0.000 0.000 0.142 0.285 0.427 0.427 0.427 0.142 0.000 0.000 0.000 0.285 0.000 0.997 0.285 0.285
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
1 2 3 4 5 6 7 8 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 -8 -7 -6 -5 -4 -3 -2 -1 0 1
1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4
0.000 0.000 0.427 0.000 0.142 0.142 0.142 1.139 0.000 0.000 0.142 0.142 0.142 0.000 0.285 0.285 0.000 0.142 0.427 0.142 0.285 0.427 0.142 0.427 0.427 0.569 0.854 0.000 0.142 0.000 0.142 0.000 0.000 0.142 0.000 0.000 0.427 0.427 1.708 0.569 1.993 0.854 2.420 0.142 0.285 0.000 0.285 0.000 0.569 0.427 0.000 0.000
291
292
2 3 4 5 6 7 8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 -5 -4 -3 -2 -1 0 1 2 3 4 5 -4 -3 -2 -1 0 1
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 7 7 7 7 7 7 7 7 7 7 7 8 8 8 8 8 8
0.142 1.424 0.142 0.427 0.142 1.424 0.427 0.000 0.142 0.142 0.997 2.278 1.424 0.854 0.569 0.142 0.569 0.142 0.285 0.142 0.712 0.142 0.000 0.285 1.993 1.566 0.569 0.285 0.569 1.708 0.712 0.285 0.427 0.142 0.285 0.142 0.427 0.712 0.569 0.142 0.285 0.569 0.427 0.997 0.285 0.285 0.142 0.569 1.708 0.569 0.000 0.142
A. Gil, J. Segura / Computer Physics Communications 136 (2001) 269–293
2 3 4 0
8 8 8 9
0.427 0.569 0.000 0.000
293