Exact distribution of the quantal content

17 The transfer of electro-chemical signals from the pre-synaptic to the post-synaptic 18 terminal of a neuronal or neuro-muscular synapse is the basic building block of 19 neuronal communication. When triggered by an action potential, the pre-synaptic 20 terminal releases neurotransmitters into the synaptic cleft through exocytotic vesi21 cle fusion. The number of synaptic vesicles that fuse, i.e., the quantal content, is 22 stochastic, and widely assumed to be binomially distributed. However, the quantal 23 content depends on the number of release-ready vesicles, a random variable that 24 depends upon a stochastic replenishment process, as well as stochastic interspike 25 intervals of action potentials. The quantal content distribution suitably accounting 26 for these two stochastic processes was not known explicitly in the literature. Here 27 we analytically obtain the exact probability distribution of the number of vesicles 28 released into the synaptic cleft, in the steady state reached during stimulation by a 29 train of action potentials. We show that this distribution is binomial, with modified 30 parameters, only when stimulated by constant frequency signals. Other forms of 31 stochastic stimulation, e.g. Poissonian action potential train, lead to distributions 32 that are non-binomial. The general formula valid for arbitrary distributions of the 33 input inter-spike interval may be employed to study neuronal transmission under 34 diverse experimental conditions. We provide exact theoretical moments which may 35 be compared with experiments to estimate the model parameters. We corroborate 36 our theoretical predictions through comparison with the quantal content distribu37 tions obtained from electrophysiological recordings from MNTB-LSO synapses of 38 juvenile mice. We also confirm our theoretically predicted frequency dependence of 39 mean quantal content by comparing with experimental data from hippocampal and 40 auditory neurons. 41

pre-synaptic neuron reach the synapse and trigger neurotransmitter-loaded vesicles The human brain is arguably the most complex natural object in the known universe.

61
Its complexity, as well as its capabilities, are a product of its roughly 86 billion [1] 62 neurons that make approximately 10 15 connections [2] with each other. The majority 63 of these connections are chemical synapses, where a pre-synaptic neuron communi-64 cates with a post-synaptic neuron through neurotransmitter release in a few nm rive exact formulae that describe the probability distribution as well as its moments 140 under these conditions. 141 For a synaptic system (Fig 1a) with a finite maximal number M of docked vesi-142 cles, sustained stimulation by APs for long enough time typically leads to a steady 143 state [7,12,20,37,38]. In such a steady state, the time-varying calcium concentration-144 dependent replenishment rate k t and release probability p r t would saturate to con-145 stant values k = k ∞ and p r = p r ∞ (see representative curves in Section 8 of supple-146 mentary information (SI) following a model of Ref. [15]). The quantal content distri-147 bution, which is the main quantity of our interest, also reaches a time-independent 148 form in this limit. In the steady state, we prove that for commonly used constant 149 ISIs, the distribution of quantal content b is indeed a binomial B(M, p rb , b) (where 150 p rb is a function of p r , k, and the frequency of stimulation). However, for a Poisson 151 input stimulation with exponentially distributed ISIs, the quantal content distribu-152 tion is non-binomial and we obtain its exact closed form. Apart from these special 153 cases, for any general stochastic stimulation with uncorrelated ISIs (t s ) drawn from 154 a distribution g(t s ) (see Fig 1b), we obtain a formal solution for the quantal content 155 distribution, in the steady state. 156 In practice, there may be systems where steady state is hard to reach due to long signed and implemented experimentally to test our predictions. Input stimulations 163 in vivo are also expected to be stochastic [40,41]. Thus, there are contexts where 164 the treatment of stochastic stimulations analyzed in this paper may be naturally 165 useful. In such situations model parameters may be estimated by direct matching 166 of exact theoretical moments of the quantal content (obtained in this paper) with 167 those obtained experimentally. Exact closed form results significantly reduce if not 168 eliminate computational cost, unlike earlier works. Furthermore, closed form results 169 are useful for subsequent theoretical development. 170 We do not explicitly consider the process of vesicle priming [42] separately and 171 assume that all processes including transport and priming that prepares vesicles for 172 fusion can be incorporated through the rate k (Fig 1a). We also do not consider 173 processes in the post-synaptic terminal after neurotransmitter release, such as post-174 synaptic receptors saturation or desensitization. Finally, we consider the distribution 175 of the quantal size to be narrow, such that the postsynaptic current may be assumed 176 proportional to the quantal content. Since the analytical tractability was the main 177 emphasis, we work with the minimal model which most earlier works in the literature 178 have used, keeping aside the above mentioned complexities for future inquiries. where n is the number of already docked vesicles and k is the docking rate per empty 203 site [20,22]. We assume that docked vesicles may be released probabilistically after 204 an AP but do not otherwise detach. We consider experimentally relevant scenar-205 ios (like in Ref. [38]) where a train of APs arrive at the pre-synaptic neuron. In 206 this study we consider the general possibility that the time intervals t s between two 207 consecutive spikes of the input AP may be random but uncorrelated, following a dis-208 tribution g(t s ). Experimentally often fixed frequency oscillatory input AP are used 209 with constant t s = T [20,21,30,38], which would correspond to g(t s ) = δ(t s −T ). The 210 number b of SVs that fuse and release neurotransmitters after the arrival of an AP is 211 a stochastic variable, referred to as the quantal content. Following [12,20,21,28,30], 212 we assume that the quantal content b after an AP follows a binomial distribution where p r is the probability of release (in the steady 214 state) and n − is the instantaneously docked vesicle number before the AP arrives. 215 We note that the quantal content b is stochastic not just because it is binomially 216 distributed, but because the number of docked vesicles n − which form its source is 217 itself stochastic. As we discussed above, the ready to release docked vesicle number 218 n grows stochastically with time, and at random intervals t s when an AP appears, 219 n = n − suddenly become n = n − −b. Thus in the steady state, the exact distribution 220 of b may be obtained by averaging the binomial over the stochastic n − , which leads 221 to a revised distribution as we would see below.

222
The quantal content distribution Q ss (b) in the steady state Schematic showing the variables we used to describe the time evolution of the number n of the vesicles inside the pre-synaptic neuron, from the value n ′ + after an action potential to n − just prior to the next action potential after time t s and then its sudden reduction to n + after a burst of size b. vesicles, resulting in a sudden decrease of n from n + + b to n + . Before the arrival 233 of APs in the pre-synaptic terminal, the instantaneous number of docked vesicles 234 n − may be assumed to be at its maximum value, M , and the quantal content b is the vesicle release and replenishment process. In this long-time steady-state limit, our goal is to obtain an explicit expression for the size distribution Q ss (b) of b, 239 summing over the likelihood of stochastic vesicle number n − , as follows: To calculate Q ss (b), we need to first obtain P ss − (n − ), the steady-state probability 241 distribution of having n = n − docked vesicles just before an AP. As mentioned 242 previously, n − and hence P ss − (n − ) is regulated by the random processes of growth of 243 n and the appearance of AP at a random t s . To study the growth of n, we introduce 244 the time-dependent conditional probability p(n, t|n ′ + ), of having n number of docked 245 vesicles at time t, given n ′ + docked vesicles at t = 0, i.e., just after the previous AP 246 (see Fig 2). In terms of this conditional probability, the distribution P ss − (n − ) in Eq.

247
(1) is given by averaging p(n, t s |n ′ + ) over both the inter-spike interval t s as well as the 248 initial number n ′ + . Note that we are using primed variables to distinguish variables 249 at the previous AP from those at the current AP.
where P ss + (n ′ + ) is the steady-state probability distribution of having n = n ′ + just after 251 an action potential. Since the distribution of ISIs is given experimentally, we need 252 to calculate the other two distributions inside the integral on the right.

253
The conditional distribution p(n, t|n ′ + ) of n between two successive action po- This equation can be solved using the generating function F (q, t|n ′ + ) = n q n p(n, t|n ′ + )

257
(see SI Section 1 for the details), and we obtain the solution as, We next need to calculate P ss + (n + ), the steady-state probability distribution of 259 having n + docked vesicles immediately after the vesicle release. This is non-trivial 260 since n + = n − − b, and since n − depends upon the n ′ + of the previous interval, the 261 result for any interval depends upon the results of the previous interval. Starting with 262 exact recursive relations between conditional probabilities P + after two successive 263 APs (see Eq. (7) in section 2 of SI for details) in similar lines as [33,34], under the 264 assumed conditions in steady state, we obtain a self-consistent equation for P ss + (n + ).

265
P ss Eq. (5) contains the unknown P ss + (n + ) on both sides. Other quantities are known, 266 including p(n + + b, t s |n ′ + ) (from Eq. (4)). In a difficult and lengthy calculation (see 267 details in SI Section 3), we show that the corresponding generation function to the 268 above distribution, namely F + (q) = ∞ n + =0 q n + P ss + (n + ) can be obtained using the 269 knowledge of F (q, t|n ′ + ), the generating function for p(n, t|n ′ + ) mentioned above. This 270 generating function leads us to the other two relevant generating functions corresponding to the distributions 272 P ss − (n − ) and Q ss (b). Since the three stochastic variables, n − , n + , and b are related to 273 each other (namely n + = n − − b), the generating functions must be too, and we can 274 show that (see details in SI Section 4) they are related as: From the knowledge of F b (q), we may finally obtain the desired distribution of 276 quantal content Q ss (b) (see SI Section 4) using the standard relationship Q ss (b) = where, Here {S (n−1) } denotes the set of all the subsets S (n−1) = {m i } = (m z , m z−1 , · · · , m 1 ) 280 with m i ∈ (1, 2, · · · , n − 1) such that m z > m z−1 > · · · > m 1 . Although the upper 281 limit of the sum is shown up to ∞, due to the combinatorial factor M n , n never Note that for an arbitrary g(t s ) there is no reason to expect that Q ss (b) in Eq.

286
(7) would reduce to a binomial distribution. As the number of subsets S (n−1) grows When t s is deterministic, we have an AP train at a fixed frequency f = 1/T , and 296 g(t s ) = δ(t s − T ). When we substitute this in Eq. (7), we find (see details in SI

297
Section 5(a)) that the expression for Q ss (b) reduces to a binomial distribution: Thus the steady state quantal content distribution corresponding to a determin-

318
The mean and CV 2 of the exact binomial distribution Eq. (9) are, We note that although the mean of the approximate binomial is ⟨n − ⟩p r = M p rb (same as the exact binomial), its CV 2 = (1−pr) M p rb is clearly different from the exact Exponentially distributed ISIs

330
When the AP train is a Poisson process, the interspike intervals are exponentially 331 distributed, g(t s ) = f 0 e −f 0 ts . We can substitute this expression in Eqs. (7) and (8) 332 to obtain (details in SI Section 5(b)) the following expression for the steady state 333 quantal content distribution: We can again calculate closed form expressions for the mean and CV 2 of the 335 distribution in Eq. (12) (see Section 5(b) in SI): Note that this distribution is no longer a binomial. Thus even in one of the 337 simplest cases of a stochastic g(t s ), a systematic treatment of the stochasticity in n − 338 leads to a quantal content distribution that is non-binomial. It is instructive to ask bution under these circumstances, it would be interesting to test it experimentally 393 by using Poissonian AP stimulations, and estimate the parameters (p r , k) using that 394 data. For that purpose, the explicit formulae for the theoretical moments from Eq.

395
(13) may be used. terminal. The steady-state distribution of the quantal content therefore depends 453 upon these coupled stochastic processes.
In this paper, we successfully perform this calculation, which yields a formula, Eq.

455
(7), for the probability distribution of the quantal contents under the assumption of 456 a reasonable kinetic replenishment model and arbitrary g(t s ). We also derive exact 457 expressions for the mean and CV 2 of this general distribution (Eq. (8)). Using  3) and (4) Rev Physiol, vol. 17, pp. 49-117, 1978.

561
[17] S. Redman, "Quantal analysis of synaptic potentials in neurons of the central 562 nervous system.," Physiological Reviews, vol. 70, no. 1, pp. 165-198, 1990.   between two successive APs is: We define a generating function for the above probability distribution: Here the limits of the above sum are n ′ + ≤ n ≤ M . Otherwise, we assume p(n, t|n ′ + ) = 29 0 if n ′ + falls outside this range. Multiplying Eq. (1) by n q n on both sides, we get 30 a PDE for the generating function F : Rearranging the terms, we have, We use the method of Lagrange characteristics to solve the Eq. (3). The result is: We do binomial expansion of Eq. (4) to get Next we can compare Eq. (5) with Eq.
(2) to get the expression of p(n, t|n ′ + ): for n ′ + ≤ n < ∞. Note that for n > M , the probability goes to zero automatically 37 in the above expression.
We may follow the same procedure for P m−1 + (n +,m−1 , t sm−1 ) and so on, and build 54 up a hierarchy of relations for any interval between two APs. The problem remains 55 hard to solve as result for every interval depends on the history of the previous 56 interval. We are particularly interested in the steady state attained after several 57 APs appear, i.e. m ≫ 1 and the output quantal contents statistically stabilise. We 58 integrate Eq. (7) over the joint probability distribution g 2 (t sm , t sm−1 ) of t sm and 59 t sm−1 . We get If the interspike intervals (ISIs) are not correlated, we would have g 2 (t sm , t sm−1 ) = 61 g(t sm )g(t sm−1 ), where g(t sm ) and g(t sm−1 ) are the normalised distributions of t sm and 62 t sm−1 respectively. Then, In the long time limit, the steady state probability distribution of vesicle number 64 after any AP is P ss Note that this is a self-consistent equation for the unknown P ss + (n + ), when we 68 know g(t s ), B(n + + b, p r , b) and p(n + + b, t s |n ′ + ) (Eq. (6)). This equation appears in 69 the main manuscript. We analytically treat it in the next section.

70
3 Generating function of the distribution P ss + (n + ) 71 in the steady state 72 We define a generating function of P ss + (n + ) as follows and solve it in this section: In Eq. (10), the limits of the summation are 0 ≤ b ≤ M − n + and 0 ≤ n ′ + ≤ n + + b.
Doing summation overñ + and using the definition of the generating function F (from 81 Eq. (2)), we get Defining h(q) = (q − qp r + p r ) and using Eq. (4) for the explicit form of function F : In the above we have used the definition of F + from Eq(11). In Eq. (14) we have a self-consistent integral equation for the function F + with two different arguments, and involving an integration over the inter-spike interval t s . Clearly, it is difficult to 86 solve the equation directly for F + (q). But as we show below, F + (q) may be obtained 87 as a series expansion around q = 1, that is where the series coefficients F (n) . The procedure of obtaining these 89 coefficients are as follows.

90
Taking derivative with respect to q on both side of Eq. (14), we have , we have Defining the quantities ψ n,m and L n as follows, where we have set ∂F + (q) + (1) according to the general definition of the n-th 96 series coefficient above.

97
Next we may take derivative of Eq. (16) with respect to q and set q = 1 and obtain an equation for F (2) + (1). We find that it is linearly coupled to F (1) + (1) as follows: We may follow the same procedure to obtain equations for the higher derivatives.

100
It is clumsy to do it manually, so we use Mathematica for this. We use symbolic 101 derivatives of Eq. (16) and automatic replacement of terms by the ψ functions using 102 the Replace syntax. The result are the equations for the next two series coefficients.

105
The explicit expressions of the four coefficients are In the above expressions, Eqs. (23) -(26), we notice the following patterns. The where the sum over {S (n−1) } runs over all such subsets as S (n−1) mentioned above.  (28) 4 Relations between the generating functions, and 124 the quantal content distribution Q ss (b) 125 Let us denote the probability distribution of having n − number of docked vesicles 126 before an AP in the steady state by P ss − (n − ). It is easy to see from similar arguments 127 as in Section 2, that P ss − (n − ) will be related to the three distributions: P ss + (n + ) (to 128 have an initial number), p (n − |t s , n + ) (for docked vesicles to build up), and g(t s ) (for 129 AP to appear). Thus We define its generating function as The ultimate object of our interest, namely the distribution of quantal content b in 132 the steady state, is related to the distribution P ss − (n − ) as: Note that the the crucial fact is that we are averaging the well-known binomial 134 distribution now, over the stochastically generated docked vesicle number n − just 135 prior to an action potential. As we would see, the resultant averaged distribution 136 would be different from the original binomial. The generating function corresponding We have four distributions p (n − , t s |n + ), P ss + (n + ), P ss − (n − ), and Q ss (b) and four cor-139 responding generating functions F , F + , F − and F b , respectively. Of these, we have 140 already calculated F (Eq. (4)) and F + (Eq. (28)). We would now proceed to relate 141 the generating functions F − and F b to F + . Once F b is known explicitly using the 142 expression of F + , we would then obtain Q ss (b) by taking suitable derivatives of F b (q) 143 with respect to q.

144
Using Eq. (29) in Eq. (30) we get, In the above we have used the definition of the function F . Substituting the explicit 146 form of F (q, t s |n + ) from Eq. (4), we have In the above we have used the definition of F + .

148
Substituting q−pr 1−pr in place of q in Eq. (14), and noticing h q−pr Comparing Eqs. (33) and (34), we immediately get, Here we have used the definition of the function F − , and the relation of F + to F − 153 from Eq. (35). Using the explicit form of F + from Eq. (28), the Eq. (36) gives, We may invert Eq. (32), to obtain the quantal content distribution Q ss (b) from 155 successive derivatives of F b (q) as follows: Note that q is present in F b (q) only through the factor (q − 1) n in Eq. (37). Using Moments:

159
(i) For the general distribution of ISIs, following Eq. (39), the mean quantal content 160 ⟨b⟩ is given by Note that  (ii) Similarly, the second moment ⟨b 2 ⟩ is given as: Where, Substituting these expressions in the Eq. (39), and recalling ϕ n,m = n m ψ n,m , we It is important to remind that {S (n−1) } is the set of all the subsets of the set 179 {1, 2, · · · , n − 1}. To proceed further, we need to use the following identity (for |x| < 1) -for proof see Section 6: Using the identity (Eq. (49)) in Eq. (48), the expression simplifies to give, Further using the identity, Hence, the steady state quantal content distribution for the deterministic g(t s ) 184 is still a binomial distribution, but with the parameters M (maximal docked vesicle 185 number) and an effective release probability p rb (different from p r ) given by: In order to compare the exact distribution (Eq. (50) .

190
Taking derivative of Eq. (35) and setting q = 1, we get If we take the derivative of Eq. (28) with respect to q at q = 1, only n = 1 term in the summation survives. Eq. (52) will become Substituting the forms of L 1 and ϕ 1,0 from Eq. (47), we get In the above we have used the expression of p rb from Eq. (51). The Eq. (53) where we have used Beta function B(a, b) = Γ(a)Γ(b)/Γ(a + b), and then Moreover using Eq. (54), we have The above Eq. (57) may be extended to get the following product: Substituting Eq. (58) into Eq. (39), we have We may simplify the expression above by using the following identity valid for any 211 arbitrary function f -for proof see Section 7: Using Eq. (60) in Eq. (59), we get Thus the steady state quantal content distribution for the spike intervals drawn 215 from a Poisson processes is not a Binomial, but instead the above distribution in Eq.
Note that 1 is added manually inside the square bracket to take care of the 241 subsets with only one element m z , as these subsets are not taken into account by LHS =1 + k + 1 1 Hence, we have shown above that the identity is true for n = k + 1. Since integer 247 k ≥ 2 is arbitrary, we have proved the identity in Eq. (49) using the method of 248 strong induction.
Let us assume that the identity is true for an arbitrary value of n = k. For, 254 n = k + 1, the set {S} of all the subsets of {1, 2, 3, · · · k} can be divided into three Clearly, the contribution of (a) in the summation term in the LHS of Eq. (60) 261 will be equal to that for n = k, that is equal to 1/ (1 − f (i)) − 1. On the other 262 hand, the contribution of (b) will be equal to that of (a) multiplied by an additional 263 factor f (k) 1−f (k) . The contribution from (c) will be simply f (k) 1−f (k) . Adding all the three 264 contributions, Hence, the identity is true for n = k +1. Thus from the principle of mathematical 266 induction, the identity is proved for any integer n > 1.  The vesicle release probability p r t and the rate of docking of vesicles k t are de- . Given an AP train, with m th AP arriving 283 at t = t m , the dynamics of CaX F and CaX D are given by and, exponentially between two APs, with a decay constant τ F (or τ D ). In Fig. 1, for the 287 parameter values given in Ref. [1] for parallel fiber, we have plotted the p r t and k t as 288 a function of t. Note that the initial values at t → 0 are p r 0 and k 0 respectively. We 289 see that p r t → p r and k t → k in the long time limit. This steady state limit is the 290 focus of the entire analysis described in this paper, in which the release probability 291 and replenishment rate has attained the steady values p r and k respectively. The code used to simulate the distribution of quantal content for constant ISIs as 294 well as Poisson train is as follows.