CALCIUM-SECRETION COUPLING IN A PRESYNAPTIC ACTIVE ZONE MODEL

A theoretical analysis of some of the relevant factors influencing the calcium time course and the strength and timing of release probabilities of vesicles evoked by an action potential in a calyx-type active zone is presented in this paper. In particular, our study focus on the comparison of cooperative vs non-cooperative calcium binding by the release site and the effect of the number of Ca2+ binding sites on the calcium sensitivity for release. Regarding the comparison of cooperative and non-cooperative kinetic schemes, our simulations show that quite different results are obtained when considering one or another: a reduction in the release probability of more than a 50% is obtained when considering the cooperative kinetic scheme. Also, a delay in the average time for release appears when using this model for the calcium sensor. Our study also shows that a non-cooperative kinetic binding scheme gives rise to a well defined average calcium level for release assuming that the same kinetic constants are considered for all the sites. Our results also suggest that the central value of the calcium sensitivity for release depends on the number of binding sites N and the dissociation constant KD with a scaling law depending on NKD . 2010 Mathematics Subject Classification. Primary: 92B05, 65C35; Secondary: 92C20.

1. Introduction.Release of neurotransmitter by small clean granules (size ∼ 45 nm of diameter) in synapses is known to depend on the Ca 2+ concentration.This dependence is interpreted by considering that N Ca 2+ ions have to bind to synaptotagmin or other proteins to produce the fusion of a vesicle with the membrane.However, the exact number of kinetic steps involved in the fusion process (the number N of "binding sites" of the secretory vesicles) is not precisely known and an extrapolation from the experimental data is usually made in order to estimate N .
In this paper, we study some of the main elements influencing granule fusion in an active zone of a presynaptic terminal model: the calyx of Held, since this is a large synapse widely used as a model to study stimulus-secretion coupling in neurotransmitter release (for recent papers see [9,11,12], for example).
The Ca 2+ dynamics in the presynaptic active zone is modeled using a fully stochastic approach [8,20,6] to model all the processes involved: given the small number of ions entering in the presynaptic active zone during an action potential, it seems appropriate to keep track of each individual ion and molecule instead of using a continuous modeling in terms of concentrations.We simulate the calcium entry corresponding to a typical action potential (AP); for this purpose, the twostate kinetic model proposed in [5] is considered for the calcium channels.For the kinetic scheme of the presynaptic vesicles, we concentrate in the description provided by the non-cooperative kinetic scheme considered in [2].This kinetic model seems to provide all the necessary ingredients for a fair description of granule fusion in this type of presynaptic terminals.However, in our study we also quantify the influence of a cooperative kinetic scheme, which have been also proposed for calyx [19], on the release probability of vesicles.Other kinetic models have been proposed in [13,21].The comparison of the cooperative and the non-cooperative kinetic schemes shows that quite different results are obtained when considering one or another: a reduction in the release probability of more than a 50 % is obtained when considering the cooperative kinetic scheme.Also, a delay in the average time for release appears for this model of the calcium sensor.As we discuss, these results are highly dependent on the particular choice of parameter values used in the kinetic models.
The influence on release of the kinetic parameters of the calcium sensors and the number of binding sites in each sensor (N ) is also analyzed in our study.The calcium sensitivity for release is seen to linearly depend on N K D (where K D is the dissociation constant parameter of the kinetic scheme) while the width of the temporal distribution of release events is narrower as N increases.A first order Magnus approximation to the solution of the problem dX dt = {c(t)A+B}X, X(0) = I, with X(t) = (X 0 , X 1 , ..., X N ) representing the probabilities of occurrence of the states X 0 , X 1 , ..., X N for a non-cooperative kinetic scheme, is obtained in order to confirm this prediction.
In order to design an "optimal" configuration which maximizes both the synaptic strength and the time efficiency of release, it is important to take into account three main effects ( supported by a qualitative analysis of the kinetic reactions involved) acting in different directions: 1.The lower the width of the temporal distribution of the release events, the more synchronous the release is in the sense that all granules will undergo exocytosis approximately at the same time.The width decreases with increasing number of binding sites.
2. However, the more binding sites the more time it will take to start exocytosis.
A delay time appears as N grows.3. Also, the more binding sites the lower the probability of release is.
The magnitude of these effects and their impact on the resulting exocytotic response, is quantified in our study for several configurations (N variable).

Methods.
2.1.Non-cooperative kinetic scheme for secretory vesicles.In the modeling of vesicle fusion in calyx of Held, we focus on the non-cooperative kinetic scheme of [3].Within this kind of kinetic schemes, it is assumed that N (N = 5) calcium ions have to be bound to the calcium sensor according to a linear chain of kinetic rate equations: A state X i describes the situation in which i binding sites are occupied by a calcium ion and X N represents the final state corresponding to a vesicle with no free binding sites.The calcium sensor has a small number of binding sites and the calcium ions in the neighborhood can bind to any of the free binding sites with equal probabilities.Conversely, any binding site occupied by a calcium ion has the same probability per unit time of releasing its bound calcium ion.When the calcium sensor has all its binding sites occupied by calcium ions, a release-promoting state (reversible) with rate constants γ and δ is considered; then, finally the kinetic chain terminates with a irreversible state, and neurotransmitters are released.The two final steps in (1) result in a delay time for fusion.
A non-cooperative kinetic scheme corresponds to a grading of the binding constants: where k + = [Ca 2+ ]k on and k − ≡ k of f , being k on the forward binding rate and k of f , the backward binding rate (or dissociation rate).Particular values for the kinetic parameters are given in Table 1.The dissociation constant(K D ) is given by

2.1.1.
A first order Magnus approximated solution to the matrix linear differential equation for X 0 , ..., X N .In matrix form, the temporal evolution of the linear kinetic subchain for X 0 , ..., X N given in (1) can be described by the following matrix problem: The cofficients of the A and B matrices are obtained from Eq. ( 2).We will take as boundary conditions: In this view, we can interpret X i as a normalized concentration or, alternatively, as the probability of occurrence of the state i.Then, mass conservation imposses To obtain an approximate solution of (3) for general N we will consider a Magnus expansion [14,1].This approximate Magnus solution will be later used in our analysis to confirm the predictions made by a microscopic (stochastic) model which accounts for discrete calcium ion binding.
The starting point of the Magnus expansion consists in writting the solution of the matrix problem X (t) = M(t)X, X(0) = I as X(t) = e Ω(t) X(0), (7) where Ω(t) is given by an infinite series Ω(t) = ∞ i=1 Ω i (t).The first two terms of the series are given by where . In general Ω m will be a m-multivariate integral involving a linear combination of nested commutators of M(t) at different times.
For our matrix problem (3), the first two terms of the Magnus series correspond to with [A, B] = AB − BA.
In order to obtain the simplest approximated solution to our problem in the Magnus scheme, we take Ω 1 .In this way, a first order Magnus approximation to the solution of the matrix linear differential equation problem will be obtained in the form X(t) ≈ e Ω1 X(0).A straighforward calculation leads to the following explicit expression for this approximation to the solution X(t) = (X 0 (t), ..., X N (t)) of (3): where The numerators in the expressions of X 0 (t), X 1 (t), ..., X N (t) in ( 10) correspond to each of the N + 1 terms, respectively, of the binomial theorem On the other hand, the denominators are all equal to < c > this way, it is easy to check the property of mass conservation It is also interesting to note that taking lim k of f → 0 in the (approximated) expression for X N (t) given in (10), we obtain the (exact) solution of a fully irreversible kinetic chain with N steps a power law which is enunciated in experimental studies by saying that the rate of neurotransmitter (or hormone) release goes as the calcium current to the power of N , being N the number of binding sites in the calcium sensor.

Cooperative kinetic scheme for secretory vesicles.
In our study, we will compare the results obtained with the non-cooperative kinetic scheme introduced in Section 2.1 against the cooperative model of [19]: The "cooperativity" is included in this model through the backward parameters k i,i−1 , which depend on a cooperativity factor b as k i,i−1 = iηb i−1 , i = 1, .., N , being i the step in the kinetic chain and η is an unbinding constant.Two additional parameters (the forward rate k on and the fusion rate γ) complete the model.Particular values of the parameters can be found in Table 1.
The cooperative kinetic scheme corresponds to the following recursive relations for the kinetic constants: The kinetic parameter k + is proportional to the forward rate k on through the calcium concentration (k + = [Ca 2+ ]k on ).As in the case of the non-cooperative scheme, the temporal evolution of the linear kinetic subchain for X 0 , ..., X N is given by Eq.( 3).The matrix A is the same as for the non-cooperative scheme (Eq.( 4)) while the matrix B is, in this case, given by However, in contrast to the non-cooperative kinetic scheme, it seems only possible to find simple expressions for the solution of the matrix problem by a first order Magnus approximation when N is small.As an example, for a linear kinetic subchain with 2 binding sites the first order Magnus approximation for the cooperative scheme gives the following density of probability of the state X 2 : where and

2.3.
A microscopic model for a presynaptic active zone.We consider that in an active zone of the calyx of Held ∼ 20 channels not evenly distributed will open during the peak of Ca 2+ influx evoked by an action potential, and that at least 2-3 vesicles are docked to membrane, ready to be fused and release their contents [4,18].The active zone is modeled as a cylinder of 0.1 μm of radius and 0.4 μm height, where a 3-D orthogonal grid maps the cylindrical domain in small cubes of Δx nm per side.The first layer of the cylindrical domain represents the cell membrane; each cube of this layer would contain a channel, a vesicle, calcium ions and/or buffer molecules.The subsequent layers, representing the cytosol, only contain buffer molecules and calcium ions.
An schematic representation of the simulation domain (active zone of the calyx of Held) considered in the simulations is shown in Figure 1.

A) B)
Calcium channel Vesicle sensor Figure 1.A) Schematic representation of the domain considered for simulations.B) Representation of the first layer of the domain, in which both calcium channels and vesicles are placed.This layer would also contain calcium ions and buffer molecules.
Ions and mobile buffer molecules move from one compartment to another compartment of the 3-D 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 and by the largest diffusion coefficient of the system.The relation between those scales is given by: Δt = (Δx) 2 /4D (17) being Δx the spatial resolution and D the largest diffusion coefficient.

2.3.1.
Input current for a single action potential.Although different types of calcium channels (P/Q-, N-and R-type) are involved in the calcium entry in the calyx of Held terminal, it seems that there are no significant differences in the time course of Ca 2+ corresponding to each of them.So, we will assume that the calcium current during a single action potential is simulated for any calcium channel by using the simple two-state model proposed in [5].The model consists of two states (Open and Closed) and transitions between states are described by two voltagedependent functions: α = α 0 exp(V/V α ) (for transitions from Closed to Open) and β = β 0 exp(−V/V β ) (for transitions from Open to Closed).We use the values for the parameters given in Table 1.
As a test of the implementation of the two-state model, Figure 1 shows the simulated current for 1000 channels (lower panel) considering a prototype action potential (upper panel).Figure 2. Simulated calcium current (lower panel) by using the two-state model proposed in [5] during an action potential (upper panel).The calcium current corresponds to a simulation with 1000 channels.
The simulated calcium influx peak is ∼ 100 pA for the average of the simulations.This seems to be in agreement with experimental data, given that the estimated number of calcium channels in the whole calyx membrane would be in the order of 25, 000 open channels during an action potential [16] and that the peak of the whole-cell current entering in the calyx terminal during an action potential seems to be of the order of 2.5 nA [5].

2.3.2.
Solving the kinetic reactions in the microscopic model.The reaction of buffers with calcium are described by first order kinetic reactions: being k i + the forward binding rate for the buffer i and k i − , the unbinding rate.These kinetic reactions are solved stochastically.This means that in terms of the number of free calcium ions N Ca , free buffer molecules N Bi and bound buffer molecules N CaB i , Eq.( 18) reads: where N B T i represents the total (free + bound) number of buffer particles of type i, N Ca T is the total number of (calcium) ions and ΔN Ca represents the variation of free calcium ions in a given cubic compartment during a time interval dt ≡ Δt/n.The effective constants K i − and K i + for unbinding and binding in each time interval dt are given by 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.This means that each Ca 2+ ion has a probability P i b = K i + N Bi of being bound by any of the N Bi free buffer molecules during the time dt = Δt/n; the corresponding probability of becoming unbound in the interval dt is P i u = K i − .The value of n in (20) is selected for every time and for each compartment of the 3-D orthogonal grid before the kinetic equations are integrated.Thus, the effective constants K i − and K i + can be different at different time intervals of the simulation algorithm, though preserving the same ratio.The selection of the appropriate number of subdivisions of the interval is made by requiring that all binding (and unbinding) probabilities are smaller enough than 1; in other words, we choose n by bisectioning the interval as many times as required to fulfill the inequalities: Numerical experiments show that taking P max = 0.1 is a safe selection that guarantees proper convergence of the method.
It is important to note that in our implementation, the calcium sensors of vesicles will be also considered as additional buffer molecules.So, the kinetic reactions described in (1) will be included in (19); in this way, all the buffer species will compete for the calcium ions according to their reaction probabilities.Given the low number of calcium ions entering in the presynaptic active zone during an action potential, this methodology (which takes into account the local variance of [Ca 2+ ]) seems more appropriated than methods using average values of [Ca 2+ ] at vesicles to predict release probabilities, as also commented in [15].

2.3.3.
The resulting stochastic algorithm.The full algorithm is executed according to the following sequence: 1. Presimulation.The number of ions and buffer molecules given by the equilibrium conditions are distributed randomly and uniformly over the cylindrical domain.2. Simulation.For each simulation step, the following actions are taken and repeated in this order until the end of the simulation: entry of ions according to the stochastic model of the calcium channel; computation of the kinetic reactions of calcium with buffers and vesicles; diffusion of Ca 2+ and mobile buffers; evaluation of the possible secretory events.
The outputs of the algorithm are concentration time courses and the total number of secretory events.The full algorithm has been extensively tested and its results proved to agree with theoretical expectations [8] at long times.As an example, the validity of the random walk computational scheme was established by testing that free diffusion of particles resulted in correct Gaussian distributions.The algorithm implementing the first-order kinetic reactions was tested by checking that starting from concentrations far from equilibrium, equilibrium was reestablished and maintained.Further details can be found in [7] .

Results and discussion.
As in [15], an endogenous fixed buffer (EFB) and a mobile buffer (ATP) with concentrations and kinetic constants given in Table 1 are assumed to be present in the medium.As we are going to compare some of our results against those given in [15], we exclude in our simulations other buffers such as parvalbumin [17].
Figure 3 shows the average calcium levels (0 − 10 nm) obtained using our microscopic model when the two buffers (the endogenous fixed buffer and ATP) are present in the medium.We show the results for an active zone with 20 calcium channels and 3 vesicles in the first layer of the simulation domain (see Figure 1B for an schematic representation).The calcium levels correspond to the average of simulated calcium concentrations from 0 to 10 nm to the cellular membrane obtained with 30 random configurations of calcium channels and vesicles.As can be seen, the calcium levels peak at ∼ 10 μM which is in agreement with estimates of [Ca 2+ ] i observed by the calcium sensors during an action potential and which are estimated using laser photolysis experiments [3].No calcium extrusion mechanisms are considered in the simulations.The simulated action potential and incoming calcium current are shown in Figure 4.The release probability of vesicles can be defined as the fraction of all readily releasable vesicles that are released during a single action potential [15].For a fixed calcium entry and under the same endogenous buffering conditions, two different type of factors have influence on the release probability of vesicles: the geometrical configuration of vesicles and calcium channels and the kinetics of the calcium sensors.We will focus our study on the influence of this last factor.The channels are modeled according to the two-state model of [5].
A first analysis concerns the particular choice of the kinetic model and the corresponding set of kinetic parameters.As commented in section 1, two standards models proposed for calyx of Held are the non-cooperative scheme of [3,2] (which is considered in our simulations) and the cooperative model of [19].Both models are described in Section 2.1 and 2.2, respectively.The comparison of both kinetic schemes with their corresponding parameter values shows that quite different results are obtained when considering one or another kinetic scheme for the secretory vesicles: in Figure 5A we show that a reduction in the release probability of more than a 50 % (from 8 % to a 3 %) is obtained when considering the cooperative kinetic scheme.This relative reduction in the release probability was also pointed out in [15].
Our results are obtained as the average of the 30 runs of our simulation algorithm considered before.For each configuration, we compute the secretory events obtained when a) the non-cooperative model or b) the cooperative model is considered as the kinetic scheme for the vesicle sensors.
Also, a delay in the average time for release of about 300 μs (see Figure 5B) appears for this model of the calcium sensor.Both results are explained because of the reduced calcium sensitivity of the cooperative scheme in comparison to the non-cooperative kinetics.This means that in order to obtain the same release probability, calcium sensors with kinetics described by the cooperative scheme should experience higher calcium levels.The results obtained are not surprising due to the particular choice of parameters for the cooperative scheme: when computing the dissociation constants ≈ 17.08 μM , which is higher than the dissociation constant for the non-cooperative scheme (K D = 10 μM ).
The differences in the results obtained with both kinetic schemes (which were proposed in connection to different sets of experimental results in calyx of Held) and which have been also mentioned in [15] and [2], seem to be related to the Figure 5. A) Release probability obtained with the noncooperative kinetic scheme of [3] and the cooperative scheme of [19].B) Average time for release obtained with the non-cooperative kinetic scheme of [3] and the cooperative scheme of [19].
maturation stage of the calyces considered in the experiments: as discussed in [22], calyces with different maturation stage were used in the experiments where the noncooperative (younger calyces) and the cooperative (more mature calyces) kinetic schemes, respectively, were proposed.In this sense, it seems that synaptic vesicles in more mature calyces show a lower calcium sensitivity of the release apparatus and, on the other hand, sense local calcium transients with higher peak concentrations during action potential-evoked phasic release, in comparison to "younger" calyces (P9-11).This would represent a tight geometrical coupling between vesicles and channels in mature calyces.
As the comparison between the non-cooperative and the cooperative kinetic scheme seems clearly influenced by the particular choice of the parameter values, it seems interesting to compare the results obtained when the same values are taken both for the dissociation constant of a non-cooperative scheme and the "effective" dissociation constant of a cooperative scheme.We will also take the same value for the k on parameter (k on = 9 × 10 7 M −1 s −1 ).In order to simplify the analysis, we will consider the linear kinetic subchain with 2 binding sites given in (15) and use the first-order Magnus approximation to the solution of the matrix problem (3).In Figures 6A,B we compare the probability of occurrence of the state X 2 for the noncooperative scheme given by (10) for N = 2 and the corresponding values obtained with the cooperative scheme ( 16) as a function of the average calcium concentration < c > t .The value of K D in the non-cooperative scheme (K non−coop D ) has been fixed 14 AMPARO GIL, VIRGINIA GONZ ÁLEZ-V ÉLEZ, JAVIER SEGURA AND LUIS MIGUEL GUTI ÉRREZ to 10 μM .The value of the effective ) in the cooperative scheme has been also fixed to 10 μM in both Figures 6A,B.Different choices of K 1 D have been considered in Figure 6A.As can be seen, the probability of occurrence of X 2 switches from being smaller than the non-cooperative scheme when K 1 D = 75, 100 μM to larger when K 1 D = 10, 25, 50 μM .This is pointing out that as the first reaction becomes more favorable (smaller K 1 D ), the cooperative scheme switches from less effective to more effective, compared with the non-cooperative scheme.As the figure shows, this switch occurs when K 1 D ≈ 50μM , which is five times larger than the non-cooperative K D value.In Figure 6B the parameter b in equation ( 16) is now varied.It should be noted that when the value of both the effective K D value and b are fixed, the value of η is fixed as well.As can be seen, the cooperative scheme is less effective than the non-cooperative scheme when a very small value of the parameter b is considered; this is not surprising because the corresponding η value is very large.As b increases, the cooperative scheme turns now to be more effective than the non-cooperative scheme, reaching the maximum of effectivity for b = 0.3.For larger values of b the density of probability of the state X 2 (t) for the cooperative scheme decreases, although it is always more effective than the non-cooperative scheme.16)).
The final part of our analysis concerns the number of binding sites (N ) of the secretory vesicles: kinetic schemes in synaptic models are proposed for a fixed number of binding sites (5, in the case of both the non-cooperative and the cooperative kinetic models); however, this number can not be exactly determined in the experiments.We have then varied the number of binding sites of the secretory vesicles in order to see which is the effect on both the synaptic timing (the temporal distribution of release events) and the synaptic strength (the release probability).In order to perform a "fair" comparison of the different configurations, we have tried to preserve the same calcium sensitivity when varying N : i.e., the different configurations have to be sensitive to a similar calcium level, being defined this level as the average calcium concentration when a exocytotic event takes places.This means that if the number of binding sites are modified, the parameters of the kinetic model should also be varied.Actually, experimental simulations shows that a scaling law of the form N K D seems to control this calcium level, as can be seen in Figure 7A, where the product N K D is kept to be constant (from the parameters of Table 1, we take this product equal to 50).This scaling law can be also tested using the approximated formula for the density of probability of the state N (X N (t)) obtained in Eq.( 10): fixing the values of t = 1 ms, which can be taken as a characteristic timing associated to an action potential, and k of f = 3 ms −1 (see Table 1, taking into account that k of f = K D k on ), Figure 8 shows X N (t = 1 ms) as a function of the average calcium concentration < c > t and for the same values of N and K D considered in Figure 7A.As can be seen, quite similar curves are obtained in all cases.As these curves represent the density of probability in the state N as a function of the average calcium, a similar sensitivity to calcium is expected for the tested configurations.
Regarding the influence of N on the synaptic strength, Figure 7B shows the release probabilities obtained for the different values of N .As expected, the highest (by far) release probability is obtained with the 2-bs configuration (∼ 26 %) while for the rest of configurations, the release probability is below 15 %.So, it is clear that the efficiency of the secretory response will be favoured with a configuration with a low number of binding sites.However, one important aspect regarding synaptic transmission is also how synchronous the signal is [10].This means that, ideally, the temporal window for the release events should be as narrower as possible.Figure 7C shows the temporal distribution of the release events as a function of the number of binding sites.Median values are plotted, as well as bars from the minimum (maximum) value to the median values.In this case, it is clear that the synchronicity of the signal will be favoured when a large number of binding sites is considered.

X
i (t) = 1 for the first order Magnus approximated solution.

Figure 3 .
Figure 3. Average calcium levels corresponding to 0−10 nm from the presynaptic membrane surface.An active zone with 20 calcium channels is considered in the simulations.

Figure 4 .
Figure 4. A) Simulated action potential.B) Incoming calcium current corresponding to 20 calcium channels in the active zone.The channels are modeled according to the two-state model of[5].
5 in each step of the cooperative kinetic chain, the values obtained are K 1 D ≈ 105.5 μM, K 2 D ≈ 52.7 μM, K 3 D ≈ 19.8 μM, K 4 D ≈ 6.6 μM and K 5 D ≈ 2 μM .This means that, approximately, the effective dissociation constant for the cooperative scheme can be estimated as K coop D

Figure 6 .K 1 D K 2 D 1 / 2
Figure 6.Comparison of the probability of occurrence of the state X 2 for the non-cooperative scheme given by(10) for N = 2 and the corresponding values obtained with the cooperative scheme (16) as a function of the average calcium concentration < c > t .The values of k on and t have been fixed to 9 × 10 7 M −1 s −1 and 1 ms, respectively.The values of the dissociation constant of the non-cooperative scheme (K non−coop D

Figure 7 .
Figure 7. A) Average calcium concentration from 0 to 10 nm when a release event takes places for different numbers of binding sites (N bs ) of the secretory vesicles: 2, 3, 4, 5, 6.The values of the dissociation constants K D have been modified in order to keep N bs K D = constant.We consider K D = 10 μM when N bs = 5, as in Table 1.B) Comparison of the release probabilities for a kinetic sensor with a different number of binding sites (N bs = 2, 3, 4, 5, 6).C) Temporal distribution of the release events as a function of the number of binding sites (N bs = 2, 3, 4, 5, 6).

Figure 8 .
Figure 8. Density of probability for the state N (X N (t) in Eq.(10)) plotted as a function of the average calcium concentration < c > t for different values of binding sites (N ) of the secretory vesicles and dissociation constant K D .

Table 1 .
Parameters used in the simulations