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.

1. 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 Ca2+ 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 Ca2+ channels (CaVs) forming ion channel complexes [2,3,4]. A prominent example of ion channel complexes is the BKCa–CaV complex: large-conductance, Ca2+-activated K+ channels (BKCa 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 BKCa channel [3,4]. Thus, the whole-cell population of BKCa channels is regulated both by global coupling, via the membrane potential, and by local coupling, via the Ca2+ nanodomains below the mouth of the CaV channels [5,6,7,8], where the local Ca2+ concentration reaches the tens of M that are required for activating the BKCa channels at physiological voltages [3,9]. Cox [10] derived a Markov chain (MC) model for characterizing the single BKCa–CaV complex with 1:1 stoichiometry, providing important insight into the open probability of BKCa channels during depolarizations and action potentials; for example, showing how inactivation of CaVs directly influence BKCa channel activity. Recently, by applying MC theory [11], reaction-diffusion models [10] and time scale analysis [12] to a more realistic BKCa–CaV complex with 1:n stoichiometry, we [13] obtained a mechanistically correct model of the BKCa current, which respects the local effects of BKCa–CaV coupling, and can be inserted in Hodgkin–Huxley-type models of whole-cell electrical activity: different CaV types and BKCa–CaV stoichiometries affect BKCa 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 Ca2+ dynamics in cardiac cells [14].

Here, we present a study on how different electrophysiological behaviors determine the downstream Ca2+-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 Ca2+ 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, Ca2+-sensing proteins, such as synaptotagmins, which trigger exocytosis upon Ca2+ binding. Therefore, the local Ca2+ concentration at the Ca2+ 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 BKCa 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 BKCa 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 BKCa–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 BKCa–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 BKCa-channels, and between CaVs and granules, and stationary approximations for characterizing the local Ca2+ 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 BKCa channels and granules.

2. Results

2.1. Whole-Cell Electrical Activity Modeling Respecting Local Control

Local coupling of ion channels is important in determining whole-cell activity. BKCa–CaV complexes provide an example of ion channel complexes, where the BKCa activity is influenced locally by associated CaVs, and differences in stoichiometry of the complex (1–4 CaV channels per BKCa 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 BKCa 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 BKCa current was modeled as a purely voltage-dependent current, neglecting Ca2+ dependency. In [13], we substituted this simplified representation of the BKCa current with a concise model that respects the local effects of BKCa–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 BKCa current , the state of every single BKCa channel does not need to be known; it is sufficient to determine the BKCa open probability over time, where the superscript n indicates the number of CaVs coupled with the single BKCa channel, and in Equation (6) has the the following form:

2.1.1. BKCa–CaV Complex with 1:1 Stoichiometry

In the case of BKCa–CaV complex with 1:1 stoichiometry, in order to compute , we devised a 6-state MC model obtained by coupling a 2-state model for the BKCa channel (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 BKCa channel activation in BKCa–CaV complex with non-inactivated CaV, denoted with , can be described by one ODE of Equation (13) (see Methods)—a quasi-steady state approximation. The BKCa channels in BKCa–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 BKCa–CaV complex with 1:1 stoichiometry, , can be expressed as with 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 Ca2+ currents are assumed to activate instantaneously: in this case, i.e., instantaneous CaV activation, , where and represent CaV activation variable and its steady-state, respectively (see Equations (8) and (9)); is given by Equation (13) with the approximations defined by Equations (16) and (17); and with b by Equation (10) (see Methods).

Figure 1

Modeling of BKCa–CaV complex with 1:1 stoichiometry. (A) Schematic representation of the BKCa model, where X and Y indicate the closed and open states of the BKCa channel. (B) Fit to the experimental data [10] of the 2-state model of panel A: steady-state BKCa open probabilities versus voltage at different Ca2+ concentrations (different markers for different Ca2+ levels as indicated in the legend) and the corresponding fit obtained by the model (dashed lines) (upper left plot); time constants (in ms) versus voltage data at given Ca2+ concentrations (see legend in each plot), and the corresponding model fits (dashed lines). (C) Scheme indicating the six states of the devised MC model for BKCa–CaV complex with 1:1 stoichiometry. , and correspond to the closed state for the BKCa channel (X) coupled with the closed (C), open (O) and inactivated (B) states for the CaV, respectively, and , and correspond to the open state for the BKCa channel (Y) coupled with the closed (C), open (O) and inactivated (B) states for the CaV, respectively. The green box indicates states with non-inactivated CaV, whereas the blue box highlights states with inactivated CaV. (D) Simulated open probabilities, in response to three different voltage steps, from −80 mV to −20 (left), 0 (middle) and 20 mV (right), for BKCa channels controlled by CaVs in complexes with 1:1 stoichiometry, obtained from different models: the original 70-state Markov chain model (gray; [10]); the 6-state Markov chain model (shown in panel C—black); the ODE model corresponding to the 6-state model (blue; Equation (12)); the simplified Hodgkin–Huxley-type model, , where is given by Equations (13)–(15) (dashed red) and the corresponding model assuming instantaneous activation with defined by Equation (13) with the approximations of Equations (16) and (17) (dash-dotted green). (E–H) Simulated BKCa currents in response to the different voltage steps (from −80 mV to from −40 mV to 40 mV in 20 mV increments for 20 ms and then back to −80 mV) obtained from the different models as in D. In D–H the average of one-thousand Monte Carlo simulations for each MC models is shown.

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 BKCa channel with a 7-state MC model for the CaV channel (see Methods and Figure S1), resulting in a 70-state MC model. The simulated open probabilities were very similar for both the models in response to different voltage steps, from −80 mV to −20, 0 and 20 mV (compare the dotted gray curve (70-state MC model) with the dash-dotted black one (6-state MC model)). Moreover, the open probability obtained by the 6-state ODE model (the solid blue curve) approximates well, the average fraction of open channels calculated from Monte Carlo simulations of the corresponding 6-state MC model, and the open-probability expression of Equation (2) (dashed red curve). Additionally, the further simplification assuming instantaneous CaV activation (dash-dotted green curve; 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 BKCa models reproduce, satisfactorily, the whole-cell BKCa currents for each step in voltage ( defined by Equation (1) with , where , with pF [31] the single-channel conductance and the number of BKCa 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 BKCa channel with (Equation (2); red plots in Figure 1G), and the corresponding model assuming instantaneous activation of the CaV currents ( given by Equation (13) with the approximations of Equations (16) and (17); green plots in Figure 1H) also work very well.

2.1.2. BKCa–CaV Complex with 1:n Stoichiometry

In the case of more than one CaV per BKCa channel, the linear buffer approximation is used to compute the Ca2+ 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 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 BKCa channel closes rapidly when all CaVs in the complex are inactivated. Hence, for the case of 1:n BKCa–CaV stoichiometry, we split the system according to the number k of non-inactivated CaVs: the BKCa activation can be described on a fast time scale by the Markov chain model of Figure 2A with states. As for the case of 1:1 stoichiometry, the dynamics of the BKCa open probability in complexes with k non-inactivated CaVs can be approximated by a single ODE (see defined by (26)). Then, the open probability of the BKCa channel coupled with n CaVs, , can be estimated by taking into account that the probability of k non-inactivated CaVs being present in a complex with n CaVs, . Then, can be expressed as

Figure 2

Multiple CaVs per BKCa–CaV complex: from local to whole-cell behavior. (A) Markov chain (MC) model for complexes with k non-inactivated CaVs. (B) Steady-state BKCa activation functions, (upper plot), and time-constants, (lower), for BKCa channels in complexes with (cyan), 2 (green) or 4 (red) CaVs, given from 1-ODE by (26) (solid) and from the simplification assuming by Equations (28) and (29) (dashed). The gray dashed curve shows the CaV activation function , for comparison. (C) Simulated BKCa open probabilities, , for the case of 1:4 BKCa–CaV stoichiometry, in response to a voltage step from –80 mV to 0 mV, obtained from the MC model of states (black), from the ODEs describing the states in panel A coupled to CaV inactivation ( defined by Equation (3) with given by Equation (18)—solid blue), from the reduced ODE model considering CaV activation kinetics ( by Equation (3) with given from 1-ODE by (26)— dashed red), and from the simplification assuming (Equation (26) for with Equations (28) and (29)—dash-dotted green). (D–F) Simulated whole-cell BKCa currents, , defined by Equation (1) in response to the different voltage steps (as in Figure 1E–H) obtained from the different ODE models used for characterizing , as in (C): given by Equation (18) (blue curves in panel D); given by Equation (26) (red curves in panel E); given by Equation (26) with assumptions of Equations (28) and (29) (green curves in panel F). Panels C–F show the average of one-thousand Monte Carlo simulations for the MC model of states (black). (G–I) Simulated electrical activity in a model of lactotrophs [26] with 1:n BKCa–CaV complexes, and with (G), (H) or (I). is described by Equation (4), where is modeled by Equation (18) (upper plots), Equation (26) (middle) and Equation (26) with assumptions of Equations (28) and (29) (lower).

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 in Equation (3) (blue curve, by Equation (18); red curve, by Equation (26); green curve, by Equation (26) with instantaneous CaV activation defined by Equations (28) and (29)) effectively approximate the Monte Carlo simulations of the full Markov Chain model with states (black curve).

The whole-cell BKCa current defined by Equation (1), with given by (3), involves the reduced model n ODEs (Equation 26) for the activation variables with . Note that if the n CaVs do not inactivate, Equation (1) reduces to which involves one ODE (Equation 26) for .

Figure 2D–F shows the simulated whole-cell BKCa currents for the case of 1:4 BKCa–CaV stoichiometry (defined by Equation (1), where , with pF [31] as in the previous section, while 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 ) in response to different voltage steps and compared with the simulated BKCa currents obtained from MC model (). In all the cases, the different ODE models for characterizing (Equation (18) (blue curves), Equation (26) (red curves) and Equation (26) with Equations (28) and (29) (green curves)) approximate very well the Monte Carlo simulations (black curve) obtained from the MC model.

2.1.3. Concise Whole-Cell Modeling Respecting Local Control

The Hodgkin–Huxley-type model of the BKCa current defined by Equation (1), that allows us to take into account the local interactions in BKCa–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 BKCa–CaV complexes with 1:1 stoichiometry, spiking electrical activity is observed as shown in Figure 2G, since insufficient BKCa current is generated. The different plots in Figure 2G (upper, middle and lower) correspond to different ODE models used for computing the open probability of the single complex defined in Equation (1) for : The 4-state ODE model ( and are not considered since CaV does not inactivate) with defined by Equation (12) and with (upper). The single ODE model for defined by Equation (13) with CaV kinetics leading to (Equation (2) since ) (middle), with simplifications for of Equation (13) assuming instantaneous CaV activation by Equations (16) and (17) (lower).

In contrast, in the case of BKCa–CaV complexes with 1:n stoichiometry with , plateau bursting appears with the number of small oscillations per burst depending on the number of CaVs per BKCa–CaV complex, as shown in Figure 2H,I for and , respectively. The different plots in Figure 2H,I (upper, middle and lower) correspond to the different ODE models used for computing , and hence, the resulting of Equation (4)—, given by Equation (18) (upper blue plots) by solving the complete ODE model of equations; , given by single ODE of Equation (26) (middle curves); , given by Equation (26) with Equations (28) and (29) by assuming instantaneous CaV activation, (lower curves).

Although the quantitative behavior is independent of the approximation for , 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 BKCa 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.

2.2. 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 Ca2+-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 Ca2+ ions bound to the Ca2+ sensor on the granule (states , with , representing the number of Ca2+ ions bound to the granule sensor) before fusing with the membrane and releasing its hormone content (state Y) (see Methods for more details). The Ca2+ concentration at the granule sensor, given by Equation (7) with being the distance from the CaV pore to the Ca2+ sensor on the granule, drives the exocytosis MC model, allowing the granule to modify its state through the Markov chain from to Y and undergo exocytosis.

Figure 3

Models for exocytosis of single granule and granule-CaV complex. (A) The 5-state model for exocytosis of a single granule adjacent to the plasma membrane, where corresponds to the state with no bound Ca2+ ions, to one, to two and to three. (B) Model of the granule-CaV complex where the granule is described by the model shown in panel A and the CaV by the model shown in Figure S1A. (C) Model of the granule coupled with n non-inactivating CaVs, where the granule is described by the model shown in panel A and the CaV by a single-channel gating with two states, closed, C, and open, O.

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 -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 -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 Ca2+ concentration at the granule sensor needed for triggering exocytosis, is equal to a basal level () when the CaV is closed or inactivated, and is equal to 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 Ca2+ levels at the granule sensor when more than one CaV is open, as performed for BKCa–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 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 Ca2+ dynamics, and hence, the CaVs states, are driven by the electrophysiological behavior, which depends on the local interactions between CaVs and BKCa 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 BKCa channels, determine the granule release by varying the number and the distance of CaVs coupled with the granules.

Spiking versus Bursting on Evoking Exocytosis

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 , 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 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 and 2, respectively), while this difference in the rate decreases when the number of CaVs increases (see the cyan and black lines for and 8, respectively). When the CaV are far from the granule (see Figure 4C,D for 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 Ca2+ entering during bursting, which becomes negligible when the number of CaVs is high and close to the granule. In this scenario, the Ca2+ concentrations at the exocytotic machinery are in saturation regimes for both the electrical patterns.

Figure 4

Probability of exocytosis for single granule coupled with n CaVs driven by spiking or bursting. (A–D). Probability of exocytosis versus time for the granule at fixed distances from n CaVs ( (blue curves), (magenta curves), (cyan curves) and (black curves), driven by spiking (solid lines) or bursting (dash-dotted lines)): nm (A); nm (B); nm (C); nm (D). (E–H) Probability of exocytosis versus total Ca2+ influx, (computed by integrating Ca2+ current over the time), for panels (A–D) respectively. Each panel shows the average of one-thousand Monte Carlo simulations, which were performed for the MC exocytosis model.

Therefore, in order to evaluate the exocytotic machinery performance with the same Ca2+ entry, we analyzed the exocytosis probability with respect to the total Ca2+ influx, . Figure 4E–H shows the granule release probability versus 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 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 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 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 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.

Figure 5

Fitting to the simulated exocytosis probabilities. (A) Trend (in logarithmic scale) of q parameter of function defined by Equation (5) with respect to the distance of the granule from n CaVs ( (blue curves), (magenta curves), (cyan curves) and (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 for the granule at from four CaVs (upper plot) and for the granule at from two CaVs. The simulated responses are the averages of one-thousand Monte Carlo simulations.

3. Discussion

In this paper, we show the role of mathematical modeling as an important tool for investigating excitable cells with focus on ion Ca2+ channels and their local interactions with BKCa 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 BKCa–CaV complexes with 1:n stoichiometry, we used a stochastic model based on Markov chain theory (of 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 BKCa 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 BKCa 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 BKCa channel determines a left shift of the BKCa 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 BKCa channel, providing theoretical insight into stochastic simulation results. Then, the concise Hodgkin–Huxley-type model of whole-cell BKCa currents, that allows taking into account local control by CaVs in BKCa–CaV complexes, can be inserted into models of cellular electrical activity, showing, in particular, how different BKCa–CaV stoichiometries cause different electrophysiological responses, including continuous spiking and bursting in pituitary cells. In [13], we also showed how BKCa–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 BKCa channels with different CaV types affects electrical activity differently in human pancreatic beta-cells, and further insight into the control by CaVs of BKCa 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 BKCa 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 ), 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 Ca2+ ions bound to the Ca2+ 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 Ca2+ 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, Ca2+ concentration. Indeed, as experimentally observed [39], global Ca2+ concentration (i.e., the bulk cytosolic Ca2+ 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 Ca2+ levels, in particular, for calculating calcium concentrations at ion channel BKCa–CaV complexes and at granules, as performed in previous works: Cox [10] found that, for computing Ca2+ levels at a BKCa 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 BKCa channel coupled with n CaVs. In order to get a deeper description of the Ca2+ levels, we can exploit compartmental modeling, as performed in our previous work [40], where we characterized the intracellular Ca2+ 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 (Ca2+ nanodomains, sub-membrane compartment, bulk cytosol and endoplasmic reticulum) involved in exocytosis.

4. Methods

4.1. Electrophysiological Model of Pituitary Cells

The original model of pituitary cells [26] included a single Ca current (), a delayed-rectifier K current (), a Ca-gated SK current () and a leak current (), in addition to the BKCa current which was modeled as a purely voltage-dependent current, neglecting Ca2+ dependency. We substituted this simplified representation of the BKCa current with our concise BKCa 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 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 BKCa currents (Equation (4)) that take into account local control in BKCa–CaV complexes with different stoichiometries, we first devised a model of a single-channel gating for the BKCa 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.

4.1.1. BKCa Channel Modeling

For describing the BKCa 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, , and , describing the transition from the open to close state and from close to open state, respectively. By assuming that, for fixed Ca2+ concentration (), BKCa activity is described by a Boltzmann function [9,27,41] and that the slope parameter of the Boltzmann function is independent of , a reasonable assumption for Ca2+ concentrations above 1 M [10,41,42], as expected in BKCa–CaV complexes [3,10], in [13], we showed that and can be expressed as a product of voltage and Ca2+-dependent terms ( and are the voltage dependent rates, and the Ca2+-dependent rates; then, , ; see [13] for mathematical expressions). Figure 1B shows the fitting to the experimental data [10], consisting of BKCa open probabilities and time constants as functions of voltage, at different Ca2+ concentrations, by using for and —the optimal parameters estimated in [13].

We also reproduced the dynamics of a more complex model for describing the BKCa 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 Ca2+ 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 Ca2+ ions and be open or closed. As shown in [10], the 10-state model is able to reproduce the characteristics and dynamics of the BKCa channel by fitting BKCa open probabilities and time constants as functions of voltage, at different Ca2+ 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 BKCa contains at least three Ca2+ binding sites (a low and two high-affinity calcium binding sites) [43,44,45] and four voltage sensor per subunit [46].

4.1.2. 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 Ca2+ channel opening rate and closing rate, respectively, and have the forms as in [47]; is a constant reverse reactivation rate and represents the Ca2+-dependent rate for channel inactivation, which is determined by the Ca2+ concentration at the Ca2+ sensor, , 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 nm) and is the single-channel Ca2+ current with the single-channel conductance and the Ca2+ 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 Ca2+ binding kinetics that are rapid relative to Ca2+ 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 Ca2+ levels at BKCa 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 Ca2+ current , 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, , where and the following equation for inactivation

Therefore, assuming instantaneous activation, , 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 Ca2+ channels not inactivated, , where o and c represent the state variables of 3-state system of Figure S1A. For 1-state model, , 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 () 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 , 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 . 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 (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 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 is the number of CaV channels, each one characterized by a single channel conductance ( pS as in [53] for CaV2.1). and are, respectively, the activation and inactivation variables of the channel: for the 3-state ODE model, , where o is the state variable corresponding to the open state of the channel; for the simplified 1-state ODE model, assuming instantaneous activation of the CaV currents, by Equation (9) and by solving the single ODE model defined by Equation (10); for the 7-state MC model, , where represents the probability to be in state O for the MC of Figure S1B. As for the CaV open probabilities, the 3-state ODE model current approximates the 7-state Markov chain model current very well for each step voltage (Figure S1F,G). Additionally, the further simplification for the 3-state ODE model assuming instantaneous activation of the CaV currents (green plots in Figure S1H) provides a good approximation of the Monte Carlo simulations.

4.1.3. BKCa Channel Open Probability for BKCa–CaV Complex with 1:1 Stoichiometry

First, we coupled the 2-state model devised for the BKCa channel with the 3-state model for CaV, resulting in a 6-state MC model for achieving a concise model of whole-cell BK current (Equation (1)) that respects the local effects of BKCa–CaV coupling. Figure 1C shows a cartoon of the model of the devised BKCa–CaV complex with 1:1 stoichiometry: , and correspond to the closed state for the BKCa channel (X) coupled with the closed (C), open (O) and inactivated (B) states for the CaV, respectively, and , and correspond to the open state for the BKCa channel (Y) coupled with the closed (C), open (O) and inactivated (B) states for the CaV, respectively. The parameters and ( and ) are functions of the calcium concentration at the BKCa channel when the associated CaV is closed, , or inactivated, (i.e., ), and to the calcium concentration at the BKCa channel when the CaV is open, . In particular, the Ca2+ levels sensed by the BKCa channel are assumed to reach steady-state immediately after CaV closure or opening [10]; in this case, is set equal to 0.2 M (background Ca2+ concentration) and is given by Equation (7) with nm representing the distance between CaV and BKCa channels [3,10]. At mV, M. Note that , since the background Ca2+ concentration is much below the levels needed for BKCa activation at physiological voltages; i.e., the probability of BKCa opening when the CaV is closed is practically zero: this approximation is supported by the fact that Ca2+ influx via CaVs is needed to open BKCa channels [54], and that the sub-membrane Ca2+ concentration of some hundreds of nM that a BKCa in a complex without open CaVs would sense, is too low to activate BKCa channels at physiological voltages [3,10].

Next, we evaluated the dynamics of the deterministic ODE model corresponding to the 6-state MC model. The model is described by the state-variables , with , representing the probabilities of the complex to be in one of the six states of the model. Then, the BKCa open probability, , is given as

The 6-state ODE model characterizing (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 () 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 BKCa channel activation in BKCa–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 BKCa activation, denoted with (): with steady-state and time constant given by

The CaV activation variable is routinely characterized in patch clamp experiments and included in models of electrical activity via the time-constant, , and the steady-state activation function, (see Equation (9)). From these quantities, and can be calculated. Note that Equation (15) makes it explicit how inherits properties of the associated Ca2+ channel type, as has been found experimentally [55,56].

In many whole-cell models (e.g., [26,28,29,30]), the Ca2+ 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 BKCa–CaV model defined by Equations (13)–(15) in the limits . In this case Equations (14) and (15) become which are completely defined from BKCa kinetics and .

By coupling the BKCa channel activation defined by Equation (13) with CaV inactivation function h, we achieved the BKCa open probability of Equation (2) for BKCa–CaV complex with 1:1 stoichiometry.

4.1.4. BKCa Channel Open Probability for BKCa–CaV Complex with 1:n Stoichiometry

BKCa channels can form ion complexes with more than one CaV with a stoichiometry of 1–4 CaVs channels per BKCa channel, as experimentally observed [3,4,23]. Here, we introduce a concise but mechanistically correct model of single BKCa–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 BKCa channel [3,10,57]: In this case, the linear buffer approximation is used to compute the Ca2+ 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 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 BKCa channel closes rapidly when all CaVs in the complex are inactivated. Hence, for the case of 1:n BKCa–CaV stoichiometry, we split the system according to the number k of non-inactivated CaVs: the BKCa activation can be described on fast time scale by the Markov chain model of Figure 2A with states, where and correspond to the states with closed and i open CaVs, with , coupled with the closed (X) and open (Y) BKCa channels, respectively. As for the case of 1:1 stoichiometry, the dynamics of the BKCa open probability can be approximated by a single ODE. In particular, we start from the complete ODE system describing the BKCa channel coupled with k non-inactivated CaVs of Figure 2A: the model consists of ODEs characterizing the state variables and , corresponding to the probability of having closed and i open CaVs, with , coupled with the closed (X) and open (Y) BKCa channels, respectively. The activation of the BKCa surrounded by k non-inactivated CaVs, denoted with , is then

By taking into account that and renaming the state variables as follows the ODE system is reduced from to equations.

Moreover, assuming the quasi-steady state approximation for , with , the ODE system of equations is reduced to the following single ODE, describing the dynamics of the BKCa activation with k non-inactivated CaVs: where and are explicit functions of V, directly or via the local Ca2+ concentration (see Equation (S36) in the Supporting Material of [13]).

By assuming instantaneous activation of CaVs, the model defined by Equation (26) can be further simplified. Indeed, in this case (i.e., ) the vertical transitions in Figure 2A are in quasi-equilibrium, and then and 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 , we obtain Equation (3) to describe the BKCa channel open probability for BKCa–CaV complex with 1:n stoichiometry.

4.2. 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 Ca2+ ions bound to the Ca2+ sensor on the granule, likely synaptotagmin [58]: in with no bound Ca2+ ions, in with one, in with two or in with three bound ions. Once it is in , 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 and its transition rate or generator matrix is given by where represents the Ca2+ binding rate, with , the Ca2+ concentration at the granule sensor, given by Equation (7) with being the distance from the CaV to the Ca2+ sensor on the granule. In this paper, the distance from the CaV to the granule means the distance from the CaV to the Ca2+ 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 Ca2+ sensor molecules, which was, therefore, included in the binding parameter . The parameter was the unbinding rate, and u was the fusion rate. The values for , and u were equal to those used in [22]. Note that, since 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 and can be described by an auxiliary variable , 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 Ca2+ 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 -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 -state MC model, as shown in Figure 3C.