Analysis, Simulation, and Optimization of Stochastic Vesicle Dynamics in Synaptic Transmission

Synaptic transmission is the mechanism of information transfer from one neuron to another (or from a neuron to a muscle or to an endocrine cell). An important step in this physiological process is the stochastic release of neurotransmitter from vesicles that fuse with the presynaptic membrane and spill their contents into the synaptic cleft. We are concerned here with the formulation, analysis, and simulation of a mathematical model that describes the stochastic docking, undocking, and release of synaptic vesicles and their effect on synaptic signal transmission. The focus of this paper is on the parameter p0, the probability of release for each docked vesicle when an action potential arrives. We study the influence of this parameter on the statistics of the release process and on the theoretical capability of the model synapse in reconstructing various desired outputs based on the timing and amount of neurotransmitter release. This theoretical capability is assessed by formulating and solving an optimal filtering problem. Methods for parameter identification are proposed and applied to simulated data. © 2019 Wiley Periodicals, Inc.


Motivation and Overview
Neurons form complex networks via synapses through which information propagates. In this article, we consider chemical synapses, in which neurotransmitters are involved (see Figures 1.1 and 1.2). Information transfer at chemical synapses occurs in three main steps [11,17,24,33,43,49,64,67,80]: (1) At the presynaptic terminal: Undocked vesicles bind to docking sites (release sites) and become docked; simultaneously, docked vesicles can become undocked without releasing their neurotransmitter through reversal of the docking process. The arrival of an action potential (nerve impulse) leads to increased membrane potential, which opens Ca 2C channels. Ca 2C stimulates the fusion of docked vesicles to the presynaptic membrane, spilling neurotransmitter molecules into the cleft. 1 (2) In the synaptic cleft: Neurotransmitter molecules bind to receptors, which are ion channels located on the postsynaptic membrane. Channels with bound neurotransmitter can open. Neurotransmitter action is terminated by enzymatic degradation, uptake into the presynaptic terminal, and diffusion out of the cleft. (3) In the postsynaptic neuron: Ionic currents flowing through the open synaptic channels displace the membrane potential. Depending on the channel type, the change may be either excitatory or inhibitory to the postsynaptic neuron.
Unlike the all-or-none action potential, synaptic transmission is graded, since the number of vesicles released by the arrival of an action potential, and also the postsynaptic response to each vesicle released, may vary. The synapse is therefore a favorite site of hormonal, pharmacologic, and neural regulation of nervous activity. Vesicle fusion and the subsequent release of neurotransmitters is stochastic and its likelihood of occurrence is a crucial factor in the regulation of signal propagation in neuronal networks [9,21,29,81]. The reliability of neurotransmitter release can be highly variable: since the 1950s, experimental data from electrophysiological, molecular, and imaging studies have demonstrated that synaptic terminals can individually set their neurotransmitter release probability dynamically through local feedback regulation [9], and also that stochastic vesicle release is the most significant source of noise in the central nervous system [2,18] It is widely believed that the synapse is the site at which learning takes place and at which memory is stored [11,61]. Work during the last half century has shown that modification of the rate of neurotransmitter release contributes to both short-term [1,23,34,47,48,79,81] and long-term changes [8,30,40,46,54,63] at synapses. Recently the rate of neurotransmitter release has also been linked to severe neurological disorders, such as Parkinson's disease [45,68] and Alzheimer's disease [53,77]. If long-lasting synaptic adaptations are indeed a mechanism by FIGURE 1.1. Docking, undocking, and release of synaptic vesicles. At a chemical synapse, one neuron influences another through the release of neurotransmitters, which are small molecules packed inside synaptic vesicles. The arrival of an action potential at the presynaptic terminal triggers the fusion of the membranes of some docked vesicles with the membrane of the presynaptic neuron, leading to the release of neurotransmitters into the synaptic cleft. The binding of one or more neurotransmitter molecules with a postsynaptic receptor triggers the opening of an ionic channel in the postsynaptic neuron, and this may either raise or lower the postsynaptic membrane potential depending on the channel type. The two black dots associated with each docked vesicle represent the cross-section of a protein ring that defines a docking site. which our experiences get translated into memories, a quantitative understanding of how various factors in synaptic transmission determine the rate of vesicle release is crucial to the understanding of the brain.
From a broader perspective, synaptic vesicle release is a form of exocytosis, a biological process through which molecules produced in a cell are secreted to the extracellular environment. Other forms of exocytosis include the release of insulin to regulate blood sugar level [74] and the release of cytotoxic granule vesicles from immune cells to attack foreign antigens [69]. Although we only consider synaptic vesicle release in this article, it is possible that the model we introduce and study herein may be more broadly applicable.

Existing Models of Synaptic Vesicle Release
The original model of synaptic vesicle release proposed by Katz [26], based on binomial statistics, is widely cited and is still being used today. Suppose there are n s independent docking sites, all of which are occupied, and that the probability that a vesicle undergoes exocytosis following the arrival of a nerve impulse is p 0 ; FIGURE 1.2. Electron microscope cross-sectional images of two synapses of cortical neurons in the mouse brain. The dark (electron dense) segment of the cell membrane corresponds to the location of docking sites, where there is a high concentration of presynaptic proteins that tether synaptic vesicles to the presynaptic membrane and mediate synaptic vesicle fusion. Docked vesicles are those vesicles located near the dark segment of the membrane. Some docked vesicles are indicated by arrows. Note that a cross-sectional image does not show vesicles that are outside the plane of the cross-section, so there may be many vesicles, including docked ones, not seen in the image. These two images are the subfigures located in the upper-left corner of figures 1A and 1B, respectively, in Wu et al. [75], used under the Creative Commons Attribution 4.0 International Public License, with arrows added to the original. then the mean of the number N of vesicles released is n s p 0 , the variance of N is n s p 0 .1 p 0 /, and the probability that k vesicles are released, Pr.N D k/, is given by Pr.N D k/ D n s Š kŠ.n s k/Š p k 0 .1 p 0 / n s k : Robinson [56] provided a set of formulae for the estimation of p 0 and n s , assuming a fixed postsynaptic response to each vesicle release event.
Katz's binomial model enjoyed great success in interpreting vesicle release data from the frog neuromuscular junction, but one crucial assumption weakens its applicability to other synapses: in his binomial model, it is assumed that prior to the arrival of each action potential all of the docking sites are filled, i.e., n s is equal to the number of docked vesicles when an action potential arrives. This assumption is not accurate in general. Several studies have reported that the number of docked vesicles prior to each action potential is variable [6,10,36,52,72]. To overcome this shortcoming, Vere-Jones [71] proposed and then carefully analyzed a model synapse with an unlimited number of docking sites, in which vesicle docking (and undocking) occurs by a homogenous Poisson process; this model was consistent with several early experimental results indicating that the numbers of vesicles released seemed to be governed by the Poisson distribution, but was nevertheless not widely accepted because of the idealized assumption of an unlimited number of docking sites. Barrett and Stevens [4,5] adopted a different approach to extend Katz's original model: they assumed that vesicle release at each docking site occurs by a Poisson process with a time-dependent rate.
In recent years, several workers have applied the extended Katz theory to neurotransmitter release at central synapses, such as those in the cerebral cortex and the hippocampus (see [15,70] for two examples). Evidence has also accumulated to indicate that in many synapses the statistics of vesicle release does not follow a Poisson distribution [51,78,80]. The Poisson assumption is only a good approximation when the rate of nerve impulse arrival is sufficiently high. Attempts to loosen the Poisson assumption led to the development of models of vesicle pool dynamics [14,38,44,55,70]. Under this framework, the arrival of each nerve impulse triggers only the probabilistic release of vesicles in the readily releasable vesicle pool (RRP) [25,59], and a set of deterministic differential equations is used to describe the replenishment of RRP from recycle and reserve vesicle pools.
A challenge in mathematical modeling of biological systems is in the design of models that accommodate not-yet-understood biological components or can account for dynamic changes in biological parameters. Some recent work in this direction include the phenomenological models proposed by Maass and Zador [29] and Bird, Wall, and Richardson [7], which allowed temporal variations of the release probability as a means to achieve short-term facilitation (and possibly a form of depression that is unrelated to depletion of docked vesicles). Another phenomenological model proposed by Fellous and Corral-Frías [16] studied the effect of heterogeneity in the initial release probability on the reliability and precision of postsynaptic output. In Merkel and Lindner [39], the authors modeled a population of facilitating and depressing synapses in which the spike trains from different presynaptic neurons are independent of each other and are Poisson. They derived analytical approximations for cross-spectra, power spectra, and the coherence function between the presynaptic and postsynaptic signals, which showed that the synaptic coherence function is largely independent of frequency. This result is consistent with the findings by de la Rocha, Nevado, and Parga [12], in which the authors found that short-term depression does not affect the efficacy of synaptic transmission for Poisson spike trains, but it can enhance synaptic transmission for more realistic temporally correlated spike trains. Using paired-cell recording in rat neocortex, Scott, Cowan, and Stricker [60] showed that short-term plasticity can not only improve but also reduce information transfer between neurons. Specifically, they showed that when information is contained in the timing of individual spikes, short-term plasticity affects information transfer relative to its impact on the probability of neurotransmitter release, whereas when information is contained only in the mean spike rate, the effect of short-term plasticity depends on the range of spike frequencies that the target network can distinguish. In Fuhrmann et al. [19], the authors analyzed the optimal frequency of presynaptic spike train for which the information content is maximal in a model synapse that includes both depression and facilitation. Another interesting line of work by Meinrenken, Borst, and Sakmann [37] and Modchang et al. [41] aims to provide a quantitative link between spatial heterogeneity in vesicle release and residual calcium dynamics in the presynaptic terminal.
Finally, we would like to refer the reader to some related work that also studied the functional roles of stochastic vesicle release in synaptic transmission. In Arleo et al. [3] and Goldman [20], an information-theoretic approach was used to examine the effect of release probability on synaptic transmission of information. Using Shannon's mutual information, the authors found that, in cerebellar granule cells and simulated model neurons, as release probability increased, the optimality of information transfer of most stimuli did not increase monotonically but instead reached a plateau at intermediate release probability levels. Their results are consistent with the results of this article. In Rosenbaum, Rubin, and Doiron [57], the authors constructed both a deterministic model and a stochastic model to study how variability in vesicle dynamics affects signal transmission. They found that the depletion of docked vesicles at higher rates of arrival of action potentials makes a stochastic synapse act as a high-pass filter, whereas a deterministic synapse that ignores the stochastic release of vesicles transfers information encoded at any frequency equally well. Note that the deterministic model considered by Rosenbaum et al. is a continuous, deterministic approximation to their stochastic model, and therefore is different from our vesicle release model proposed in this article (see Section 2). In our model, vesicle dynamics is stochastic even when the release probability is 100% because vesicle docking (and undocking) is always governed by a stochastic process; thus the number of vesicles released at a spike with a 100% release probability is equal to the number of vesicles that happened to dock during the immediately preceding interspike interval.
Manwani and Koch [31] also adopted an optimization point of view concerning stochastic vesicle release. In their model, each synapse has a single docking site that is always filled, both the probability of release per docked vesicle p 0 and the postsynaptic response are predetermined based on experimental observation of the properties of cortical synapses. They found that single stochastic synapses cannot transmit presynaptic spike density S.t / reliably, but redundancy obtained using a small number of multiple synapses leads to a significant improvement in the reconstruction of S.t /. Our work is very different, since (1) we consider a model synapse with an arbitrary number of docking sites with docking and undocking processes, (2) we do not prescribe a fixed p 0 but instead study how different values of p 0 affect the fidelity of synaptic transmission, and (3) we do not presume the biophysical details of a synapse's postsynaptic response to the presynaptic vesicle release events but instead use optimal linear filter theory to measure the capability of our idealized postsynaptic neuron in the estimation of any desired signal derived from the presynaptic spike density S.t /, given that we know the statistics of the signal ensemble from which S.t / has been drawn. Our theory of synaptic vesicle docking, undocking, and release has allowed us to show that an appropriately chosen value of p 0 improves synaptic transmission by lowering the mean square error in the estimation of various desired signals, and how the best choice of p 0 is determined by various factors.

Our Prior Work
In a recent paper [76] we considered an idealized model synapse with an unlimited number of vesicle docking sites and no undocking, in which we assumed that vesicle docking occurs by a homogeneous Poisson process with mean rate˛0, that presynaptic action potentials arrive by a stochastic process with mean rate s.t / > 0, and that each vesicle that is docked has a probability p 0 to be released upon the arrival of each action potential, independently of other docked vesicles. In this idealized case, we found that the expected rate of vesicle release r.t / is governed by d dt This implies that a stimulus at any steady level s.t / D const leads eventually to the same vesicle release rate equal to˛0, the mean vesicle docking rate. This complete insensitivity to the absolute level of stimulation is a consequence of the assumption that there is an unlimited number of docking sites available and is consistent with several experimental observations [1,22,32,81] (see also [47] for a review). We also examined the theoretical capability of a synapse in the estimation of desired signals using information from the stochastic vesicle release events under the framework of optimal linear filter theory: we found that a small p 0 , such as 0.1, reduces the error in the reconstruction of the input, or in the reconstruction of the time derivative of the input in comparison to a larger value of p 0 such as p 0 D 0:5 or 1. This implies that the probabilistic nature of synaptic vesicle release can play a beneficial role in synaptic transmission.

Main Results
In this paper, we consider a model of stochastic vesicle release that is characterized by four parameters: the number of docking sites, n s ; the rate (i.e., probability per unit time) of vesicle docking at each empty site,˛; the rate of undocking for each filled site,ˇ; and the probability of release, p 0 , when an action potential arrives, of each vesicle that is docked at that time. The input to our model synapse is a sequence of action potential arrival times, and the output of the model presynaptic terminal is a sequence of random nonnegative integers, each of which is the number of vesicles released by the corresponding action potential. Conditioning on the sequence of action potential arrival times, we derive and solve a recursion relation for the expected numbers of vesicles released, and also a correlation function that partially characterizes the statistics of vesicle release. Then we adopt the point of view that the action potentials themselves are generated by a stochastic process and are carrying information about an underlying continuous signal, and we ask to what extent that signal can be reconstructed by linear filtering of the time series of numbers of vesicles released. We address this question both analytically and numerically. In the analytic case, we make simplifying assumptions that are not needed when the problem is tackled numerically. In both cases, we focus on the choice of the parameter p 0 , and we find that the quality of the best signal reconstruction that can be done depends on this choice. Roughly speaking, the result is that p 0 should be equal to 1 when the effective number of docking sites is small, but p 0 should be small when the effective number of docking sites is large. The latter case is interesting, since it implies that randomness in vesicle release can be helpful for signal preservation during synaptic transmission. The terminology "effective number of docking sites" in the foregoing refers to the influence of the undocking process in setting an upper bound that is smaller than n s on the expected number of docked vesicles. The optimal choice of p 0 is also influenced by other parameters such as the rate of arrival of action potentials. We conclude by showing how the parameters of the model can be identified from experimental data, and also how the model can be tested experimentally.

A Model of Synaptic Vesicle Release for Spiking Neurons
Throughout this paper, we make the following assumptions: (1) The input to a synapse is a sequence of action potential arrival times Later, we will use the capital letter : : : T k : : : when we consider action potential arrival times that are generated by a stochastic process. When we use the lowercase : : : t k : : : , we are assuming that these action potential arrival times are given. (2) The synapse has some number n s of equivalent vesicle release sites (or docking sites). Any particular site may be occupied or unoccupied by a synaptic vesicle.
(3) Between action potential arrival times, every unoccupied site has a probability per unit time˛of becoming occupied, and every occupied site has a probability per unit timeˇof becoming unoccupied. Thus, between action potential arrival times, each site obeys the reaction scheme 0˛ * ) ˇ1 in which 0 denotes an unoccupied site and 1 denotes an occupied site. The changes that occur at one site are independent of those occurring at any other site. Note that the reaction with rate constantˇdoes not involve the release of neurotransmitter into the synaptic cleft. Instead it is simply the undocking of a previously docked vesicle, without neurotransmitter release [43]. (4) At each action potential arrival time, every site that is occupied immediately before the action potential arrival time has the possibility of releasing the contents of its vesicle into the synaptic cleft and thereby becoming an unoccupied site. The probability that such release occurs at any particular site is denoted by p 0 , and the decision whether to release the vesicle or not is made independently for each site. p 0 is also known as the vesicle fusion probability. Let D.t / be the number of docked vesicles at time t, and let N k be the number of vesicles released by the arrival of the k th action potential. At any given time t between action potential arrival times, D.t / changes in steps of˙1, and the probability per unit time that D.t / increases by 1 is and then, of course, We regard the sequence .N k ; T k / as the output of the synaptic vesicle release process (i.e., the output of the presynaptic terminal).
3 Synaptic Vesicle Release with a Finite Number of Docking Sites 3.1 S N k , the Expected Number of Vesicles Released at Each Action Potential Conditioned on the Action Potential Arrival Times The first result of this paper is a recursion formula for the expected number of vesicles released at each action potential (spike) time given that we know the action potential arrival times.
Conditioned on the action potential arrival times t k , we denote by . / the expected value of . / conditioned on these action potential arrival times. We have 1 e .˛Cˇ/.t t k 1 / /; which holds for t 2 .t k 1 ; t k /: In particular, at t D t k , we have where we used (3.2).
Finally, multiply both sides by p 0 and make use of (3.1) to obtain the recursion formula then the expected number of vesicles released at each action potential, conditioned on the action potential arrival times ft k g, is given by the recurrence We call n s the effective number of docking sites.
We can use equation (3.8) to express S N k in terms of S N i for any i < k. Multiplying both sides of (3.8) by the summation factor e t k =.1 p 0 / k , we obtain This gives, for any i < k, .1 p 0 / k j e .t k t j / .1 e .t j t j 1 / /: (3.10)

The Autocovariance of N k Conditioned on the Action Potential Arrival Times
Denote by ' i k the autocovariance of N k : We first consider a synapse with a single docking site and possibly with undocking allowed. This should yield results that are useful in the general case, since the individual sites are independent of each other, provided that we condition on the arrival times of the action potentials t k .
With a single site, D s .t / D 0 or 1, and likewise, N s k D 0 or 1, in which the superscript s is the index of the docking site (s D 1 in this case). Thus, we can re-interpret S D s as the probability that D s .t / D 1 and N s k as the probability that N s k D 1: The autocovariance of the sequence of random variables N s k is defined by Since the definition of ' s i k is symmetrical under interchange of .i; k/, there is no need to consider separately the case i > k.
Note that Thus, to obtain ' s i k we only need to find Pr.N s k D 1 j N s i D 1/. Equations (3.1), (3.2), and (3.3) for x D and S N k now become N s k D p 0 D s t k ; (3.14) 16) in which, as before, D˛Cˇ, and we have now introduced the notation p DCˇ: Note that p is the same as n s in the special case that n s D 1 (see (3.7)), but we have introduced a special notation for this to avoid confusion when we later consider n s > 1, and also to emphasize that p is a probability. Specifically, p is the probability that a site is occupied if there have not been any action potentials for a long time. In particular, if there is no undocking, thenˇD 0 and p D 1.
Substituting (3.16) into (3.15), we get a recursion relation for D s t C k (instead of the recurrence for x D t k derived previously) by This differs from (3.4) by the appearance of .1 p 0 / in both terms instead of just the first term. To solve the above recurrence we multiply both sides by the summation factor e t k =.1 p 0 / k to obtain It follows that for any i < k Now suppose it is known that a vesicle was released by the i th action potential, that is, N s i D 1: Then it must be the case that D s .t C i / D 0: Conditioning on this, (3.18) gives, for i < k, .1 p 0 / k j e .t k t j / .1 e .t j t j 1 / /: Combining (3.14) and (3.15), we have Recall that N s k is the probability that N s k D 1, and in this case we have done the calculation conditioned on N s i D 1. Therefore, we have the result that for i < k Moreover, by letting i ! 1 in (3.19), we find Substituting (3.19) and (3.20) into (3.11) and (3.12), we obtain the following result, valid for all .i; k/: where ı i k is the Kronecker delta.
Note that since all of the sites are statistically identical, the superscript s refers to any one of the sites, and therefore it also serves as a convenient label for a singlesite quantity. The results involving this superscript do not actually depend on s (unless we are concerned with a random variable, such as N s k , which does depend on s, even though N s k does not). Now we are ready to scale up to the case of n s independent sites. We have where N s k is the number of vesicles (0 or 1) released at the s th site by the k th action potential. Note that the whole sequences N s 1 k and N s 2 k are independent for s 1 ¤ s 2 , provided that we condition on the spike times. Also, these sequences have the same probability distribution. It follows that since N s k is independent of s. Also, and finally Substituting (3.21) into (3.23) we obtain the formula of the autocovariance of N k , conditioned on the action potential arrival times ft k g, valid for all .i; k/ where ı i k is the Kronecker delta.

Example: A Regular Spike Train
Here we consider the special case of a regular spike train. Results from this example will be used later for parameter estimation in Section 6. We would also like to refer the reader to Matveev and Wang [35], in which the authors derived analytic results for a regular spike train under a different hypothesis that at most one vesicle can be released per spike.
Suppose t k t k 1 D t ; then (3.8) has the steady-state solution given by This shows that N k ! x N geometrically with ratio This observation allows us to construct a closed-form solution of the recurrence (3.8) for the special case of a sudden change in the rate of arrival of action potentials. Conditioned on the action potential arrival times ft k g, where the expected number of vesicles released at the time of the k th action potential is  We can re-express the foregoing in terms of the rate of arrival of action potentials and the rate at which vesicles are releasing their content of neurotransmitter by making the definitions so that x R.s/ is the steady-state rate of release of vesicles when the rate of arrival of action potentials is constant and equal to s.
In terms of these variables, (3.26) becomes where n s D˛n s is the mean rate of vesicle docking of the entire synapse when all docking sites are empty. That is, for sufficiently fast arrival of action potentials, the steady-state rate of vesicle release does not depend on the rate of arrival of action potentials, since it is limited by the maximum rate at which docking can occur. This shows that in the regime of fast arrival of action potentials, a synapse with a finite number of docking sites behaves in a similar manner to a model synapse with an unlimited number of docking sites [76]. On the other hand, as s ! 0, we have x R.s/ p 0 n s s: That is, for sufficiently slow arrival of action potentials, the steady-state rate of vesicle release of a synapse with a finite number of docking sites is proportional to the rate of arrival of action potentials. Furthermore, in the special case of a regular spike train ft k g, where for all k, the autocovariance of N k , given by (3.24), simplifies to where ı i k is the Kronecker delta function, and x N .t / is defined in equation (3.27). Figure 3.2 illustrates ' i k as a function of .k i /. Note that the height of the central peak x N .
x N / 2 =n s can be larger or smaller than the amplitude of the negative tails, which we may define by extrapolation to jk i j D 0, and is simply . x N / 2 =n s (in magnitude). The ratio of these two magnitudes is which is greater than 1 if x N < 1 2 n s and smaller than 1 if x N > 1 2 n s . (Note that r cannot be negative, since the mean number of vesicles released by an action potential cannot exceed the number of vesicle docking sites.) The above formula for r can be used as a check for parameter fitting (discussed in Section 6), so it is important to note that the value of r can be determined from data on ' i k . The numerator is simply ' kk , and the denominator is obtained by fitting the negative tails of ' i k to a geometric sequence.
Note that if n s is large, the negative tail of the autocovariance will be undetectable. This means that, in the limit n s ! 1, the random variables N i and N k are uncorrelated for i ¤ k. This observation is consistent with the limiting case discussed below of an unlimited number of docking sites, since N i and N k are in that case independent for i ¤ k if we condition on the spike times. 2. An illustration of the autocovariance ' i k . We plot two versions of ' i k : the theoretical result, as predicted by equation (3.30), is plotted as blue circles; the stochastic simulation result is plotted as red crosses. The negative tail, obtained by fitting the values of ' i k over i ¤ k to an exponential function using the least squares method, is plotted in the red dashed curve. Model parameters are t D 0:1 sec (interspike interval), p 0 D 0:5, n s D 100,˛0 D 1000 sec 1 ,ˇD 0.
To obtain the stochastic simulation result, we first compute N k numerically by simulating the stochastic vesicle release process for an evenly spaced spike train that lasts 100 sec. Once N k is obtained, we compute its autocovariance ' i k . To improve the accuracy of ' i k , we repeat the above process 1000 times and then take the average of ' i k over those 1000 sample paths. See the online supplement for the algorithms used in the stochastic simulation of N k . Recall that N k is the number of synaptic vesicles released by the arrival of the k th action potential. The autocovariance of N k in (3.24) shows that, as n s ! 1, the random variables N i and N k become uncorrelated for i ¤ k. This suggests that N k are independent in a model synapse with an unlimited number of docking sites; this is indeed true, as proven below for arbitrary spike trains.
Let P D .m; t / D Pr.D.t / D m/ for m D 0; 1; 2; : : : Between action potentials, i.e., on a time interval .t k 1 ; t k /, the process governing D.t / is described by the diagram in Figure 4.1 corresponding to the equation where the factor OEm ¤ 0 is 1 if the statement "m ¤ 0" is true and is 0 if "m ¤ 0" is false.
We look for a solution in which P D .m; t / is given by a Poisson distribution with some unknown mean D .t /: and raising m by 1 in this gives in which we omit the factor OEmC1 ¤ 0, since this is true for every m D 0; 1; 2; : : : . Combining (4.3) and (4.4), Thus, every term of (4.1) contains the factor OEm ¤ 0P D .m 1; t / P D .m; t /; Since D is the expected value of D, we have D Á x D. Thus, (4.5) is the limiting case of (3.3) obtained by taking the limits˛! 0 and n s ! 1 while keeping˛n s Á˛0 constant.
The above shows that if D is Poisson immediately after any action potential, it remains Poisson up to the time of the next action potential. But we also know that for every k the random variables N k and D t C k are obtained from the random variable D t k by binomial splitting, that is, we have This shows that if D.t k / is Poisson, then N k and D.t C k / are Poisson and moreover they are independent random variables. Since D.t C k / is the only possible link between N k and the whole future of the process, it follows that the value of N k has no influence at all upon that future, i.e., that all of the N k are independent. Thus, conditioned on the spike times t k , if the process starts with a Poisson-distributed number of docked vesicles (e.g., 0), then all of the N k are Poisson-distributed and independent. The expected value of N k conditioned on ft k g is obtained by letting !ˇand n s !˛0=ˇin (3.8). The result is the following theorem: In a model synapse with an unlimited number of docking sites obtained by letting n s ! 1 while keeping˛n s Á˛0 constant, let˛0 be the probability per unit time that the number of docked vesicles increases by 1 and letb e the probability per unit time of undocking of a given docked vesicle. Also, let p 0 be the probability of release of a given docked vesicle when an action potential arrives. Finally, let N k be the number of vesicles released by the k th action potential. Then, conditioned on the action potential arrival times ft k g, if the process starts with a Poisson-distributed number of docked vesicles (such as 0), then all of the N k are independent and Poisson-distributed with mean given by the following recurrence: It is surprising that the N k are independent because it may seem that N k should depend on D.t k /, the number of vesicles docked right before the arrival of the k th action potential, which in turn should depend on N k 1 . As we have seen, however, the independence of the N k follows from the Poisson nature of the numbers of docked vesicles, and from the behavior of a Poisson random variable under binomial splitting. We emphasize that the independence of the N k only holds in the limit of an unlimited number of docking sites.
Since the statistics of a Poisson-distributed random variable are determined completely by its mean, Theorem 4.1 provides a computationally efficient way to generate the time series of vesicle release events without the need to simulate the vesicle docking, undocking, and release dynamics, which would require keeping track of the number of docked vesicles. Note, however, that this shortcut is not applicable to a case in which there is only a finite number of docking sites.
In a model synapse with an unlimited number of docking sites and no undocking, the recurrence for N k can be simplified further by lettingˇ! 0 in (4.7): COROLLARY 4.2. In a model synapse with an unlimited number of docking sites obtained by letting n s ! 1 while keeping˛n s Á˛0 constant, suppose that vesicle docking occurs by a homogeneous Poisson process with mean rate˛0, and that there is no undocking .ˇD 0/. Conditioned on the action potential arrival times ft k g, all of the N k are independent and Poisson-distributed with mean given by the recurrence

Optimal Filtering Problem for Stochastic Vesicle Docking, Undocking, and Release
Let P 1 be a stationary stochastic process that generates the presynaptic spike density S.t / 0, which is then used to generate a sequence The optimal filtering problem for stochastic vesicle docking, undocking, and release. The input to the synapse is the sequence of action potential arrival times T k that encodes some continuous signal Q.t /. The stationary stochastic process P 1 generates the presynaptic spike density S.t /, which in turn generates T k and Q.t /. The process P 2 of stochastic vesicle docking, undocking, and release generates the output of the presynaptic terminal, which is the sequence .N k ; T k /, in which N k is the number of vesicles released at the time T k . The rate of vesicle release is P k N k ı.t T k /. This is filtered by convolving it with the function h.t / to produce the reconstructed signal The error in the reconstruction at time t is e.t /. Our definition of error ignores mean values; this is indicated by the capacitor symbol in the path to e.t /. The optimal filtering problem is to choose the impulse response function h. / of the filter to minimize E e 2 .t / .
of ordered action potential arrival times (see Appendix A for an example of such a process). A desired signal Q.t / with mean zero is also generated from S.t /; depending on the function of the synapse, Q.t / can be S.t / itself or some other signal derived from S.t /. Let P 2 be another stochastic process that takes < T k < T kC1 < as input and generates the output of the presynaptic terminal, which is the sequence .N k ; T k / , in which the nonzero integer N k is the number of vesicles released at the time T k . We assume that the random variables N k conditioned on T k have the following properties: First, the N k are independent of each other and also of Q.t /. Next, each of the N k is Poisson-distributed with mean N k . Thus, we are considering here the case of an unlimited number of docking sites, as in Theorem 4.1. Finally, each of the N k depends in a deterministic manner on the T l for l Ä k. We emphasize that these properties hold only when we condition on the T k . Without this conditioning the N k are certainly not independent of each other and also not independent of Q.t /.
In fact, we shall discuss now a procedure for estimating Q.t / from the N k .
We use the notation EOE to denote the expectation over both stochastic processes, but when E is applied to any function of Q or the T k , only the first process P 1 is involved, so in such a case E denotes the expectation over P 1 . We use the notation . / to denote conditional expectation of . / over the process P 2 , with the sequence T k regarded as known.
The optimal filtering problem for stochastic vesicle docking, undocking, and release may now be stated as follows (see Figure 4.2). Let Note that this reconstruction makes use not only of the number N k of vesicles released by the k th action potential but also of the time T k at which that release occurred. We seek h.t / to minimize Recall that by hypothesis. Thus, we are trying to find an impulse response h.t / of the filter such that R.t / approximates Q.t / the best, but our definition of error ignores mean values.
We might, for example, seek to reconstruct the rate S.t /, defined later, at which the stochastic process P 1 generates action potentials. In that case, we would let Q.t / D S.t / EOES.t /. Alternatively, we might be more interested in detecting changes in this rate, in which case we would set Q.t / D dS.t /=dt EOEdS.t /=dt .
In general, Q. / can be any stochastic function of time generated by the process P 1 with the only restrictions being that EOEQ.t / D 0 for all t, and that the joint statistics of Q. / and T k are stationary, i.e., that they are the same for all as the joint statistics of Q . / and It is important to emphasize that, although we speak of reconstructing the signal Q.t /, we are not claiming that a reconstruction of the form (4.8) is actually done by the postsynaptic neuron, or indeed anywhere in the brain. The purpose of considering the reconstruction problem is to make quantifiable the notion of how much information (not in the sense of information theory, but speaking more generally) about the relevant signal Q.t / is transmitted by the synapse. One way to measure this is to attempt a reconstruction (4.8), make it optimal by choosing h to minimize the mean square error, and then use the optimal mean square error as a measure of infidelity of the synapse. Since we are not actually building a device to do the reconstruction, nor are we claiming that such a device exists in the brain, there is no reason to restrict the impulse response h of the filter to be causal [50]; that is, we allow h.t / to be nonzero for all t including t < 0.
Optimal linear filter theory [28,73] states that the optimal acausal filter h. / can be computed by (4.14) and y h; y ' RR ; y ' RQ are the Fourier transforms of their respective functions. y ' RR is known as the auto power spectral density of R, and y ' RQ is known as the cross power spectral density of R and Q.
Therefore, to find the optimal filter h. /, we need to evaluate the expressions on the right-hand sides of equations 4.12 and 4.13. An important tool in evaluating expectations over both processes P 1 and P 2 is the identity that That is, we can evaluate the expectation over both processes by first evaluating the expectation over P 2 with the output of P 1 regarded as known, and then evaluating the expectation of the result over P 1 .
Under the regime of small, band-limited signals, a closed-form expression for the impulse response h. / of the optimal filter can be found in terms of the crossspectral density function of the desired signal and the presynaptic spike times.
We have and therefore (4.16) The last term on the right-hand side of (4.16) is unaffected by h and therefore plays no role in the optimization process. Application of (4.15) to the middle term gives To evaluate the first term in (4.16) , according to (4.15) we add and subtract x R.t / and proceed as follows: In the last term of the foregoing, all quantities have already been averaged over the stochastic process P 2 , so the expectation is only over the stochastic process P 1 . In the middle term, the factor . x R.t / EOE x R.t // has this same property, so that when we average over P 2 the first factor gives x R.t / x R.t / D 0. Thus the middle term vanishes. Hence, the application of (4.15) to (4.18) gives The first term on the right-hand side of (4.19) is the variance of the stochastic process P 2 conditioned on the outcome of P 1 , averaged over P 1 . The second term is the variance of the mean of the outcome of P 2 resulting from the fluctuations produced by P 1 . Substituting (4.17) and (4.19) into (4.16), we get (4.20) From (4.8), we have and therefore Now recall the assumption that the N k , conditioned on T k , are independent and Poisson-distributed. (This assumption holds for docking, undocking, and release models with an unlimited number of docking sites with or without undocking.) This gives Thus (4.20) becomes (4.21) in which x R.t / is given by (4.2). Assume that the sequence of action potential arrival times T k is a perturbation of a sequence of equally spaced times where is a given constant (the unperturbed period of the spike train) and is a small parameter. Then N k can be written as in which x N . / is the mean number of vesicles released by each spike when the spike train is perfectly regular with constant interspike interval .
Note that the leading terms in (4.22) and (4.2) are nonrandom; the random variables are T .1/ k and N k .1/ . We assume that these have mean zero: Equation (4.23) can always be made true by a shift in the origin of time, and it will be shown later that (4.24) follows from (4.23).
Now we make use of the expansions (4.22) and (4.2) to evaluate some of the expressions that appear in (4.21). Substituting (4.22) and (4.2) into (4.2), we get in which x N . / is nonrandom; a formula for it will be derived later.
From (4.25) and (4.2) we get In contrast to (4.26), which has the leading factor 2 in every term, we have since the leading-order term here is actually nonrandom; its expected value is simply its value.
Although we do not yet know how h will depend on , we can see already that the expression given by (4.26) is O. 2 / in relation to the expression given by (4.27). Therefore, we shall neglect (4.26) in the following.
Another term in (4.21) that we need to evaluate is At this point, we need to use the relationship between N k and T k that comes from our model of synaptic vesicle docking, undocking, and release from Theorem 4.1. Making use of the expansions (4.22) and (4.2), (4.7) becomes in which we used the expansion The terms of (4.30) that do not involve give This defines the function x N . /, which up to now has been used without giving it a specific definition. Note that x N . / can also be obtained from (3.27) and (3.6)-(3.7) by letting n s ! 1,˛! 0, while keeping˛n s Á˛0.
The first-order terms of (4.30) give It follows that (4.23) ) (4.24) as claimed above.
Another consequence of (4.36) is that Because P 1 is a stationary stochastic process, EOEQ.t /T and since this holds for all t , we have Putting everything together, we may now write the following formula for the mean square error so that ' QQ .0/ D EOEQ 2 .t /: Recall that P 1 , which generates Q.t /, is a stationary stochastic process.
It is easy to check that the right-hand side of (4.40) is a periodic function of t with period . A priori, it is not clear that we can minimize EOEe 2 .t / separately for each t (although we shall see later this can be done for band-limited signals), so we seek, at first, to minimize the average mean square error: After integration by parts in the term involving h 0 , this becomes and with this choice we find that the optimal average mean square error is given by where y f .!/ is the Fourier transform of f .t /. Our next task is to express this optimal average mean square error as a function of the parameters of stochastic vesicle docking, undocking, and release. To do so, we first evaluate y f .!/ and then j y f .!/j 2 . From (4.43), we have (4.51) With ! and ! 0 both in . ! 0 ; ! 0 / where ! 0 satisfies (4.50), the only value of n in the sum that gives a nonzero result is n D 0, so the whole sum over n becomes ı..! C ! 0 / / ı.! C ! 0 /= . Thus (4.51) becomes Note that this result is independent of t. It follows (see (4.40)), that the optimal choice of h in this band-limited case actually makes the mean square error be independent of t . Since this choice also minimizes the average over t of the mean square error, it must also minimize To see this, let e 0 .t / be the error when h is chosen optimally as described above. Then EOEe 2 0 .t / is independent of t, so its maximum is equal to its mean, which, in turn, is optimal, so we have Thus, (4.55) becomeš where is a given constant (the unperturbed period of the spike train) and is a small parameter. Suppose the stochastic process P 1 that generates both Q.t / and the sequence T k is band-limited in the sense that y ' QT .!/ is supported on some interval . ! 0 ; ! 0 / with ! 0 < ; Then the impulse response h.t / of the filter that minimizes the mean square error (4.9), to lowest order in , has Fourier transform y h.!/ given by in which x N . / is the mean number of vesicles released by each spike when the spike train is perfectly regular with constant interspike interval , and The corresponding minimal mean square error, to lowest order in , is Note that the band-limited assumption (4.63) needs to be made on y ' QT . If we just assume that the auto-power spectral density y ' QQ is band-limited, this tells us nothing about the relationship between Q and the T k , which is of the concern here (see the Appendix for an example).
In Figure 4.3 we plot the analytical estimate of the mean square error (4.64) as a function of p 0 (blue curve) for a small-signal example with D 0:05, in which the desired signal Q.t / is the presynaptic spike density S.t / generated by a smoothed dichotomous jump process. We also plot a numerically evaluated mean square error (red crosses) for comparison. See Section 5 for descriptions of the signal S.t / and numerical procedures.

Optimal p 0 for Synaptic Transmission
To make sense of the above result in the context of how p 0 affects the fidelity of synaptic transmission, we first consider the special case of no undocking (ˇD 0). We have x N . / D˛0 ; D 1 p 0 ; D p 0˛0 : Then (4.64) becomes, to lowest order in , Recall that Â 0 < =2. We now choose Â 0 more specifically as the largest value with the property that for all Â 2 . Â 0 ; Â 0 /. This is clearly satisfied for Â D 0. For all other Â 2 . =2; =2/ it is equivalent to .tan Â Â / 2 < Â 2 :  Thus, we have shown that in a model synapse with no undocking (ˇD 0) and with an unlimited number of docking sites obtained by letting n s ! 1 while keeping˛n s Á˛0 constant, the best choice of p 0 for minimizing the mean square error (4.9) is p 0 D 0, and with this choice the minimal mean square error, to lowest order in , is provided that the assumptions made in Theorem 4.3 hold, and that the desired signal Q.t / is band-limited with tan.! 0 =2/ < ! 0 : (4.71) As p 0 ! 0, the buildup of docked vesicles compensates for the low probability of release of each docked vesicle. The above result shows that in a model synapse with an unlimited number of vesicle docking sites, a smaller probability of release per docked vesicle when an action potential arrives leads to a smaller mean square error in the reconstruction of various desired signals. As p 0 ! 0 the mean square error approaches the value given by (4.70).
For the more general case of an unlimited number of docking sites with a nonzero undocking rate, we consider the limit of smallˇ . When the undocking ratě > 0, the quantities that we need are We are interested in small (but nonzero)ˇ , so we use the following asymptotic versions of these formulae: x Substituting these results into (4.64), we get Note the singular nature of the factor p 0 =.ˇ C p 0 / with regard to the limitš ! 0 and p 0 ! 0. We have On the other hand, Consequently, we can safely neglectˇ within the integral, but not in the factor p 0 =.ˇ C p 0 /. This gives the simpler formula We would like to maximize the term in (4.73) that is subtracted from ' QQ .0/. Whenˇ D 0, we know from the earlier analysis in this section that its maximum occurs at p 0 D 0, but this cannot continue to be true forˇ > 0, since then the term in question is 0 when p 0 D 0. Nevertheless, we expect that the optimal p 0 will approach 0 asˇ ! 0.
Accordingly, we do a Taylor series approximation of the integral with respect to p 0 . We already know its value when p 0 D 0 is and we also have From this we see that the next term in the Taylor series expansion of the integral is where Putting everything together, we now have, to lowest order in andˇ (but with fixed), Our task is to choose p 0 to maximize and we are especially interested in the behavior of the optimal p 0 asˇ ! 0. The equation that determines the optimal p 0 is .ˇ /I 0 3.ˇ /I 2 p 2 0 2I 2 p 3 0 D 0; in which the second term is negligible compared to the first. Therefore, we have shown that, in a model synapse with undocking (ˇ> 0) and with an unlimited number of docking sites obtained by letting n s ! 1 while keeping˛n s Á˛0 constant, the optimal p 0 is given asymptotically by provided that the assumptions made in Theorem 4.3 hold.
A nonzero undocking rate prevents the unlimited accumulation of docked vesicles, so the above result suggests that, in a synapse with a finite number of docking sites, the best choice of p 0 should be some nonzero number. The exact optimal value of p 0 would depend on the parameters of vesicle docking and the statistics of the signal ensemble. Section 5 provides several numerical examples of the optimal filtering of stochastic vesicle release where the optimal p 0 is a nonzero number under various biologically relevant scenarios.

Simulation of Stochastic Vesicle Dynamics and Its Optimal Filtering
This section illustrates the effect of the probabilistic release of synaptic vesicles on the fidelity of synaptic transmission by means of numerical simulations under various biologically relevant parameter regimes. In this section, we no longer require the signal to be of small amplitude, and also we do not assume an unlimited number of docking sites. The numerical methods we use here for the optimal filtering problem are also not restricted to band-limited signals, but we do use the bandlimited assumption in constructing our examples. The rationale for this is that the band-limited assumption relates the rate at which action potentials are generated to the frequency content of the underlying signal in a manner that is reminiscent of the Shannon sampling theorem, so we expect that the spike train will represent the signal well only if such a restriction is imposed.
In all of the following examples, we assume that the desired signal Q.t / is either the presynaptic spike density S.t / or its derivative dS.t /=dt . S.t / is a continuous random function generated by a dichotomous jump process z S .t / (i.e., a telegraph process with two discrete levels), with smoothing, as follows: Let 0 < s 1 < s 2 , and suppose z S .t 0 / D s 1 . We define 12 to be the probability per unit time that S.t / jumps up from s 1 to s 2 , and 21 to be the probability per unit time that S.t / jumps down from s 2 to s 1 . Finally, we remove those frequency components of S.t / that are higher than the mean frequency of jump events (see Algorithm 1 in Appendix B). Note that this smoothed signal S.t / now satisfies the band-limited condition in (4.71), and the condition in (4.63) assumed in Theorem 4.3.
We assume that the presynaptic spike trains T k are generated by the following deterministic integrate-and-fire model with no leakage current over the time interval OEt 0 ; t end (see Algorithm 2 in Appendix B).
Z T k 0 S.t /dt D k; k D 1; 2; : : : : Note that in our recent paper [76] we used the "faithful copy neuron" approach [27,66], which allows the user to generate a spike train with a user-prescribed interspike interval distribution such that the expected spike rate is equal to S.t / at any given time t . Here, however, we use a simple integrate-and-fire model for our numerical examples to remove unnecessary distraction from the stochastic processes of greatest interest: the stochastic docking, undocking, and release of vesicles. A spike train generated by this deterministic integrate-and-fire model is, of course, still stochastic because its spike density S.t / is a random process.
For the special case of an integrate-and-fire neuron, we should remark that the process P 1 is not quite stationary, even if S.t / is generated by a stationary process, since the phase of the spikes is initially fixed by the fact that we start the integration at a particular time. Nevertheless, for most reasonable choices of the process that generates S.t /, the whole process P 1 becomes stationary as t ! 1, and this is true even if S.t / also has some specified initial condition at t D 0 as in our examples. Once we obtain the action potential arrival times T k from the presynaptic spike density function S.t /, we simulate numerically the stochastic docking, undocking, and release of vesicles in our model synapse with a finite number of docking sites at different values of p 0 (see Algorithm 3 in Appendix B).
In the limiting case of an unlimited number of docking sites (possibly with undocking allowed), Theorem 4.1 provides an efficient method for generating independent, Poisson-distributed vesicle release events (see Algorithm 4 in Appendix B).
We then compute the impulse response h.t / of the optimal filter and the mean square error EOEe 2 .t / by averaging over the ensemble of the presynaptic spike density S.t / using a sufficiently large number of samples of S.t / together with the spike train generated by each of them (see Algorithms 5 and 6 in Appendix B). Note that once we know h, the impulse response of the optimal filter for the estimation of S , then the impulse response of the optimal filter for the estimation of the time derivative of S is dh=dt.
The two panels in Figure 5.1 plot the mean square error EOEe 2 .t / as a function of p 0 in the estimation of the presynaptic spike density S.t / and its derivative for various numbers of docking sites. A smaller value of EOEe 2 .t / indicates a more accurate estimation of the desired signal. Figure 5.1 shows that the optimal p 0 for a synapse with 1000 docking sites (without undocking) is p 0 D 0:06 and 0:10 in the estimation of S.t / and dS.t /=dt , respectively. The mean square error curve for the 1000 docking site case (blue curve) shows that an increase in p 0 from its optimal value results in a significantly larger mean square error in the estimation of the two desired signals. In a synapse with 100 docking sites, the optimal p 0 is 0:42 for both estimation tasks. However, the mean square error curve in the case of 100 docking sites (red curve) is less steep compared with the case of 1000 docking sites; any value of p 0 in the range of 0:3-0:5 produces satisfactory estimations for both tasks. As we decrease the number of docking sites further to 10 docking sites or a single docking site, the optimal p 0 for estimating both desired signals becomes 1. Despite the fact that both cases (i.e., 10 sites and a single site) have the same optimal p 0 D 1, their mean square error curves exhibit different shapes and different slopes. A model synapse with 10 docking sites is robust to deviations in p 0 when p 0 is close to its optimal value 1 since the mean square error curve (yellow curve) is relatively flat and concave up. This is in contrast to the case of a single docking site, which has a steeper and concave down mean square error curve (purple curve): a reduction of p 0 from 1 to 0.9 leads to a 55 times increase in the mean square error in the estimation of S.t / and a 159 times increase in the estimation of dS.t /=dt.
The two panels in Figure 5.2 plot the optimal p 0 as a function of the number of docking sites (n s ) in the estimation of the presynaptic spike density S.t / and its derivative. In both estimation tasks, an initial increase in the number of docking sites (from 1 to about 200) leads to a sharp decrease in the optimal p 0 (from 1 to about 0.2). As the number of docking sites increases further, the optimal p 0 continues its approach to 0. However, the rate of approach is faster in the task of estimating S.t / compared with the task of estimating dS.t /=dt.
Recall that p 0 D 0 is optimal in a model synapse with an unlimited number of docking sites, as predicted by our analysis in Section 4.3. To see this, Figure 5.3 plots the mean square error EOEe 2 .t / as a function of p 0 in the estimation of S.t / and its derivative for three cases: an unlimited number of docking sites without undocking, an unlimited number of docking sites with undocking, and 500 docking sites without undocking. In a model synapse with an unlimited number of docking sites (blue curve), the mean square error in the estimation of both S.t / and its derivative increases monotonically as p 0 increases, with p 0 D 0 being the optimal. In contrast, any model synapse with a finite number of docking sites has a nonzero optimal p 0 as shown earlier in Figure 5 undocking rate (red curve) also has a nonzero optimal p 0 , as predicted by our analysis in Section 4.3. This is not surprising because undocking prevents the unlimited accumulation of docked vesicles, and thus it makes a synapse with an unlimited number of docking sites (red curve) behave like a synapse with a finite number of docking sites (yellow curve). Figure 5.4 illustrates the performance of the optimal filter in a model synapse with 100 docking sites for different values of p 0 D 1, 0:3, and 0:01. The left column shows the filtered output R.t / from equation (4.8) in the estimation of S.t /, obtained by applying the optimal filter to the vesicle release events N k . The middle column shows the instantaneous rate of vesicle release (unfiltered output), defined [27] by N k =.T k T k 1 /, where N k is the number of vesicles released at the k th action potential and T k is the time of the occurrence of that action potential. The right column shows the filtered output R.t / from equation (4.8) in the estimation of dS.t /=dt. Overall, p 0 D 0:3 leads to significantly better performance than p 0 D 0:01 and slightly better performance than p 0 D 1 in the estimation of S.t / and dS.t /=dt from the rate of vesicle release.

Determination of Model Parameters
Consider the inverse problem of estimating model parameters. We first note that if we measure and n s , then n s can be any integer such that n s n s : Once n s has been chosen,˛andˇare then determined by It is interesting to note that by considering the mean behavior, it is impossible to distinguish models with the same . ; n s / but different .˛;ˇ/. Such models, however, produce different statistics. We first consider the determination of n s by analyzing the observed time series N k in a regular spike train with equal interspike intervals (see Section 3.3). We have It may appear that the above formula for n s is impractical to use when n s x N , since in that case ' kk = x N may be close to 1 and the difference 1 ' kk = x N may not be measured accurately. In practice, however, this formula works sufficiently well in the parameter regime of our interest; it gives accurate, and many times exact, results for a wide range of parameters, e.g., p 0 from 0:01 to 1, n s from 1 to 1000, and˛0 D 1000 with or without undocking, especially in the regime of n s = x N < 100. An alternative method for estimating n s is based on fitting the negative tail of the autocovariance '. Recall that ' kk in (3.30) gives the height of the central peak of ', and that . x N / 2 =n s , the coefficient that multiplies the geometric sequence in ' i k , measures the magnitude of the negative tail in '. The ratio r of the height of the central peak divided by the magnitude of the negative tail is given by equation (3.31). Thus, if we measure x N and also the ratio r, we can find n s very simply as n s D x N .r C 1/: (6.4) Unfortunately, this method based on "tail-fitting" turns out to be more sensitive to noise compared to the first method based on formula (6.3). In practice, we can use (6.3) for the estimation of n s , but then compare the experimental autocovariance to the theoretical one as a check that this theory is applicable (see Figure 3.2).
We have now discussed how to determine the parameter n s from the autocovariance, or if n s is so large that its value is of no importance we can determine that fact from the absence of a negative tail in the autocovariance. Note that the determination of n s can be done with spike trains of different periods t to check that the result is independent of t. We shall next discuss how to determine the other parameters of the model.
The parameters n s and are best determined from depletion experiments [13,15]: Give a rapid train of action potentials of sufficient number to release all docked vesicles, wait a time of duration T , and repeat. See Figure 6.1 for an illustration.
During each of the bursts, count the number of vesicles released. Since the mean number docked at the end of each burst is (approximately) 0, the mean number and this is the average of the number that will be released in each burst. By measuring this for at least two values of T , we can identify n s and . If more than two values of T are used, then we also get a check on the theory.
For example, if we have n 1 D n s .1 e T 1 /; (6.5) n 2 D n s .1 e T 2 /; (6.6) where .n 1 ; T 1 / and .n 2 ; T 2 / are the measured data, then must satisfy Note that Since Â=.e Â 1/ is a decreasing function, we have the following for > 0: Without loss of generality, suppose T 1 < T 2 . Then there exists a unique solution to (6.7) provided that If this is not the case, the data are not consistent with our model. Once has been found, e.g., with the MATLAB ® function fzero, we can solve (6.5) or (6.6) for n s (and of course we should get the same answer either way).
We would like to offer two remarks here. First, the above method for determining and n s is insensitive to the value of p 0 . This is good because p 0 may be different under different conditions as a result of facilitation [9,65,67,81]. Even if p 0 varies during the depletion burst, this is perfectly all right, since the burst is designed to deplete and count all of the vesicles that were docked prior to the burst, and as long as enough action potentials are applied for this to occur, it does not matter how many vesicles are released by each action potential during the burst. In fact, facilitation improves the accuracy of the result, since it makes it more likely that the depletion burst will be 100% effective. Second, the proposed method uses only the mean number of vesicles released by a depletion burst, not the statistics of that number. We could try to extract more information from the statistics, but we do not need to do so. Note, however, that we did need to use the statistics of the process to determine the actual number of sites n s as opposed to the effective number of sites n s . Recall that n s was determined above from the correlation function ' i k .
Once and n s have been determined by the depletion experiment, we can use the average number of vesicles released per spike in a regular (i.e., periodic) spike train with period t to determine p 0 . The relevant equation is (3.27), which defines x N .t /. Solving this equation for p 0 we obtain Note that p 0 may depend on t because of facilitation. The above procedure and formula still work for each t and can also be used to determine p 0 .t /, i.e., p 0 as a function of t .
We have now discussed the determination of all of the parameters of the model. In particular, we can determine n s ; n s ; ; p 0 from mean data only. Finally,˛andˇare given by (6.1) and (6.2). Note that the depletion experiment provides yet another test of our theory: we can compare the experimentally observed distribution of the number of vesicles released in depletion experiments for any interval T to our theoretical prediction discussed below. The depletion experiment measures how many vesicles are docked Since the sites are independent, the number of vesicles docked at time T , D.T /, will have a binomial distribution with parameters n s and P .T /. Note that we have a prediction for every T . Note also that this new test has the advantage, as we have discussed already, of not requiring p 0 to be constant, and not depending at all on the value of p 0 . Ifˇ= 0, P .T / ! 1 as T ! 1, so the result is simply that n s vesicles are released every time T is sufficiently large. In several experimental studies, however, this is not what has been observed [15,42].
If n s is large, and recall that˛D˛0=n s , then P .T / behaves like so in this case we expect a Poisson distribution with mean˛0 .1 e ˇT /.
Together with comparing the experimental autocovariance of N k to the theoretical one as discussed earlier in this section, we have now provided two separate checks on our theory of vesicle docking, undocking, and release.
Note also that our method of parameter estimation is related to "quantal analysis," also known as "mean-variance analysis" or "multiple probability fluctuation analysis," which is a general method that uses mean, variance, and covariance from vesicle release statistics to recover synaptic parameters (see [58,62] for reviews and analyses of quantal analysis). The method presented here is designed for our model of synaptic vesicle docking, undocking, and release, whereas the quantal analysis method does not explicitly model vesicle docking and undocking processes but instead requires the user to prescribe the probability of vesicle occupancy (and its adaptation over time) at each docking site. In particular, our method of synaptic parameter estimation allows for vesicle undocking, which is a constant source of noise that leads to fluctuation in the number of docked vesicles even in the absence of stimuli, a phenomenon observed at single synaptic release sites in hippocampal slices [15,43].

Discussion
The first topic that we would like to discuss is the role of undocking in our model. As a starting point for this discussion, consider the case of an unlimited number of docking sites with no undocking. This model, which was the focus of our previous paper [76], has some beautiful mathematical properties. In particular, the mean rate of release of synaptic vesicles is completely insensitive to the mean rate of arrival of action potentials, so that the model synapse responds only in a transient way to changes in the rate of arrival of action potentials. This is an idealized version of the typical behavior of much of the nervous system-in which there is a strong response to sudden changes, but mean values are largely ignored. Also, the numbers of vesicles released by the different action potentials, conditioned on the action potential arrival times, are independent and Poisson-distributed. Finally, the optimal choice of p 0 , as measured by the ability of a linear filter to reconstruct the signal that is being carried by the incident spike train, is 0. This last result is paradoxical, and points to the singular nature of the model. The reason that a low value of p 0 is helpful is that it allows for a buildup in the number of docked vesicles. With an unlimited number of docking sites, and with no undocking, there is no upper limit to this buildup.
The above pathology can be cured in either (or both) of two ways. One way is to assume that there is a finite number of docking sites, and the other way is to allow for undocking. Either of these features is sufficient to prevent the unlimited buildup of docked vesicles, and the two choices have very similar consequences insofar as the mean behavior of the model is concerned, but they have different statistical consequences, and this should allow them to be distinguished experimentally, although we emphasize that it is not a question of choosing between mutually exclusive alternatives, since both may be operating at the same time.
If we retain the assumption of an unlimited number of docking sites and allow for undocking, then some of the beautiful mathematical properties mentioned above are retained. In particular, it is still the case that the numbers of vesicles released, conditioned on the action potential arrival times, are independent and Poisson distributed. This facilitates mathematical analysis and is fundamental to theory developed in this paper on the optimal choice of p 0 . With undocking, however, it is no longer the case that the optimal choice of p 0 is equal to 0. Instead, it is O..ˇ / 1=3 /, whereˇis the undocking rate and is the unperturbed time between action potentials in the incident spike train.
When the number of docking sites is finite, the numbers of vesicles released are no longer conditionally independent, but the consequences for the optimal filtering problem are similar to those of undocking, since either process prevents the unlimited buildup of docked vesicles and thereby shifts the optimal p 0 towards larger values. Our demonstration of this is numerical, however, since we do not have a theoretical way to evaluate the optimal p 0 when the number of docking sites is finite.
More generally, we have shown that there is a whole class of models that are equivalent with regard to their mean behavior but differ with regard to their numbers of docking sites and their rates of undocking. These models can only be distinguished by considering their fluctuations. It would not be surprising, given this equivalence in the mean, that some phenomena that have been attributed to a limitation in the number of docking sites are actually caused by undocking or vice versa.
The second topic for discussion is the relationship between synaptic failure and the probability of vesicle release. A quantity that is often measured in experiments is p r D Pr.at least one vesicle is released by an action potential/; so that 1 p r is the probability of synaptic failure. The parameter we used throughout this article is p 0 , the probability of release per docked vesicle upon the arrival of an action potential. These two probabilities are related by where D.t / is the number of docked vesicles at time t . Hence, p r is not a constant but a dynamic variable. Our proposed method of parameter identification estimates p 0 directly. Moreover, as discussed in Section 6, even if p 0 depends on the period t of the regular spike train that is used in the proposed identification process, because of facilitation, our estimation formulas still work for each t and can be used to determine p 0 as a function of t .
Finally, we would like to point out some limitations of our model. First, our model does not take into account facilitation. Facilitation could be modeled by allowing the parameter p 0 to depend on the rate of the spike train. We do not consider this effect in this article except that we have designed a parameter identification method that is independent of p 0 , so that if other parameters remain constant all parameters will be successfully recovered despite the presence of facilitation. Second, our model only considers activity-dependent synchronous release, and there is no partial release, or variation in the content of a vesicle. Third, our model also assumes that p 0 is the same for all sites, and random variables associated with different sites are independent, conditioned on the action potential arrival times. Note also that, except for docked vesicles, our model does not consider vesicle pool dynamics. Last, the optimal filter considered in this article is acausal in the sense that a causality constraint is not imposed on h.t /, which can take nonzero values for all t. We did so because we aimed to obtain a theoretical lower bound on the mean square error in the linear estimation of desired signals using the presynaptic vesicle release events. From the perspective of this paper, the observer is allowed to record the vesicle release events and process them at leisure to reconstruct the desired signal. It may also be of interest to consider the causal case in which the requirement is imposed that h.t / D 0 for t < 0.
Let be such that s 0 D 1; and look for a solution T k of the above recursion relation with the property that Expanding in powers of we find Then, since s 0 D 1, we get

Appendix B Algorithms
Algorithm 1 Generate a smoothed dichotomous signal S.t / to be used as the presynaptic spike density Input: Scalars s 1 , s 2 , 12 , 21 , t 0 , t end , n t , and t , in which n t is the number of time steps for the discretization of the time interval OEt 0 ; t end , t D .t end t 0 /=n t , and the rest of the parameters are defined in the main text. Output: Array S D S 0 ; S 1 ; : : : ; S n t , the presynaptic spike density function S.t / in the form of a vector. Each entry in the array S corresponds to a time node in the vector t 0 ; t 0 C t; : : : ; t end , that is, S j D S.t 0 C jt /.
1: Initialize an array variable S D S 0 ; S 1 ; : : : ; S n t with length n t C 1 2: S 0 WD s 1 3: Initialize a scalar variable j WD 1 4: Draw a new random number u uniformly distributed on OE0; 1 5: while j < n t do j WD j C 1 18: Compute the discrete Fourier transform of S and store it in the array y S D y S .!/ , in which y S .!/ denotes the Fourier coefficient corresponding to the frequency ! 19: y S .!/ WD 0 for all frequencies ! that satisfy j!j > . 12 C 21 /=2 20: Compute the inverse discrete Fourier transform of y S and store it in the array S 21: Return S Algorithm 2 Generate action potential arrival times (i.e., spike times) for an integrate-and-fire neuron with presynaptic spike density S.t / Input: Array S generated by Algorithm 1 and scalars t 0 , t end , and t . Output: Array T D T 1 ; T 2 ; : : : ; T k , whose entries are spike times. 1: Initialize a scalar variable k WD 0 2: Initialize a scalar constant T 0 WD t 0 3: Initialize an adaptive-length vector variable T D T 1 ; T 2 ; : : : 4: while T k < t end do 5: k WD k C 1 6: Solve R T k T k 1 S.t /t D 1 for T k , in which S.t / is the linear interpolation of S over evenly spaced time nodes OEt 0 ; t 0 C t; : : : ; t end 7: Return T Algorithm 3 Simulation of stochastic vesicle docking, undocking, and release in a model neuron with a finite number of docking sites Input: Array T generated by Algorithm 2 and scalars t 0 ; n s ,˛, andˇas defined in the main text. Output: Array N D N 1 ; N 2 ; : : : ; N k , in which each entry is the number of vesicles released at the spike time given by the corresponding entry in T. if "site occupied" D 1 then if v < p 0 then 12: N j WD N j C 1 x WD x C .T j T j 1 /.1 p 0 / l j x WD x C .1 p 0 / j 1 e ˇ.T l T l j /

14:
N l WD a Poisson-distributed random number with mean˛0 p 0 .1 p 0 x/ 15: Return N Algorithm 5 Compute the impulse response h.t / of the optimal filter of the rate of stochastic vesicle release in the estimation of S.t / Input: Scalars s 1 , s 2 , 12 , 21 , t 0 , t end , n s ,˛(or˛0 in the case of an unlimited number of docking sites),ˇ, p 0 , n t , t , and n sample , in which n t is the number of time steps for the discretization of the time interval OEt 0 ; t end , t D .t end t 0 /=n t , n sample is the number of simulated sample paths of the stochastic process S.t /, and the rest of the parameters are defined in the main text. Output: Array h D h 0 ; h 1 ; : : : ; h n t =2 , the optimal impulse response h.t / in the form of a vector (assuming n t is even). Each entry in the array h corresponds to a time node in the vector t 0 ; t 0 C t; : : : ; t end =2 , that is, h j D h.t 0 C jt /. Generate the array S (presynaptic spike density ) using Algorithm 1 5: Generate the array T (action potential arrival times) using Algorithm 2 6: Generate the array N (vesicle release numbers) using Algorithm 3 (or Algorithm 4 in the case of an unlimited number of docking sites) 7: Construct the function N.t / D P k j D1 N k g.t T j /, the linear interpolation of N over time nodes T, in which g.t / is a single square pulse centered at t D 0 with width t and height 1=t Evaluate N.t / on evenly spaced time nodes on the second half of the time interval, and store the obtained values in the array N WD N n t =2 ; N n t =2C1 ; : : : ; N n t , in which N j D N.t 0 C jt / 9: S WD S n t =2 ; S n t =2C1 ; : : : ; S n t , in which S j D S.t 0 Cjt / (i.e., remove the first half of the entries in S) 10: S WD S 1 n t =2C1 P n t j Dn t =2 S j (i.e., subtract the mean from itself) 11: N WD N 1 n t =2C1 P n t j Dn t =2 N j (i.e., subtract the mean from itself) 12: Compute the array S NS WD cross power spectral density between N and S 13: Compute the array S NS WD auto power spectral density of N 14: Compute the array y h WD S NS =S N N (component-wise division) 15: Compute the array h sample WD inverse discrete Fourier transform of y h 16: h WD h C h sample 17: h WD h=n sample 18: Return h Algorithm 6 Compute the mean square error in the estimation of S.t / using the optimal filter applied to the rate of stochastic vesicle release Input: Scalars s 1 , s 2 , 12 , 21 , t 0 , t end , n s ,˛(or˛0 in the case of an unlimited number of docking sites),ˇ, p 0 , n t , and t , and array h D h 0 ; h 1 ; : : : ; h n t =2 (see Algorithm 5 for parameter descriptions).
Output: Mean square error e 2 .
1: Initialize a scalar variable e 2 WD 0 2: for l D 1 to n sample do 3: Generate the array S (presynaptic spike density ) using Algorithm 1 4: Generate the array T (action potential arrival times) using Algorithm 2 5: Generate the array N (vesicle release numbers) using Algorithm 3 (or Algorithm 4 in the case of an unlimited number of docking sites) 6: Construct the function R.t / D P k j D1 N j h.t T j / 7: Evaluate R.t / on evenly spaced time nodes on the second half of the time interval, and store the obtained values in an array R WD R n t =2 ; R n t =2C1 ; : : : ; R n t , in which R j D R.t 0 C jt / 8: S WD S n t =2 ; S n t =2C1 ; : : : ; S n t , in which S j D S.t 0 Cjt / (i.e., remove the first half of the entries in S)

9:
S WD S 1 n t =2C1 P n t j Dn t =2 S j (i.e., subtract the mean from itself) 10: R WD R 1 n t =2C1 P n t j Dn t =2 R j (i.e., subtract the mean from itself) 11: Compute the scalar variable e 2;sample WD 1 n t =2C1 P n t j Dn t =2 .R j S j / 2 12: e 2 WD e 2 C e 2;sample 13: e 2 WD e 2 =n sample 14: Return e 2