Hybrid Markov-mass action law model for cell activation by rare binding events: Application to calcium induced vesicular release at neuronal synapses

Binding of molecules, ions or proteins to small target sites is a generic step of cell activation. This process relies on rare stochastic events where a particle located in a large bulk has to find small and often hidden targets. We present here a hybrid discrete-continuum model that takes into account a stochastic regime governed by rare events and a continuous regime in the bulk. The rare discrete binding events are modeled by a Markov chain for the encounter of small targets by few Brownian particles, for which the arrival time is Poissonian. The large ensemble of particles is described by mass action laws. We use this novel model to predict the time distribution of vesicular release at neuronal synapses. Vesicular release is triggered by the binding of few calcium ions that can originate either from the synaptic bulk or from the entry through calcium channels. We report here that the distribution of release time is bimodal although it is triggered by a single fast action potential. While the first peak follows a stimulation, the second corresponds to the random arrival over much longer time of ions located in the synaptic terminal to small binding vesicular targets. To conclude, the present multiscale stochastic modeling approach allows studying cellular events based on integrating discrete molecular events over several time scales.

Cellular activation is described by the binding of few diffusing molecules to specific small and hidden targets. What is the time scale of a cell response, triggered by rare molecular events? Such a process requires exploring microdomains at a nanometer precision. Examples are vesicular fusion triggered by the binding of calcium ions, detection of a morphogen concentration by a growth cone, molecular events underlying the induction of synaptic plasticity at neuronal synapses or activating cellular check points where molecular events are transformed into a cellular decision.
Traditional chemical kinetics 1 are based on mass-action laws or reaction-diffusion equations, but this is a deficient description of stochastic chemical reactions in micro-domains, where only a small number of molecules are involved [2][3][4] . Several models were developed to describe stochastic reactions, such as Markov reaction-diffusion equations based on the joint probability distribution for the concentration of bound and unbound reactants, leading to coupled systems of reaction-diffusion 5 . Another direction involves the Narrow Escape Theory 5 to coarse-grain the binding and unbinding reactions to the time scale of diffusing reactant in and out of a binding site. This approach reduces the classical reaction-diffusion approach to a Markovian jump process description of the stochastic dynamics for the binding and unbinding of molecules. The Markovian approximation 1,6 can be used for example to obtain predictions for the rate of molecular dynamics underlying spindle assembly checkpoint during cell division 7 or the probability that a messenger RNA escapes degradation through binding a certain number of microRNAs. The Markovian approach can also be used to compute the mean time the number of bound molecules reaches a threshold, called the mean time to threshold theory (MTT), for characterizing the stability of activation 6 .
However, none of the approaches described above have been used to model the interaction between a large bath of particles modeled by a continuous dynamic and a stochastic dynamic induced by rare events where few particles bind to small targets to trigger a response.
The goal of this report is to introduce a novel hybrid model and to provide an application to synaptic physiology for computing the probability of vesicular release following calcium influx. Calcium ions induce vesicular release when few of them accumulate underneath a vesicle. This accumulation can follow the direct influx at a short time scale or can be due to ions originating from the pre-synaptic terminal with a longer time scale. In all cases, triggering a vesicle release is a rare molecular event described as the accumulation of few ions at a single protein location 8 . However, due to instrumental limitations, it is not yet possible to measure calcium dynamics at a nanometer level at neuronal synapses. Yet, synchronous and asynchronous synaptic release might depend on calcium dynamics and we provide here several predictions using the present hybrid stochastic model of chemical reactions, that we compare to Gillespie simulations. The proposed mechanism of asynchronous vesicular release that we found here clarifies the role of calcium ions during short-term synaptic plasticity 9 .

Results
Hybrid Markov-mass action model. Reactant molecules R are modeled as Brownian particles, entering a domain Ω through localized point sources with an inward flux J(t) at each single source. The entering Brownian particles can reach any target site to activate them. After particles enter, they can either hit a small target or reach the bulk, where they are lost in the undifferentiated state consisting in many particles. In the bulk, the reactant particles R can either interact with a substrate S, leave the bulk with a rate k es or bind a target site with rate k T , as shown in Fig. 1. The arrival of T particles at a single target is the event that we are interested in, as it models a cellular activation from a molecular event. As we shall see here, the Brownian reactant particles and the small targets can be described by a coarse-grained model, where the mass action description (for the continuum limit) is coupled to a finite state Markov chain (Fig. 1).
For the application to vesicular release described in the second part, the bulk represents the pre-synaptic terminal ( Fig. 2A), and N T the number of small target sites (red ribbon in Fig. 2B) hidden underneath a ball representing a vesicle.
The flux J(t) of particles enters the domain at a source point x. The entering particles can reach a small target with a splitting probability p s (x) or go to the bulk with probability 1 − p s . A particle entered at point x that reaches one of the N T target sites, will reach specifically site i with probability q(x, i), for i = 1..N T (see the Methods section for the description). Thus the flux of ions arriving at a target site i is the sum of the conditional fluxes over the ensemble of point source locations : We recall that the splitting probability p s (x) for a Brownian particle to reach the small ribbon below a ball ( Fig. 2B) before entering the bulk, when balls are distributed on a square lattice, has been estimated previously using conformal mapping methods and asymptotic approximations 10 : where r is the size of the ball, H is half the distance between two targets, A = 9.8, ρ(x) is the distance between the point source and the closest target site, and ε is the height of the ribbon representing the target site. This probability accounts for the particular geometry of the target and depends on the relative distance between the targets and the source points 10 . and are then divided between the target sites, with probability p s (Eq. 2) and the bulk with probability 1 − p s . In the bulk, particles R react with the substrate S, modeled by the classical mass action equations, with forward rate k 0 , and backward rate k −1 . The escape rate from the bulk is k es . Particles located in the bulk can also bind to the targets with constant rate k T . The overall scheme represents the coupled mass action equations (in the bulk) with a Markov chain and the ensemble models the arrival of T particles to a single target. So far, we discussed the direct arrival of a particle from a point source. We now consider the arrival of a Brownian particle coming from the bulk to one of the small targets. This arrival to a small ribbon (red in Fig. 2B) is Poissonian with a rate (computed in ref. 10 T where D is the particle diffusion coefficient. We shall now present the set of equations that describe the arrival of a few number of particles to the target either coming directly from the point source or from the bulk. For each target site i and for a given distribution P of source points, the probability Pr i {k, t} that k particles are bound at time t is computed from a Markov model: once a particle binds to a target, it cannot unbind, and thus the transition probability from the states k − 1 to k (bound particles), between t and t + Δ t, is due to the fluxes of particles from the point sources and from the bulk: where λ dis i represents the arrival rate of particles to the target: in the right-end side of the rate equation 5 represents the cumulative flux of particles reaching target site i summed over all point sources. It represents the amount of particles reaching target i directly after their entry through a source point. The second term k T N f (t) represents the fraction of particles originating from the bulk that will bind to the target sites. Here N f (t) is the number of free particles at time t in the bulk, and k T is the arrival rate of particles located in the bulk. The rate k T is computed using the Mean First Passage Time theory 10 , and when the target size is small compared to the size of the bulk, the binding events are rare and thus can be modeled as Poissonian rate which is the reciprocal of the mean first time for a particle to reach the target 5,10 . These two terms account for the targets geometry, as well as their relative positions with respect to the source points .
We shall now described the model of activation. A target is considered to be activated when there are exactly T bound particles. We obtain for each target 1 ≤ i ≤ N T , the Markov chain for the probabilities  (A) A Brownian particle (green) entering the presynaptic terminal either reaches a target (blue) located near the source points (orange) or reaches the bulk. In the bulk, the particle can bind to a substrate or leave the domain. (B) Schematic description of the arrival of a particle to a target site. The Brownian particles enter through source points. They can reach the small ribbon (red) located underneath a target (blue). The arrival of T particles to the ribbon is needed to activate the target.
Scientific RepoRts | 6:35506 | DOI: 10.1038/srep35506 and the normalization condition is We can now introduce the time τ T i to a threshold T for target i, when the distribution of point source is . This time corresponds to the first time that there are exactly T particles at target i. It is also the last passage time of a particle when T − 1 particles are already at the binding sites. The distribution of time τ T i is called the activation time for the target i and is given by T i T i hence the probability density function τ f T i of the activation time is defined by To conclude this section, we have presented here a Markov chain that counts exactly T binding events. We note that the particles originate from two sources: either from the bulk or a transient influx. We shall introduce in the next section a new set of equations to couple the Markov chain and the density of free particles N f (t) that can further interact with a varying population of substrates S.

The mass action law equations for reactant particles interacting with a substrate in the bulk.
The reactant particles in the bulk can interact with a substrate according to the chemical reaction where k −1 (resp. k 0 ) is the backward (resp. forward) rate. We determine the mass-action equations for the number of free N f (t), and bound N b (t) particles. The total number of substrate sites S tot is fixed, and at time t, the number of available sites is S tot − N b (t). The particles can escape the domain Ω at a Poissonnian rate k es , and bind to a target site with a rate k T . We note that the probability that target i is free, i.e. not activated, is given by In summary, N f and N b satisfy the mass action equations: The first two terms in the first equation represent the classical unbinding and binding of molecules to buffers. The third term represents the fraction of particles directly entering the bulk from the initial influx. The fourth term corresponds to particles that leave the bulk or bind to a free target (we assume that the escape rate k es is fast compared to the binding rates). Finally the last term accounts for the release into the bulk of the T particles bound to a target, following the target activation. Indeed, after activation, the T bound particles are released into the bulk, leading to an increase of the free particles N f by a jump event of T particles (Eq. 10). The second equation is the classical mass action law for the number of bound buffers. Interestingly, the system of equations 6-14 can be decoupled and resolved. We use a direct integration between t 0 and t for the probability p t ( )  We shall now provide an application to neurobiology about calcium signaling at single neuronal synapses. There are indeed 10 11 neurons in a human brain, and each has a mean of 10 3 synapses. Synapses are thought to code part of the memory, mainly by modulating the neuronal response through time variations in the release of vesicles. We will study now some aspects of the stochastic release of a vesicle from the accumulation of T = 5 calcium ions to key proteins. This example illustrates the theory presented above in the context of synaptic activation from rare molecular binding events.
Dynamics of vesicular release at neuronal synapse. The pre-synaptic terminal of a neuronal cell contains a large amount of vesicles that can be randomly released following an action potential, characterized by a release probability 11 . However, the release mechanism is not well understood and computing the release probability remains a challenge [12][13][14][15] . The mechanism is described as follows: after an action potential invades the pre-synaptic terminal, voltage-gated calcium channels (VGCCs) open, leading to a calcium influx into the pre-synaptic domain 8,16,17 . When several diffusing calcium ions have succeeded to find small molecular targets, which form a complex of molecular machinery (such as synaptotagmin and any other key molecules) located underneath a vesicle, these molecules are changing their conformation, activating the fusion of the vesicle with the cell membrane and thus triggering the release of neurotransmitters 11 . This is a key step of neuronal communication.
At the same time, calcium ions can bind and unbind to buffer molecules located in the bulk of the pre-synaptic terminal, or exit at the end of the terminal. The success of vesicular release depends on rare events: the arrival of diffusing calcium ions to the small target molecules [18][19][20] . Ions are either coming directly from the initial entering calcium flux or from ions located the bulk. Vesicular release also depend on the relative position between vesicles and the VGCCs and on their organization on the surface 10 . To investigate the time distribution of vesicular release triggered by an accumulation of T = 5 calcium ions at a ribbon (Fig. 2B), we solved numerically the system of ODE Eqs 6-14 using Matlab (see the Methods section for the description).
The initial influx of calcium ions through a VGCC J(t) is simulated using a simplified Hodgkin-Huxley model for calcium current 21 : the membrane voltage depolarizes following an action potential, modeled by a transient current I app (t). The induced calcium current lasts 2.75 ms, with a maximal value of I Ca,max = 36.2 nA (Fig. 3A). The total charge Q = 0.025 fC corresponds to an entry of 80 calcium ions (Methods section). We consider that vesicles are located on a square lattice and that channels are uniformly distributed. The probabilities q(x, i) are estimated numerically using Brownian simulations (Methods section). All parameters used for the simulations are described in Table 1.
To confirm the result of this novel Markov-mass action model (red curves in Fig. 3B), we propose an alternative approach using a Gillespie algorithm based on Poissonian rates. The Gillespie's simulations track rare binding events to small targets, which is well approximated by single exponential distribution the rate of which is the mean first binding time 22 . Due to the small size of the binding sites compared to the volume of the domain (Table 1), the motion of calcium ions in the bulk -the binding on the SNARE Complex, the binding on buffers an the escape from the domain-can all be approximated by a rate process, computed from the Narrow Escape Theory 10,23 . This approximations allows us to use the Gillespie algorithm (blue traces in Fig. 3B), using the rates presented in Table 1 (see the Methods section for a description). We found that the two models quite different in nature agree together (Fig. 3B).
When there are no buffers in the bulk (Fig. 3B), the time distribution of vesicular release is bimodal with two consecutive stages: a short time period following the immediate entrance of calcium ions. In that case, vesicular release is triggered by ions that are directly coming from the channel influx. The second regime characterized by a broader peak is induced by the random arrival of calcium ions from the bulk to the target sites, until T sites per vesicle are occupied. The large release time distribution is due to the long time scale of arrival for ions that could be bound and need to diffuse to the small targets. Using formula 3, we obtain for a volume of 1 μm 3 , D = 20 μm 2 /s and ε = 0.001 μm, a mean time τ = = = . . This time scale is much longer than the initial calcium transient from the channels and explains the long residual effect of calcium on vesicular release.
To further clarify the origin of the residual (second) peak, we tested numerically the effect of considering buffers in the synaptic bulk at various concentrations. Their role is to buffer the calcium ions in the bulk. We found that increasing the number of buffers abolishes the second peak (Fig. 3B-D), confirming that it is due to the arrival of free calcium ions located in the bulk.
We conclude this report with the following prediction: the diffusion time course in the bulk (which is modeled as a ball of radius 600 nm), where target binding sites are located underneath vesicles of size 20 nm leads to a residual vesicular release that peaks around 60 ms. However, in the case of a high buffer concentration, because the mean unbinding time to buffer is 1/k −1 = 2 s, at the time scale of tens of ms most ions are still buffered and the second peak disappears (Fig. 3C,D) (the escape rate is 1/k es = 0.7 s). However, we can still observe rare events in the distribution tail due to the arrival of free ions from the bulk.

Discussion and Conclusion
We presented here a hybrid model that connect a Markov process to a mass action model called the Markov-mass action model. The aim is to study the coupling between a continuum ensemble of particles that participate into a small number of rare events, representing the binding of a small and finite tiny subset of these Brownian molecules to small targets. It was already noticed that discrete-state stochastic models for particles are not well approximated by continuum equations 24,25 . The present model is also quite different from previous numerical simulation efforts to generate a finite number of Brownian trajectories from a continuum ensemble 26,27 .
To illustrate the applicability of this method, we model vesicular release at synapses from calcium dynamics, where we found a first peak at short time scale (less than 10 ms) induced by the direct entrance of calcium ions, followed by a second peak, which is due to the return of calcium from the bulk to the small hidden targets, a process that also depends on the interaction with buffer molecules and can last hundreds of milliseconds. The emergence of these two phases can possibly explain asynchronous quantal release, without introducing any additional time constant at a molecular level. The present approach could also be used to describe short-term plasticity at a molecular level 20,28 and to model the noise inherent to synaptic dynamics 29 .

Methods
Influx of particles. The initial influx of particles through a channel following an action potential is computed using a simplified Hodgkin-Huxley model. The membrane depolarization following an action potential is modeled using an applied current I app = 50 mV during 1 ms. The equations are summarized below The initial values V 0 , n 0 , m 0 , h 0 are the equilibrium solutions of the ODE. They are summarized with all the parameters in Table 2. The induced calcium current lasts 2.75 ms, with a peak value I Ca,max = 36.2 nA (Fig. 3A).
Numerical estimation of the splitting probability between targets. In this section, we describe the numerical simulations that we used to estimate the probability to distributed an ions among one of the N T target sites (vesicles), distributed on a square lattice. There exactly N c source points (channels) distributed in a square S of length 2H near the targets 10 (Fig. 4A,B and Table 1).
After particles either enter the bulk or bind to one of the targets. In the latter case, using Brownian simulations we determined the target. For that, we subdivide each square S into 8 sub-triangles (s i ) i=1...8 (Fig. 4B). Using symmetry we restrict our analysis to triangle s 1 . For a point source x ∈ s 1 (Fig. 4B), we estimate the probability q(x, i 1 ) that the entering particle binds to the closest target i 1 , (blue in Fig. 4B), to the second closest i 2 , q(x, i 2 ) (green or magenta) and so on, with respect to the distance dist(x) to the closest vesicle. The normalization identity is The exact numerical procedure is the following: To compute the distribution q(x, i), we further decompose the triangle s 1 into blocks parameterized by polar coordinates (ρ k , θ k ) k . We generated 200 runs for each block k and fitted the results using a routine procedure in Matlab (Fig. 4, the color code in B and C are the same). We neglect the probability to bind to the 5th-closest target and anyone further away. This numerical result allows to compute the probability q(x, i) to reach specifically the target site i, when starting from point x. Finally, the fraction of particles entering through a channel positioned at x, that reaches a target i, is p s (x)q(x, i), and the fraction of flux coming from a channel located at x and arriving to target i is J i (x, t) = J(t)p s (x)q(x, i).
Mass action model. We solve numerically the coupled system of ODE Eqs 6-14 for different numbers of buffer sites and for a uniform distribution of channels using Matlab and Monte Carlo simulations. For a given number of buffer sites and a given positioning of the N c channels, we solved the system using a Matlab solver. To get the results for a uniform distribution of channels (Fig. 3B-D Stochastic simulations based on the Gillespie algorithm. The stochastic model is a Gillespie algorithm based on Poissonian rates. The arrival time of a Brownian particle to a small hole is well approximated by a Poisson process and the rate is the reciprocal of the mean first passage time. In our model, the arrival of a Brownian particle to a small hole can represent binding to buffers, escape from the domain and binding to a target site. We shall now present the rates used in simulations.
The buffers are positioned at independent sites, modeled by small spheres. The mean first arrival time τ B of a Brownian particle to a small spherical target of radius r B in the domain Ω, with diffusion coefficients D for the particle and D B for the buffer, is given by 23 The release rate to a buffer site is Poissonian with backward rate k −1 ( Table 1). The domain is described in Fig. 2 and consists of a head, of volume |Ω h | and a cylindrical neck with radius r neck and height l neck ( Table 1). The escape at the end of the domain is Poissonian, with rate the reciprocal of the mean escape time 23 , given by where D is the diffusion coefficient. Finally, the arrival of a particle to a small ribbon (colored in red in Fig. 2B) is Poissonian with rate given in Eq. 3. These arrival times are exponentially distributed with rates equal to the reciprocal of the binding time X X when there are N f independent Brownian particles in the domain, the rate constant for the first binding event is The conditional probability that a binding event T b occurs between time t and t + Δ t, when there are N f particles is  After entry through channels, particles either bind to the targets with probability p s Eq. 2, or go to the bulk with probability 1 − p s . Using a Gillespie algorithm and the rates described above, we compute the time course of the number of particles in the bulk N f , and the time distribution of vesicular release (see also Table 1).  A). They enter the bulk when they reach the upper face (gray square in A) of the rectangular cuboid around a sphere with lower face S and height 2r. Each square S can be subdivided into triangular subunits s i (B), where Brownian simulations are used to estimate the probability that a particle reaches the closest target (blue in B,C), the second closest (green or magenta), and so on. The length is H = 70 nm (Table 1) and we show in (C) the results for θ = π/4. The color code in (B,C) is identical.