Abstract
The exocytotic process in neurons and neuroendocrine cells consists of a sequence of reactions between well defined proteins. In the present study, we have created for the first time a comprehensive kinetic model that demonstrates the dynamics of interactions between key synaptic proteins that are associated with exocytosis. The interactions between the synaptic proteins were transformed into differential rate equations that, after their integration over time, reconstructed the experimental signal. The model can perfectly reconstruct the kinetics of exocytosis, the calcium-dependent priming and fusion processes, and the effects of genetic manipulation of synaptic proteins. The model suggests that fusion occurs from two parallel pathways and assigns precise, non-identical synaptic protein complexes to the two pathways. In addition, it provides a unique opportunity to study the dynamics of intermediate protein complexes during the fusion process, a possibility that is hidden in most experimental systems. We thus developed a novel approach that allows detailed characterization of the temporal relationship between synaptic protein complexes. This model provides an excellent platform for prediction and quantification of the effects of protein manipulations on exocytosis and opens new avenues for experimental investigation of exocytosis.
Introduction
Exocytosis is thought to be mediated by a sequence of interactions between cytosolic, vesicular, and plasma membrane proteins and can be described as a sequential process of vesicle maturation between various vesicular pools. Accordingly, the kinetics of exocytosis are well described by fusion of vesicles from specific vesicular pools (Neher, 1998; Sudhof, 2000; Rettig and Neher, 2002; Rosenmund et al., 2003). In the last decades, the functions of specific proteins in this process, the way they interact with each other, and the effect of calcium on exocytosis have been studied intensively (Littleton and Bellen, 1995; Fernandez-Chacon and Sudhof, 1999; Jahn and Sudhof, 1999; Sudhof, 2000; Rettig and Neher, 2002; Richmond and Broadie, 2002; Rosenmund et al., 2003). However, no attempt was ever made to assemble this knowledge into a rigorously detailed mechanism. In the present study, we have created for the first time a comprehensive kinetic model describing the dynamics of interactions between key synaptic proteins that are associated with exocytosis. The model accurately reconstructs the kinetics of exocytosis under various experimental conditions. The model significantly extends the macroscopic formalism in which the dynamics are presented as a sum of exponential processes, by replacing them with a detailed reaction mechanism. Thus, instead of time constants and amplitudes, we quantitated the exocytosis process by the concentrations of the protein and the rate constants of their interactions. This allows a description of the releasable vesicular pools as defined protein complexes.
Materials and Methods
The exocytotic process is a sum of many sequential events, and its modeling was based on an accepted order of reactions between defined proteins, the concentrations and rates of reactions of which were expressed by first- and second-order rate constants. Accordingly, the exocytotic process was transformed into a set of binary reactions that describes the reactions between synaptosome-associated protein of 25 kDa (SNAP)-25, syntaxin, vesicle-associated membrane protein (VAMP)/synaptobrevin, Munc13, complexin, synaptotagmin, and Ca2+ ions. The sequence of events and the parameters controlling the reactions are given in Scheme I and Table 1, respectively. The reaction steps were written as a linear sequence of chemical equilibria and transformed into a set of chemical rate reactions in which the product of one step is the reactant of the next one. Once the sequence of events was established, the model was converted into a set of coupled rate equations using standard kinetic formalism, in which the velocity of the reaction is a product of the concentrations of reactants times the rate constant of the reaction. The differential rate equations corresponding with the sequence of events are given below:
The sequence of events used in the model to simulate the exocytotic process is shown. All reaction steps are arranged as a linear sequence up to step 5, in which the flux is divided into two parallel pathways. Each pathway operates with a different type of calcium sensor. Please note that reaction 0 is not part of the linear sequence. The parameters are listed in Table 1.
Values for the protein concentrations and the rate and equilibrium constants that reconstruct the experimental protocols presented in this study
Reaction 1: d[binary SNARE complex]/dt = k1 × [syntaxin] × [SNAP-25] - k-1 × [binary SNARE complex] + k-2 × [SNARE] - k2 × [binary SNARE complex] × [VAMP].
Reaction 2: d[SNARE]/dt = k2 × [binary SNARE complex] × [VAMP] - k-2 × [SNARE] + k-3 × [SNARE′] - k3 × [SNARE] × [Munc13′].
Reaction 3: d[SNARE′]/dt = k3 × [SNARE] × [Munc13′] - k-3 × [SNARE′] + k-4 × [SNARE⋆] - k4 × [SNARE′] × [Ca2+].
Reaction 4: d[SNARE⋆]/dt = k4 × [SNARE′] × [Ca2+] - k-4 × [SNARE⋆] + k-5 × [SNARE#] - k5 × [SNARE⋆] × [complexin].
Reaction 5: d[SNARE#]/dt = k5 × [SNARE⋆] × [complexin] - k-5 × [SNARE#] + k-6 × [RC-I] - k6 × [SNARE#] × [synaptotagmin-I] + k-11 × [RC-II] - k11 × [SNARE#] × [synaptotagmin⋆].
Reaction 6: d[RC-I]/dt = k6 × [SNARE#] × [synaptotagmin-I] - k-6 × [RC-I] + k-7 × [RC-I-1Ca] - k7 × [RC-I] × [Ca2+].
Reaction 7: d[RC-I-1Ca]/dt = k7 × [RC-I] × [Ca2+] - k-7 × [RCI-1Ca] + k-8 × [RC-I-2Ca] - k8 × [RC-I-1Ca] × [Ca2+].
Reaction 8: d[RC-I-2Ca]/dt = k8 × [RC-I-1Ca] × [Ca2+] - k-8 × [RC-I-2Ca] + k-9 × [RC-I-3Ca] - k9 × [RC-I-2Ca] × [Ca2+].
Reaction 9: d[RC-I-3Ca]/dt = k9 × [RC-I-2Ca] × [Ca2+] - k-9 × [RC-I-3Ca] - k10 × [RC-I-3Ca].
Reaction 10: d[fusion-RC-I]/dt = k10 × [RC-I-3Ca].
Reaction 11: d[RC-II]/dt = k11[SNARE#] × [synaptotagmin⋆] - k-11[RC-II] + k-12[RC-II-1Ca] - k12[RC-II] × [Ca2+].
Reaction 12: d[RC-II-1Ca]/dt = k12[RC-II] × [Ca2+] - k-12[RC-II-1Ca] + k-13[RC-II-2Ca] - k13[RC-II-1Ca] × [Ca2+].
Reaction 13: d[RC-II-2Ca]/dt = k13 × [RC-II-1Ca] × [Ca2+] - k-13 × [RC-II-2Ca] + k-14 × [RC-II-3Ca] - k14 × [RC-II-2Ca] × [Ca2+].
Reaction 14: d[RC-II-3Ca]/dt = k14 × [RC-II-2Ca] × [Ca2+] - k-14 × [RC-II-3Ca] - k15 × [RC-II-3Ca].
Reaction 15: d[fusion-RC-II]/dt = k15 × [RC-II-3Ca].
Reaction 0: d[Munc13′]/dt = k0 × [Munc13] × [Ca2+] - k-0 × [Munc13′].
All of the reactions are reversible, except for the fusion events (equations 10 and 15), which are irreversible. The upstream events controlling the regulation of monomeric SNAP receptors (SNAREs) and the downstream reactions of endocytosis were not included in the model. In addition, recycling and regeneration of vesicles were ignored at this stage of the model. Please note that reaction 0 is not a part of the mainstream process and progresses independently. Finally, the equilibrium constant for each reaction (Table 1) is given by the ratio of the partial rate constants (Keq = k-1/k1), which means that the forward and backward rate constants are coupled through the equilibrium constant.
Each one of the differential rate equations corresponds with one step of the reaction scheme and includes the influx from the preceding reaction and the flux of matter into the next step. For example, reaction 1 states that the velocity of the formation of the binary SNARE complex (d[binary SNARE complex]/dt) is the sum of the rate of its formation by the reaction between syntaxin and SNAP-25 (k1 × [syntaxin] × [SNAP-25]) and the dissociation of SNARE (k-2 × [SNARE]), minus the rates of its dissociation (-k-1 × [binary SNARE complex]) and its consumption through the interaction with the VAMP to form the SNARE complex (-k2 × [binary SNARE complex] × [VAMP]).
The equations comply with the detailed balance principle and the conservation of mass law. For example, the amount of syntaxin, at any time during the progression of the reaction ([syntaxin]t), is given by its initial concentration ([syntaxin]0) minus the sum of all downstream intermediates that contain syntaxin as a component:
[syntaxin]t = [syntaxin]0 - [binary SNARE complex]t - [SNARE]t - [SNARE′]t - [SNARE⋆]t - [SNARE#]t - [RC-I]t - [RC-I-1Ca]t - [RC-I-2Ca]t - [RC-I-3Ca]t - [fusion-RC-I] - [RC-II]t - [RC-II-1Ca]t - [RC-II-2Ca]t - [RC-II-3Ca]t - [fusion-RC-II].
Rate constants evaluation. The integration of the differential rate equations given above would reconstruct the experimental observations only if the model is supplied with the appropriate protein concentrations and rate constants. Thus, the main task was to search for these values by treating both the concentrations and the rate constants as adjustable parameters. For the rate constants, the search was performed in the range between an upper limit (as approximated by the Debye-Smoluchowski equation; see Appendix) and six orders of magnitude below. For example, the upper limit for the rate constant for the reaction between free diffusing Ca2+ ions with proteins was set to be as high as 109 m-1sec-1. In contrast, the reaction rate between the cytosolic protein Munc13 with a SNARE complex is slower because of the low diffusion coefficient of the protein, and the upper limit for the rate constant was set to ∼108m-1sec-1. The rate of reaction between the membrane proteins syntaxin and SNAP-25, having a two-dimensional diffusion coefficient of ∼10-9cm2sec-1, was set with an upper limit of ∼107 m-1sec-1. The interaction between VAMP, which resides on a docked vesicle, and the binary complex, which resides on the plasma membrane, was calculated according to the two-dimensional diffusion coefficient of the binary complex (∼10-9cm2sec-1). To account for the existence of several copies of VAMP on the same vesicle, which increase the probability of its encounter with the binary SNARE complex, the upper limit of the rate constant was estimated according to Berg and Purcell (1977) to be 3 × 107 m-1sec-1 (see Appendix). The upper limit for the rate constant of the reaction between SNARE′ and synaptotagmin was estimated as ∼108 m-1sec-1, as calculated for the diffusion of proteins on the surface of a vesicle using the theoretical expressions of Hardt (1979) for diffusion controlled reaction on a two-dimensional space.
Previous measurements, performed with various experimental systems, had defined the value of 19 of the rate constants (k1, Keq1, k-2, k5, k-5, k7, k-7, k8, k-8, k9, k-9, k10, k12, k-12, k13, k-13, k14, k-14, k15) (Table 1), and the published values were taken as a first estimation. These rate constants were further refined by a search within a narrow range of one order of magnitude. Table 1 lists the most suitable values of all adjustable parameters that were investigated in this study. It should be stressed that, in all cases, the search yielded rate constants that were equal or smaller than the upper estimation, meaning that they are theoretically acceptable.
The concentrations evaluation. The reactant concentrations are given in molar units, which in the case of some reactants calls for algebraic transformations as detailed below. For soluble proteins, or Ca2+ ions, the molar units are readily applicable. In contrast, when VAMP is present in more than a single copy on a vesicle and each VAMP can initiate an interaction with the binary complex, the formal concentration of the protein is inadequate for kinetic formulation. To account for that, we based our calculations on the molar concentration of the vesicles, and the fusion process was reconstructed in discrete “quanta,” each of them representing the fusion of one vesicle. The correlation between the number of identical copies of the reactive sites (VAMP molecules) on the vesicle and the rate of its interaction with the binary complex was treated by the theoretical model of Berg and Purcell (1977) and is discussed in Appendix. The molar concentration of vesicles was calculated to be 10 nm, assuming 1000 docked vesicles in a chromaffin cell (Ashery et al., 2000, and references therein) at a distance of up to 200 nm from the membrane (in a volume of 1.6 × 10-13 liter). Accordingly, the concentrations of VAMP, synaptotagmin-I, and synaptotagmin⋆, which reside on the vesicle, were taken as identical to the vesicle concentration (10 nm). The concentration of syntaxin and SNAP-25, which are the most abundant proteins in the nervous system, was varied in the range of 0.1-100 μm (Lang et al., 2001). The concentration of Munc13, which is known to be low in chromaffin cells (Ashery et al., 2000), was searched for in the range of 0.1-10 nm. The concentration of complexin was varied in the range of 0.1-1 μm.
In the experiments that were simulated in this study, both in vitro or in vivo, the Ca2+ ions were associated with an infinite reservoir that was hardly depleted during the exocytotic process. Accordingly, in the model the Ca2+ up-take by the proteins was not considered, and the Ca2+ concentration was kept constant, as set by the experimental conditions.
The simulation procedure. Precise derivation of the equilibrium equations led to a set of 16 coupled, nonlinear ordinary differential equations. Such a system (Lindfield and Penny, 1995) is readily subjected to numeric integration, using the Runge Kutta methods as embedded in the Matlab program (version 6.5R13), using the ODE23S routine (Lindfield and Penny, 1995).
To fit the experimental results, the program was fed with the initial concentrations of the proteins, and the rate constants and the parameters were assigned to the reaction steps by a thorough search that was performed in the parameter space (between the ranges that were described above). The adjustable parameters were varied until the deviations between the reconstructed and experimental curves were minimal, as shown in the figures below. The equations were propagated in time, and the concentrations of the intermediates and the final products were saved as vectors. To quantitatively reconstruct the experimental signal, the fusion vector (the sum of reactions 10 and 15; see above), was multiplied by the capacitance of a single vesicle (1.25 fF, for mouse chromaffin cells) (Moser and Neher, 1997). Thus, the fusion process was reconstructed in discrete quanta, each representing the fusion of one vesicle. Please note that we used the same set of parameters for all the experiments described in this study.
In the present study, the search was performed “manually” in which the operator used his kinetic intuition and biochemical knowledge to fit the curves. Recently (our unpublished results), the search procedure was adapted for a computerized search (using genetic algorithm analysis) (Moscovitch et al., 2004). The results of the “unbiased” genetic algorithm search were subjected to statistical analysis, which indicated that the solution is the global minimum in the multi-dimensional parameter space. The results presented in this study are located within the global minimum. Full details of the statistical study will be published elsewhere (our unpublished results).
Results
The model system
The process of exocytosis is presented as a linear array of interactions between synaptic proteins with a single bifurcation, as detailed in Scheme 1. It starts with the formation of the binary SNARE complex, between SNAP-25 and syntaxin I, which is the first intermediate protein complex in the model (Lang et al., 2001; Fasshauer and Margittai, 2004; Rickman et al., 2004). The complex later binds synaptobrevin/VAMP II to form the mature ternary SNARE complex (Sutton et al., 1998; Weber et al., 1998; Margittai et al., 2003). At this step, the SNARE complex is sequentially activated, first by Munc13-1 and then by calcium, and forms the SNARE⋆ complex. A novel step in which Munc13-1 is preactivated by Ca2+ was added to the flow chart and accounts for the finding that Munc13 proteins have several calcium-binding domains (Brose et al., 1995; X. Z. Xu et al., 1998; Rhee et al., 2002). The binding of calcium to SNARE and to Munc13-1 represents the calcium-dependent priming processes in chromaffin cells (von Ruden and Neher, 1993; Voets, 2000). The next step in the maturation of the SNARE complex is the binding of complexin to SNARE⋆ to form the SNARE# complex (Reim et al., 2001; Chen et al., 2002; Marz and Hanson, 2002; Pabst et al., 2002). From that point, the process splits into two parallel branches. Path I interacts with synaptotagmin-I (Davis et al., 1999; Gerona et al., 2000; Chapman, 2002; Zhang et al., 2002; Rickman and Davletov, 2003) to form the first releasable complex (RC-I), whereas path II interacts with a second putative calcium sensor, synaptotagmin⋆, and forms the second releasable complex (RC-II). When the intracellular Ca2+ concentration ([Ca2+]i) increases, each path reacts successively with three calcium ions, and the vesicles fuse with the cell membrane (Voets et al., 1999; Voets, 2000; Sorensen et al., 2002, 2003). The bifurcation was introduced to account for the fact that, in chromaffin cells derived from synaptotagmin-I knock-out (SytIKO) mice, the fusion process still continued, albeit at a slower rate (Voets et al., 2001). At present, we do not know whether the fusion in the SytIKO cells represents the function of another synaptotagmin isoform (Sudhof, 2002) or a different calcium sensor protein. For this reason, in the second branch we used the general notation “synaptotagmin⋆.”
Reconstruction of calcium-triggered exocytosis
The first experimental protocol that we simulated was the flash-photolysis protocol in chromaffin cells. In this experiment, the cells are dialyzed with intracellular solution containing 280 nm Ca2+ for 2 min. After this prepulse equilibration, a flash of a UV light elevates [Ca2+]i to 10-50 μm for 5 sec, and secretion is measured continuously by membrane capacitance measurement (T. Xu et al., 1998; Ashery et al., 2000).
Chromaffin cells respond to the flash stimulation with a typical biphasic capacitance increase, composed of a phase of fast secretion (the exocytotic burst), followed by a slower, sustained phase of secretion (T. Xu et al., 1998; Voets, 2000). According to the macroscopic analysis, the increase in the measured capacitance is described by a sum of three exponential reactions, each representing the fusion of a different “pool” of vesicles (Neher, 1998; Ashery et al., 2000; Rettig and Neher, 2002). The vesicles in the readily releasable pool (RRP) fuse rapidly with the cell membrane with a fast time constant (∼20 msec) and contribute to the initial fast burst component. This is followed by a slower fusion of vesicles from the slowly releasable pool (SRP) that occurs with a slower time constant (∼200 msec) (Voets et al., 1999). The linear sustained components represent the priming and subsequent fusion of vesicles from the unprimed pool (Ashery et al., 2000).
The response of chromaffin cells derived from SytIKO mice to a flash stimulation is characterized by release dynamics that are reminiscent of the SRP. In these cells, the rapid burst component is missing, whereas the slow burst component is as in the wild-type (WT) cells (Voets et al., 2001).
The simulation program reconstructed the very same protocol. The intermediate protein complexes were allowed to propagate for 10 min at low [Ca2+] (50 nm) to simulate resting conditions. After this simulation period, the final concentrations of the intermediates were fed as the initial concentrations for 2 min simulation with an elevated Ca2+ concentration (280 nm) to simulate the prepulse equilibrium period. The concentrations of the intermediate complexes, as calculated at the end of the prepulse, were used as the initial values for the final simulation in which [Ca2+] was set to 30 μm (i.e., the Ca2+ pulse).
The reconstruction of the release kinetics (Fig. 1), as measured for the WT and the SytIKO cells, accurately reproduces both in amplitude and shape the secretion process with minimal deviations from the experimental signals. Moreover, the simulated curves reconstruct with high accuracy also the burst components of WT and the SytIKO (Fig. 1, inset). The values of the adjustable parameters needed for the reconstruction of the experimental data are given in Table 1.
Reconstruction of the secretion dynamics from WT and SytIKO chromaffin cells. The main frame presents the average increase in membrane capacitance recorded from chromaffin cells after a flash-photolysis experiment (black line) and the reconstruction curve obtained by the model (gray line). The inset depicts the experimental results (black lines) and simulations (gray lines) at a higher time resolution. The top curves correspond with measurements performed with WT cells, whereas the bottom curves correspond with measurements performed with the SytIKO cells. Please note that the reconstructed curves fall within the margins of the experimental noise. The y-axis indicates the fusion process in molar units of vesicles.
The reconstruction of the dynamics for the WT cells was performed according to the bifurcated model. To obtain the reconstructed curve for SytIKO cells, the only parameter that was changed was the concentration of synaptotagmin-I, which was set to zero. Thus, the model can reconstruct basic exocytosis profiles and deletion of synaptotagmin-I.
It is our experience from previous studies that the reconstruction of a single measurement can be attained by more than one combination of adjustable parameters (Moscovitch et al., 2004). To gain confidence in the solution, we tested the ability of the model to reconstruct signals that examine different aspects of the exocytotic process and were recorded under a variety of experimental protocols: refilling of the RRP, calcium-dependent priming and fusion steps, and the overexpression of Munc13-1. Despite the complexity of the different experimental procedure, the model reconstructs successfully these experiments with the same set of adjustable parameters given in Table 1.
Refilling process and the rapid recovery of the exocytotic burst
Experimentally, the response to a second stimulus that is delivered 2 min after the first stimulation displays similar release dynamics to the first response (Ashery et al., 2000). This scenario was successfully reproduced by the model, with no need to modify any of the adjustable parameters (data not shown). Thus, the model can reproduce not only the release dynamics but also the refilling process after the pulse.
It had been demonstrated that after the depletion of the RRP, a fast refilling occurs (Smith et al., 1998; Voets et al., 1999; Dinkelacker et al., 2000). In flash experiments, the initial refilling of the RRP was mirrored by a decrease in the slow burst component (SRP), and therefore it was suggested that the two pools of vesicles are found in dynamic equilibrium (Voets et al., 1999). In the present model, the equilibrium between the fast burst and slow burst components occurs at the bifurcation step, in which both synaptotagmin-I and synaptotagmin⋆ interact reversibly with the SNARE# state. To test this scenario, we simulated four successive short stimulations that deplete the fast burst component (mainly RC-I). We then applied a calcium pulse at different time intervals after the short stimulations and analyzed the kinetics of the response to the calcium pulse. Analysis of the recovery kinetics of the RC-I complex shows that the initial recovery of RC-I is mirrored by a decrease in the size of RC-II, and after 10 sec, both pathways (RC-I and RC-II) are refilled in parallel, as shown in Figure 2. This behavior reconstructs the recovery of the RRP and SRP pools after depletion of RRP (Voets et al., 1999). Thus, we suggest that the two branches of the model interconvert via the bifurcation step (SNARE#) in a similar way as was shown for the RRP and the SRP.
Reconstruction of the recovery of RRP. Reconstruction of the recovery process of the RRP in chromaffin cells after four short stimulations (100 msec each; data not shown). At given time intervals after the last stimulation, the cells were challenged by a Ca2+ pulse, and the fusion signal was analyzed for the contribution of the RC-I and the RC-II components (for details, see Voets et al., 1999). The results of the model are given as the molar concentrations of the RC-I pathway (triangles) and the RC-II pathway (circles) at given time intervals after the last short stimulation. Please note that the initial recovery of RC-I is mirrored by a decrease in the size of RC-II.
Ca2+ dependence of the rate constants for fusion
The release rate is known to depend steeply on [Ca2+], yielding a third-fourth power relationship (Dodge and Rahamimoff, 1967; Heinemann et al., 1994; Voets, 2000). We have simulated exocytosis under increasing concentrations of Ca2+ (between 10 and 50 μm) and examined the relationship between the simulated rate of release and [Ca2+]. Each curve was fitted with two exponents for WT simulation, as is usually done with experimental data, and with one exponent for SytIKO simulation (Voets et al., 2001). Analysis of the WT curves resulted in two rate constants, different from each other by a factor of 10 and with a similar dependence on calcium (Fig. 3, filled circles and filled squares). These curves reconstruct the calcium dependence of the experimental RRP and SRP (Voets, 2000). According to the model, these slopes correspond to the sequential binding of three calcium ions to the final, releasable complexes.
The effect of the Ca2+ concentration on the rate of fusion. Signals were generated by the model, at increasing concentrations of Ca2+ pulse, and analyzed by the curve-fitting procedure (Igor Pro; Wavemetrics, Lake Oswego, OR) to determine the rate constants of the simulated fast and slow burst components. The results are presented on a log-log scale and demonstrate the calcium dependence for the fast (circles) and the slow (filled squares) burst components for WT cells. Repeating the same calculations for the SytIKO cells yielded only one component, and its rate constants (open squares) fall close to the slow rate constants of the WT cells. If we used a linear sequence without the bifurcation step in the model, the rate constants of the slow burst component were insensitive to the Ca2+ concentration (triangles).
The analysis of the curves obtained for the SytIKO simulations revealed that the rate constants fall next to the lower curve of the WT cells (Fig. 3, open squares), which corresponds to the rate constants of the SRP. These results support the general hypothesis that vesicles in chromaffin cells can fuse from both the SRP and the RRP (Voets et al., 1999; Voets, 2000; Sorensen et al., 2002, 2003) and that fusion in SytIKO represents fusion from the SRP. In addition, it indicates that the macroscopic SRP corresponds mostly with the fusion of vesicles with a Ca2+ sensor that is not synaptotagmin-I.
Reconstruction of the same experiments while eliminating the bifurcation (i.e., assuming that only synaptotagmin-I is operating as the Ca2+ sensor in WT cells) failed to reproduce the dependence of the slow burst component on [Ca2+] (Fig. 3, open triangles). Thus, the bifurcation is an essential element in the reconstruction of the calcium-dependent release dynamics of the WT cells and the secretion measured with the SytIKO cells.
The main difference between the synaptotagmin-I and synaptotagmin⋆ pathways is the rate of fusion: RC-II fuses at a rate that is only 10% of that of RC-I (Table 1, reaction 10 and 15). Thus, the kinetic analysis suggests that the main requirement for the observed calcium dependence of the fast and slow components of release is a considerably slower fusion rate in pathway II.
Modulation of priming by prepulse Ca2+ and Munc13-1 overexpression
Elevation of prepulse [Ca2+] has been shown to modulate the size of the exocytotic burst in chromaffin cells, a phenomenon that is known as calcium-dependent priming (von Ruden and Neher, 1993; Voets, 2000). To simulate the priming, [Ca2+] during the prepulse equilibrium was adjusted to increasing concentrations (50-1500 nm) for 50 sec, and the size of the exocytotic burst was measured after a Ca2+ pulse to 30 μm. Similar to the experimental data (Fig. 4, open circles), the simulation (Fig. 4, filled squares) demonstrates a gradual increase in secretion when calcium is elevated at the low [Ca2+] range (<0.6 μm), representing calcium-dependent priming, and a decrease in secretion at high prepulse [Ca2+](>0.6 μm).
Modulation of priming by Ca2+. The figure depicts the dependence of the burst amplitude, as measured 1 sec after the Ca2+ pulse, on the prepulse calcium. The results of the simulation (squares) show a similar pattern to the experimental data (circles) (Voets, 2000). The y-axis indicates the fusion process in molar units of vesicles.
One of the advantages of the mathematical model is its ability to predict the effect of modulation of the level of specific proteins. To test this possibility, we compared the amount of exocytosis under different levels of Munc13-1, a presynaptic protein that facilitates priming in chromaffin cells. Overexpression of Munc13-1 in chromaffin cells led to a significant increase of exocytosis in response to a stimulation (Fig. 5) (inset from Ashery et al., 2000). When the model was simulated with an increased Munc13-1 level, the outcome of the simulation reproduced the experimental observations without a need to modulate any of the other parameters (Fig. 5). Thus, in practice, the model “predicted” the effect of the Munc13-1 overexpression and can be used to predict the outcome of manipulations of other components of the fusion machinery.
The effect of Munc13-1 overexpression on exocytosis. The main frame depicts the release dynamics as calculated in the model for normal Munc13-1 (0.28 nm; bottom curve) and for high Munc13-1 (2.8 nm; top curve) concentrations. The y-axis indicates the fusion process in molar units of vesicles. The inset depicts the experimental observation of overexpression of Munc13-1 in chromaffin cells (Ashery et al., 2000).
The dynamics of the intermediates during the fusion process
Another advantage of the mathematical model is that, during the reconstruction of the fusion dynamics, it generates the corresponding dynamics of all the intermediate complexes in the process. An example of the dynamic changes in the intermediate complexes after a calcium pulse is shown in Figure 6. This kind of analysis enables us to study the contribution of each intermediate to the different phases of fusion and predicts the different scenarios emerging from deletion manipulation of the fusion machinery components. During the pre-equilibrium period, ∼120 vesicles have matured to the releasable complexes RC-I and RC-II and to the SNARE# complex (Fig. 6e-g). Please note that the ordinates in Figure 6 indicate the number of vesicles that had matured to a given stage in the process. Careful examination of the kinetics of the intermediates revealed that the initial burst, that the macroscopic analysis equates with the RRP, is related to the depletion of the RC-I complex through its intermediates (Fig. 6g,i,k,m,p). The slow release phase is mostly associated with the release of the RC-II complex through its intermediates (Fig. 6e,h,j,l,n) and might represent the macroscopic SRP. In addition, the maturation of the SNARE# complex during the stimulation slightly contributes to the slow burst component in the model. The sustained component is controlled by the upstream reactions that cause a steady consumption of the SNARE and SNARE# populations (Fig. 6a,c,d,f). In the present model, the SNARE intermediate that represents the number of vesicles with preassembled SNARE complexes consists of ∼840 vesicles. Because the precise sequence of events that regulates the formation of the SNARE complexes (modulation by tomosyn or Munc18, formation of monomeric SNARE by soluble N-ethylmaleimide-sensitive factor) is still an open problem, these steps are not simulated. Therefore, all the upstream interactions are represented by the formation of the SNARE complex (step 2 in Scheme I), leading to an overestimation of the actual number of the SNARE complexes. We assume that in vivo, these regulatory steps down-regulate the number of preassembled SNARE complexes. In the future, we intend to implement these upstream steps into the model that will then improve our resolution of the dynamic interactions in early stages of exocytosis.
Detailed dynamics of the individual intermediate protein complexes during the fusion reaction. The variation in the concentrations of the intermediate protein complexes is given as the number of vesicles that matured to the state indicated in each frame. The concentration of Munc13-1 is given in nanomoles. The simulation represents the last second of the prepulse period (120 sec at 280 nm Ca2+; gray lines), and at t = 1, the Ca2+ concentration was increased to 30 μm (the calcium pulse). The dynamics of the intermediate complexes are shown during the 5sec stimulation period (black lines). The frames that appear in a-o are in the same order and locations as the reactions appear in Scheme 1, starting from the formation of the SNARE complex (step 2).
Discussion
In the present study, we applied classical chemical dynamics for describing the complex intracellular process of vesicle exocytosis. The advantage of this formalism is the replacement of a macroscopic description, based on curve fitting, by a precise mechanism based on ordered reactions between proteins, in which the rate constants of these reactions are suitable for quantitative interpretation. The model is capable of accurate reconstruction of various experimental protocols and of reproduction of the calcium-dependent priming process. Moreover, when the simulated curves are analyzed by the curve-fitting program used for analysis of experimental data, the dependence of the time constants on the Ca2+ concentration match those derived from real measurements. Finally, the model can reproduce the effect of selective mutations on exocytosis, such as overexpression of Munc13-1 or deletion of synaptotagmin-I.
This new method of analysis provides a novel platform to study the molecular mechanism of exocytosis and appears to bear a high predictive power, which brought about the following conclusions: (1) In WT cells, the secretion proceeds through two parallel pathways. The rapid one utilizes synaptotagmin-I as the Ca2+ sensor and corresponds mostly to the macroscopic RRP. The second pathway utilizes another calcium sensor (synaptotagmin⋆) and might be related to the SRP. In SytIKO cells, only the synaptotagmin⋆ pathway operates. (2) The kinetic analysis indicates that the main difference between the pathways is in the final step of the fusion process. (3) The so-called “priming step” is composed of sequential steps of protein-protein interactions, consisting of the association of the SNARE complex, its initial activation by Munc13-1 and by Ca2+ ions, and finally the binding of complexin and the calcium sensors.
The model supports the hypothesis that fusion occurs from two populations of primed vesicles (Voets et al., 1999; Voets, 2000; Sorensen et al., 2002, 2003) and provides information regarding their molecular nature. The model suggests that the fusion rate from RC-II is much slower than that of RC-I, similar to the differences between the fusion rates of the SRP and the RRP (Voets, 2000; Sorensen et al., 2003). The slower rate of fusion of RC-II can be attributable to different binding kinetics of the RC-II complex to the plasma membrane or some intrinsic properties that slow down its fusion. Additional experiments will be needed in the future to clarify this phenomenon.
The accurate reconstruction of several experimental signals, measured under a variety of initial conditions, implies that the whole sequence of events is properly modeled. Thus, the model can now shed light on the intermediate dynamics of the process, events that until now were concealed in a “black box.” Examination of the dynamic changes of the intermediate protein complexes allowed us to correlate between the various kinetic components of secretion and the protein complexes that might contribute to these components. Furthermore, the model provides a tool to examine the effect of various manipulations and to predict their outcome. However, it must be stressed that the model is open to refinements while still retaining its accuracy and predictive power and can be expanded to include more reaction steps that, at present, are experimentally unresolved.
The model can also accommodate subtle alternations in the sequence of events; for example, the precise location of the bifurcation step, an essential mechanistic element of the present model, is not fully established. Presently, the complexin binding precedes the interaction with the Ca2+ sensor. Yet, the binding of complexin after the splitting point is, model-wise, kinetically feasible. The presently published information does not permit defining the precise order of events at this step. Once more data are available, additional refinement can be made.
The model combines 30 rate constants for reconstruction of the sequential interactions between seven synaptic proteins and calcium ions. Of these, 19 rate constants and three protein initial concentrations have already been independently measured or calculated in previous studies (the published parameters are marked in Table 1). Thus, the number of unknown parameters that had to be determined for fitting the experimental curves is only 15. A systematic search over this multidimensional parameter space (our unpublished observations) using a genetic algorithm (Moscovitch et al., 2004) has indicated that the fitting of the experimental observations operated as a very selective filter so that only one combination survives. Thus, for the presented sequence of events, only one combination of parameters describes all experimental manipulations that are discussed in this study. This combination is presented in Table 1.
The rate constants, as given in Table 1, are the best-fitting values derived both by manual fitting to the experimental results and independently searched by a thorough genetic algorithm analysis (our unpublished observations). It should be pointed out that the rate constants associated with the formation of the SNARE complex and the binding of Munc13-1 vary over a wide range, whereas the rate constants of the later steps are much more accurate, and their range of variation is less than one order of magnitude (data not shown). These initial steps will be further characterized in the future as more information becomes available.
In this model, we neither included upstream steps of vesicle docking and mobilization nor those steps that are related to the regulation of monomeric SNARE proteins. In addition, the scenario described in this study was optimized to describe exocytosis of large, dense core vesicles from chromaffin cells, and secretion from other systems may be different. However, because the molecular mechanisms of exocytosis are highly conserved in evolution, the model can serve as a framework for optimization of other systems. We assume that the fitting of secretion kinetics from other systems will mainly demand modulation or addition of regulatory steps to the sequence.
Finally, we want to stress that the present model consists of the minimal steps needed to reconstruct the fusion process. All the steps that were implemented in the model are essential and consistent with the present knowledge of the proteins involved and their sequence of reaction. It should be noted that most of the rate constants described here are slower than the estimated diffusion-controlled values. In these cases, we assume that the actual reaction step consists of more than a single step. The hidden steps can either be a conformational change of the protein, a reaction of the indicated complex with another (still unidentified) protein(s) or lipid(s), or even oligomerization of the complex/protein. These equations can therefore be used as a basis for additional refinement when new components and kinetic features are identified. Thus, the model serves as a dynamic platform that can gain additional accuracy without losing the fundamental features that are already implemented in the basic reactions.
Appendix 1
The estimation of the upper limit of a diffusion-controlled reaction (kmax) is given by the following Debye-Smoluchowski equation: 
which correlated the rate constant with the radius of encounter (∑ Rj,i) and the diffusion coefficients of the reactants ∑ Dj.i, where Nav is Avogadro number.
The diffusion coefficients of soluble proteins were approximated by the Einstein Stokes equation. The diffusion coefficient of protein in membrane was taken as 10-9 cm2sec-1.
The multiplicity of identical sites on the same vesicle increases the probability of a fruitful encounter of a ligand with any of the identical sites. Thus, a vesicle that carries a few VAMP molecules on its surface (Walch-Solimena et al., 1995) will have a higher probability of encounter with the binary SNARE complex than one carrying only one copy of VAMP. This problem was treated by Berg and Purcell (1977). According to their formulation, the rate of encounter k(N) between the free reactant with any of the N identical sites (each having a radius r) located on a sphere with a radius R is given below, where kmax is given by the Debye-Smoluchowski equation: 
Thus, 10 copies of VAMP II (Walch-Solimena et al., 1995), each having an effective gyration radius of ∼10 nm, will lead to a rate of encounter that is 25% of the maximal rate, the one predicted for a vesicle that is fully covered by VAMP II molecules. The advantage of this treatment is that it avoids the necessity of defining the number of identical copies of a reactant on a vesicle, because the number of copies is incorporated into the adjustable parameter (i.e., the rate constant).
Footnotes
This work was supported by grants from The Israel Science Foundation (424/02-16.6), the Minerva Junior Research Group, and the Ministry of Absorption to U.A.; and from The Israeli Science Foundation (472/01-2) and the American Israel Binational Science Foundation (2002129) to M.G. We thank Dr. Thomas Voets for providing the raw data for some experiments and Prof. Erwin Neher and Ofer Yizhar, Uri Nili, and Adam Grundland for comments on previous versions of this manuscript.
Correspondence should be addressed to either of the following at the above address: Uri Ashery, E-mail: uria{at}post.tau.ac.il; or Esther Nachliel, E-mail: Eti{at}hemi.tau.ac.il.
Copyright © 2004 Society for Neuroscience 0270-6474/04/248838-09$15.00/0
