A first-passage approach to diffusion-influenced reversible binding: insights into nanoscale signaling at the presynapse
AA first-passage approach to diffusion-influencedreversible binding for synaptic vesicle fusion
Maria Reva , David A. DiGregorio , and Denis S. Grebenkov Institut Pasteur, Unit of Synapse and Circuit Dynamics, CNRS UMR 3571, Paris, France Laboratoire de Physique de la Mati `ere Condens ´ee (UMR 7643) , CNRS – Ecole Polytechnique, IP Paris, 91128Palaiseau, France * senior authors ** correspondence: [email protected] and [email protected] ABSTRACT
Signal transduction driving synaptic transmission between neurons is governed by a cascade of stochastic reaction-diffusionevents that lead to calcium-induced vesicle release of neurotransmitter. As experimental measurements of such systemsare challenging due their nanometer and sub-millisecond scale, numerical simulations remain the principal tool for studyingcalcium dependent synaptic vesicle fusion, despite limitations of time-consuming calculations. In this paper we develop ananalytical solution to explore dynamical stochastic reaction-diffusion problems, based on the concept of first-passage times.This is the first analytical model that accounts simultaneously for relevant statistical features of calcium ion diffusion, buffering,and its binding/unbinding reaction with vesicular sensor. In particular, unbinding kinetics are shown to have a major impacton the calcium sensor’s occupancy probability on a millisecond scale and therefore cannot be neglected. Using Monte-Carlosimulations we validated our analytical solution for calcium influx through the voltage-gated calcium channel and instantaneousinflux. Overall we present a fast and rigorous analytical tool to study simplified reaction-diffusion systems that allows asystematic exploration of the biophysical parameters at the molecular scale, while correctly accounting for the statistical natureof molecular interactions within cells, that can also serve as a building block for more general cell signalling simulators.
Introduction
Intracellular transport of molecules is crucial for the normal function and growth of living cells . Many intracellularsignalling cascades, such as those mediated by calcium ions ( Ca + ), are generated by biochemical reactions and diffusion,which are often hard to accurately measure, particularly on the submillisecond and submicron temporal and spatial scales .Mathematical modeling can be used to interpret and predict features of intracellular signalling that are not yet observable. Oneparticularly interesting signaling process is the ability of neurons in the brain to communicate to each other by transformingelectrical into chemical signals and then back to electrical signals at specialized junctions called synapses. Electrical impulses,or action potentials (APs) in the presynaptic neuron recruit voltage gated calcium channels (VGCC) that mediate Ca + fluxacross the membrane followed by diffusion and binding to buffer molecules throughout the presynaptic terminal. When free Ca + reach their target, a Ca + sensor protein tethered to neurotransmitter containing synaptic vesicles (SVs), rapid fusion withthe plasma membrane is triggered, resulting in the release of neurotransmitter into the synaptic cleft. The neurotransmittermolecules diffuse towards receptors on the postsynaptic cell which then initiate an electrical signal. Ca + entry and diffusion tothe Ca + sensor are thought to occur within tens of nanometers on a millisecond timescale . Therefore, numerical simulationsof chemical reactions and Ca + diffusion have been essential for understanding the spatial-temporal dynamics of the calciumion concentration [ Ca + ] driving synaptic vesicle fusion. Ca + buffer molecules play an important role in shaping the spatial and temporal dynamics of intracellular unbound [ Ca + ]gradients generated by the opening of calcium channels , which can in turn influence SV fusion probability over time. Bufferscompete with SV fusion sensors for free Ca + . Buffer/sensor competition generally does not affect the unbinding SV rate.However, the occupancy of the sensor, which ultimately determines the probability of SV fusion, depends on both forward andbackward rate constants.In general, an analytical solution of Ca + reaction-diffusion equations is not possible. Therefore, deterministic andstochastic numerical simulations have been workhorses to study this problem. However, both strategies are time-consumingand suffer from inaccuracies under certain parameter regimes. The finite element methods do not account for the stochasticopening of VGCCs or fluctuations in Ca + flux. Moreover, it has been shown that such fluctuations should be simulatedexplicitly in order to accurately predict vesicle fusion probability . This motivates the use of stochastic simulators that arenotoriously time consuming and inaccurate for simulating low probability events associated with low numbers of molecules. a r X i v : . [ phy s i c s . b i o - ph ] A p r umerous methods and software suites have been developed over the past decades to simulate stochastic reaction-diffusionproblems. These can be generally divided in two groups: particle based and lattice based. In the particle based methods eachparticle is treated individually that makes computation time prohibitively expensive when the concentrations of particles arehigh and/or there are numerous species. The lattice based methods divide space in voxels and treat a number of moleculeswithin this voxel, rather than individuals . This approach can speed up the simulations but at the price of reduced spatialand temporal resolutions. In this context, analytical models can provide new insights and intuition on the complex system andhelp in building improved simulation techniques with higher speed and accuracy.The linearized buffer approximation (LBA ) yields an approximate analytical solution of a diffusion-reaction problem.This method allows one to compute the concentration of calcium in a chosen location without simulations. This can be extendedto multiple buffers and has provided important intuition about their spatial-temporal impact on intracellular Ca + , and potentialeffect on the probability of SV fusion. However, this approach has a significant limitation, as it can only be applied in the steadystate conditions, thus not suitable for the brief calcium entries followed by the action potential . Another recent multi-scaleapproach is based on the narrow escape problem of searching for a hidden target by a single calcium ion. An analyticalsolution of this problem was found and coupled with a Markovian jump process that models buffering and calcium influx . Inspite of its advantages, this hybrid method does not account for the Ca + sensor’s binding and unbinding kinetics, which canbe of crucial importance for the vesicle release dynamics, as shown below. Recent first-passage approaches have been usedto account for the finite backward rate constant of binding for multiple particles , but do not consider competing bindingpartners for diffusants and thus do not provide the intuitive power of LBA.Here we propose a probabilistic diffusion-influenced reversible calcium binding model that overcomes the aforementioneddeficiencies and accounts for the forward and backward binding rate constants of the Ca + sensor, as well as competing bindingpartners (fixed endogenous (EFB) and mobile buffers). This novel analytical model simulates a point source Ca + entry,reaction with buffer, diffusion and binding to a Ca + -binding protein that mediates SV fusion. The solution allows us to studythe effect of binding reaction rate constants on the occupancy probability of the sensor by Ca + at all temporal scales withoutcomputational costs. We confirm the necessity of taking into account the unbinding kinetics in the simulations of vesicle releaseprobabilities. We also demonstrate the validity of the analytical solution by Monte Carlo simulations and study the effectsof the sensor’s kinetics and geometrical properties of the synapse on the probability of the single site occupancy. Moreover,an extension to multiple calcium ions and its limitations are discussed. To our knowledge, this is the first analytical solutionfor a stochastic reaction-diffusion problem that accounts simultaneously for target binding/unbinding kinetics in the presenceof competing buffer species, and accurately predicts target occupancies following stochastic flux from an arbitrary spatialarrangement of ion channels. This approach is applicable to a wide range of diffusion-reaction problems within cells beyondthat of fast synaptic transmission. Results
Impact of unbinding kinetics on the single vesicle release probability and time course
Reversible first order chemical reactions are described by forward ( k on ) and backward ( k off ) rate constants. However, inthe case of Ca + diffusion and binding to a Ca + sensor for SV fusion it has been argued that the first passage time at thetarget is the dominant physical process influencing the probability of SV fusion over time, and thus approximations without k off might be sufficient. We tested the importance of Ca + sensor k off for AP-evoked SV fusion using a simple finite-elementsreaction-diffusion simulation. Spatio-temporal profiles of free [ Ca + ] were simulated for sensor distances of 10 and 50 nmfrom the perimeter of a VGCC cluster (perimeter model , see Fig.1A) in the presence of Ca + buffers (ATP or EFB).We modeled the probability of SV fusion within the five binding site kinetic model of the Ca + sensor , and comparedto a model in which k off s were set to zero. For sensor-to-channel distances (coupling distance, CD) as short as 10 nm, thetime course of SV release within the first millisecond is hardly different with and without a k off (Fig.1.C, blue lines), while therelease probability is increased 2.4 times (Fig.1.B, blue lines). However, in the case of 50 nm CD (which is physiological atsome synapses ), setting the k off to zero increased the single vesicle release probability by 7-fold (Fig.1.B, green lines),and increased the half-width of the time course of fusion probability by 61% (Fig.1.C, green lines). These simulations showthat for the shortest CDs, first passage time models that only consider k on could qualitatively reproduce the time evolution insensor occupancy, but not the final probability that a SV would fuse. However, for longer CDs both the time course and fusionprobability were altered in the absence of unbinding. Thus, a reversible Ca + binding reaction (finite k off ) must be consideredfor such simulations, particularly since the estimated distances range from 10 nm to 100 nm across synapses . In the followingsections, we investigate analytical solutions that describe, specifically, Ca + diffusion and consider explicitly the binding andunbinding kinetics of the sensor. CA
50 nmVGCCSensor (50 nm)Sensor (10 nm) P v r e l ea s e r a t e
50 nm no k off ( Cnt)10 nm no k off ( Cnt) 50 nm no k off ( Cnt)10 nm no k off ( Cnt)
Figure 1. A. Diagram of the active zone arrangement, with a cluster of VGCC (hollow circles) and two positions of SVsensors at 10 and 50 nm. B. Vesicle release probability for two scenarios: control (with unbinding reaction,"Cnt", dashed lines)and without unbinding reactions in the sensor kinetics (solid lines), and two CDs: 10 nm (blue) and 50 nm (green). C. Releaserates (corresponding to panel B), normalized to peak amplitude.
Analytical model of reaction-diffusion in synapse
We developed an analytical model of Ca + diffusion based on a modified first passage time process, in which the target Ca + sensor binding occupancy is modeled as a first order reaction with reversible kinetics (both k on and k off ) . This probabilisticdiffusion-influenced reversible calcium binding model is described in detail in the Methods Section (see also Fig. 2). In brief,we placed a single Ca + sensor, with the values of k on and k off taken from models agreeing with experimental data , at thecenter of the circular surface of a half sphere. The hindering effect of the synaptic vesicle was not considered, since it wasshown not to influence sensor occupancy . The dynamics of Ca + ions is described as switching diffusion between free andbuffer-bound states . In summary our model has the following parameters: the size of the sensor ρ , the distance betweenorigin of the simulation domain and calcium channel r , the radius of the simulation domain R , k on and k off of the Ca + sensor,the exchange rates k i and k i (product of concentration and forward rate constant of the buffer) for binding/unbinding to i -thbuffer, and diffusion coefficients of free Ca + ( D ) and those bound to buffers ( D i - diffusion coefficient of i -th buffer, Table 1). Figure 2. A. Schematic diagram of an axonal bouton (presynaptic terminal) containing a release site (active zone). Inset:Idealized active zone scheme showing VGCC clusters and their tens of nanometers proximity to the Ca + target sensor for SVfusion. B. Geometric representation of the simulation domain. The simulation compartment is a reflecting half sphere of radius R , target is a partially absorbing half sphere of radius ρ , the point source of Ca + entry is located on the membrane (horizontalsurface) at distance r − ρ from the sensor (target).In order to derive a set of equations describing the solution, we took a two-step approach. First, we found the probabilitydistribution of the first-passage time of a Ca + ion to a simplified (single binding site) sensor in the presence of competingbinding partners for a single Ca + entering the bouton. Second, a renewal technique allowed us to incorporate unbindingkinetics on the sensor and to relate the distribution of the first-binding time to the occupancy probability ( P ( t , r ) ) for a single Ca + ion started at some position r , to be on the sensor at time t (see Methods). The derived analytical solution (19) describeshow the occupancy probability evolves with time t and approaches its steady-state limit that we also found explicitly as P ∞ = (cid:18) + k off π ( R − ρ ) N A k on (cid:18) + M ∑ j = k j k j (cid:19)(cid:19) − , (1) arameter Notation Value Unit ReferenceGeometrical parameters Simulation domain radius R
300 nm adapted from Sensor’s radius ρ Coupling distance (CD) r − ρ
15 nmCalcium diffusion coefficient D µ m ms − Ca + sensor Forward rate constant k on ·
127 mM − ms − Backward rate constant k off − EFB ( i = ) Diffusion coefficient D µ m ms − Backward rate constant k
10 ms − Forward rate constant k on ,
100 mM − ms − Total concentration c Binding rate k
400 ms − ATP buffer ( i = ) Diffusion coefficient D µ m ms − Reverse rate constant k
10 ms − Forward rate constant k on ,
100 mM − ms − Total concentration c Binding rate k
20 ms − EGTA buffer ( i = ) Diffusion coefficient D µ m ms − Backward rate constant k − Forward rate constant k on , − ms − Total concentration c
10 mM Binding rate k
105 ms − Table 1.
Biophysical parameters of our diffusion-reaction model of the synapse.where N A is the Avogadro number. The sophisticated form of Eq. (19) highlights the complexity of the time-dependentdiffusion-reaction problem as compared to its steady-state limit, to which most former works were dedicated. While ourfirst-passage approach is applicable to any number of co-existing buffers, we focused on two cases of no buffer and singlebuffer, for which explicit analytic formulas for P ( t , r ) were provided. We explored the accuracy of these formulas and theirassumptions using Monte Carlo simulations (see Methods). By solving analytically the governing reaction-diffusion equations,we provide an efficient framework for studying and modeling the dynamics of calcium diffusion and binding a target, inparticular electrical impulse driven, calcium-ion-mediated SV fusion. Calculating the single ion occupancy probability for a simple Ca + sensor for SV fusion For our simplified vesicle fusion model we assumed that SV fusion probability was directly determined by the Ca + occupancy of target sensor, P ( t , r ) . Using our analytical solution, it was possible to calculate the P ( t , r ) across seven ordersof magnitude in time scales, from submicroseconds to seconds, using different model parameters. For the idealized case ofinstantaneous binding and no unbinding ( k on = ∞ , k off = Ca + ion that hits the sensor remains bound forever. As aconsequence, P ( t , r ) is equal to the probability of the first passage time to the sensor. As expected, this probability monotonicallyincreases with time and approaches 1 after one second (Fig. 3. A, black solid line), consistent with a pure diffusion-limitedreaction. The analytical solution is in excellent agreement with MC simulations realized with the same model parameters (Fig.3. A, black dashed line). When using a finite forward rate constant ( k on = ·
127 mM − ms − , k off = P ( t , r ) was reduced(Fig. 3. A, blue solid line), in perfect agreement with MC simulations (Fig. 3. A, blue dashed line). The use of both a finiteforward and backward rate constants generated a biphasic occupancy curve: the P ( t , r ) increased to a maximum value andfollowed by a decreased to a steady state value, as expected physiologically. The rising phase of this curve matched that ofthe curve when k off was set to zero (Fig.3. A, green solid line), thereby delineating the time scale where only ion binding isdominant. MC simulations reproduced the analytical solutions (Fig. 3. A, green dashed line), despite the inherent fluctuationsdue to a limited number of MC trials. hen increasing the size of the bouton (simulation volume), R =500 nm, the peak of P ( t , r ) was not altered, but thesteady-state limit P ∞ was decreased (1 · − for R =
300 nm, 3 · − for R =
500 nm, see Eq. (1)), consistent with alterationin the steady-state, volume-averaged [ Ca + ] (Fig. 3. B). For smaller bouton sizes ( R =100 nm), the peak of P ( t , r ) was increased(2 · − ). The rising phases, however, were identical for all tested radii, suggesting that diffusion determines the initial timecourse of SV fusion, provided that the bouton volume is much larger than the CD. CD = 5 nmCD = 15 nmCD = 25 nmCD = 45 nmCD = 95 nm P ( t, r) s10 -5 -1 -4 -6 -3 -2 -7 P ( t, r) No bufferATP ( MC)EGTAEFB ( MC)
No bufferATPEGTAEFB × -4 × -4 × -4 -3 -4 -5 -2 -1 P ( t, r) -3 -4 -5 -2 -1 s10 -5 -1 -4 -6 -3 -2 -7 k on = 5*12700 mM -1 ms -1 , k off = 1570 ms -1 k on = ∞ k off = 0k on = 5*127 mM -1 ms -1 , k off = 15.7 ms -1 k on = 5*0.127 mM -1 ms -1 , k off = 0.157 ms -1 P ( t, r) -3 -4 -5 -2 -1 s10 -5 -1 -4 -6 -3 -2 -7 R= 100 nmR= 300 nmR= 500 nm= 5*127 mM -1 ms -1 ,k off = 15.7 ms -1 ( MC)= 5*127 mM -1 ms -1 , k off = 0 ( MC)k on k on k on = ∞ k off = 0 ( MC) -5 -1 s s s10 -4 -6 -3 -2 -7 -5 -1 -4 -6 -3 -2 -7 -5 -1 -4 -6 -3 -2 -7 FDBEC A P ( t, r) P ( t, r) Figure 3.
Calculation of P ( t , r ) . A. P ( t , r ) computed using analytical solution (solid lines) and MC simulations (dashed lines)for the case of instant binding (black lines); in the presence of binding kinetics (blue lines); in the presence of both referentbinding and unbinding kinetics (green lines). B. Analytically computed P ( t , r ) for different sizes of the domain ( R ): 100 nm(blue), 300 nm (red) and 500 nm (green). C. P ( t , r ) computed using analytical solution for instant binding (black line), fastreaction rate constants (blue line), referent (green line) and slow reaction rate constants (brown line). D. Analytically computed P ( t , r ) for CDs, r − ρ , varying from 5 to 95 nm. E. P ( t , r ) computed using analytical solution (solid lines) and MC simulations(dashed lines) in the absence (blue) or presence of the following buffers: ATP (green), EGTA (gray) and EFB (red), at the CDof 15 nm. F. Similar to E, but at the CD of 45 nm. he occupancy probability was then computed for different pairs of forward and backward rate constants, each chosen suchthat equilibrium dissociation constant remains constant ( K D = k on / k off = µ M). In all cases, fast ( k on = · − ms − , k off = − ), the reference ( k on = ·
127 mM − ms − , k off = . − ), and slow ( k on = · .
127 mM − ms − , k off = .
157 ms − ) rate constants, the P ( t , r ) reaches the same equilibrium (Fig. 3. C). However, within the first few hundreds ofmicroseconds, faster rate constants enable a rapid capturing of Ca + as well as faster unbinding, thus generating biphastic P ( t , r ) that is larger, briefer and earlier (Fig. 3. C, blue line). Interestingly even with a forward rate-constant of greater than10 M − s − , the diffusion-limited case is not matched on the sub-microsecond time scale (Fig. 3. C, black line), due to thefast k off . Thus, the k off can strongly influence target occupancy and the regime of diffusion-limited binding.The distance between [ Ca + ] sources (VGCCs) and the Ca + sensor of SVs can vary in nature, influencing the kinetics andprobability of SV fusion . We, therefore, examined the effect of varying this VGCC-SV coupling distance (CD = 5-95 nm)on P ( t , r ) (Fig. 3. D). The peak of P ( t , r ) decreased from 0.027 at 5 nm to 0.001 at 95 nm. As expected from finite elementssimulations , longer CDs resulted in smaller peak fusion probabilities and longer times to peak fusion rates (6.1 µ s at 5 nmto 47.5 µ s at 95 nm).Finally, we explored the effect of various Ca + buffers on P ( t , r ) . Ca + buffers critically shape the spatio-temporalprofile of intracellular Ca + . We considered ATP, a naturally occurring low-affinity, fast and mobile endogenous calciumbuffer , non-specific low-affinity endogenous fixed buffers (EFB) , and the mobile exogenous buffer EGTA. EGTA is awell-characterized buffer, with slow forward Ca + binding rate constant, that has been used to infer VGCC-SV CDs throughcompetition with the Ca + sensor, thus producing an observed inhibition of synaptic transmission proportional to the CD .Because of its slow k on , large concentrations of EGTA are needed to interfere with the sensor for SV fusion (typically between1 mM and 10 mM ). For particle-based simulations this can be computationally prohibitive. The effect of all three buffers havebeen studied extensively, and thus have well-characterized binding rate constants , see Table 1. Using our analytical approachwe could rapidly calculate P ( t , r ) for a CD of 15 nm in the presence of 0.2 mM ATP and 10 mM EGTA. Both buffers onlyslightly decreased the peak amplitude of P ( t , r ) from 0.012 to 0.01 and shifted the time of its peak from 10 µ s to 8.5 µ s. Onother hand, high concentration of EFB (4 mM) had more prominent effect, decreasing the peak probability of being boundto (7 · − ) and shifting its time to 4 µ s (Fig. 3. E). These results are consistent with the lack of effect of ATP being due toits low concentration, the lack of effect of EGTA being due to its slow forward rate constant. At a CD of 45 nm, the P ( t , r ) peak was decreased from 1 · − to 3 · − (EGTA), 4 · − (EFB) and 6 · − (ATP); the time of peak was shifted from4.7 µ s to 2.1 µ s (EGTA), 9 µ s (EFB) and 4.5 µ s (ATP) (Fig. 3. F). The steady-state Ca + occupancy is dramatically reducedby the large concentration of the high affinity buffer, EGTA. These differential affects of EGTA on the peak occupancy for CDof 15 and 45 nm, as well as on the steady-state occupancy, are very similar to previous analytical results , finite elementssimulations , and MC simulations , due largely to the slow forward rate constant. The analytical solution was verified withMC simulations for ATP and EFBs (Fig. 3. E), but not for EGTA as the large number of molecules associated with 10 mMEGTA was too time-consuming for MC simulations. These results validate the use of the analytical approach to provide anintuitive understanding of stochastic reaction and diffusion across a wide range timescales and for large numbers of molecules,conditions that are prohibitive when using particle-based simulators. Temporal regime in which reaction-diffusion models must consider reversible kinetics
Equipped with an analytical solution, we reexamined the importance of k off in dictating P ( r , t ) . Figure 4 shows P ( r , t ) calculated for different k off s and for different CDs: 15 nm (Fig. 4. A), 45 nm (Fig. 4. B) and 95 nm (Fig. 4. C). The hightemporal resolution of the simulations show that there is a characteristic time window ( , t c ) at which the sensor occupancy isindependent of the k off . The upper limit t c of this characteristic time window was defined as the time point when two P ( t , r ) curves for different unbinding kinetics start to deviate (blue and green dots, Fig. 4). This characteristic time increases as k off decreases. For CDs less than 50 nm, the characteristic time window was less than 10 microseconds for physiological rateconstants (Fig. 4. A.B). The simplified first passage time approach confirms finite element simulations (Fig. 1) showing thatbackward rate constants in the physiological range can influence Ca + sensor’s occupancy for physiological source to targetdistances, and therefore must be modeled explicitly for accurate predictions of Ca + -dependent SV fusion. Calculation of Ca + sensor occupancy probability for multiple calcium ions Thus far we considered the case when only a single Ca + enters the presynaptic volume. We next explored performance ofour analytical solution for Ca + fluxes. During action potential-induced opening of a single channels we estimate approximately200 Ca + ions enter at a single VGCC driving release over a Gaussian-like time course (half-width 0.3 ms) . Knowing P ( r , t ) , for the single Ca + ion, we used Eq. (25) to approximate the probability, denoted as P N ( r , t ) , that at least one Ca + isbound to the sensor at time t following an instantaneous influx of N Ca + ions (see Methods). For an instantaneous flux of 200ions, P N ( t , r ) increased to nearly 1 for a CD of 15 nm, in contrast to the very low probabilities ( < .
01) in the single ion case. Asimilar peak was estimated using MC simulations, but the time course of analytical P N ( t , r ) was broader than that computedusing MC simulations (dashed line on Fig. 5. A; shaded region is standard error of the mean (SEM)). This difference between -5 -1 -4 -3 -2 -6 -4 -2 k off = 15.7 ms -1 k off = 0.157 ms -1 k off = 0 C P ( t, r) -6 -4 -2 -5 -4 -3 -2 -1 -5 -1 s s10 -4 -6 -3 -2 -7 -5 -1 -4 -6 -3 -2 k off = 15.7 ms -1 k off = 0.157 ms -1 k off = 0 k off = 15.7 ms -1 k off = 0.157 ms -1 k off = 0 BA P ( t, r) P ( t, r) Figure 4.
Influence of backward rate constant on P ( t , r ) . A. P ( t , r ) computed using analytical solution with fixed k on = ·
127 mM − ms − , and the following unbinding rates: k off = k off = .
157 ms − (red line) and k off = . − (blue line), for the CD of 15 nm. The departure between P ( t , r ) ’s is depicted by dots and vertical dashed lines.Blue dot: t c (cid:39) . µ s, green dot: t c (cid:39) µ s. B. Similar to A, for CD of 45 nm. Blue dot: t c (cid:39) . µ s, green dot: t c (cid:39) µ s. C. Similar to A, for CD of 95 nm. Blue dot: t c (cid:39) . µ s, green dot: t c (cid:39) µ s.the analytical solution and the MC simulation was smaller for the longer CD of 45 nm (Fig. 5. A, blue lines). This discrepancyidentifies a shortcoming of the analytical approximation for multiple Ca + ions, which does not account for binding exclusionwhen the target experiences another ion while already bound. In other words, the discrepancy between the MC and analyticalapproximation can be attributed to the saturation of the single binding site, which was not taken into account in the analyticalsolution for the Ca + influx. Since the probability that two ions might interact with the target is decreased in the presence of thecompeting buffer molecules, we tested whether the presence of physiological concentrations of EFBs influences the differencebetween the analytical approximation and MC simulations. Indeed, the presence of EFB decreases P N ( t , r ) for both CDs (Fig.5. B), as well as the discrepancy between analytical and simulated curves (Fig. 5, black lines). The error in the time courseestimate was still present for shorter CDs, but for the longer CD, the two curves are indistinguishable (Fig. 5, green lines). Thebetter accuracy in the presence of EFB can be attributed to lower binding probabilities experienced by the target sensor, whichis consistent with lower [ Ca + ]. BA P N ( t, r) P N ( t, r) Analyt. CD =15 nmMC CD =15 nm ( SEM)Analyt. CD =45 nmMC CD =45 nm ( SEM) Analyt. CD =15 nmMC CD =15 nm ( SEM)Analyt. CD =45 nmMC CD =45 nm ( SEM) BA no EFB 4 mM EFB -8 -6 -4 -2 s -8 -6 -4 -2 s Figure 5.
Occupancy probability ( P N ( t , r ) ) for an instantaneous influx of many Ca + ions ( N = A. In the absence of fixed buffer, forCD of 15 nm (red lines) and CD of 45nm (blue lines). B. Similarly in the presence of EFB (4 mM), CD of 15 nm (black lines)and CD of 45 nm (green lines).Because both the presence of a competing Ca + buffer and increasing the distance to the target would be expected to reducefirst passage time probabilities, we next explored how the number of released Ca + ions and unbinding kinetics influencedthe discrepancy between MC and analytical P N ( t , r ) solutions. The strength of applied instantaneous Ca + influx was varied(50, 100, and 200 ions) and the difference between MC and analytical curves was quantified by the Mean Absolute difference(MAE) and full width at half maximum (FWHM) error. We saw that with decreasing number of released ions, the dissimilaritydecreases for all values of CDs, reflected in the values of MAE and FWHM error (Fig. 6). However, for long CDs we notice adecrease in the FWHM error with increasing number of released ions (from 81% for 50 ions to 53% for 200 ions) (Fig 6. C), his is due to reduced trial variability in MC simulations arising from a higher P N ( t , r ) . Moreover, altering P N ( t , r ) by adjusting k off (slow (0 .
157 ms − ) and fast (1570 ms − )) was also consistent with the primary source of error being due to high occupancy(see Figs. 1 and 2 in the SI).In summary, these simulations show that the analytical solution exhibits deviations under conditions that generate peak of P N ( t , r ) greater than 0.5, in which case errors in the half-width of the P N ( t , r ) could be larger than 25%. In contrast, when thetarget occupancy is less than 50%, the analytical solution yields accurate results even for a large number of Ca + . -8 -6 -4 -2 × -3 × -3 × -3 × -3 CD =15 nm
A B C i on s CD =45 nm CD =95 nm -8 -6 -4 -2 -8 -6 -4 -2 × -3 × -3 × -3 × -2 Analyt.MC i on s -8 -6 -4 -2 -8 -6 -4 -2 -8 -6 -4 -2 -8 -6 -4 -2 -8 -6 -4 -2 -8 -6 -4 -2 i on s ms ms s s ss s ss s s msms ms msms ms ms P N ( t, r) P N ( t, r) P N ( t, r) . % ( F W H M ) . ( M AE ) Figure 6. P N ( t , r ) for various numbers of ions and differing coupling distances (CD). P N ( t , r ) for 50 (1st row), 100 (2nd row)and 200 (3rd row) simultaneously released ions for CD of 15 nm ( A ), 45 nm ( B ) and 95 nm ( C ). Main plots represent semi-logscale, while linear scale plots are on insets. Black and green lines show respectively analytical and MC results. The black andblue inset text on each plot represent FWHM error and MAE between analytical and MC results correspondingly. Analytical solution to SV fusion sensor occupancy for Ca + entry via stochastic gating of calcium channels Thus far we considered instantaneous entry of Ca + ions, however, in neurons Ca + fluxes arise from a gradual openingof VGCCs (Fig. 7. A) throughout the duration of an electrical impulse or action potential. Moreover, it is also known thataccurate estimates of the occupancy probability must consider the stochastic nature of VGCC opening, particularly as comparedto deterministic approximations of mean open channel probability . Here, we studied the analytical solution, denoted as P AP ( r , t ) , for Ca + fluxes generated from stochastic opening of VGCCs (see Methods). VGGC openings and associated Ca + fluxes were obtained from MC simulations. The VGCC model was constrained by experimental estimates of single channelopen probability, single channel conductance and duration of the current (see Methods and Section V.A of the SI). For each a + ion entry time the P ( t , r ) ’s were calculated, and the final occupancy probability P AP ( r , t ) was the average over all theaforementioned P ( t , r ) ’s calculated at a given time instance t (see Eq. (27) in Methods). SensorVGCCAP
10 us m V
15 nm B A ICa ms P AP ( t, r) Figure 7. Ca + entry through a stochastic VGCC. A. Opening of the single VGCC is driven by an AP (black line) triggering Ca + influx (gray line). Released Ca + ions diffuse towards a sensor at distance of 15 nm. B. The P AP ( r , t ) computed usinganalytical approximation (red) and MC (blue).For a single VGCC located 15 nm from the sensor (Fig. 7. A), in the presence of EFB, the calculated P AP ( t , r ) was similarto that from MC simulations (Fig. 7. B), and much smaller than the peak occupancies predicted from an instantaneous Ca + flux. Thus our analytical method is capable of describing how a simplified SV fusion sensor could be driven by stochasticallygated channels. Moreover, the peak occupancies are two orders of magnitude smaller than the calculation for an instantaneousentry of 200 ions (Fig. 6. A), suggesting that the errors due to multiple Ca + binding are minimal. The extension to multiplechannels is simple via the principle of superposition , provided the total sensor occupancy remains less than 50%. Discussion and Conclusion
In this work, we introduced an analytical framework for computing a Ca + binding occupancy of a Ca + sensor for SVfusion following reactions and diffusion from a VGCC source. The probabilistic diffusion-influenced reversible calcium bindingmodel is based primarily on a modified first passage method . The main novelty of this approach is its ability to account forbinding and unbinding kinetics of the sensor in the presence of competing Ca + buffers. The new method was validated usingMonte Carlo-based particle reaction-diffusion simulation.We demonstrated that this method has several advantages over traditional deterministic and particle-based simulations: (i)the analytical nature of the solution facilitates exploration of diffusion and binding over a wide range of timescales (nanosecondsto seconds), a range that is often inaccessible, tedious or inaccurate when using numerical methods; (ii) explicit calculation ofprobabilities facilitates more accurate simulation of biological stochasticity, particularly for rare events; and (iii) our mean-fieldapproach circumvents the limitations of tracking the diffusion of thousands of molecules that represent millimolar concentrationsfound in biology.Significant variations of the sensor occupancy probability for physiologically relevant parameters, and over timescalesfrom nanosecond to seconds, can be easily calculated by our analytical solution without increasing the computation time orgenerating integration errors from small time steps. The peak of P ( t , r ) was usually at few tens of microseconds, a timescalecomparable to experimental measurements of SV fusion times, but a thousand times smaller than the mean first-passage time,which gives a measure of the average rate occurrence of a stochastic event (in our setting, the sensor binding). This is a strikingbiological example of a process in which the mean first-passage time is misleading while the commonly employed reduction ofthe first-passage time distribution to its mean is not permitted. Similarly, the focus on the steady-state limit can be justified insome applications, but understanding of vesicle release requires attention to its dynamics over a very broad range of time scales.In this light, analytical models offer unique opportunities for understanding generic features of diffusion-influenced reactions inbiological systems.Strictly speaking, our analytical solution describes the occupancy probability for a single sensor that can bind an unlimitednumber of Ca + . This, however, is not compatible with the finite number of binding sites for most biochemical reactions.Therefore the fact that each binding site can be occupied only by a single Ca + at any given moment, highlights a limitation ofthe analytical approach. At the same time, we showed that for P N ( t , r ) less than 0.5 the discrepancy between MC and analyticalsolution simulations is minimal. In particular high P N ( t , r ) are achieved with large molecular fluxes and short CDs (<20 nm). he difference can be reduced by adding competing Ca + buffers, as they effectively lower the peak occupancy of the sensor.Future mathematical studies could provide statistical approaches to account for molecular occlusion that occurs when P ( t , r ) ishigh (>0.5). Nevertheless, the current model is accurate for particular conditions such as low Ca + fluxes (corresponding to0.3 ms single channel currents of 0.3 pA), long CDs (>40 nm) and in the presence of the EFB. Also, the method was shownto be accurate for the case of Ca + fluxes via single stochastic VGCC, even at 15 nm. This seems likely due to the lowerinstantaneous flux that occurs for time-distributed influx during an action potential.Our method could be extended to include multiple SV fusion sensors, through the definition of the spatial domains aroundeach sensor, as proposed earlier . The size of this domain can be given by the distance where an ion can contribute to thebinding to the particular sensor. Accounting for the saturation of the single binding site of the sensor is another importantperspective which can be explored by adapting recently proposed models for the dynamics of impatient particles . Inparticular, Lawley and Madrid suggested to model the distribution of the first-passage time onto the target by a mono-exponential function, in which case the number of bound particles can be described by a Markov birth-death process, for whichthe first-passage time statistics are well known . However, the accuracy of such an approximation remains questionable in oursetting, especially for re-binding steps when the particle unbinds from the sensor and diffuses in the bulk until the next binding.Once the bound Ca + occlusion problem is correctly implemented, it would then be feasible to investigate the cooperativebehavior generated by multiple binding sites .Finally, the time efficiency of our method is superior to the existing MC alternatives. The duration of the MC simulationsdepends on the complexity of the system. In our case, brute force MC simulations with 4mM of fixed buffer took several hoursto complete and other simulations with higher buffer concentrations were not feasible. In contrast, our analytical method allowsone to achieve same results as MC in a fraction of minute, without limitations on molecular concentrations. This opens doorsfor exploration of complex diffusion-reaction systems that were previously out of reach. Methods
Theoretical model (i)
Geometric settings . A Ca + ion (or multiple ions) was injected on the membrane on the distance r from the centerof a synaptic terminal (bouton) with radius R, and allowed to diffuse throughout. At the center, we placed a single Ca + sensor of hemispherical shape and radius ρ , a value analogous to an interaction radius (Fig. 2. A). The outer boundary of thebouton (a hemi-spere of radius R ) is modeled as reflecting, i.e., the flux of Ca + ions at this boundary is zero. Note that moreelaborate partially reflecting boundary could also be considered to account for Ca + ions escaping far from the synaptic boutonmembrane but we strick to the reflecting condition here. In summary, we consider the active zone of the shape (Fig. 2. B) Ω = { xxx = ( x , y , z ) ∈ R : ρ < | xxx | < R , z > } . (2)Importantly, we neglect the presence of the synaptic vesicle whose reflecting boundary might hinder the motion of Ca + ions;in fact, it has been shown by Monte Carlo simulations that in physiological conditions, the synaptic vesicle does not influencethe single vesicle release probability .(ii) Ca + ions are modeled as independent point-like diffusing particles that undergo Brownian motion with diffusioncoefficient D in the region Ω between the boundaries of sensor and active zone; in particular, the charge of Ca + ions isignored due to bulk screening of electrostatic interactions. A Ca + ion source (VGCC) was set at a fixed distance r − ρ fromthe sensor (we discuss below how to deal with multiple sources).(iii) Buffers are modeled as co-existing continuous homogeneous reactive media that can bind, transport, and release Ca + ions; their functioning is assumed to be in a linear regime (i.e., low occupancy), i.e. the exchange between the free Ca + state(denoted by index 0) and the bound state with the i -th buffer (denoted by index i ) occurs through the standard first-order kinetics,with the exchange rates k i and k i ; the Ca + ion in a bound state diffuses with the diffusion coefficient D i but cannot bindto the sensor. Under the assumption of a homogeneous reactive medium, the “binding rate” can be expressed as k i = k on , i c i ,where k on , i is the conventional binding constant and c i is the concentration of the i -th buffer.(iv) Sensor kinetics . We consider a sensor with a single binding site, its kinetics is determined by k on and k off binding rateconstants. In most cases (unless stated otherwise), we used reaction rate constants identical to the reaction rate constant of thefirst binding site from the 5 state sensor model (see also Section V.B of the SI). When the Ca + ion reaches the surface of thesensor it can be reflected from it or bind to it, the random choice of either depending on the sensor binding constant k on33,34 .When bound, the Ca + ion remains in this state for a random time τ distributed by an exponential law Φ ( t ) = P { τ > t } = exp ( − k off t ) , (3)1 / k off being the mean waiting time before unbinding reaction. After the unbinding from the sensor, the Ca + ion resumes itsdiffusion in the intrasynaptic region Ω until the next binding event. ensor occupancy probability We are interested in computing the so-called occupancy probability P ( t , xxx ) that a particle (here, a Ca + ion), started froma fixed point xxx at time 0, is at the bound state on the sensor at a later time t . Due to the sensor kinetics, the particle canundergo numerous binding/unbinding events up to time t . To account for these events, we introduce an auxiliary probabilitydensity ψ n ( t , xxx ) of the n -th binding at time t . These densities can be obtained via recurrent functional relations. In fact, theindependence between the time spent in the bound state on the sensor, and the time of a bulk excursion after unbinding, implies ψ n ( t , xxx ) = t (cid:90) dt t (cid:90) t dt ψ n − ( t , xxx ) φ ( t − t ) ψ ( t − t ) , (4)where φ ( t ) = k off exp ( − k off t ) is the probability density of the exponential waiting time in the sensor-bound state, and ψ ( t ) isthe probability density of re-binding at time t after the release at time 0. This is a standard renewal relation, which states that,after the ( n − ) -th binding of the particle at some time t (with the density ψ n − ( t , xxx ) ), the particle remains bound during time t − t and unbinds at time t (with the density φ ( t − t ) ), diffuses in the bulk during time t − t and re-binds at time t (withthe density ψ ( t − t ) ). Since the intermediate binding/unbinding events may occur at any times between 0 and t , one has tointegrate over t and t . The integral relation (4) is reduced to a product in the Laplace space, i.e.,˜ ψ n ( p , xxx ) = ˜ ψ n − ( p , xxx ) ˜ φ ( p ) ˜ ψ ( p ) = ˜ ψ ( p , xxx ) (cid:2) ˜ φ ( p ) ˜ ψ ( p ) (cid:3) n − , (5)where tilde denotes Laplace-transformed quantities, e.g.,˜ ψ n ( p , xxx ) = ∞ (cid:90) dt e − pt ψ n ( t , xxx ) . (6)The probability of a particle to be in the bound state at time t can be expressed as follows P ( t , xxx ) = ∞ ∑ n = t (cid:90) dt (cid:48) ψ n ( t (cid:48) , xxx ) Φ ( t − t (cid:48) ) . (7)In this infinite sum, the n -th term is the probability that after the n -th binding at time t (cid:48) (with the density ψ n ( t (cid:48) , xxx ) ), the particleremains at the bound state for time t − t (cid:48) (with the probability Φ ( t − t (cid:48) ) given by Eq. (3)). This relation simply reflects thefact that the particle, which is at the bound state at time t , has experienced either 1, or 2, . . . or n , or . . . binding events. In theLaplace space, we get˜ P ( p , xxx ) = ∞ ∑ n = ˜ ψ n ( p , xxx ) ˜ Φ ( p ) = ˜ ψ ( p , xxx ) (cid:18) − ˜ φ ( p ) ˜ ψ ( p ) (cid:19) − ˜ Φ ( p ) . (8)The three factors in the product have a clear interpretation: the first arrival and binding onto the sensor, multiple re-bindingevents on the sensor, and waiting after the last re-binding. For exponential waiting times, one has ˜ Φ ( p ) = / ( p + k off ) and˜ φ ( p ) = k off / ( p + k off ) , and thus we finally get˜ P ( p , xxx ) = ˜ ψ ( p , xxx ) (cid:18) p + k off ( − ˜ ψ ( p )) (cid:19) − . (9)The inverse Laplace transform of Eq. (9) allows one to return to the time domain to get P ( t , xxx ) . In this way, the probability P ( t , xxx ) is reduced to the analysis of the “elementary” diffusion step – the binding to the sensor – that determines both ˜ ψ ( p , xxx ) and ˜ ψ ( p ) . We emphasize that Eq. (9), written in terms of survival probabilities, is well known for describing reversible kineticsin chemical physics, see and references therein. Distribution of the first-binding time
We start by noting that the symmetry of the considered domain Ω allows one to effectively remove the synaptic boutonmembrane at z = Ω by a simpler spherical layer Ω = { xxx = ( x , y , z ) ∈ R : ρ < | xxx | < R } . (10) n other words, the distribution of first-passage times computed in Ω is identical to that computed in Ω . The advantage of thelatter domain is that it is rotation invariant so that the problem can be reduced to one-dimensional radial part, as discussedbelow. When there is no buffer, the computation of the first-passage time to a target is rather standard but technicallyinvolved in the case of a spherical layer . Accounting for buffers presents one of the major challenges and originalities of thiswork. Note that our approach generalizes some earlier results for two-channel diffusion .We investigate the model with M distinct buffers by using an ( M + ) -state switching diffusion model: the Ca + ion can beeither in a free state (0) or in a buffer-bound state ( i ), with i = , . . . , M . Given that the buffers are modeled as continuous andhomogeneous media, a transition from the state i to the state j happens spontaneously, with a given rate k i j (see Section I of theSI for a formal definition of the model). A general scheme for studying first-passage times for switching diffusions was recentlydeveloped in . We introduce ( M + ) survival probabilities S i ( t , xxx ) for a Ca + ion started at xxx in the state i to be unbound fromthe sensor until time t . These probabilities satisfy ( M + ) coupled backward Fokker-Planck (or Kolmogorov) equations : ∂ t S i = D i ∆ S i + M ∑ j = k i j ( S j − S i ) ( i = , , . . . , M ) , (11)subject to the initial condition: S i ( , xxx ) =
1. Here D i is the diffusion coefficient of the Ca + ion in the state i , ∆ is the Laplaceoperator, and we set k ii = Ca + ion exchange between bound states: k i j = ( ≤ i , j ≤ M ) . (12)In other words, any exchange between the states i and j occurs through the free state 0. The last term in Eqs. (11) describestransitions between states i and j .Equations (11) should be completed by boundary conditions at the inner sphere at | xxx | = ρ (the sensor) and the outer sphereat | xxx | = R (the frontier of the active zone). The outer reflecting boundary simply confines the Ca + ions within the active zone,i.e., it ensures that there is no flux of Ca + ions across this boundary: − D i ∂ n S i ( t , xxx ) = ( | xxx | = R , i = , , . . . , M ) , (13)where ∂ n is the normal derivative directed outwards the domain. Since the Ca + ions in bound states cannot bind to the sensor,the same Neumann boundary condition is imposed at the inner sphere: − D i ∂ n S i ( t , xxx ) = ( | xxx | = ρ , i = , . . . , M ) . (14)Finally, the calcuim ions in the free state can bind to the sensor that implies the Robin boundary condition − D ∂ n S ( t , xxx ) = k on S ( t , xxx ) N A ( πρ ) ( | xxx | = ρ ) . (15)It is obtained by equating the net diffusive flux at the sensor (left-hand side) to the reactive flux (right-hand side) controlled bythe reaction constant k on , where N A is the Avogadro number, and 4 πρ is the surface area of the sensor . We emphasize thatthe presence of buffers has two effects: change in the diffusion coefficient and impossibility of a buffer-bound Ca + ion to bindto the sensor. Since the calcuim ions are released in the free state, we are interested exclusively in S ( t , xxx ) . However, findingthis probability requires solving the coupled system of equations for all S i . Note that 1 − S ( t , xxx ) describes the fraction of Ca + ions that have been bound to the sensor up to time t . This is the cumulative probability distribution for the first-binding time. Inparticular, its probability density reads ψ ( t , xxx ) = ∂ t (cid:0) − S ( t , xxx ) (cid:1) = − ∂ t S ( t , xxx ) . (16)Due to the rotational invariance of the problem, the probabilities S i ( t , xxx ) and the probability density ψ ( t , xxx ) , written in sphericalcoordinates, depend only on the radial coordinate r = | xxx | . From now on, we replace xxx by r .The solution of the system (11) of coupled partial differential equations is detailed in Section II of the SI. In a nutshell, theLaplace transform reduces these equations to a system of ordinary differential equations with respect to the radial coordinate r that is then solved by standard methods. Once the solution is found, one gets from Eq. (16)˜ ψ ( p , r ) = − p ˜ S ( p , r ) . (17)In addition, as a bulk excursion after the unbinding event starts at the sensor surface, r = ρ , one has˜ ψ ( p ) = ˜ ψ ( p , ρ ) . (18) s a consequence, the knowledge of ˜ S ( p , r ) yields both dynamical characteristics, ˜ ψ ( p , r ) and ˜ ψ ( p ) , that determine theLaplace-transformed probability ˜ P ( p , r ) according to Eq. (9).The last step for getting P ( t , r ) in time domain requires the inverse Laplace transform of ˜ P ( p , r ) . This is performed bydetermining the poles of this function and applying the residue theorem. When the poles are simple, the occupancy probabilityadmits the following exact representation: P ( t , r ) = P ∞ + ∞ ∑ n = exp (cid:0) − α n D t / ρ (cid:1) M ∑ j = b ( j ) n u ( α ( j ) n , r ) , (19)where u ( α , r ) = ρ sin (cid:0) α R − r ρ (cid:1) − R α cos (cid:0) α R − r ρ (cid:1) r , (20)and the steady-state limit P ∞ is given by Eq. (1). The coefficients b ( j ) n and α ( j ) n are determined by exact but complicatedformulas provided in the SI, whereas α n are found as strictly positive solutions of some trigonometric equation (provided in theSI). For instance, in the simplest case of no buffer ( M = α ( ) n = α n , b ( ) n = µ sin ( α n β ) (cid:0) α n w + w (cid:1) + α n cos ( α n β ) (cid:0) α n w + w (cid:1) , (21)where w = ( + β ) + β ( β + µ ( + β )) , (22a) w = ( + µ − λ ( + β )) − λ β , (22b) w = β ( + β ) , (22c) w = β ( + µ − λ ( + β )) − ( β + µ ( + β )) , (22d)with dimensionless parameters β = ( R − ρ ) / ρ , λ = k off ρ / D , µ = k on / ( πρ D N A ) , (23)and α n are strictly positive solutions of the trigonometric equationsin ( α n β ) = (cid:2) α n ( β + µ ( + β )) − λ β (cid:3) α n cos ( α n β ) α n ( + β ) + α n ( + µ − λ ( + β )) − λ ( n = , , . . . ) . (24)For a single buffer ( M = M >
1, numerical computations based on ouranalytical solution remain fast and accurate.The exact solution (19) is the main analytical result of the paper. Although this solution may look cumbersome and involvessome numerical steps (truncation of the infinite series, numerical computation of the coefficients, etc), its explicit form allowsfor both analytical and numerical investigation of the occupancy probability P ( t , r ) . Instant calcium influx If N independent ions are released simultaneously from the same fixed position r , then the single site occupancy probabilitycan be computed as a probability of at least one out of N ions being bound to the sensor, according to the formula: P N ( t , r ) = − ( − P ( t , r )) N , (25)where P ( t , r ) is the occupancy probability of single binding site by a single ion.In the MC simulations P N ( t , r ) was computed as the probability of finding a single ion bound to the sensor, given that N ions were released at the same time from the same position. alcium influx through single VGCC Single VGCC was modeled as a three-step Hodgkin-Huxley process (see also Section V.A of the SI), the model andparameters are taken Ref. . This model reproduces the single channel characteristics that were measured previously: openingprobability 0.3, maximum of the single channel current 0.3 pA, full width half maximum of the single channel current 250 µ s.Simulation of Ca + influx through the channel was done using MC tool, with 1000 trials. For the j -th trial, we storedthe random times t ( j ) , . . . , t ( j ) N when N ions “entered” the system. Given the instances of ion appearance we calculated theprobability of sensor occupancy for this trial using Poisson binomial distribution at each point: P ( j ) trial ( t , r ) = − N ∏ i = ( − P ( t + t ( j ) i , r )) . (26)Then these probabilities were averaged among the trials: P AP ( t , r ) = ∑ j = P ( j ) trial ( t , r ) . (27)More generally, if the ions entered from different VGCC channels, one could use P ( t + t ( j ) i , r ( j ) i ) with the appropriate location r ( j ) i of the source of the i -th ion in the j -th trial. In this way, one can easily implement sophisticated spatio-temporal characteristicsof the Ca + ions release. Stochastic simulations
For verification of analytical results we use particle-based stochastic numerical simulations (MCell software ). In MCelldiffusion of individual molecules is modeled using Brownian dynamics, while chemical reactions occur due to the collisionof molecules and follow Poisson distribution. All the parameters for simulations are identical to the parameters of analyticalsolution. The presynaptic domain of radius 300 nm and sensors are modeled as spheres, intersected by a reflecting plane inthe origin (Fig. 2.B). The sensor is located in the origin of the volume and has a radius of 5 nm. Depending on the context ofthe simulation, the Ca + input, number of calcium channels and the distance between the sensor and calcium channels weremanipulated; for instance, for the computation of the occupancy probability by a single ion, it was released at time 0 from asingle source. Each time the particle hits the sensor was recorded. The first-passage time distribution was computed based onthe recorded times.To compute occupancy probabilities we stored the time instances of the reaction between Ca + ion and the sensor, then thenumber of binding events at each time instance was divided by the total number of trials. The interaction range between twoparticles was set to 5 nm, the time step was chosen to be 5 ns. Deterministic simulations
The release rates were simulated using a 5-state model of Ca + triggered vesicle fusion . The model was integratedusing forward Euler scheme in a customary Matlab routine. The input Ca + transients are results of the spatial deterministicsimulations for the channel-vesicle arrangement as in Fig. 1, provided by Yukihiro Nakamura. Acknowledgements
This study was supported by the Centre National de la Recherche Scientifique, Fondation pour la Recherche Medicale(Equipe FRM), Agence Nationale de la Recherche (ANR-2010-BLANC-1411, ANR-13-BSV4-0016, ANR-17-CE16-0019,and ANR-17-CE16-0026), and Ile de France (Domaine d’Interet Majeur (DIM) MALINF: DIM120121). M.R. was supportedby the Pasteur Paris University (PPU) doctoral program. The laboratory of D.D. is a member of the Bio-Psy Laboratory ofExcellence. We would like to thank Yukihiro Nakamura for sharing simulated calcium transients, Jean-Baptiste Masson andNelson Rebola for the comments on the manuscript and discussions.
Author contributions statement
M.R, D.D. and D.G. conceived the study, analysed results and wrote manuscript, D.G. made analytical derivations, M.R.performed simulations.
Additional information
The authors declare no competing interests. eferences Alberts, B. et al.
Specialized tissues, stem cells, and tissue renewal.
Molecular biology of the cell . 5th ed. New York:Garland Science, Taylor & Francis Group (2008). Berridge, M. J., Bootman, M. D. & Roderick, H. L. Calcium: calcium signalling: dynamics, homeostasis and remodelling.
Nat. Rev. Mol. Cell Biol. , 517 (2003). Eggermann, E., Bucurenciu, I., Goswami, S. P. & Jonas, P. Nanodomain coupling between Ca2+ channels and sensors ofexocytosis at fast mammalian synapses.
Nat. Rev. Neurosci. , 7–21 (2012). Nakamura, Y. et al.
Nanoscale distribution of presynaptic Ca2+ channels and its impact on vesicular release duringdevelopment.
Neuron , 145–158 (2015). Roberts, W. M. Localization of calcium signals by a mobile calcium buffer in frog saccular hair cells.
J. Neurosci. ,3246–3262 (1994). Matveev, V., Zucker, R. S. & Sherman, A. Facilitation through buffer saturation: constraints on endogenous bufferingproperties.
Biophys. J. , 2691–2709 (2004). Dittrich, M. et al.
An excess-calcium-binding-site model predicts neurotransmitter release at the neuromuscular junction.
Biophys. J. , 2751–2763 (2013). Andrews, S. S., Addy, N. J., Brent, R. & Arkin, A. P. Detailed simulations of cell biology with smoldyn 2.1.
PLoS Comput.Biol. , e1000705 (2010). Modchang, C. et al.
A comparison of deterministic and stochastic simulations of neuronal vesicle release models.
Phys.Biol. , 026008 (2010). Blackwell, K. T. An efficient stochastic diffusion algorithm for modeling second messengers in dendrites and spines.
J.Neurosci. Methods , 142–153 (2006).
Chen, W. & De Schutter, E. Parallel steps: large scale stochastic spatial reaction-diffusion simulation with high performancecomputers.
Front. Neuroinformatics , 13 (2017). Naraghi, M. & Neher, E. Linearized buffered Ca2+ diffusion in microdomains and its implications for calculation of[Ca2+] at the mouth of a calcium channel.
J. Neurosci. , 6961–6973 (1997). Nakamura, Y., Reva, M. & DiGregorio, D. A. Variations in Ca2+ influx can alter ca2+-chelator-based estimates of ca2+channel-synaptic vesicle coupling distance.
J. Neurosci.
Holcman, D. & Schuss, Z. The narrow escape problem.
SIAM Rev. , 213–257 (2014). Metzler, R., Oshanin, G. & Redner, S. E.
First-Passage Phenomena and Their Applications (World Scientific Press, 2014).
Grebenkov, D. S. First passage times for multiple particles with reversible target-binding kinetics.
J. Chem. Phys. ,134112 (2017).
Guerrier, C. & Holcman, D. Hybrid markov-mass action law model for cell activation by rare binding events: Applicationto calcium induced vesicular release at neuronal synapses.
Sci. Reports (2016). Grebenkov, D. S. First passage times for multiple particles with reversible target-binding kinetics.
J. Chem. Phys. ,134112 (2017).
Lawley, S. & Madrid, J. First passage time distribution of multiple impatient particles with reversible binding.
J. Chem.Phys. , 214113 (2019).
Rebola, N. et al.
Distinct nanoscale calcium channel and synaptic vesicle topographies contribute to the diversity ofsynaptic function.
Neuron , 693–710 (2019).
Wang, L.-Y., Neher, E. & Taschenberger, H. Synaptic vesicles in mature calyx of held synapses sense higher nanodomaincalcium concentrations during action potential-evoked glutamate release.
J. Neurosci. , 14450–14458 (2008). Vyleta, N. P. & Jonas, P. Loose coupling between ca2+ channels and release sensors at a plastic hippocampal synapse.
Science , 665–670 (2014).
Yin, G. & Zhu, C.
Hybrid Switching Diffusions: Properties and Applications (Springer, New York, 2010).
Grebenkov, D. S. A unifying approach to first-passage time distributions in diffusing diffusivity and switching diffusionmodels.
J. Phys. A: Math. Theor. , 174001 (2019). Klafter, J. & Sokolov, I. M.
First steps in random walks: from tools to applications (Oxford University Press, 2011). Shahrezaei, V. & Delaney, K. R. Consequences of molecular-level Ca2+ channel and synaptic vesicle colocalization for theCa2+ microdomain and neurotransmitter exocytosis: a monte carlo study.
Biophys. J. , 2352–2364 (2004). Araç, D. et al.
Close membrane-membrane proximity induced by Ca2+-dependent multivalent binding of synaptotagmin-1to phospholipids.
Nat. Struct. & Mol. Biol. , 209 (2006). Allbritton, N. L., Meyer, T. & Stryer, L. Range of messenger action of calcium ion and inositol 1, 4, 5-trisphosphate.
Science , 1812–1815 (1992).
Xu, T., Naraghi, M., Kang, H. & Neher, E. Kinetic studies of Ca2+ binding and Ca2+ clearance in the cytosol of adrenalchromaffin cells.
Biophys. J. , 532–545 (1997). Nägerl, U. V., Novo, D., Mody, I. & Vergara, J. L. Binding kinetics of calbindin-d 28k determined by flash photolysis ofcaged Ca2+.
Biophys. J. , 3009–3018 (2000). Bucurenciu, I., Kulik, A., Schwaller, B., Frotscher, M. & Jonas, P. Nanodomain coupling between Ca2+ channels andCa2+ sensors promotes fast and efficient transmitter release at a cortical gabaergic synapse.
Neuron , 536–545 (2008). Grebenkov, D. S. Reversible reactions controlled by surface diffusion on a sphere.
J. Chem. Phys. , 154103 (2019).
Grebenkov, D. S., Filoche, M. & Sapoval, B. Spectral properties of the brownian self-transport operator.
Eur. Phys. J. B , 221–231 (2003). Grebenkov, D. S. Imperfect diffusion-controlled reactions. In Lindenberg, K., Metzler, R. & Oshanin, G. (eds.)
ChemicalKinetics: Beyond the Textbook , chap. 8, 191–219 (World Scientific, 2019).
Agmon, N. & Szabo, A. Theory of reversible diffusion-influenced reactions.
J. Chem. Phys. , 5270 (1990). Prüstel, T. & Tachiya, M. Reversible diffusion-influenced reactions of an isolated pair on some two dimensional surfaces.
J. Chem. Phys. , 194103 (2013).
Redner, S.
A Guide to First Passage Processes (Cambridge University press, Cambridge, 2001).
Crank, J.
The Mathematics of Diffusion (Clarendon, Oxford, 1975), 2 edn.
Carslaw, H. S. & Jaeger, J. C.
Conduction of Heat in Solids (Clarendon, Oxford, 1975), 2 edn.
Grebenkov, D. S., Metzler, R. & Oshanin, G. Strong defocusing of molecular reaction times results from an interplay ofgeometry and reaction control.
Commun. Chem. , 96 (2018). Godec, A. & Metzler, R. First passage time statistics for two-channel diffusion.
J. Phys. A: Math. Theor. , 084001(2017). Shoup, D. & Szabo, A. Role of diffusion in ligand binding to macromolecules and cell-bound receptors.
Biophys. J. , 33(1982). Lauffenburger, D. A. & Linderman, J.
Receptors: Models for Binding, Trafficking, and Signaling (Oxford University Press,Oxford, 1993).
Hodgkin, A. L. & Huxley, A. F. A quantitative description of membrane current and its application to conduction andexcitation in nerve.
J. Physiol. , 500–544 (1952).
Kerr, R. A. et al.
Fast Monte Carlo simulation methods for biological reaction-diffusion systems in solution and on surfaces.
SIAM J. on Sci. Comput. , 3126–3149 (2008). G. Yin and C. Zhu
Hybrid Switching Diffusions: Properties and Applications (Springer, New York, 2010).
G. Yin and C. Zhu, “Properties of solutions of stochastic differential equations with continuous-state-dependent switching”,J. Diff. Eq. , 2409-2439 (2010).
N. A. Baran, G. Yin, and C. Zhu, “Feynman-Kac formula for switching diffusions: connections of systems of partialdifferential equations and stochastic differential equations”, Adv. Diff. Eq. 2013:315 (2013).
D. S. Grebenkov, “Time-averaged mean square displacement for switching diffusion”, Phys. Rev. E , 032133 (2019). M. Freidlin,
Functional Integration and Partial Differential Equations , Annals of Mathematics Studies (Princeton UniversityPress, Princeton, New Jersey, 1985).
D. S. Grebenkov, “Partially Reflected Brownian Motion: A Stochastic Approach to Transport Phenomena”, in “Focus onProbability Theory”, Ed. L. R. Velle, pp. 135-169 (Nova Science Publishers, 2006).
D. S. Grebenkov, “Residence times and other functionals of reflected Brownian motion”, Phys. Rev. E , 041139 (2007). A. Singer, Z. Schuss, A. Osipov, and D. Holcman, “Partially Reflected Diffusion”, SIAM J. Appl. Math. , 844 (2008). D. S. Grebenkov, R. Metzler, and G. Oshanin, “Strong defocusing of molecular reaction times results from an interplay ofgeometry and reaction control”, Commun. Chem. , 96 (2018). J. Sheng, L. He, H. Zheng, L. Xue, F. Luo, W. Shin, T. Sun, T. Kuner, D. T. Yue, and L.-G. Wu, “Calcium-channel numbercritically influences synaptic strength and plasticity at the active zone”, Nature Neurosci. , 998-1006 (2012). upplemental Information SI1 Formal definition of the switching diffusion model
We reproduce here a formal definition of the ( M + ) -state switching diffusion model following Ref. . We consider atwo-component process ( XXX t , ν t ) , in which XXX t is the diffusion process in R , and ν t is the pure jump process with the states at { , , . . . , M } . When there is no boundary, the process is defined by a standard stochastic equation dXXX t = (cid:112) D ν t I dWWW t , ( XXX , ν ) = ( xxx , i ) , (S1)where WWW t is the standard Wiener process in R , I is the identity matrix, and D i is the diffusion coefficient at the state i . Thejump process is defined for any i (cid:54) = j by P { ν t + dt = j | ν t = i , XXX s , ν s , s ≤ t } = k i j dt + o ( dt ) , (S2)where k i j is the rate of transition from the state i to the state j . The propagator p ( xxx , i , t | xxx , i , ) is the probability density for theprocess to be in (the vicinity of) the point xxx in the state i at time t when stated from the point xxx in the state i . The propagatorsatisfies ( M + ) coupled forward Fokker-Planck equations ∂ t p ( xxx , i , t | xxx , i , ) = D i ∆ p ( xxx , i , t | xxx , i , ) + M ∑ j = (cid:2) k ji p ( xxx , j , t | xxx , i , ) − k i j p ( xxx , i , t | xxx , i , ) (cid:3) , (S3)subject to the initial condition p ( xxx , i , | xxx , i , ) = δ i , i δ ( xxx − xxx ) . Some properties of the propagator were discussed in (seealso the references therein).In turn, for a given smooth function f , the expectation of a functional f ( XXX t , ν t ) given that the process has started at xxx and i , u ( xxx , i , t ) = E { f ( XXX t , ν t ) | X = xxx , ν = i } , (S4)satisfies the ( M + ) coupled backward Fokker-Planck (or Kolmogorov) equations for each i , ∂ t u ( xxx , i , t ) = D i ∆ u ( xxx , i , t ) + M ∑ j = k i j (cid:0) u ( xxx , j , t ) − u ( xxx , i , t ) (cid:1) , (S5)subject to the initial condition u ( xxx , i , ) = f ( xxx , i ) (strictly speaking, this is a terminal condition but as the rates k i j do not dependon time, one can recast it as the initial condition).In the presence of a (partially) reflecting boundary, the diffusion component of the process is modified in a standard way(via the Skorokhod equation) , whereas the forward and backward Fokker-Planck equations need to be completed by theassociated boundary conditions, see . Setting f =
1, one can interpret u ( xxx , i , t ) as the probability for a particle started at xxx inthe state i to survive up to time t . SI2 General analytical solution
In this section, we present the derivation of the analytical solution for a general case with M buffers. Two particular cases(without buffer and with one buffer) will be detailed in Sections SI3 and SI4. SI2.1 Survival probabilities
We aim to find the survival probabilities S i ( t , xxx ) satisfying Eqs. (S5) with f = Ω = { xxx ∈ R : ρ < | xxx | < R } (S6)between two concentric spheres of radii ρ and R . The rotation symmetry of this domain implies that S i ( t , xxx ) depend only onthe radial coordinate r = | xxx | so that we can drop the dependence on angular coordinates and write S i ( t , r ) . Equations (S5) aresubject to the initial condition S i ( t = , r ) = , (S7)and have to be completed by boundary conditions (see the main text) ρ (cid:0) ∂ r S i ( t , r ) (cid:1) r = ρ = µ i S ( t , ρ ) , (S8a) (cid:0) ∂ r S i ( t , r ) (cid:1) r = R = , (S8b) t the inner and outer spheres, respectively, where µ = µ = k on πρ D N A , µ i = ( i = , . . . , M ) (S9)are dimensionless reactivities, with N A being the Avogadro number, and k on the on-rate binding constant.Introducing the Laplace-transformed survival probabilities (denoted by tilde),˜ S i ( p , r ) = ∞ (cid:90) dt e − pt S i ( t , r ) , (S10)one can rewrite the above equations as ( p + k i − D i ∆ ) ˜ S i − M ∑ j = k i j ˜ S j = ( ρ < r < R ) , (S11a) ∂ r ˜ S i = ( r = R ) , (S11b) µ i ˜ S i − ρ ∂ r ˜ S i = ( r = ρ ) , (S11c)where ∆ = ∂ r + ( / r ) ∂ r is the radial part of the Laplace operator, and k i = M ∑ j = k i j . (S12)As the rate k ii is undefined, we set k ii = S i ( p , r ) = a i + M ∑ j = b i j v ( δ j , r ) , (S13)where a i and b i j are unknown coefficients, and v ( δ , r ) = ρ r (cid:18) sinh ( δ ( R − r ) / ρ ) − ( + β ) δ cosh ( δ ( R − r ) / ρ ) (cid:19) , (S14)with β = ( R − ρ ) / ρ , (S15)and δ j are unknown factors. In fact, the function v ( δ , r ) is a linear combination of two independent solutions e δ r / r and e − δ r / r of the equation ∆ u − δ u =
0, and the chosen form (S14) ensures the Neumann boundary condition at the outer sphere for any δ : (cid:0) ∂ r v ( δ , r ) (cid:1) r = R = . (S16)Substituting Eq. (S13) into Eq. (S11), we get for i = , . . . , M ( p + k i ) (cid:18) a i + M ∑ j = b i j v ( δ j , r ) (cid:19) − D i ρ M ∑ j = b i j δ j v ( δ j , r ) − M ∑ (cid:96) = k i (cid:96) (cid:18) a (cid:96) + M ∑ j = b (cid:96) j v ( δ j , r ) (cid:19) = M + r ∈ ( ρ , R ) that implies M + i = , . . . , M : ( p + k i ) a i − M ∑ (cid:96) = k i (cid:96) a (cid:96) = (cid:0) p + k i − ( D i / ρ ) δ j (cid:1) b i j − M ∑ (cid:96) = k i (cid:96) b (cid:96) j = . (S19) he first set (S18) of M + a i is uncoupled from the rest and can be solved separately. Inverting theunderlying matrix, W = γ − k − k · · · − k M − k γ · · · − k γ · · · · · · · · · · · · · · · · · ·− k M · · · γ M (S20)(with γ i = p + k i ) and applying to the vector ( , , . . . , ) † , one gets a = (cid:18) + M ∑ i = k i p + k i (cid:19)(cid:18) p + k − M ∑ i = k i k i p + k i (cid:19) − = p . (S21)The other a i can also be found but their contribution will be canceled by µ i = i > b i j , in which δ j are some parameters. One can note that,for each j , there are M + j . In other words, we can decouple these equationsinto M + M + δ instead of δ j for one block. The equations in each block arehomogeneous, so that there is either none, or infinitely many solutions. For the existence of solutions, the determinant of theunderlying matrix in front of coefficients b i j should be zero. This matrix has precisely the same form as W in Eq. (S20), butwith γ i = p + k i − ( D i / ρ ) δ . The determinant of this matrix as a function of z = δ is the polynomial of degree ( M + ) H ( z ) = γ · · · γ M (cid:18) γ − M ∑ i = k i k i γ i (cid:19) . (S22)The M + z i , determine the unknown δ i : δ i = √ z i (here one can use either of two values ±√ z i , thefinal results remaining unchanged).For each j , the set (S19) of equations on b i j has infinitely many solutions. One can express b i j (for i = , . . . , M ) in terms of b j as b i j = k i p + k i − ( D i / ρ ) δ j b j . (S23)The remaining M + b j are determined by the ( M + ) boundary conditions at the inner sphere: (cid:0) µ i ˜ S i ( p , r ) − ρ ∂ r ˜ S i ( p , r ) (cid:1) r = ρ = ( i = , . . . , M ) (S24)that implies ( M + ) linear relations M ∑ j = b i j c i j = a i µ i ( i = , . . . , M ) , (S25)where c i j = (cid:18) ρ (cid:0) ∂ r v ( δ j , r ) (cid:1) r = ρ − µ i v ( δ j , ρ ) (cid:19) = sinh ( β δ j ) (cid:0) ( + β ) δ j − − µ i (cid:1) + δ j cosh ( β δ j ) (cid:0) β + µ i ( + β ) (cid:1) . (S26)Substituting Eqs. (S23) into these relations, one gets M + M + b j : M ∑ j = C i j b j = a i µ i ( i = , . . . , M ) , (S27)with C i j = c i j × ( i = ) , k i p + k i − ( D i / ρ ) δ j ( i > ) . (S28) nverting the matrix C , one obtains b j and thus fully determines ˜ S i ( p , r ) . Given that µ i = i > b j can formally bewritten as b j = µ f j ( p ) p f ( p ) , (S29)with f ( p ) = det ( C ) , f i j ( p ) = ( − ) i + j C i j , (S30)where C i j is the ( i , j ) minor of C , i.e., the determinant of the M × M matrix that results from deleting row i and column j of C .We get thus˜ S ( p , r ) = p (cid:18) + w ( p , r ) f ( p ) (cid:19) , (S31)where w ( p , r ) = µ M ∑ j = f j ( p ) v ( δ j , r ) . (S32)This is the exact analytic solution of the problem in the Laplace domain. In order to get the solution in time domain, one needsto compute the poles of ˜ S ( p , r ) which are given by zeros of the function f ( p ) considered in the complex plane ( p ∈ C ).The survival probability ˜ S ( p , r ) also determines the probability density of the first binding time, ˜ ψ ( p , r ) = − p ˜ S ( p , r ) ,from which˜ ψ ( p , r ) = − w ( p , r ) f ( p ) . (S33)In the general case k i > ψ ( , r ) = ψ ( t , r ) (we omit the related asymptotic analysis of the minors f i j ( p ) and of f ( p ) as p →
0; see the example for one buffer in Sec. SI4). As a consequence, p = S ( p , r ) , and S ( t , r ) vanishes in the long time limit. In turn, if k i = i , the calcium ions can be trapped forever by that buffer, and ˜ S ( t , r ) reaches a nonzero limit (the fraction of such trapped ions). In this specific case, ˜ ψ ( , r ) <
1, i.e., the normalization of ψ ( t , r ) does not hold. In practice, even if k i is very small, it is nonzero, and this pathologic situation does not occur. Note also that˜ S ( p , r ) determines the moments of the first binding times; in particular, the mean time is simply (cid:104) T (cid:105) = ∞ (cid:90) dt t ψ ( t , r ) = ∞ (cid:90) dt S ( t , r ) = ˜ S ( , r ) , (S34)where we integrated by parts and used that ψ ( t , r ) = − ∂ t S ( t , r ) and S ( ∞ , r ) = SI2.2 Occupancy probability
As discussed in the Methods Section, the probability density of the first binding times determines the occupancy probability P ( t , r ) in the Laplace domain as˜ P ( p , r ) = ˜ ψ ( p , r ) ˜ Q ( p ) , (S35)where ˜ Q ( p ) = (cid:18) p + k off ( − ˜ ψ ( p , ρ )) (cid:19) − . (S36)Substituting Eq. (S33) into this equation yields˜ Q ( p ) = (cid:18) p + k off + pk off M ∑ j = b j v ( δ j , ρ ) (cid:19) − . (S37) ext, substituting this expression and Eq. (S75) into Eq. (S35), we get explicitly˜ P ( p , r ) = − w ( p , r ) F ( p ) , (S38)with F ( p ) = ( p + k off ) f ( p ) + k off w ( p , ρ ) . (S39)The poles of ˜ P ( p , r ) are given by zeros of the function F ( p ) : F ( p n ) = ( n = , , . . . ) . (S40)One can invert the Laplace transform by using the residue theorem. In particular, if the poles are simple, one gets P ( t , r ) = ∞ ∑ n = ¯ b n w ( p n , r ) exp ( p n t ) , (S41)where ¯ b n = − p → p n ∂ p F ( p ) , (S42)in which the derivative can be computed by using ∂∂ p v ( δ , r ) = − (cid:18) cosh ( β δ ) + β ( + β ) δ sinh ( β δ ) (cid:19) ∂ δ∂ p (S43)and ∂∂ p c i j = ∂ δ j ∂ p (cid:26) cosh ( β δ j ) (cid:0) µ i + β ( + β ) δ j (cid:1) + δ j sinh ( β δ j ) (cid:0) µ i β ( + β ) + ( β + β + ) (cid:1)(cid:27) . (S44)It can checked that p = p n <
0. As a consequence, as t → ∞ , the probability P ( t , r ) approaches a stationary value P ∞ , which is independent of the starting point r and given by the residue at p = P ( t , r ) = P ∞ + ∞ ∑ n = exp ( p n t ) M ∑ j = b jn v ( δ j ( p n ) , r ) , (S45)where b jn = µ f j ( p n ) ¯ b n . (S46)Setting α n = ρ (cid:112) − p n / D , α ( j ) n = − i δ j ( p n ) , b ( j ) n = ib jn , one can rewrite the occupancy probability in a more conventional form: P ( t , r ) = P ∞ + ∞ ∑ n = exp (cid:0) − α n D t / ρ (cid:1) M ∑ j = b ( j ) n u ( α ( j ) n , r ) , (S47)where u ( δ , r ) = ρ r (cid:18) sin ( δ ( R − r ) / ρ ) − ( + β ) δ cos ( δ ( R − r ) / ρ ) (cid:19) . (S48)We note that the functions ˜ S ( p , r ) and ˜ P ( p , r ) involve complicated combinations of roots (e.g., square roots, see below)emerging from the zeros of the polynomial H ( z ) in Eq. (S22). As a consequence, the use of the residue theorem for evaluating he inverse Laplace transform of these functions is not straightforward as one needs to introduce cuts in the complex plane toproperly deal with such multivariate functions. In addition, the application of Eq. (S41) relies on the assumption of simple poles.In this paper, we do not provide rigorous mathematical analysis of both statements. In turn, we have checked the correctnessand the accuracy of the derived formulas in time domain by comparison with the numerical inversion of the Laplace transform(not shown).In summary, the analytic solution requires three numerical steps: (i) computation of δ j as the zeros of Eq. (S22); (ii)inversion of the matrix C in Eq. (S28), from which f j ( p ) , f ( p ) and thus b j are found; and (iii) finding the zeros of f ( p ) (for getting S ( t , r ) ) or of F ( p ) (for getting P ( p , r ) ) for the inversion of the Laplace transform. We emphasize that δ j and b j depend on p , i.e. one needs to perform the first two steps for all values of p at which ˜ S ( p , r ) has to be found. In practice, thenumber of buffers, M , is not large so that these numerical steps can be done very rapidly and with any accuracy. We will discussthe cases M = M = SI2.3 Steady-state limit P ∞ As time t goes to infinity, the occupancy probability P ( t , r ) from Eq. (S47) approaches the steady-state limit P ∞ , whichis determined by the residue of ˜ P ( p , r ) at the pole p =
0. Even though all the formulas determining ˜ P ( p , r ) are given, thecomputation of this residue is technically involved, see the related analysis below for the particular cases of no buffer and onebuffer. For this reason, we prefer to rely here on qualitative physical arguments that allow us to get the exact form of P ∞ withouttedious computations.In the steady-state, the system reaches an equilibrium between the free state, the buffer-bound states, and the sensor-boundstate. Moreover, as the binding/unbinding kinetics on the sensor occurs only through the free state, one can separate thekinetics with the sensor and the kinetics with the buffers. The equilibrium kinetics with the sensor can be understood asa two-state switching model, governed by two exchange rates: k off describes the transition from the sensor-bound stateto the free state, whereas an effective rate k , s = k on c characterizes the opposite transition, where c is the equilibrated(homogeneous) concentration of calcium ions. If p is the equilibrium fraction of calcium ions in the free state, then theconventional concentration (in M = mol/liter) reads c = p / ( N A V ) , where V = π ( R − ρ ) / P ∞ = k , s / ( k , s + k off ) or, equivalently, P ∞ = + k off N A Vk on p . (S49)The fraction p of calcium ions in the free state can be determined from the equilibrium between the free state andbuffer-bound states. For this purpose, we only consider the dynamics of the ( M + ) -state switching model governed by thetransition matrix W from Eq. (S20) with γ i = k i (i.e., at p = W † that corresponds to the eigenvalue 0: ( , k / k , k / k , . . . , k M / k M ) † . After normalization to 1, the probability offinding the calcium ion in the free state (i.e., the fraction of calcium ions in this state) is p = (cid:18) + M ∑ j = k j k j (cid:19) − . (S50)We get therefore P ∞ = (cid:18) + k off N A Vk on (cid:18) + M ∑ j = k j k j (cid:19)(cid:19) − . (S51)The same expression for P ∞ is retrieved for cases M = M = SI3 No buffer case
The survival probability for no buffer case is well known (see and references therein). For illustrative purposes, weretrieve this survival probability from our general approach. This step is also needed for finding the occupancy probability P ( t , r ) .When there is no buffer ( M = H = γ = p − ( D / ρ ) δ =
0, from which δ = ρ (cid:112) p / D .The matrix C consists of a single element C = c , from which f ( p ) = f ( p ) = c , and thus b = µ / ( pc ) . TheLaplace-transformed survival probability becomes then˜ S ( p , r ) = p + µ v ( δ , r ) p c ( p ) , (S52) here c ( p ) is given by Eq. (S26). Setting δ = i ˆ α , one can rewrite the equation f ( p ) = S ( p , r ) in atrigonometric formsin ( ˆ αβ ) = ( β + ( + β ) µ ) ˆ α + µ + ( + β ) ˆ α cos ( ˆ αβ ) , (S53)which has infinitely many nonnegative solutions denoted as ˆ α n , enumerated by n = , , , . . . (we use hat symbol here todistinguish the quantities determining S ( t , r ) from similar quantities determining P ( t , r ) below). The poles are ˆ p n = − D ˆ α n / ρ .Note that the pole corresponding to ˆ α = − / p that precisely compensates the term a = / p , andthus it will be excluded. The inverse Laplace transform is then obtained by the residue theorem: S ( t , r ) = ∞ ∑ n = ˆ b n u ( ˆ α n , r ) exp (cid:0) − ˆ α n D t / ρ (cid:1) , (S54)where ˆ b n = µ ˆ α n (cid:18) cos ( ˆ α n β ) (cid:2) µ − β ( β + ) ˆ α n (cid:3) − ˆ α n sin ( ˆ α n β ) (cid:2) β ( β + ) µ + ( β + β + ) (cid:3)(cid:19) − (S55)and u ( δ , r ) is given by Eq. (S48). The derivative with respect to t yields the probability density of the first-binding time ψ ( t , r ) = D ρ ∞ ∑ n = ˆ α n ˆ b n u ( ˆ α n , r ) exp (cid:0) − ˆ α n D t / ρ (cid:1) . (S56) Occupancy probability P ( t , r ) The computation of the probability P ( t , r ) follows the same lines. Setting δ = i α in Eq. (S39), one gets F = iD ρ (cid:26) sin ( αβ ) (cid:18) ( α − λ )(( + β ) α + + µ ) + λ µ (cid:19) − α cos ( αβ ) (cid:18) ( α − λ )( β + µ ( + β )) + λ µ ( + β ) (cid:19)(cid:27) , (S57)from which the equation on α readssin ( αβ ) = (cid:2) α ( β + µ ( + β )) − λ β (cid:3) α cos ( αβ ) α ( + β ) + α ( + µ − λ ( + β )) − λ , (S58)where λ = k off ρ / D . This equation has infinitely many positive zeros that we denote as α n , with n = , , . . . (the zero α = p n = − D α n / ρ . Since w ( p n , r ) = i µ u α n ( r ) with u ( δ , r ) given by Eq. (S48), we obtain by the residue theorem P ( t , r ) = + λ ( β + ) − µ + ∞ ∑ n = b n u ( α n , r ) e − α n D t / ρ , (S59)where the first term comes from the residue at p =
0, and b n = − i µ lim p → p n (cid:0) ∂ p F ( p ) (cid:1) = i µ ρ D α n (cid:0) ∂ α F ( α ) (cid:1) α = α n . (S60)Recalling the definition of dimensionless parameters λ , µ and β , one easily checks that the first term in Eq. (S59) coincideswith the steady-state limit P ∞ in Eq. (S51).Taking the derivative of Eq. (S57) with respect to α , one gets an explicit formula for b n : b n = µ sin ( α n β ) (cid:0) α n w + w (cid:1) + α n cos ( α n β ) (cid:0) α n w + w (cid:1) , (S61)with w = ( + β ) + β ( β + µ ( + β )) , w = ( + µ − λ ( + β )) − λ β , w = β ( + β ) , w = β ( + µ − λ ( + β )) − ( β + µ ( + β )) . (S62) imiting cases In the limit k off = λ = Q ( p ) = / p according to Eq. (S36). In this case,˜ P ( p , r ) = ˜ ψ ( p , r ) p = − ˜ S ( p , r ) p , and thus P ( t , r ) = − S ( t , r ) , as expected. One can also check that the solutions α n coincide with ˆ α n .In turn, in the limit of perfectly adsorbing sensor (i.e., with infinitely fast binding kinetics: k on = µ = ∞ ), Eq. (S53) isreduced tosin ( ˆ α n β ) = ( + β ) ˆ α n cos ( ˆ α n β ) , (S63)and the survival probability becomes S ( t , r ) = ρ r ∞ ∑ n = exp (cid:0) − ˆ α n D t / ρ (cid:1) × sin (cid:0) ˆ α n R − r ρ (cid:1) − ( + β ) ˆ α n cos (cid:0) ˆ α n R − r ρ (cid:1) ˆ α n (cid:0) cos ( ˆ α n β ) − β ( + β ) ˆ α n sin ( ˆ α n β ) (cid:1) . (S64)The probability density ψ ( t , r ) = − ∂ t S ( t , r ) is obtained by taking the derivative with respect to t . Note that in this limit,the unbinding events are effectively suppressed as a particle that unbinds from such a sensor immediately re-binds. As aconsequence, one gets again P ( t , r ) = − S ( t , r ) . Other results
Mean first-binding time
The mean first-binding time reads (cid:104) T (cid:105) r = ˜ S ( , r ) = r ρ ( R − ρ ) / µ + ρ R ( r − ρ ) − r ρ ( r − ρ ) r ρ D , (S65)while the mean excursion time (at r = ρ ) is (cid:104) T (cid:105) ρ = R − ρ D ρ µ = ρ V µ D A . (S66)where V is the volume of the domain Ω and A is the area of the sensor. As a consequence, the mean first-binding time, whichis essentially proportional to the volume of the active zone, is a useless characteristics in this situation. In turn, the mode (i.e.,the position of the density maximum, i.e., the most probable value) is representative. Asymptotic analysis of the smallest eigenvalue
The long-time behavior of ψ ( t , r ) , Q ( t ) , and the probability P ( t , r ) , is determined by the smallest absolute value of thepole | p | of the underlying Laplace-transformed quantity. Let us first consider the density ψ ( t , r ) , for which the smallest | p | is determined by ˆ α . Denoting x = ˆ α β and assuming that x →
0, one can use the Taylor expansion of Eq. (S53) to determinethe asymptotic behavior of ˆ α for large β . In the lowest order in 1 / β , we getˆ α (cid:39) µ ( + µ ) β (cid:39) µρ ( + µ ) R . (S67)According to Eq. (S54), the above relation determines the slowest decay rate of the survival probability, ρ / ( D ˆ α ) , which isclose to the mean time (S66) when ρ (cid:28) R . Short-time asymptotic behavior
The short-time asymptotic behavior corresponds to the limit p → ∞ . In this limit, Eq. (S38) becomes in the leading order in1 / p : ˜ P ( p , r ) (cid:39) µ √ D exp (cid:0) − ( r − ρ ) (cid:112) p / D (cid:1) r p / , (S68)from which the short-time asymptotic behavior follows for r > ρ P ( t , r ) (cid:39) ( D t ) / µ √ π r ( r − ρ ) exp (cid:0) − ( r − ρ ) / ( D t ) (cid:1) . (S69)This asymptotic behavior is applicable at times as short as t (cid:28) ( r − ρ ) / ( D ) . In turn, for r = ρ , Eq. (S68) yields P ( t , ρ ) (cid:39) √ D µ √ πρ t / ( t → ) . (S70) I4 One buffer case
For a single buffer ( M = H = (cid:0) p + k − ( D / ρ ) z (cid:1)(cid:0) p + k − ( D / ρ ) z (cid:1) − k k , (S71)and its two zeros determine δ and δ : δ = ρ D D (cid:18) D ( p + k ) + D ( p + k ) − (cid:113) ( D ( p + k ) − D ( p + k )) + D D k k (cid:19) , (S72) δ = ρ D D (cid:18) D ( p + k ) + D ( p + k ) + (cid:113) ( D ( p + k ) − D ( p + k )) + D D k k (cid:19) . (S73)Getting f ( p ) = C C − C C , f ( p ) = C , f ( p ) = − C (S74)from the 2 × C , one obtains the coefficients b j b = C µ p ( C C − C C ) , b = − C µ p ( C C − C C ) , (S75)where the elements C i j are given explicitly by Eq. (S28). We obtain thus˜ S ( p , r ) = p + b v ( δ , r ) + b v ( δ , r ) . (S76)In order to invert the Laplace transform, one needs to determine the poles of ˜ S ( p , r ) that are given by the zeros ˆ p n of thefunction f ( p ) . There are infinitely many zeros and they are nonpositive: ˆ p n ≤
0. To compute the residues, one needs thederivative of f ( p ) with respect to p , which can be evaluated by using Eq. (S44) and ∂ δ j ∂ p = ρ δ j p + k + k − ( D + D ) δ j / ρ D ( p + k ) + D ( p + k ) − D D δ j / ρ . (S77)Finally, we proceed to check that the two zeros of f ( p ) , p = p = − ( k + k ) , are not the poles of ˜ S ( p , r ) , and thusexcluded from the analysis.(i) In the limit p →
0, we get δ (cid:39) ρ k + k D k + D k p + O ( p ) , (S78a) δ (cid:39) ρ ( D k + D k ) D D + O ( p ) , (S78b) v ( δ , ρ ) (cid:39) − δ − δ ( β / + β / ) + O ( δ ) , (S78c) C (cid:39) µδ + δ (cid:18) β + β + β + µ (cid:18) β + β (cid:19)(cid:19) + O ( δ ) , C (cid:39) δ (cid:18) β + β + β (cid:19) + O ( δ ) , (S78d)whereas C and C approach constants. We get thus b (cid:39) µ p C µδ C − O ( δ ) = δ p + O ( p − / ) , b (cid:39) − µ p O ( δ ) µδ C − O ( δ ) = O ( ) . (S79) ince v ( δ , r ) = − δ + O ( δ ) , the singularities from a = / p and b v ( δ , r ) cancel each other so that p = S ( p , r ) .(ii) Setting p = − ( k + k ) + ε , one has δ = − D k + D k D D + O ( ε ) , δ = k + k D k + D k ε + O ( ε ) . (S80)As a consequence, we get C (cid:39) µδ = O ( ε / ) , C (cid:39) δ = O ( ε / ) , (S81)whereas a , C and C approach constants. We obtain then b = µ p C C C − C C = O ( ) , b = − µ p C C C − C C (cid:39) µ C p (cid:39) δ p = O ( ε − / ) . (S82)Since v ( δ , r ) (cid:39) − δ , the term b v ( δ , r ) has no singularity so that p = − ( k + k ) is not a pole of ˜ S ( p , r ) .We conclude that S ( t , r ) = ∞ ∑ n = (cid:18) ˆ b n v (cid:0) δ ( ˆ p n ) , r (cid:1) + ˆ b n v (cid:0) δ ( ˆ p n ) , r (cid:1)(cid:19) exp ( ˆ p n t ) , (S83)where v ( δ , r ) is given by Eq. (S14), andˆ b n = µ C ( ˆ p n ) ˆ p n f (cid:48) ( ˆ p n ) , ˆ b n = − µ C ( ˆ p n ) ˆ p n f (cid:48) ( ˆ p n ) . (S84)The derivative with respect to t yields ψ ( t , r ) = ∞ ∑ n = (cid:18) ˆ b n v (cid:0) δ ( ˆ p n ) , r (cid:1) + ˆ b n v (cid:0) δ ( ˆ p n ) , r (cid:1)(cid:19) × | ˆ p n | exp ( ˆ p n t ) . (S85) Probability P ( t , r ) Similarly, the inversion of ˜ P ( p , r ) involves the zeros p n of F ( p ) from Eq. (S39) that can be written explicitly as: F ( p ) = ( p + k off ) f ( p ) + k off µ (cid:0) C v ( δ , ρ ) − C v ( δ , ρ ) (cid:1) , (S86)with f ( p ) from Eq. (S74). As previously, one can show that the zero p = − ( k + k ) is not a pole of ˜ P ( p , r ) . In turn, p = p → F ( p ) (cid:39) δ p (cid:26) µ C ( ) + λ D ( k + k ) D k + D k (cid:18) β + β + β (cid:19) × (cid:0) C ( ) − C ( ) − µ v ( δ ( ) , ρ ) (cid:1)(cid:27) , where C ( ) , C ( ) and δ ( ) denote the values of these functions evaluated at p =
0. In turn, w ( p , r ) = µ (cid:0) C v ( δ , r ) − C v ( δ , r ) (cid:1) (cid:39) − µ C ( ) δ + O ( δ ) , so that the residue at p = P ∞ = (cid:26) + λ D ( k + k ) D k + D k ( + β ) − × C ( ) − C ( ) − µ v ( δ ( ) , ρ ) µ C ( ) (cid:27) − , (S87)which is independent of the starting point r . After simplifications, we have P ∞ = (cid:18) + λ ( + β ) − µ ( + k / k ) (cid:19) − . (S88) his expression coincides with Eq. (S51).We get thus P ( t , r ) = P ∞ + ∞ ∑ n = (cid:18) b n v ( δ ( p n ) , r ) + b n v ( δ ( p n ) , r ) (cid:19) e p n t , (S89)with b n = − µ C ( p n ) F (cid:48) ( p n ) , b n = µ C ( p n ) F (cid:48) ( p n ) . (S90) One fixed buffer
For the fixed buffer ( D → δ = ρ D (cid:18) p + k − k k p + k (cid:19) , δ → ∞ . (S91)As a consequence, one needs to treat this case separately to avoid diverging terms.The last relation in Eqs. (S91) implies that c i (cid:39) ( + β ) δ e βδ → ∞ ( i = , ) . In addition, we have C = c , C = c k p + k , C = c k p + k , so that in the limit D →
0, we get b = µ p c , b = , (S92)given that C C = c c k p + k − ( D / ρ ) δ = c c p + k − ( D / ρ ) δ k → . (S93)We conclude that˜ S ( p , r ) = p + µ v ( δ , r ) p c , (S94)i.e., we retrieved the solution (S52) for the case without buffer, in which δ = ρ (cid:112) p / D is replaced by δ = ρ (cid:112) p (cid:48) / D , where p (cid:48) = p + k − k k p + k . (S95)The fixed buffer is expected to slow down the arrival onto the sensor because of binding calcium ions and thus stopping theirdiffusion. In particular, one can notice this effect in an increase of the mean first-binding time to the sensor, given by ˜ S ( , r ) .Noting that p (cid:48) = p =
0, one finds that the mean first-binding time without buffer, (cid:104) T nb (cid:105) , is multiplied bythe factor ( + k / k ) in the presence of a fixed buffer: (cid:104) T (cid:105) = ˜ S ( , r ) = (cid:104) T nb (cid:105) (cid:18) + k k (cid:19) . (S96)The relation to the former solution without buffer allows one to easily invert the Laplace transform. In fact, the formerpoles of ˜ S ( p , r ) were ˆ p n = − D ˆ α n / ρ . Inverting the relation (S95), one can see that each former pole ˆ p n splits in two newpoles ˆ p n , = − λ n , and ˆ p n , = − λ n , , with λ n , = σ n − (cid:112) σ n − k D ˆ α n / ρ , (S97a) λ n , = σ n + (cid:112) σ n − k D ˆ α n / ρ , (S97b) ith σ n = D ˆ α n / ρ + k + k . As a consequence, the inverse Laplace transform of Eq. (S94) becomes S ( t , r ) = ∞ ∑ n = ˆ b n u ( ˆ α n , r ) (cid:18) c n , e − λ n , t + c n , e − λ n , t (cid:19) , (S98)where u ( δ , r ) is given by Eq. (S48), the coefficients ˆ b n are given by Eq. (S55), and the weights c n , = D ˆ α n / ρ λ n , + k k ( λ n , − k ) , (S99a) c n , = D ˆ α n / ρ λ n , + k k ( λ n , − k ) , (S99b)appear from the change of variables: d p (cid:48) / d p = + k k / ( p + k ) , see Eq. (S95), and from the factor 1 / p in the secondterm of Eq. (S94). Note that if k = k =
0, one has λ n , = λ n , = D ˆ α n / ρ , and one retrieves Eq. (S54).Note also that λ n , → k and λ n , → D ˆ α n / ρ as n → ∞ and thus c n , → c n , →
1. In other words, the exchangekinetics does not affect the high-frequency eigenmodes.Substituting Eq. (S98) into (S35, S36), we get˜ P ( p , r ) = − µ v ( δ , r )( p + k off ) c + k off µ v ( δ , ρ ) , (S100)so that one needs to find zeros of the denominator of this expression. As in the former case for ˜ S ( p , r ) , one can expect twosequences of zeros: p n , → − k and p n , → − ∞ . In fact, when p → − k + p (cid:48) from Eq. (S95) diverges to − ∞ , so thatthere are infinitely many zeros accumulating towards − k . This accumulation requires a more subtle numerical procedure tocalculate zeros. SI5 Models used for Monte Carlo simulations
SI5.1 Calcium channel model
We describe a VGCC by a 3-state Hodgkin and Huxley gating model so that the calcium release was modeled accordingto: C α ( V ( t )) (cid:29) β ( V ( t )) C α ( V ( t )) (cid:29) β ( V ( t )) O k ( t ) → Ca + , (S101)with two closed states C , C and one open state O of the VGCC. Here α ( V ( t )) and β ( V ( t )) are voltage dependent rates,computed as α ( V ( t )) = exp ( V ( t ) / . ) , β ( V ( t )) = .
14 exp ( − V ( t ) / ) , (S102)for a given AP waveform V ( t ) in mV. The dynamics starts from the close state C . The parameters in these rates were adjustedsuch that the resulting single channel open probability, current duration, and peak match experimentally observed quantities .The calcium ions are released from the open channel with the rate: k ( t ) = g e ( V ( t ) − V rev ) , (S103)where g = . , e is the elementary charge, and V rev = −
45 mV is thereversal potential . SI5.2 Sensor kinetics model
The synaptic vesicle release is often described by the following 5-state sensor kinetics model: V k on (cid:29) k off b V k on (cid:29) k off b V k on (cid:29) k off b V k on (cid:29) k off b V k on (cid:29) k off b V γ → F , (S104)where V i denote the binding states of the sensor (i.e., the sensor with i calcium ions bound, and V meaning the unbound state),and F is the fused state of the vesicle. The parameters are: k on =
127 mM − ms − , k off = . − , b = . γ = − .Except for the deterministic simulations shown in Fig. 1, we consider only first step of this model, namely, between V and V ,using the k on constant 5 times higher than the above value 127 mM − ms − . I6 Supplementary Figures -8 -6 -4 -2 A B C
Analyt.MCs P N ( t, r) P N ( t, r) . % ( F W H M ) . ( M AE ) . % .
004 31 . % . . % . -8 -6 -4 -2 . % . . % . -8 -6 -4 -2 -8 -6 -4 -2 -8 -6 -4 -2 s -8 -6 -4 -2 ss ss Supplementary Figure S1. P ( t , r ) with slow unbinding kinetics ( k on = .
157 mM − ms − ) for 50 ( A ), 100 ( B ) and 200( C ) simultaneously released ions for CD of 15 nm (top row) and 45 nm (bottom row). Black and green lines show respectivelyanalytical and MC results. The black and blue inset text on each plot represent FWHM error and MAE between analytical andMC results correspondingly. -8 -6 -4 -2 A B C
Analyt.MCs P N ( t, r) P N ( t, r) . % ( F W H M ) . ( M AE ) . % ~0 1 . % . . % ~0 -8 -6 -4 -2 . % . . % . -8 -6 -4 -2 -8 -6 -4 -2 -8 -6 -4 -2 s -8 -6 -4 -2 ss ss Supplementary Figure S2. P ( t , r ) with fast unbinding kinetics ( k on = − ms − ) for 50 ( A ), 100 ( B ) and 200 ( C )simultaneously released ions for CD of 15 nm (top row) and 45 nm (bottom row). Black and green lines show respectivelyanalytical and MC results. The black and blue inset text on each plot represent FWHM error and MAE between analytical andMC results correspondingly.)simultaneously released ions for CD of 15 nm (top row) and 45 nm (bottom row). Black and green lines show respectivelyanalytical and MC results. The black and blue inset text on each plot represent FWHM error and MAE between analytical andMC results correspondingly.