From Local to Global Modeling for Characterizing Calcium Dynamics and Their Effects on Electrical Activity and Exocytosis in Excitable Cells

Electrical activity in neurons and other excitable cells is a result of complex interactions between the system of ion channels, involving both global coupling (e.g., via voltage or bulk cytosolic Ca2+ concentration) of the channels, and local coupling in ion channel complexes (e.g., via local Ca2+ concentration surrounding Ca2+ channels (CaVs), the so-called Ca2+ nanodomains). We recently devised a model of large-conductance BKCa potassium currents, and hence BKCa–CaV complexes controlled locally by CaVs via Ca2+ nanodomains. We showed how different CaV types and BKCa–CaV stoichiometries affect whole-cell electrical behavior. Ca2+ nanodomains are also important for triggering exocytosis of hormone-containing granules, and in this regard, we implemented a strategy to characterize the local interactions between granules and CaVs. In this study, we coupled electrical and exocytosis models respecting the local effects via Ca2+ nanodomains. By simulating scenarios with BKCa–CaV complexes with different stoichiometries in pituitary cells, we achieved two main electrophysiological responses (continuous spiking or bursting) and investigated their effects on the downstream exocytosis process. By varying the number and distance of CaVs coupled with the granules, we found that bursting promotes exocytosis with faster rates than spiking. However, by normalizing to Ca2+ influx, we found that bursting is only slightly more efficient than spiking when CaVs are far away from granules, whereas no difference in efficiency between bursting and spiking is observed with close granule-CaV coupling.


Introduction
Mathematical modeling has played an important role in characterizing the electrical properties of neurons and other excitable cells. In this field, Hodgkin and Huxley were the first to propose a mathematical model for explaining ionic mechanisms underlying the generation and propagation of action potentials (APs) giant squid axons [1]. In the Hodgkin-Huxley model and most of its descendants, the system of ion channels is coupled globally via the membrane potential or the bulk cytosolic Ca 2+ concentration. However, some ion channels are colocalized, implying that the activity of one channel may affect the other via local control: there is increasing evidence for direct coupling between certain ion channels and Ca 2+ channels (CaVs) forming ion channel complexes [2][3][4]. A prominent example of ion channel complexes is the BK Ca -CaV complex: large-conductance, Ca 2+ -activated K + channels (BK Ca channels), ubiquitously expressed in excitable cells determining the electrical behavior [3], form ion channel complexes with CaVs, with a stoichiometry of 1-4 CaVs per BK Ca channel [3,4]. Thus, the whole-cell population of BK Ca channels is regulated both by global coupling, via the membrane potential, and by local coupling, via the Ca 2+ nanodomains below the mouth of the CaV channels [5][6][7][8], where the local Ca 2+ concentration reaches the tens of µM that are required for activating the BK Ca channels at physiological voltages [3,9]. Cox [10] derived a Markov chain (MC) model for characterizing the single BK Ca -CaV complex with 1:1 stoichiometry, providing important insight into the open probability of BK Ca channels during depolarizations and action potentials; for example, showing how inactivation of CaVs directly influence BK Ca channel activity. Recently, by applying MC theory [11], reaction-diffusion models [10] and time scale analysis [12] to a more realistic BK Ca -CaV complex with 1:n stoichiometry, we [13] obtained a mechanistically correct model of the BK Ca current, which respects the local effects of BK Ca -CaV coupling, and can be inserted in Hodgkin-Huxley-type models of whole-cell electrical activity: different CaV types and BK Ca -CaV stoichiometries affect BK Ca channel activity and the resulting whole-cell electrical activity in neurons and other excitable cells. This kind of local-global modeling is similar to previous work on Ca 2+ dynamics in cardiac cells [14].
Here, we present a study on how different electrophysiological behaviors determine the downstream Ca 2+ -regulated exocytosis process in endocrine cells, by which hormones are released from cells. In most endocrine cells, hormones are contained in granules that, in response to a series of cellular mechanisms culminating with an increase in the intracellular Ca 2+ levels, fuse with the cell membrane allowing the release of their content (i.e., hormones) into the extracellular environment. The main mechanisms regulating hormone exocytosis are shared with exocytosis of synaptic vesicles underlying neurotransmitter release in neurons [15,16]. Different proteins mediate the process; in particular, the soluble N-ethylmaleimide sensitive factor attachment receptor proteins (SNAREs) [17]: SNAP and syntaxin, which are located in the cell membrane and VAMP, also called synaptobrevin, inserted into the vesicle/granule membrane. The v-SNAREs (v for vesicle) contained in the granule can form, with t-SNAREs (t for target) inserted in the cell membrane, the so-called SNARE complex [15], driving fusion of the two membranes, which-in the case of endocrine cells-allows the hormone molecules contained in the granule to exit the granule and enter the blood stream. SNARE complexes interact with other proteins, notably, Ca 2+ -sensing proteins, such as synaptotagmins, which trigger exocytosis upon Ca 2+ binding. Therefore, the local Ca 2+ concentration at the Ca 2+ sensor of the exocytotic machinery is a key factor determining the probability (rate) of exocytosis of the secretory granule [18][19][20]. Experimental evidence of local coupling between single granules and CaVs showed that the exocytosis rate of a single granule increases significantly when it is close to CaVs [21]. Recently, we developed a method [22] for characterizing the local interactions between the single granule and the surrounding CaVs by exploiting a strategy that is similar to the methodology devised for describing the BK Ca current [13]: using absorbing MC models allows for achieving analytic results for the expected exocytosis rate of a single granule, showing how coupling different numbers of CaVs at difference distances with the granule determines different responses.
In the following, we use the devised BK Ca current model [13] as an example for analyzing the behavior (i.e., electrical activity) of a biologically realistic and important system (i.e., an endocrine pituitary cell) composed of units (i.e., ion channels) that interact both globally and locally: varying the number of CaVs per BK Ca -CaV complex results in different electrophysiological responses of pituitary cells, in particular, continuous spiking and so-called pseudo-plateau bursting, the latter characterized by few small oscillations riding on a depolarized plateau. Then, we couple electrical activity and exocytosis by combining the devised models [13,22] in order to investigate how local BK Ca -CaV coupling via continuous spiking or bursting in pituitary cells determine the downstream exocytotic response, by varying the distance and number of CaVs coupled with the single granules. In particular, we exploit MC models for describing the local coupling between CaVs and BK Ca -channels, and between CaVs and granules, and stationary approximations for characterizing the local Ca 2+ levels, which allow us to couple electrical activity and exocytosis in a straightforward manner. Moreover, we reduce model complexity to achieve simplified ordinary differential equation (ODE) models that respect the local control by CaVs of BK Ca channels and granules.

Whole-Cell Electrical Activity Modeling Respecting Local Control
Local coupling of ion channels is important in determining whole-cell activity. BK Ca -CaV complexes provide an example of ion channel complexes, where the BK Ca activity is influenced locally by associated CaVs, and differences in stoichiometry of the complex (1-4 CaV channels per BK Ca channel [3,4,23]) can affect not only single channel activity but also whole-cell behavior significantly [13].
A model of pituitary cells has been used to study the role of BK Ca channels in the generation of so-called plateau bursting, which consists of a few small oscillations riding on a depolarized plateau, and is important for secretion [24,25]. In the original model [26] the BK Ca current was modeled as a purely voltage-dependent current, neglecting Ca 2+ dependency. In [13], we substituted this simplified representation of the BK Ca current with a concise model that respects the local effects of BK Ca -CaV coupling and can be inserted in the Hodgkin-Huxley-type model of electrical activity of pituitary cells (see Methods, Equation (6)). For characterizing the whole-cell BK Ca current I BK , the state of every single BK Ca channel does not need to be known; it is sufficient to determine the BK Ca open probability p (n) Y over time, where the superscript n indicates the number of CaVs coupled with the single BK Ca channel, and I BK in Equation (6) has the the following form:  Figure 1A), whose parameters are set by fitting the available experimental data [10] (see Methods and Figure 1B) and a 3-state model for the CaV channel ( Figure S1C). By assuming the corresponding 6-state ODE model and using time-scale analysis (in particular, CaV inactivation is slow compared to other processes), we split the model in two submodels, one including states with non-inactivated CaV (green box in Figure 1C) and the other with inactivated CaV (blue box in Figure 1C). Moreover, the submodel describing the BK Ca channel activation in BK Ca -CaV complex with non-inactivated CaV, denoted with m (1) BK , can be described by one ODE of Equation (13) (see Methods)-a quasi-steady state approximation. The BK Ca channels in BK Ca -CaV complexes exhibit inactivation because of inactivation of the associated CaVs, and with approximately identical dynamics, as found in experiments [27] and Monte Carlo simulations (see Figure 1; [10]). Then, the open probability for BK Ca -CaV complex with 1:1 stoichiometry, p (1) Y , can be expressed as with m BK given by Equation (13) and h represents the inactivation function of the CaVs (see Methods). We make a further simplification assuming instantaneous CaV activation, since in many whole-cell models (e.g., [26,[28][29][30]), the Ca 2+ currents are assumed to activate instantaneously: in this case, i.e., instantaneous CaV activation, m CaV = m CaV,∞ , where m CaV and m CaV,∞ represent CaV activation variable and its steady-state, respectively (see Equations (8) and (9)); m (1) BK is given by Equation (13) with the approximations defined by Equations (16) and (17); and h = 1 − b with b by Equation (10) (see Methods).       As performed in [13], we compared the devised 6-state MC model of the BK-CaV complex with the first stochastic model of the gating of this complex devised by Cox [10], obtained by coupling the 10-state MC model for the BK Ca channel with a 7-state MC model for the CaV channel (see Methods and Figure   BK given by Equation (13) with the approximations of Equations (16) and (17)) approximates the full system decently, except for the initial phase before CaV activation reaches equilibrium. Finally, Figure 1E-H shows how the simplified BK Ca models reproduce, satisfactorily, the whole-cell BK Ca currents I BK for each step in voltage (I BK defined by Equation (1) with n = 1, where g BK = N BKḡBK , withḡ BK = 100 pF [31] the single-channel conductance and N BK = 1000 the number of BK Ca channels): the simplified 6-state Markov chain model (black plots in Figure 1F-H) and the corresponding 6-state ODE model (blue plots in Figure 1F) approximate the 70-state Markov chain model current very well ( Figure 1E; [10]); the simplified Hodgkin-Huxley-type model current for the BK Ca channel with p (1) (2); red plots in Figure 1G), and the corresponding model assuming instantaneous activation of the CaV currents (m (1) BK given by Equation (13) with the approximations of Equations (16) and (17); green plots in Figure 1H) also work very well.

BK Ca -CaV Complex with 1:n Stoichiometry
In the case of more than one CaV per BK Ca channel, the linear buffer approximation is used to compute the Ca 2+ profile from n channels by superimposing n nanodomains found for single, isolated CaVs. Then, the MC model of Figure 1C can be extended to a model with 3 × n × 2 states. However, as discussed previously, we assumed that, on a fast timescale, the fraction h of non-inactivated CaVs is constant, and note that the BK Ca channel closes rapidly when all CaVs in the complex are inactivated. Hence, for the case of 1:n BK Ca -CaV stoichiometry, we split the system according to the number k of non-inactivated CaVs: the BK Ca activation can be described on a fast time scale by the Markov chain model of Figure (26)). Then, the open probability of the BK Ca channel coupled with n CaVs, p (n) Y , can be estimated by taking into account that the probability of k non-inactivated CaVs being present in a complex with n CaVs, As performed in [13], we compared the simulated open probabilities obtained from the different models in response to a voltage step from −80 to 0 mV, and the results reported in Figure 2C show how the different ODE models exploited for computing m which involves one ODE (Equation 26) for m (n) BK . Figure 2D-F shows the simulated whole-cell BK Ca currents I BK for the case of 1:4 BK Ca -CaV stoichiometry (defined by Equation (1), where g BK = N BKḡBK , withḡ BK = 100 pF [31] as in the previous section, while N BK is lower and equal to 700 in order to reproduce the experimental data; i.e., the maximum amplitude values reported in [10,32]) obtained from the different ODE models (exploited for computing m   (26) with Equations (28) and (29) (green curves)) approximate very well the Monte Carlo simulations (black curve) obtained from the 3 × 4 × 2 MC model.

Concise Whole-Cell Modeling Respecting Local Control
The Hodgkin-Huxley-type model of the BK Ca current defined by Equation (1), that allows us to take into account the local interactions in BK Ca -CaV complexes, can be exploited to investigate how the stoichiometry of the complex affects whole-cell electrical behavior of pituitary cells.
In the case of BK Ca -CaV complexes with 1:1 stoichiometry, spiking electrical activity is observed as shown in Figure 2G, since insufficient BK Ca current is generated. The different plots in Figure 2G   In contrast, in the case of BK Ca -CaV complexes with 1:n stoichiometry with n > 1, plateau bursting appears with the number of small oscillations per burst depending on the number of CaVs per BK Ca -CaV complex, as shown in Figure 2H,I for n = 2 and n = 4, respectively. The different plots in Figure  Although the quantitative behavior is independent of the approximation for m (n) BK , minor qualitative differences are present. The approximation given by Equation (26) reproduces very well, the behavior obtained from the complete ODE model for the activation of the BK Ca channel surrounded by n non-inactivated CaVs defined by Equation (18) (Figure 2G-I, upper and middle panels), whereas the further simplification given by Equations (28) and (29) produces smaller and more spikes per burst (lower panels). Nonetheless, considering parameter uncertainties and experimental variations, even Equations (28) and (29) produce reliable results.

Coupling Electrical Activity with Exocytosis
We investigated how cellular electrical activity regulates the downstream exocytosis process. We modeled the single granule containing the hormones released from the cell by Ca 2+ -regulated exocytosis process with the 5-state MC model of Figure 3A: the granule can be in one of four states depending on the number of Ca 2+ ions bound to the Ca 2+ sensor on the granule (states G i , with i = 0, . . . 3, representing the number of Ca 2+ ions bound to the granule sensor) before fusing with the membrane and releasing its hormone content (state Y) (see Methods for more details). The Ca 2+ concentration at the granule sensor, Ca CaV (r G ) given by Equation (7) with r = r G being the distance from the CaV pore to the Ca 2+ sensor on the granule, drives the exocytosis MC model, allowing the granule to modify its state through the Markov chain from G 0 to Y and undergo exocytosis.
In the following, in order to characterize the local interactions between granules and CaVs, we analyzed the property of the secretory complex obtained by coupling a single granule with one or more CaVs: we coupled the granule, described by the 5-state MC model of Figure 3A, with n CaVs, each one described by the 3-state MC model of Figure S1A, and obtained the (5 × n × 3)-state MC model, as developed in [22]. For the electrical activity model in pituitary cells defined by Equation (6), CaVs do not inactivate, and the exocytosis can be described by the 5 (n + 1)-state MC model shown in Figure 3C. Using this model, the granule state depends on the CaVs states. In particular, for the case with one CaV, the Ca 2+ concentration at the granule sensor needed for triggering exocytosis, is equal to a basal level (Ca c ) when the CaV is closed or inactivated, and is equal to Ca CaV (r G ) defined by Equation (7) when the CaVs is open; for a more general case with n CaVs, the linear buffer approximation [7] is used to summarize Ca 2+ levels at the granule sensor when more than one CaV is open, as performed for BK Ca -CaV complex. In [22], we used phase-type distribution results for Markov chains [11] for estimating the expected exocytosis rate (the release probability) of a single granule. We found that the distance r G is a major factor in determining the exocytosis rate, as recently demonstrated and quantified explicitly [21]. Furthermore, and in agreement with the experiments [21], the results showed that increasing the number of CaVs coupled with the granule determines a much higher rise of the exocytosis rate, which in the case of inactivating CaVs is more pronounced when the granule is close to CaVs (about 10 nm), whereas for non-inactivating CaVs the highest relative increase in rate is obtained when the granule is far from CaVs (about 50 nm), suggesting that it is not necessary that the granule is very close to CaVs for triggering exocytosis. However, in [22] we did not take into account the coupling between electrical activity and exocytosis, assuming constant values for the membrane potential. Here, the Ca 2+ dynamics, and hence, the CaVs states, are driven by the electrophysiological behavior, which depends on the local interactions between CaVs and BK Ca channels, as shown in the previous section. In particular, we studied how the two typical electrophysiological behaviors (spiking or bursting), observed in pituitary cells due to different local coupling between CaVs and BK Ca channels, determine the granule release by varying the number and the distance of CaVs coupled with the granules.   Figure 4A-D shows the time evolution of the exocytosis probability of a single granule coupled with different numbers of CaVs placed at difference distances. In each panel, for the granule coupled with n CaVs placed at fixed distance r G , we assumed that CaVs dynamics are driven by continuous spiking in the form of Figure 2G (upper plot) or bursting, as in Figure 2H (upper plot). When the CaVs are close to the granule (see Figure 4A,B for r G = 10 and 20 nm, respectively), the bursting pattern (dash-dotted lines) evokes release at a higher rate than the spiking pattern (solid lines) for a reduced number of CaVs coupled to the granule (the blue and magenta lines for n = 1 and 2, respectively), while this difference in the rate decreases when the number of CaVs increases (see the cyan and black lines for n = 4 and 8, respectively). When the CaV are far from the granule (see Figure 4C,D for r G = 50 and 100 nm, respectively), bursting promotes exocytosis with faster rates than spiking even for a higher number of CaVs coupled to the granule. This difference in secretion rate between spiking and bursting is mainly due to larger amount of Ca 2+ entering during bursting, which becomes negligible when the number of CaVs is high and close to the granule. In this scenario, the Ca 2+ concentrations at the exocytotic machinery are in saturation regimes for both the electrical patterns.

Spiking versus Bursting on Evoking Exocytosis
Therefore, in order to evaluate the exocytotic machinery performance with the same Ca 2+ entry, we analyzed the exocytosis probability with respect to the total Ca 2+ influx, Q Ca . Figure 4E-H shows the granule release probability versus Q Ca for the same cases simulated in Figure 4A-D. From this analysis, it is seems that the efficiency in evoking exocytosis between spiking and bursting is similar. This finding is also confirmed by fitting the simulated responses using an exponential function f e defined as with the aim to estimate the parameter q, whose value can provide insight into the efficiency in evoking granule release (a high value of q means high efficiency). Figure 5A shows the trend of q with respect to the distance r G of the granule from different numbers of coupled CaVs (different colors) for the two electrical patterns, spiking (solid lines) and bursting (dash-dotted lines). When the granule is close to CaVs, independently of their number, bursting and spiking have a similar effect on exocytosis, whereas, when the granule is far from CaVs, in the case of few CaVs (one or two), bursting is slightly more efficient than spiking. For all the cases, the fitting approximates the simulated data well, as shown for two cases in Figure 5B: the upper plot, reporting the simulated exocytosis probabilities of a granule at distance r G = 20 nm from four CaVs for spiking and bursting patterns, and the corresponding fitting shows how there is no virtually difference between the two electrical patterns in evoking exocytosis; the lower plot, displaying the simulated exocytosis probabilities of a granule at distance r G = 100 nm from two CaVs for spiking and bursting patterns and the corresponding fitting, shows how there is a small difference between spiking and bursting, with the latter resulting more efficient in granule release.       (5) with respect to the distance r G of the granule from n CaVs (n = 1 (blue curves), n = 2 (magenta curves), n = 4 (cyan curves) and n = 8 (black curves)), driven by spiking (solid lines) or bursting (dash-dotted lines). Note that for each estimate the 95% confidence interval is very limited. (B) Fit to the to the simulated exocytosis probabilities versus Q Ca for the granule at r G = 20 from four CaVs (upper plot) and for the granule at r G = 100 from two CaVs. The simulated responses are the averages of one-thousand Monte Carlo simulations.

Discussion
In this paper, we show the role of mathematical modeling as an important tool for investigating excitable cells with focus on ion Ca 2+ channels and their local interactions with BK Ca potassium channels in influencing electrical activity of pituitary cells and with hormone-containing granules determining the granule release by exocytosis process. Therefore, whole-cell models have to be consistent with the local mechanisms operating at molecular levels, and hence, we show how to exploit all the available information in a coherent and structural way by using mathematical modeling in order to handle the complexity of cellular electrophysiology.
In order to handle the local interactions in BK Ca -CaV complexes with 1:n stoichiometry, we used a stochastic model based on Markov chain theory (of 3 × n × 2 states; see Figure 2), as a starting point for analyzing the single complex dynamics. However, the fluctuations resulting from stochastic gating kinetics observed at the molecular level tend to become negligible as a system's size approaches the whole-cell scale, where, for describing the electrical behavior, it is not necessary to know the state of each single complex, but it is sufficient to know the open probability of BK Ca channel population. Therefore, we used the corresponding deterministic model, and by exploiting time-scale analysis and quasi-steady state approximation, we achieved a concise, deterministic Hodgkin-Huxley-type model of BK Ca currents defined by Equations (1) and (3). This approach allowed us not only to reduce model complexity and computational costs of the stochastic model, but also to achieve an explicit interpretation of the parameters and their effects on whole-cell behavior through the direct read of the formula of the simplified deterministic model: We showed that increasing the number of CaVs coupled with BK Ca channel determines a left shift of the BK Ca activation curve, since the probability of at least one CaV being open is greater with more channels in the complex (see Equations (28) and (29) assuming instantaneous activation of CaVs). We also derived in [13] an analytic expression for the time to first opening of a BK Ca channel, providing theoretical insight into stochastic simulation results. Then, the concise Hodgkin-Huxley-type model of whole-cell BK Ca currents, that allows taking into account local control by CaVs in BK Ca -CaV complexes, can be inserted into models of cellular electrical activity, showing, in particular, how different BK Ca -CaV stoichiometries cause different electrophysiological responses, including continuous spiking and bursting in pituitary cells. In [13], we also showed how BK Ca -CaV stoichiometry controls the fast after-hyperpolarization (fAHP) in a model of hypothalamic neurosecretory cells; i.e., the undershoot seen after an action potential controlling firing frequency and transmitter release [3,33]. Moreover, coupling BK Ca channels with different CaV types affects electrical activity differently in human pancreatic beta-cells, and further insight into the control by CaVs of BK Ca channels, whose block stimulates insulin secretion in human [34] and mouse [35] beta-cells, may lead to a better understanding of beta-cell function and how it becomes disturbed in diabetes. Furthermore, this approach should be relatively straightforward to apply to other ion channel complexes; e.g., the CaV3-Kv4 complex [2].
Recently, we [22] exploited the methodology devised for modeling the BK Ca currents for handling the local interactions between granules and CaVs, and specifically, by using phase-type distribution results for Markov chains [11], we obtained analytic results for the expected exocytosis rate of a granule coupled with different numbers of CaVs placed at different distances. We also exploited a quasi-steady state approximation for the corresponding ODE model of the 5-state MC model for exocytosis of a single granule adjacent to the plasma membrane, as shown in Figure 3A, in order to reduce the model complexity, especially in the case of the granule coupled with n CaVs. However, in this case, the quasi-stead approximation only works for the final state of the chain (i.e., state G 3 ), before the granule can fuse with the cell membrane (i.e., state Y), since its dynamics are the fastest. The sequence of the different states of the MC for the granule before fusing, according to the number of Ca 2+ ions bound to the Ca 2+ sensor on the granule, allows us to reproduce the delayed exocytosis with respect to a raise in the calcium concentration, as observed by flash-release experiments [36][37][38], and in this study, we used the complete MC models of Figure 3 for describing granule exocytosis. Hence, by modeling the local Ca 2+ dynamics, we coupled exocytosis and electrical models with the aim to investigate how the two main electrophysiological behaviors in pituitary cells, continuous spiking and bursting, affect the downstream exocytosis process. From our results, we found that, surprisingly, bursting is only slightly more efficient than spiking, and only when CaVs are few and far from the granule. These differences to previous findings [25] can be explained by the fact that bursting has an important role in the resupply of the primed granule pool, which depends on the global, rather than local, Ca 2+ concentration. Indeed, as experimentally observed [39], global Ca 2+ concentration (i.e., the bulk cytosolic Ca 2+ concentration) is higher during bursting than spiking. We did not model the resupply, since we assumed that the granule was already primed for exocytosis.
For our aims, we used steady-state reaction-diffusion equations for characterizing Ca 2+ levels, in particular, for calculating calcium concentrations at ion channel BK Ca -CaV complexes and at granules, as performed in previous works: Cox [10] found that, for computing Ca 2+ levels at a BK Ca channel coupled with one CaV, steady-state solutions of reaction-diffusion equations solved by CalC approximate the numerical solutions of these equations well, and we confirmed these results in [13], assuming more realistic cases with BK Ca channel coupled with n CaVs. In order to get a deeper description of the Ca 2+ levels, we can exploit compartmental modeling, as performed in our previous work [40], where we characterized the intracellular Ca 2+ dynamics in glucagon-secreting pancreatic alpha-cells. Compartmental modeling allows us to couple electrical activity and exocytosis, providing a deeper knowledge of the relative contributions of the various sub-cellular compartments (Ca 2+ nanodomains, sub-membrane compartment, bulk cytosol and endoplasmic reticulum) involved in exocytosis.

Electrophysiological Model of Pituitary cells
The original model of pituitary cells [26] included a single Ca 2+ current (I Ca ), a delayed-rectifier K + current (I K ), a Ca 2+ -gated SK current (I SK ) and a leak current (I leak ), in addition to the BK Ca current which was modeled as a purely voltage-dependent current, neglecting Ca 2+ dependency. We substituted this simplified representation of the BK Ca current with our concise BK Ca model controlled by CaVs in complexes. Then, the dynamics of membrane potential V are described by the following ODE: where C is the membrane capacitance and I BK is defined by Equation (4) since CaVs do not inactivate.
All the other currents are modeled as in [26]. In order to achieve a Hodgkin-Huxley-type model of BK Ca currents (Equation (4)) that take into account local control in BK Ca -CaV complexes with different stoichiometries, we first devised a model of a single-channel gating for the BK Ca channel and then we coupled the model with different number of CaVs, each one described by a standard 3-state model, which can be reduced to two states, since CaVs do not inactivate, as reported in the following.

BK Ca Channel Modeling
For describing the BK Ca channel, we used a single-channel gating with two states, as shown in Figure 1A, where X corresponds to the closed state and Y to the open state. The dynamics of the channel are determined by voltage and calcium-dependent rates, k − , and k + , describing the transition from the open to close state and from close to open state, respectively. By assuming that, for fixed Ca 2+ concentration (Ca), BK Ca activity is described by a Boltzmann function [9,27,41] and that the slope parameter of the Boltzmann function is independent of Ca, a reasonable assumption for Ca 2+ concentrations above 1 µM [10,41,42], as expected in BK Ca -CaV complexes [3,10], in [13], we showed that k − and k + can be expressed as a product of voltage and Ca 2+ -dependent terms (w − (V) and w + (V) are the voltage dependent rates, f − (Ca) and f + (Ca) the Ca 2+ -dependent rates; then, k − = w − f − , k + = w + f + ; see [13] for mathematical expressions). Figure 1B shows the fitting to the experimental data [10], consisting of BK Ca open probabilities and time constants as functions of voltage, at different Ca 2+ concentrations, by using for k − and k + -the optimal parameters estimated in [13].
We also reproduced the dynamics of a more complex model for describing the BK Ca channel proposed by Cox's lab [10,42]. They assumed that each alpha subunit (four subunits overall), four of which compose the tetrameric structure of the channel [43], has a single Ca 2+ binding site and a single voltage sensor, and through simplifications, they obtained a 10-state model, where each state can have from zero to four bound Ca 2+ ions and be open or closed. As shown in [10], the 10-state model is able to reproduce the characteristics and dynamics of the BK Ca channel by fitting BK Ca open probabilities and time constants as functions of voltage, at different Ca 2+ concentrations. Note that, although the model is complex, it represents an empirical model: individual rate constants are not likely correspond to any real calcium binding events or movements of the channel's voltage sensors, since real BK Ca contains at least three Ca 2+ binding sites (a low and two high-affinity calcium binding sites) [43][44][45] and four voltage sensor per subunit [46].

CaV Modeling
We described CaV by using the 3-state ODE model of Figure S1A (see Equations in [13]), where C corresponds to the closed state, O to the open state and B to the inactivated (blocked) state of the calcium channel [47]. α and β represent the voltage-dependent Ca 2+ channel opening rate and closing rate, respectively, and have the forms as in [47]; γ is a constant reverse reactivation rate and δ represents the Ca 2+ -dependent rate for channel inactivation, which is determined by the Ca 2+ concentration at the Ca 2+ sensor, Ca CaV (r), having the following form by using reaction-diffusion theory [7,10,48]: r represents the distance of the sensor from the channel pore (in this case r = 7 nm) and i Ca max =ḡ Ca (V − V Ca ) is the single-channel Ca 2+ current withḡ Ca the single-channel conductance and V Ca the Ca 2+ reversal potential. The formula defined by Equation (7), called excess buffer approximation (EBA), is based on the assumption that the buffer is unsaturable [48], while another common formula, called rapid buffer approximation (RBA), is valid for buffers that are saturated near an open channel and have Ca 2+ binding kinetics that are rapid relative to Ca 2+ diffusion [49,50].
Here, we used EBA, as was done by Cox in his work [10], from which we took the data for our study. Cox found that EBA well approximates the solutions of reaction-diffusion equations solved by the simulation software CalC [51] for computing Ca 2+ levels at BK Ca channel with the buffer conditions specified in his work, and we confirmed these results in [13] with more realistic cases with BK channel coupled with n CaVs. Note that for the single-channel Ca 2+ current i Ca max , we use the same formalism used by Cox (Ohm's law), although the use of the Goldman-Hodgkin-Katz equation may be more appropriate. As shown in [47], the processes of activation and inactivation can be approximately separated in time, since activation is much faster than inactivation. In particular, we achieve the following model for the CaV activation variable, m CaV , where and the following equation for inactivation Therefore, assuming instantaneous activation, m CaV = m CaV,∞ , the 3-state system can be approximated by 1-state ODE model described by Equation (10). Note that we define with h the CaV inactivation function, representing the fraction of Ca 2+ channels not inactivated, h = o + c, where o and c represent the state variables of 3-state system of Figure S1A. For 1-state model, h = 1 − b, with b given by Equation (10).
We also reproduced a more complex model for describing CaVs, proposed by Cox [10], based on the model of Boland and Bean [52], where the channel is described by a 7-state Markov chain ( Figure S1B), with the states C i (i = 0, . . . , 4) used for representing the movements of four voltage sensors. Here, the transitions labeled α and β represent the movements of four voltage sensors and are voltage-dependent, having the forms as in [47]. When all four subunits have moved to an activated position, i.e., the channel is in state C 4 , the complex can undergo a final voltage-independent step to the open state O with a constant rate ; ζ represents the reverse rate from O to C 4 . The inactivation rate δ and the reverse reactivation rate γ have the form as for the 3-state model. Figure S1C shows the fitting to the data [10] (the red circles in the plots; i.e., the mean values of peak open probabilities (upper plot) and time constants (lower plot)) of the 3-state ODE model (the blue curves) and 7-state Markov chain (the gray curves), both satisfactorily reproducing activation curves and times. For the 3-state ODE model, we used the estimated parameters reported in [13]. Moreover, Figure S1D shows the simulated CaV open probabilities in response to three different voltage steps from −80 mV to −20 (left), 0 (middle) and 20 mV (right), obtained from the 7-state Markov chain model (gray curves), the 3-state ODE model (blue curves) and the corresponding 3-state Markov chain model (dash-dotted black curves), and the 1-state ODE model defined by Equation (10) assuming instantaneous activation m CaV = m CaV,∞ (green curves). The CaV open probabilities for the 3-state ODE model and the corresponding 3-state MC model in response to the different voltage steps approximate the 7-state MC model very well. Note that in response to voltage step from −80 to −20 mV (see left plot in Figure S1D), the 7-state MC model shows an initial delay for the CaV open probability (see dash-dotted gray plot) due to an overestimation of the time constant of the model (compare the gray dash-dotted curve with the red data in the lower plot in Figure 1C for V = −20 mV). Finally, Figure S1F-H shows the whole-cell CaV currents of the different models with different step voltages ( Figure S1E). The whole-cell CaV current is defined by using the Hodgkin-Huxley formalism: where N CaV = 1000 is the number of CaV channels, each one characterized by a single channel conductanceḡ Ca (ḡ Ca = 2.8 pS as in [53] for In particular, the Ca 2+ levels sensed by the BK Ca channel are assumed to reach steady-state immediately after CaV closure or opening [10]; in this case, Ca c is set equal to 0.2 µM (background Ca 2+ concentration) and Ca o is given by Equation (7) with r = 13 nm representing the distance between CaV and BK Ca channels [3,10]. At V = 0 mV, Ca o ≈ 19 µM. Note that k + c ≈ 0, since the background Ca 2+ concentration Ca c is much below the levels needed for BK Ca activation at physiological voltages; i.e., the probability of BK Ca opening when the CaV is closed is practically zero: this approximation is supported by the fact that Ca 2+ influx via CaVs is needed to open BK Ca channels [54], and that the sub-membrane Ca 2+ concentration of some hundreds of nM that a BK Ca in a complex without open CaVs would sense, is too low to activate BK Ca channels at physiological voltages [3,10].
Next, we evaluated the dynamics of the deterministic ODE model corresponding to the 6-state MC model.
The 6-state ODE model characterizing p Y (which can be described by a system of five ODEs because the probabilities sum to 1) can be further simplified by applying timescale analysis. Indeed, the inactivation and reactivation of CaV are slower than (de-)activation; thus, on a fast timescale, the fraction of non-inactivated CaVs (h = 1 − (p BX + p BY )) is assumed to be constant, and the model can be split into two submodels with, respectively, four and two states (indicated by green and blue boxes in Figure 1C). Moreover, the 4-state ODE model of the corresponding reduced 4-state MC (green box in Figure 2C) describing the BK Ca channel activation in BK Ca -CaV complex with non-inactivated CaV can be further simplified. Indeed, by exploiting quasi-steady state approximation, we obtained a 1-ODE model for describing the BK Ca activation, denoted with m with steady-state and time constant given by The CaV activation variable m CaV is routinely characterized in patch clamp experiments and included in models of electrical activity via the time-constant, τ CaV , and the steady-state activation function, m CaV,∞ (see Equation (9)). From these quantities, α = m CaV,∞ /τ CaV and β = 1/τ CaV − α can be calculated. Note that Equation (15) makes it explicit how m BK,∞ inherits properties of the associated Ca 2+ channel type, as has been found experimentally [55,56].
In many whole-cell models (e.g., [26,[28][29][30]), the Ca 2+ currents are assumed to activate instantaneously, which precludes calculation of α and β. Implicitly, such models assume that CaV gating is infinitely faster than the kinetics of other channels in the model. In our setting, this assumption corresponds to investigating the BK Ca -CaV model defined by Equations (13)- (15) in the limits α, β → ∞. In this case Equations (14) and (15) become  [3,4,23]. Here, we introduce a concise but mechanistically correct model of single BK Ca -CaV complexes with 1:n stoichiometry developed in [13] to be inserted into a whole-cell model of electrical activity.
We extend the 6-state MC model for the complex with 1:1 stoichiometry to incorporate different stoichiometries assuming that n CaVs are all located 13 nm from the BK Ca channel [3,10,57]: In this case, the linear buffer approximation is used to compute the Ca 2+ profile from n channels by superimposing n nanodomains found for single, isolated CaVs. Then, the MC model of Figure 1C can be extended to a model with 3 × n × 2 states. However, as discussed in the previous section, we assume that, on a fast timescale, the fraction h of non-inactivated CaVs is constant, and note that the BK Ca channel closes rapidly when all CaVs in the complex are inactivated. Hence, for the case of 1:n BK Ca -CaV stoichiometry, we split the system according to the number k of non-inactivated CaVs: the BK Ca activation can be described on fast time scale by the Markov chain model of Figure 2A By taking into account that and renaming the state variables as follows . . .
the ODE system is reduced from 2(k + 1) to (k + 1) equations. Moreover, assuming the quasi-steady state approximation for p Y i , with i = 0, . . . , k − 1, the ODE system of (k + 1) equations is reduced to the following single ODE, describing the dynamics of the BK Ca activation with k non-inactivated CaVs: where m By assuming instantaneous activation of CaVs, the model defined by Equation (26) can be further simplified. Indeed, in this case (i.e., m CaV = m CaV,∞ ) the vertical transitions in Figure 2A are in quasi-equilibrium, and then and m (k) BK follows Equation (26) with By taking into account that the probability that k non-inactivated CaVs are present in a complex with n CaVs is ( n k )h k (1 − h) n−k , we obtain Equation (3) to describe the BK Ca channel open probability for BK Ca -CaV complex with 1:n stoichiometry.

Exocytosis Model
We modeled the single granule, adjacent to the plasma membrane and primed for exocytosis, as performed in [22], by the Markov chain shown in Figure 3A. The granule can be in one of four different states depending on the number of Ca 2+ ions bound to the Ca 2+ sensor on the granule, likely synaptotagmin [58]: in G 0 with no bound Ca 2+ ions, in G 1 with one, in G 2 with two or in G 3 with three bound ions. Once it is in G 3 , the granule can fuse with the membrane and release its hormone content, assuming the final state Y [25,59]. Therefore, the model takes values in the state space S = {G 0 , G 1 , G 2 , G 3 , Y} and its transition rate or generator matrix M G is given by where represents the Ca 2+ binding rate, with Ca CaV (r G ), the Ca 2+ concentration at the granule sensor, given by Equation (7) with r = r G being the distance from the CaV to the Ca 2+ sensor on the granule. In this paper, the distance from the CaV to the granule means the distance from the CaV to the Ca 2+ sensor on the granule, which will be of the order of tens of nm. For comparison, secretory granules have diameters on the order 100-500 nm [60][61][62][63]. We assumed a constant number of Ca 2+ sensor molecules, which was, therefore, included in the binding parameter k Ca . The parameter k − was the unbinding rate, and u was the fusion rate. The values for k − , k + and u were equal to those used in [22]. Note that, since G 3 dynamics are fastest (the value of u is much higher than those of the other parameters), the 5-state MC model can be reduced to a 4-state MC model, where the dynamics of states G 2 and G 3 can be described by an auxiliary variable G 23 , using quasi-steady state approximation for the corresponding ODE model, as performed in [22]. Here, we exploit the complete sequence of the 5-state MC model in order to reproduce the delayed exocytosis with respect to a raise in the Ca 2+ concentration, as observed by flash-release experiments [36][37][38], although the 4-state MC model does not modify the results substantially. Instead, a further state-reduction can significantly modify the exocytosis probability.
By coupling the 5-state exocytosis model with the 3-state model of Figure S1A for CaV, we obtained the 15-state MC model of Figure 3B. For the general case, where the granule is coupled with n CaVs, we obtain a (5 × n × 3)-state MC model (see [22] for model details and mathematical description). For the electrical activity model in pituitary cells described by Equation (6), CaVs do not inactivate, and then, each CaV can be described by single-channel gating with two states, closed, C, and open, O, and the granule release by 5 (n + 1)-state MC model, as shown in Figure 3C.