Discrete-State Stochastic Models of Calcium-Regulated Calcium Influx and Subspace Dynamics Are Not Well-Approximated by ODEs That Neglect Concentration Fluctuations

Abstract

Cardiac myocyte calcium signaling is often modeled using deterministic ordinary differential equations (ODEs) and mass-action kinetics. However, spatially restricted “domains” associated with calcium influx are small enough (e.g., 10⁻¹⁷ liters) that local signaling may involve 1–100 calcium ions. Is it appropriate to model the dynamics of subspace calcium using deterministic ODEs or, alternatively, do we require stochastic descriptions that account for the fundamentally discrete nature of these local calcium signals? To address this question, we constructed a minimal Markov model of a calcium-regulated calcium channel and associated subspace. We compared the expected value of fluctuating subspace calcium concentration (a result that accounts for the small subspace volume) with the corresponding deterministic model (an approximation that assumes large system size). When subspace calcium did not regulate calcium influx, the deterministic and stochastic descriptions agreed. However, when calcium binding altered channel activity in the model, the continuous deterministic description often deviated significantly from the discrete stochastic model, unless the subspace volume is unrealistically large and/or the kinetics of the calcium binding are sufficiently fast. This principle was also demonstrated using a physiologically realistic model of calmodulin regulation of L-type calcium channels introduced by Yue and coworkers.

1. Introduction

Concentration changes of physiological ions and other chemical species (such as kinases, phosphatases, and various modulators of cellular activity) influence and regulate cellular responses [1]. These dynamics are often modeled using systems of deterministic ordinary differential equations (ODEs) that assume chemical species concentrations are nonnegative real-valued quantities (i.e., the state-space is continuous). In such descriptions, the rate of change of the concentration of each species is usually specified under the assumption of mass-action kinetics, that is, the rate of a reaction is proportional to the product of reactant concentrations. However, under physiological conditions the concentrations of chemical species are often quite low and, in some cases, restricted subspaces in which these species are contained are very small. For example, L-type calcium channels in cardiac myocytes are typically clustered in small “diadic subspaces” that have a volume of ~10⁻¹⁷ liters, with approximately 20,000 diadic subspaces per cell [2, 3]. Resting calcium concentration in the diad is typically 0.1 micromolar, a value that corresponds to an average of 0.6 calcium ions per subspace [4]. Because only whole numbers of calcium ions can be present in a subspace at any given time, the question arises: is it appropriate to use deterministic ODEs to model subspace calcium dynamics?

Previous studies have compared discrete-state (stochastic) and continuous-state (deterministic) models in the analysis of biological and chemical systems, including models of biochemical networks, enzyme kinetics, and population dynamics [5–21]. These studies have shown that in the “large-system limit" (i.e., a large “copy number” of each chemical species), the solution of discrete and continuous models are equivalent [12]. However, for a small system, concentration values obtained from a continuous deterministic model (an approximation that neglects concentration fluctuations) can significantly deviate from the expected value obtained from the discrete stochastic model. When chemical reactions are higher than first order, there is no guarantee that the deterministic mass-action formulation will agree with, or be a good approximation to, the expected value of species concentrations obtained from a chemical master equation that accounts for discrete system states and concentration fluctuations [5]. An excellent study by Goutsias discusses the relationship between the discrete and continuous formulations for general biochemical systems [22] (for theoretical context, see [23]).

Because of recent interest in the physiological relevance of spatially localized control of voltage- and calcium-regulated calcium influx and sarcoplasmic reticulum calcium release in cardiac myocytes [24–26], we sought to determine precisely when the conventional deterministic formulation of these processes are a valid approximation. When is it appropriate to model the dynamics of subspace calcium using deterministic ODEs? When does one require a stochastic description that accounts for the fundamentally discrete nature of calcium-regulated calcium influx?

To answer this question, we constructed and analyzed a minimal Markov model of a calcium-regulated calcium channel and associated subspace. We compared the expected steady-state subspace calcium concentration in this stochastic model (a result that accounts for the small subspace volume) with the result obtained using the corresponding deterministic ODE model (an approximation that assumes large system size). Section 2.1 introduces our model formulation and shows the agreement between deterministic and stochastic descriptions when subspace calcium does not regulate calcium influx. However, when calcium binding regulates channel activity (through either activation or inactivation), the deterministic and stochastic descriptions often disagree (Sections 2.2 and 2.3). In general, the effect of concentration fluctuations in a spatially restricted calcium domain with a calcium-regulated calcium influx pathway (e.g., a stochastically gating L-type calcium channel) is only well-approximated by the deterministic description when the subspace volume is sufficiently (unphysiologically) large or the kinetics of calcium binding to the calcium-regulated channel are sufficiently fast. This principle was also demonstrated using a physiologically realistic model of calmodulin regulation of L-type calcium channels produced by Yue and coworkers (Section 2.4).

2. Methods and Results

2.1. Calcium Influx and Subspace Calcium Concentration Fluctuations

We begin with the case of a single calcium channel that is associated with a spatially restricted subspace but not regulated by subspace calcium (Figure 1 and Section 2.1). The description of the model in the absence of calcium regulation simplifies the initial presentation of the model and allows us to illustrate general properties of subspace calcium concentration fluctuations. Subsequently, we present a more complete model formulation that includes calcium-regulated calcium influx (Figure 3 and Section 2.2). For simplicity, we neglect the presence of endogenous calcium binding proteins and assume a constant flux of calcium, denoted by 𝛼, into the subspace. The subspace calcium concentration is passively coupled with relaxation rate 𝛽=0.01 ms⁻¹ to the constant bulk concentration of 𝑐_∞=0.1 μM [27]. These assumptions lead to the following deterministic model of subspace calcium dynamics: 𝑑𝑐𝑑𝑡=𝛼−𝛽𝑐−𝑐_∞, (1) where the influx rate 𝛼 has units of concentration per time (e.g., μM/ms) and the calcium concentration 𝑐 is a continuous real-valued quantity.

Figure 1 

Diagram of the components and fluxes in a minimal subspace model. Calcium influx 𝛼 (in units of μM/s) leads to increased calcium concentration 𝑐 (units of μM) in a diadic subspace of volume 𝑣 (liters). Subspace calcium moves to the bulk passively via diffusion at rate 𝛽 (given by 0.01 ms−1). Bulk calcium at the concentration 𝑐∞=0.1 μM returns to the subspace at the same rate. The equilibration time of subspace calcium is 𝜏=1/𝛽=100ms [27].

2.1.1. Stochastic Model

In the corresponding stochastic description of calcium influx into a diadic subspace, the state variable is the number of calcium ions in the subspace (a discrete quantity that we will denote by C∈{0,1,…,∞}, where the caret (hat) indicates that C is a dimensionless number of molecules rather than concentration, and the capitalization indicates a random variable. The fluctuating subspace calcium concentration (also a random variable, denoted by C) depends on both C and the subspace volume (𝑣), that is, CC=𝑣. (2) Using this relationship, it is straightforward to derive the transition rates between the discrete states of the stochastic model that are consistent with (1). The resulting state-transition diagram for the stochastic model is 0⇌^𝛼_𝛽1⇌^𝛼_2𝛽2⋯𝑛−1⇌^𝛼_𝑛𝛽𝑛⇌^𝛼_()𝛽𝑛+1𝑛+1⋯, (3) where the index that labels states, 𝑛∈{0,1,…,∞}, ranges over all possible numbers of calcium ions in the subspace and the constant 𝛼 is proportional to the subspace volume, that is, 𝛼=𝑣𝛼+𝛽𝑐_∞. (4)

2.1.2. Master Equation and Steady-State Probability Distribution

If we write 𝑝_𝑛(𝑡)=Pr{C(𝑡)=𝑛}, the equations for the dynamics of the probability of each state in (3), that is, the chemical master equation for the number of calcium ions in the subspace, is given by 𝑑𝑝₀𝑑𝑡=−𝛼𝑝₀+𝛽𝑝₁,𝑑𝑝_𝑛𝑑𝑡=−𝛼+𝑛𝛽𝑝_𝑛+𝛼𝑝_𝑛−1+(𝑛+1)𝛽𝑝_𝑛+1,𝑛=1,2,…. (5) Note that the correspondence between the rate constants in the deterministic (1) and stochastic (3)–(5) models is established by substituting 𝑐=𝑐/̂𝑣 in (1) to find the rate of change of the number of calcium ions in the deterministic model, that is, 𝑑𝑐̂=𝑑𝑡𝛼−𝛽𝑐.̂ (6) This equation indicates that 𝑐̂ increases at rate 𝛼 (due to influx and diffusion from the bulk), a value that is independent of the number of calcium ions in the subspace. At the same time, 𝑐̂ decreases at rate 𝛽𝑐̂, a value that is proportional to 𝑐̂ because each ion has an opportunity to diffuse into the bulk. Consequently, the transition rates leading out of state C=𝑛 in the stochastic model are given by 𝛼 for the C=𝑛 to 𝑛+1 transitions and 𝛽𝑛 for the C=𝑛 to 𝑛−1 transitions.

To find the steady-state probability distribution of C, we set the left hand sides of (5) to zero to obtain 𝑛𝛽𝑝_𝑛=𝛼𝑝_𝑛−1,𝑛=1,2,… (7) from which it follows that {𝑝_𝑛} is a Poisson distribution with parameter 𝜆=𝛼/𝛽, that is, 𝑝_𝑛=𝑒^−𝜆𝜆^𝑛.𝑛! (8)

2.1.3. Analysis of Concentration Fluctuations

To see how the subspace calcium concentration fluctuations predicted by this minimal model depend on the parameters 𝛼, 𝛽, 𝑐_∞, and 𝑣, recall that the mean and variance of the Poisson distribution (8) is equal to the parameter 𝜆 and, consequently, the steady-state expected number of calcium ions in the subspace is given by EC=^∞_𝑛=0𝑛𝑝_𝑛=𝜆=𝛼𝛽𝛼=𝑣𝛽+𝑐_∞=𝑣𝑐_∗, (9) where the last equality defines 𝑐_∗ as follows: 𝑐_∗=𝛼𝛽+𝑐_∞. (10) Using (9) and the fact that C=C/𝑣 implies E[C]=E[C]/𝑣, we can identify 𝑐_∗ as the expected subspace calcium concentration: E[C]=𝑐_∗. (11) Similarly, the steady-state variance of the number of calcium ions in the subspace is VarC=^∞_𝑛=0(𝑛−E[C])²𝑝_𝑛=𝑣𝑐_∗, (12) and Var[C]=Var[C]/𝑣² implies that the variance of the subspace calcium concentration is [C]=𝑐Var^∗𝑣. (13) Note that the coefficient of variation of C and C are identical and inversely proportional to subspace volume, that is, CV[C]=(Var[C])^1/2/E[C]=1/√𝑣𝑐_∗ and, similarly, CV[C]=Var[C]E[C]=VarC/𝑣²EC=1/𝑣√𝑣𝑐_∗. (14) This is a well-known principle from statistical physics: fluctuation amplitudes scale with the reciprocal of the square root of system size (the subspace volume 𝑣).

Figure 2 illustrates fluctuation amplitudes in the minimal subspace model by plotting the steady-state probability distribution of C and C (left and right columns, resp.). In the first row, using subspace volume of 𝑣=𝑣₀=10⁻¹⁷ liters and influx rate of 𝛼=0.049 μM/ms, the expected calcium concentration is E[C]=𝑐_∗=𝛼/𝛽+𝑐_∞=5 μM, and the expected number of subspace calcium ions is E[C]=𝑣₀𝑐_∗=30. In both cases the coefficient of variation is 1/√30=0.18 (the spread of the distributions as illustrated is due to the different 𝑥-axis scales). The following rows of Figure 2 show that in a subspace three or ten times larger (𝑣=3𝑣₀ or 10𝑣₀), the coefficient of variation drops to 0.11 and 0.058, respectively, when the calcium influx rate is scaled to result in the same expected calcium concentration (𝑐_∗ fixed, see (14)). As might be expected, concentration fluctuations in the stochastic model are more pronounced for small volumes and become negligible for large volumes, because CV[C]=1/√𝑣𝑐_∗→0 as 𝑣→∞ for fixed 𝑐_∗.

Figure 2 

Steady-state probability distribution of the number of calcium ions (C, left column) and subspace calcium concentration (C, right column) for subspace volume of 𝑣0=10−17 liters and subspaces that are 3 and 10 times larger. Parameters: 𝛼=0.049 μM/ms, 𝛽=0.01 ms−1, 𝑐∞=0.1 μM; the steady-state expected subspace calcium concentration is E[C]=𝑐∗=5 μM.

Figure 3 

Diagram of the components and fluxes in a subspace model that includes calcium-regulated calcium influx. A single calcium channel (with two calcium binding sites) is associated with a subspace of volume 𝑣. The calcium influx rate is 𝛼0 and 𝛼1 when calcium is unbound and bound, respectively, and the transition rates between these states are 𝑘+𝑐2 and 𝑘−, where 𝑐 is the subspace calcium concentration. Subspace calcium is passively coupled at rate 𝛽 to the bulk cytosol with constant concentration 𝑐∞.

Most importantly, the deterministic and stochastic descriptions of this minimal subspace model agree in the following sense: the expected value of the fluctuating calcium concentration in the stochastic model E[C]=𝑐_∗=𝛼/𝛽+𝑐_∞ is equal to the steady-state of the deterministic ODE that neglects concentration fluctuations (found by setting the left hand side of (1) to zero). Readers familiar with fluctuations in biochemical models will understand that this agreement is a consequence of the fact that the minimal subspace model involves three elementary reactions, all of which are zeroth or first order (see arrows in Figure 1).

2.1.4. Moment Calculation

The numerical results presented above can be obtained analytically by considering the dynamics of the moments of the number of calcium ions in the subspace, defined as 𝜇_𝑞=^∞_𝑛=0𝑛^𝑞𝑝_𝑛. (15) By conservation of probability, the zeroth moment 𝜇₀=1 and the first moment is the expected number of calcium ions in the subspace (9), 𝜇₁=EC. (16) The second moment 𝜇₂ is related to the variance of the number of calcium ions via VarC=𝜇₂−𝜇₁². (17) By differentiating (15) with respect to time and substituting for the time derivatives using the master equation (5), it can be shown that the zeroth moment is constant (𝑑𝜇₀/𝑑𝑡=0) and, furthermore, 𝑑𝜇₁=𝑑𝑡𝛼−𝛽𝜇₁,𝑑𝜇₂=𝑑𝑡𝛼+2𝛼+𝛽𝜇₁−2𝛽𝜇₂, (18) where we have used 𝜇₀=1. Setting the left hand sides of these equations to zero, we see that steady-state first and second moment are 𝜇₁=𝛼/𝛽 and 𝜇₂=𝛼/𝛽+(𝛼/𝛽)², consistent with (11) and (17).

2.2. Stochastic Subspace Model with Calcium Regulation

This section augments the subspace model presented above to include calcium regulation of a calcium channel (see Figure 3). We assume that calcium binding instantaneously modifies the conductance of the channel, that is, the rate of calcium influx into the domain is 𝛼₀ when the channel is calcium-free and 𝛼₁ when the channel is calcium-bound. We further assume the channel has two binding sites for calcium and, for simplicity, approximate rapid sequential binding of calcium ions with instantaneous binding. Thus, the transitions between the two distinct states of the subspace (the so-called “stochastic functional unit" or “calcium release unit") occur at rates 𝑘⁺𝑐² and 𝑘⁻, respectively, (Figure 3, curved arrows). Note that the rate constant 𝑘⁻ has units of ms⁻¹, 𝑘⁺ has units of μM⁻² ms⁻¹, and the dissociation constant for calcium binding, denoted by 𝜅, has units of μM and is given by 𝜅²=𝑘⁻/𝑘⁺.

2.2.1. Stochastic Model

Let us denote the states of the stochastic system by (𝑛,0) and (𝑛,1), where 𝑛∈{0,1,…,∞} and the second element of the ordered pairs, either 0 or 1, indicates calcium-free and bound channel, respectively. With a little thought we can sketch the following state-transition diagram for the stochastic subspace model with calcium influx,(0,0)⇌^𝛼0_𝛽(1,0)⇌^𝛼0_2𝛽(2,0)⇌^𝛼0_3𝛽(3,0)⇌^𝛼0_4𝛽()⋯24,0𝑘⁺⥯𝑘⁻6𝑘⁺⥯𝑘⁻12𝑘⁺⥯𝑘⁻(0,1)⇌^𝛼1_𝛽(1,1)⇌^𝛼1_2𝛽(2,1)⋯, (19)

where 𝑘⁺=𝑘⁺/𝑣². The rate of calcium binding to the channel in the stochastic model is inversely proportional to the square of the volume, because of the concentration dependence of the association reaction (𝑘⁺𝑐²=𝑘⁺𝑐²̂/𝑣²). The downward transitions between states (𝑛,0) and (𝑛−2,1) include the combinatorial coefficient, 𝑛(𝑛−1), double the number of ways that two indistinguishable calcium ions can be chosen from the 𝑛 ions in the subspace. This factor of two is required so that the microscopic propensity 𝑘⁺ agrees with the macroscopic rate 𝑘⁺𝑐² for large 𝑛 and 𝑣 with 𝑐=𝑛/𝑣 fixed, that is, 𝑛(𝑛−1)𝑘⁺()𝑘=𝑛𝑛−1⁺𝑣²=𝑘⁺𝑐²−𝑐𝑣, (20) an expression that approaches 𝑘⁺𝑐² as 𝑣→∞ [28].

2.2.2. Master Equation

Let us write 𝑝⁰_𝑛(𝑡) to indicate the probability that at time 𝑡 the channel is calcium-free and C=𝑛. Similarly, 𝑝¹_𝑛(𝑡) is the probability that C(𝑡)=𝑛 and the channel is calcium bound. Reading off the transition rates from the state-transition diagram (19), we write the following master equation for the calcium-regulated channel and subspace: 𝑑𝑝⁰_𝑛𝑑𝑡=−𝛼₀+𝑛𝛽+𝑛(𝑛−1)𝑘⁺𝑝⁰_𝑛+𝛼₀𝑝⁰_𝑛−1+(𝑛+1)𝛽𝑝⁰_𝑛+1+𝑘⁻𝑝¹_𝑛−2,𝑑𝑝¹_𝑛𝑑𝑡=−𝛼₁+𝑛𝛽+𝑘⁻𝑝¹_𝑛+𝛼₁𝑝¹_𝑛−1+(𝑛+1)𝛽𝑝¹_𝑛+1+(𝑛+2)(𝑛+1)𝑘⁺𝑝⁰_𝑛+2. (21) Similar to the approach described in the previous section, we define the moments of the number of calcium ions in the subspace jointly distributed with the state of the channel, as follows: 𝜇_𝑞^0/1=^∞_𝑛=0𝑛^𝑞𝑝_𝑛^0/1, (22) where the superscript 0/1 indicates either index occurring on both the left and right hand sides of the equality. Note that the zeroth moments sum to unity by conservation of probability (𝜇⁰₀+𝜇¹₀=1). The expected number of calcium ions in the subspace conditioned on the channel being calcium free or bound, respectively, is given by E^0/1C=∑^∞_𝑛=0𝑛𝑝_𝑛^0/1∑^∞_𝑛=0𝑝_𝑛^0/1=𝜇₁^0/1𝜇₀^0/1. (23) Similarly the second moments 𝜇₂^0/1 are related to the conditional variances via Var^0/1C=𝜇₂^0/1𝜇₀^0/1−𝜇₁^0/1𝜇₀^0/12. (24)

2.2.3. Moment Calculation

By differentiating (22) with respect to time and substituting for the time derivatives using (21), it can be shown that the time-derivatives of the zeroth-moments, 𝜇⁰₀ and 𝜇¹₀—that is, the probability of the channel being in the calcium free or bound state—are given by 𝑑𝜇⁰₀𝑑𝑡=−𝑘⁺𝜇⁰₂+𝑘⁺𝜇⁰₁+𝑘⁻𝜇¹₀,𝑑𝜇¹₀=𝑑𝑡𝑘⁺𝜇⁰₂−𝑘⁺𝜇⁰₁−𝑘⁻𝜇¹₀, (25)(26) where we note that 𝑑𝜇⁰₀/𝑑𝑡+𝑑𝜇¹₀/𝑑𝑡=0 and 𝜇⁰₀+𝜇¹₀=1. In the same way, the equations for the first moments, 𝜇⁰₁ and 𝜇¹₁, are found to be 𝑑𝜇⁰₁=𝑑𝑡𝛼₀𝜇⁰₀−𝛽𝜇⁰₁−𝑘⁺𝜇⁰₃+𝑘⁺𝜇⁰₂+𝑘⁻𝜇¹₁+2𝑘⁻𝜇¹₀,𝑑𝜇¹₁=𝑑𝑡𝛼₁𝜇¹₀−𝛽𝜇¹₁+𝑘⁺𝜇⁰₃−3𝑘⁺𝜇⁰₂+2𝑘⁺𝜇⁰₁−𝑘⁻𝜇¹₁. (27) Setting the left hand side of (25) to zero, we find that the steady-state probability of a calcium-bound channel is 𝜇¹₀=𝜇⁰₂−𝜇⁰₁𝜅²𝑣²=E⁰C²−E⁰[C]/𝑣𝜅²+E⁰C²−E⁰[C]/𝑣, (28) where in the second equality we have used 𝜇⁰₂=𝜇⁰₀E⁰[C²]=𝜇⁰₀𝑣²E⁰[C²].

Note that as the volume increases (𝑣→∞), E⁰[C]/𝑣 becomes negligible compared to E⁰[C²], while E⁰[C²]→(E⁰[C])² as the conditional variance goes to zero (Var⁰[C]→0). Thus, in the large system limit, the probability that the channel is in the calcium-bound state is given by 𝜇¹₀=E⁰[C]²𝜅²+E⁰[C]². (29) In the case of a calcium-activated channel, 𝜇¹₀ is the open probability.

2.2.4. Analysis of Concentration Fluctuations

The moment analysis in the previous section suggests that the expected calcium concentration in the subspace given by E[C]=𝜇⁰₀E⁰[C]+𝜇¹₀E¹[C] (30) and the probability that a calcium-activated channel is open, 𝑝_open=𝜇¹₀, may depend on the subspace volume. In order to analyze the effect of small system size and concentration fluctuations at steady-state, we integrated (21) and determined the probability distributions {𝑝_𝑛^0/1} for various model parameters.

Figures 4(a) and 4(b) show the probability distribution for 𝑣=𝑣₀ and 8𝑣₀ using a representative set of parameters (see caption). In these calculations, the channel is closed when calcium-free and open when calcium-bound, that is, 𝛼₀=𝑣𝛽𝑐_∞<𝑣𝛽𝑐_∗=𝑣𝛼+𝛽𝑐_∞=𝛼₁. (31) For this reason, Figures 4(a) and 4(b) show a conditional expectation for the calcium concentration (vertical dotted lines) that is greater when the channel is calcium bound (E¹[C]>E⁰[C]). Note that the eight-fold increase in system size leads to a significant increase in the channel open probability, that is, 𝑝_open=𝜇¹₀ = 0.23 and 0.78 for 𝑣=𝑣₀ and 8𝑣₀, respectively. Thus, the open probability of the channel is significantly influenced by the subspace volume, in spite of the fact that the calcium influx rate is scaled so that in the absence of calcium-regulation there is no effect of volume (𝛼 constant, as in Section 2.1). Comparison of Figures 4(a) and 4(b) also shows a qualitative change in the probability distribution of the subspace calcium concentration (unimodal when 𝑣=𝑣₀, bimodal when 𝑣=8𝑣₀).

(a)

(b)

(c)

Figure 4 

Subspace volume-dependence of calcium fluctuations and open probability of a calcium-activated channel. Steady-state probability distribution for 𝑣 = 𝑣0 (a) and 8𝑣0 (b) for the calcium-activated channel (𝜅=2 μM, 𝑘+=0.05 μM−2 ms−1) [29, 30]. (c) Steady-state E[C] and 𝑝open=𝜇10 for integer multiples of the unitary volume 𝑣0 and different rate constants for calcium binding 𝑘+ (0.005 to 0.15 μM−2 ms−1) with 𝜅 fixed.

Figure 4(c) shows the expected calcium concentration E[C] (30) and open probability (𝑝_open=𝜇¹₀) for the calcium-activated channel as a function of subspace volume 𝑣 and different rate constants for calcium binding 𝑘⁺ (fixed dissociation constant 𝜅). Both E[C] and 𝑝_open increase with subspace volume 𝑣, that is, the restricted volume of a physiological subspace leads to an open probability and expected calcium concentration that is less than predicted in the corresponding (approximate) continuous description.

Both the open probability and expected calcium concentration asymptotically approach values in a range that are easily precalculated. For example, E[C] and 𝑝_open must be greater than the values obtained by assuming channel gating is extremely slow, in which case E[C]≈𝑐_∞ and 𝑝_open≈𝑐²_∞/(𝜅²+𝑐²_∞), because the transition from the calcium-free to -bound channel usually occurs with a subspace that is equilibrated with bulk calcium. In addition, E[C] and 𝑝_open are always less than the values obtained under the assumption of rapid channel binding, values given by simultaneous solution of 𝑝_open=𝑐²/(𝜅²+𝑐²) and 𝑐=𝑝_open𝛼/𝛽+𝑐_∞. These fast and slow system limits are indicated in Figure 4(c) by red and blue horizontal lines, respectively.

In order to further characterize the effect of subspace volume on the calcium-regulated channel and subspace dynamics, we defined the small system deviation Δ as E[C][C]Δ=−E_∞E[C]_∞, (32) where E[C] is calculated using a system volume of 𝑣=𝑣₀ and E[C]_∞ is the same quantity calculated in the large system size limit (𝑣→∞, numerically estimated using 𝑣⩾10𝑣₀). Figure 5 shows the small system deviation as a function of unitary subspace volume 𝑣₀ and influx parameter 𝑐_∗ for a calcium-activated channel. In all cases, Δ was negative, meaning that E[C] for 𝑣=𝑣₀ was suppressed below the large system size limit and increased with volume (cf. Figure 4(c)). For small 𝑐_∗, Δ was near zero. For intermediate values of 𝑐_∗ (5–10 𝜇M), the suppression was quite large (Δ≈−60%). As 𝑣₀ increased, Δ becomes less negative and approaches zero. In general, as 𝑐_∗ increased above this range, the suppression ultimately becomes negligible.

(a)

(b)

Figure 5 

Percentage small system deviation (Δ, (32)) as a function of unitary subspace volume 𝑣0 and influx parameter 𝑐∗ for a single calcium-activated channel (𝜅=2 μM, 𝑘+=0.005 μM−2 ms−1) and calcium-inactivated channel (𝜅=0.63 μM, 𝑘+=0.05 μM−2 ms−1).

For comparison, Figure 5(b) shows the small system deviation for a calcium-inactivated channel. In general, the magnitude of Δ for the calcium-inactivated channel was smaller than the calcium-activated channel. For small 𝑐_∗ (1–3 μM), the magnitude of Δ increased with 𝑐_∗, while above this range Δ was essentially independent of 𝑐_∗. As with the calcium-activated channel, the magnitude of Δ decreased and approached zero as 𝑣₀ increased.

2.3. Calcium Regulation of Multiple Channels

The previous section analyzed the effect of subspace volume when the influx pathway involves calcium regulation of a single channel. In this section, we assume that the total number of channels increases with subspace volume (see Figure 6). As before, we assume that calcium binding instantaneously modifies the calcium channel conductance, that is, the rate of calcium influx into the domain is determined by 𝛼₀ when all channels are calcium-free and 𝛼₁ when all channels are calcium-bound.

Figure 6 

Illustration of two possible volume scalings. For the single channel volume scaling, calcium influx 𝛼 increases proportional to the increase in 𝑣, but the single channel has only two conductance levels, 𝛼0 and 𝛼1, depending on whether calcium is free or bound. In the alternative scaling, the number of channels increases in proportion to the volume 𝑣, and when there are many channels the calcium influx rate may take many values between 𝛼0 and 𝛼1.

2.3.1. Deterministic Model

Assuming as before that two free calcium ions C bind to channel B to form the complex C₂B, we can write the following kinetic scheme: 2C+B⇌^𝑘+_𝑘−C₂B. (33) The deterministic ODE system that applies in the case of a large subspace volume is 𝑑𝑐𝑑𝑡=𝛼₀𝑏𝑏_𝑡+𝛼₁𝑏_𝑡−𝑏𝑏_𝑡−𝛽𝑐−𝑐_∞−𝑘⁺𝑐²𝑏+𝑘⁻𝑏_𝑡−𝑏,𝑑𝑏𝑑𝑡=−𝑘⁺𝑐²𝑏+𝑘⁻𝑏_𝑡−𝑏, (34) where we have written 𝑐=[C] and 𝑏=[B]. Because the total (calcium-free plus-bound) concentration of channels, 𝑏_𝑡=[B]+[C₂B], is a constant determined by initial conditions, we have eliminated the equation for [C₂B]. At steady-state the channels will be in equilibrium with subspace calcium, that is, 𝑏/𝑏_𝑡=𝜅²/𝜅²+𝑐².Thus, in the case of a calcium-activated channel (𝛼₀=0,𝛼₁>0), the steady-state calcium concentration satisfies 0=𝛼₁𝑐²𝜅²+𝑐²−𝛽𝑐−𝑐_∞, (35) while in the case of a calcium-inactivated channel (𝛼₀>0, 𝛼₁=0), 0=𝛼₀𝜅²𝜅²+𝑐²−𝛽𝑐−𝑐_∞. (36) Figure 7 shows bifurcation diagrams for the steady-state calcium concentration in both cases. For the calcium-activated channel there is a range of 𝜅 that leads to bistability (Figure 7(a)), while no bistable regime exists for a calcium-inactivated channel (Figure 7(b)).

(a)

(b)

Figure 7 

Bifurcation diagram showing the steady-state calcium concentration 𝑐 as a function of dissociation constant 𝜅 in the deterministic ODE model for a subspace containing multiple calcium-activated (a) and calcium-inactivated (b) channels. Other parameters as in Figure 2.

2.3.2. Stochastic Model

Following the notation developed in the previous section, we write 𝑝^𝑚_𝑛=P{C=𝑛,C₂B=𝑚,𝑡}=P{C=𝑛,B=𝑏̂_𝑡−𝑚} for 𝑛∈{0,1,…,∞} and 𝑚∈{0,1,…,𝑏̂_𝑡} and, where 𝑏̂_𝑡 is the total number of channels (for integer ℓ, 𝑏̂_𝑡=ℓ when 𝑣=ℓ𝑣₀). The state-transition diagram for the Markov process (not shown) is analogous to (19) but with 𝑏̂_𝑡+1 rows as opposed to two. The master equation for the dynamics of the calcium channel and subspace calcium concentration is 𝑑𝑝^𝑚_𝑛𝑑𝑡=−𝛼_𝑚+𝑛𝛽+𝑚𝑘⁻+𝑛(𝑛−1)𝑏̂_𝑡−𝑚𝑘⁺𝑝^𝑚_𝑛+𝛼_𝑚𝑝^𝑚_𝑛−1+(𝑛+1)𝛽𝑝^𝑚_𝑛+1+(𝑚+1)𝑘⁻𝑝^𝑚+1_𝑛−2+(𝑛+2)(𝑛+1)𝑏̂_𝑡−𝑚+1𝑘⁺𝑝^𝑚−1_𝑛+2, (37) where 𝛼_𝑚=𝑣(𝛼_𝑚+𝛽𝑐_∞) and 𝛼_𝑚=𝛼₀𝑏̂_𝑡−𝑚𝑏̂_𝑡+𝛼₁𝑚𝑏̂_𝑡. (38) In (38), it is understood that terms in the master equation involving negative indices (i.e., 𝑛 or 𝑚<0) evaluate to zero.

2.3.3. Concentration Fluctuations

Figure 8(A) shows the steady-state probability distribution for 𝑣=𝑣₀, 2𝑣₀ and 4𝑣₀ for a calcium-activated channel with dissociation constant chosen so that the deterministic system is monostable (𝜅=0.45 μM). For 𝑣=𝑣₀, there is one channel and two channel states (closed and open). For the closed channel, the distribution of calcium concentration is Poisson-like with conditional mean near 𝑐_∞, while for the open channel, the conditional mean is near 𝑐_∗. For 𝑣=2𝑣₀ and 4𝑣₀, there are two or four channels and thus three or five system states, each corresponding to a particular number of free versus bound channels (38). While the conditional expectation of the calcium concentration is always between 𝑐_∞ and 𝑐_∗, these distributions deviate from Poisson.

Figure 8 

Subspace volume-dependence of concentration fluctuation and channel open probability for multiple calcium-activated channels. (A) Steady-state probability distribution for 𝑣 = 𝑣0 (a), 2𝑣0 (b), and 4𝑣0 (c) (𝜅=0.45 μM, 𝑘+=5×10−4 μM−2 ms−1). For each panel, the dashed black line denotes the conditional expected concentration (E𝑚[C]). The steady-state probability distribution is shown for each possible number of closed (𝑁𝐶) and open (𝑁𝑂) channels. (B) (a) Steady-state E[C] and 𝑝open for the monostable system as a function of 𝑣 for rate constants of calcium binding (𝑘+=5×10−5 to 5×10−3 μM−2 ms−1). The fast/large and slow system limits are shown in red and blue, respectively. (b) Steady-state E[C] and 𝑝open for the bistable system (𝜅=2 μM) as a function of 𝑣 using 𝑘+= 0.005 to 0.015 μM−2 ms−1. In the bistable system, the larger of the two stable equilibrium (large system limit) is shown in red. The smaller equilibrium is approximately equal to the slow system limit (shown in blue).

Figure 8(B)a shows E[C] and 𝑝_open for subspace volumes 𝑣 given by different discrete multiples of the unitary volume 𝑣₀. Using parameters that lead to a monostable deterministic ODE system, we find, similar to the case of the single channel (Figure 4), a significant deviation between the expected calcium concentration and open probability for a small subspace as compared to the large system limit (35). As expected, both E[C] and 𝑝_open approached the fast/large system limit as 𝑣 increased. This also occurs for fixed 𝑣 with increasing 𝑘⁺, that is, the rate constant for calcium binding. For fixed 𝜅, smaller values of 𝑘⁺ can cause Δ to approach −100%, that is, the small volume associated with a diadic subspace can almost completely suppress the open probability of a calcium-activated channel. When parameters are chosen so that the deterministic ODE system is bistable, the dependence of E[C] and 𝑝_open is more complex (Figure 8(B)b). Interestingly, the small system deviation in this case is often a biphasic function of system volume.

Figure 9 shows analogous results for calcium-inactivated calcium influx. As with the calcium-activated channel, E[C] and 𝑝_open were suppressed below the fast/large system limit (Figure 9(B)). Δ is often negative, but became negligible as 𝑘⁺ increased. Similarly, as 𝑣 increased, both E[C] and 𝑝_open approached the fast/large system limit.

Figure 9 

Subspace volume-dependence of concentration fluctuation and channel open probability for multiple calcium-inactivated channels. (A) Steady-state probability distribution for 𝑣 = (a) 𝑣0, (b) 2𝑣0, and (c) 4𝑣0 (𝜅=0.63 μM, 𝑘+=0.005 μM−2 ms−1). (B) Steady-state E[C] and 𝑝open as a function of 𝑣 (𝑘+=0.0015 to 0.15 μM−2 ms−1).

Figure 10 summarizes the dependence of the small system deviation (Δ, (32)) on the unitary subspace volume (𝑣₀) and calcium influx parameter (𝑐_∗) for the scaling that involves multiple calcium-activated and -inactivated channels. In all cases, Δ was negative, meaning that E[C] was suppressed compared with the large system values predicted by the deterministic ODE model. Up to 80% suppression was observed for the calcium-activated channel, but for the calcium-inactivated channel the maximum suppression was 20%. In both cases, the largest suppression (most negative Δ) occurs when 𝑝_open is small (i.e., small 𝑐_∗ for the calcium-activated channel and large 𝑐_∗ for the calcium-inactivated channel). In general, as the unitary volume 𝑣₀ is increased, there is less suppression compared to the large system size limit.

(a)

(b)

Figure 10 

Percentage small system deviation (Δ) as a function of unitary subspace volume 𝑣0 and influx parameter 𝑐∗ for multiple (a) calcium-activated channels (𝜅=0.45 μM, 𝑘+=0.005 μM−2 ms−1) and (b) calcium-inactivated channels (𝜅=0.63 μM, 𝑘+=0.05 μM−2 ms−1).

2.4. The Effect of Domain Size in a Model of Calmodulin-Mediated Channel Regulation

In the previous sections, we demonstrated that the expected steady-state subspace concentration determined using a minimal model of a calcium-activated or -inactivated channel was volume-dependent and could greatly differ from the steady-state concentration computed from deterministic ODEs. In this section, we show similar results for a state-of-the-art model of calmodulin-mediated calcium regulation.

Both the N-lobe and C-lobe of calmodulin have two binding sites for calcium. Depending on the calcium channel type (L, N, or P/Q), calcium binding to the C-lobe has been shown to be responsible for either activation or inactivation of the channel, while N-lobe binding appears to be primarily responsible for channel inactivation [32]. Yue and colleagues demonstrated that the C-lobe responds primarily to the local subspace calcium concentration, while the N-lobe responds to the global or bulk concentration [31]. Tadross et al. developed a 4-state model for calmodulin regulation of the calcium channel (see Figure 11(A)) that includes states for the calmodulin regulator lobe (either the C-lobe or N-lobe) bound to a preassociation site that does not alter channel activity (state 1), unbound (state 2), bound to two calcium ions (state 3), or bound to two calcium ions and an effector site that does alter channel activity (state 4) [31]. Tadross et al. demonstrated that depending on the model parameters, in particular the ratio of the transition rates between states, the calmodulin regulation was sensitive to either local or global calcium levels.

Figure 11 

Calmodulin regulation of the calcium channel at steady-state. (A) State diagram of calmodulin regulation of a calcium channel (modified from [31]). (B) Steady-state probability distribution for the calmodulin-activated channel using (a) “slow CaM" and (b) “SQS" parameters [31].

Using this published model as a starting point, we formulated the corresponding discrete Markov model. The elementary reactions for calmodulin-mediated regulation of the channel are 𝑆⇌₁^𝛾−_𝛿−𝑆⇌₂^𝑘+𝑐2_𝑘−𝑆⇌₃^𝛾+_𝛿+𝑆₄, (39) where states 𝑆₁ and 𝑆₂ are calcium-free, states 𝑆₃ and 𝑆₄ are calcium-bound, and state 𝑆₄ determines the fraction of channels activated (or inactivated) by calmodulin. When it is assumed that a single calmodulin molecule is colocalized with the calcium channel (as in Section 2.2), the master equation takes the following form: 𝑑𝑝¹_𝑛𝑑𝑡=−𝛼₀+𝑛𝛽+𝛾⁻𝑝¹_𝑛+𝛼₀𝑝¹_𝑛−1+(𝑛+1)𝛽𝑝¹_𝑛+1+𝛿⁻𝑝²_𝑛,𝑑𝑝²_𝑛𝑑𝑡=−𝛼₀+𝑛𝛽+𝛿⁻+𝑛(𝑛−1)𝑘⁺𝑝²_𝑛+𝛼₀𝑝²_𝑛−1+(𝑛+1)𝛽𝑝²_𝑛+1+𝛾⁻𝑝¹_𝑛+𝑘⁻𝑝³_𝑛−2,𝑑𝑝³_𝑛𝑑𝑡=−𝛼₀+𝑛𝛽+𝛾+𝑘⁻𝑝³_𝑛+𝛼₀𝑝³_𝑛−1+(𝑛+1)𝛽𝑝³_𝑛+1+𝛿⁺𝑝⁴_𝑛+(𝑛+2)(𝑛+1)𝑘⁺𝑝²_𝑛+2,𝑑𝑝⁴_𝑛𝑑𝑡=−𝛼₁+𝑛𝛽+𝛿⁺𝑝⁴_𝑛+𝛼₁𝑝⁴_𝑛−1+(𝑛+1)𝛽𝑝⁴_𝑛+1+𝛾⁺𝑝³_𝑛, (40) where for a calmodulin-activated channel 𝛼₀ and 𝛼₁ are given by (31).

Figure 11(B) shows the steady-state probability distribution numerically calculated from these equations. Using a parameter set referred to as “slow CaM,” Tadross et al. showed that calmodulin was primarily sensitive to the local subspace calcium level (representing the C-lobe) when calcium binding to calmodulin was slow [31]. With “slow CaM” parameters, we found that calmodulin bound to the effector site (𝑆₄) had the greatest steady-state probability. Because the calmodulin binding was slow, each conditional distribution had their respective largest peaks near the slow limit (𝑐_∞ for states 𝑆₁ and 𝑆₂, 𝑐_∗ for 𝑆₃ and 𝑆₄) (Figure 11(B)a). Using an alternate parameter set referred to as “SQS," Tadross et al. showed that calmodulin was primarily sensitive to the global calcium level (representing the N-lobe), when calcium binding to calmodulin was fast. Similar to the slow CaM case, state 𝑆₄ had the greatest steady-state probability using the SQS parameters. Due to the fast binding kinetics, the conditional distributions were more similar than in the slow CaM case (Figure 11(B)b). For both parameter sets, the calcium concentration distribution for the low occupancy states (𝑆₁ and 𝑆₂) were bimodel, with peaks near 𝑐_∞ and 𝑐_∗.

Figure 12 shows the small system size suppression Δ for both slow CaM and SQS parameter sets assuming a single calmodulin-regulated channel. As in our simplified model (Figure 5), Δ was quite large in magnitude for some conditions (up to 30% suppression). For the calmodulin-activated channel, the dependence of Δ on 𝑣₀ and 𝑐_∗ was similar to our simplified model (cf. Figure 5(a)), decreasing in magnitude and approaching 0 as 𝑣₀ or 𝑐_∗ increased.

Figure 12 

Small system deviation for a single calmodulin-activated and -inactivated channel using “slow CaM” and “SQS” parameters (see text).

The parameter space for the calmodulin-inactivated channel differed somewhat from our simplified model (Figure 5(b)). For both the slow CaM and SQS parameters, Δ decreased as 𝑐_∗ or 𝑣₀ increased. Additionally, for both the calmodulin-activated and -inactivated channels, Δ had greater dependence on 𝑣₀ for the SQS parameters, which is consistent with calmodulin being more sensitive to the bulk concentration (since increasing 𝑣₀ greatly influences the number of ions entering from the bulk).

We also calculated the small system deviation Δ for the case of multiple calmodulin-regulated channels (Figure 13). For the calmodulin-activated channel, results were similar to our simplified model (Figure 10(a)), in particular Δ approached 0 as both 𝑐_∗ and 𝑣₀ increased. The magnitude of Δ was smaller for the SQS parameters compared with the slow CaM parameters, consistent with faster kinetics approaching the large system limit and calmodulin being less sensitive to the local calcium concentration. The parameter space for multiple calmodulin-inactivated channels also differed somewhat from our simplified model (Figure 10(b)). In general, the magnitude of Δ decreased as 𝑐_∗ increased. However, in contrast with the parameter space using SQS parameters, Δ was fairly insensitive to 𝑣₀ using slow CaM parameters, which is consistent with calmodulin being, in this case, less sensitive to the bulk calcium concentration.

Figure 13 

Small system deviation for multiple calmodulin-activated and -inactivated channels, using “slow CaM” and “SQS” parameters.

3. Discussion

We developed a minimal model of a calcium-regulated channel in a small subspace and formulated a Markov model in which each possible discrete state is represented. For small subspace volumes, we found that the value predicted by a continuous-state, deterministic ODE model often deviated from the expected steady-state calcium concentration in the discrete-state, stochastic model. We analyzed how this deviation depends on channel binding kinetics, subspace volume, and calcium influx rate. We demonstrated that the deterministic description also deviated from the stochastic model in a physiologically realistic model of calmodulin-mediated calcium channel regulation.

3.1. Physiological Implications

Many studies have modeled the influence of signaling proteins on intracellular and transmembrane ion channel/receptor kinetics, such as calcium/calmodulin-dependent kinase II phosphorylation [33] or beta-adrenergic signaling [34] in cardiac myocytes and glutamate receptor activation in neurons [35]. Many of these signaling interactions occur in small volumes (e.g., the cardiac dyad [4] and neuronal synapse [36]) and include binding interactions with species present in low concentration (calcium and glutamate, resp.). In cardiac myocytes, the local calcium concentration can greatly influence the whole cell response through calcium-induced calcium release, the sodium-calcium exchanger current (which can trigger activation of an action potential), and a host of intracellular signaling pathways [37]. We found that a stochastic model that accounts for the discrete nature of such interactions may deviate from the corresponding deterministic ODE model. Under certain conditions, the small system deviation is negligible, in particular for the case of a large calcium influx rate (Figure 10). During a cardiac action potential, many L-type calcium channels are synchronously opened, and thus the calcium concentration rapidly increases from the micro- to millimolar range. Similarly, following neuronal firing the glutamate concentration in the synaptic cleft can increase several orders of magnitude [35]. In these situations, the deviation of species concentrations from that suggested by deterministic ODE models may not be physiologically significant. However, during resting conditions, the deviation may be significant, and concentration fluctuations due to the small subspace volume could influence channel dynamics (Sections 2.2–2.4). It has been shown that stochastic openings in calcium release channels in the dyadic subspace of cardiac myocytes during diastole can lead to spontaneous calcium release and arrhythmias during heart failure [38]. Our findings demonstrate that a discrete model of the subspace concentration may be important in this physiological context, because it is likely that fluctuations due to the small number of calcium ions play a significant role in generating spontaneous calcium release events.

In addition to demonstrating that a discrete/stochastic model of calcium-regulated calcium influx often deviates from a continuous/deterministic description, we analyzed how subspace volume and concentration fluctuations influence channel dynamics. Because calmodulin effectively colocalizes with the L-type calcium channel [39], the results associated with the “single channel” volume scaling (Figure 5 and Section 2.2) are most relevant. Such colocalization is ubiquitous; many regulators have been shown to colocalize with channels or receptors, including phospholamban with calcium ATPase in the sarcoplasmic reticulum membrane [40], G-protein receptor kinases with G-protein receptors on the cell membrane [41], and Bax with voltage-dependent ion channels in the mitochondrial membrane [42]. Additionally, the volume of diadic subspaces can be greatly altered during pathophysiological conditions. For example, the L-type calcium channels and ryanodine receptors localization in the dyad is disrupted during heart failure and the subspace volume in which these channels reside is much greater in heart failure than during physiological conditions [43]. Our findings show that when a small number of molecules are present in the subspace (small 𝑣₀ and 𝑐_∗), subspace volume can greatly influence the steady-state properties of stochastically gating channels (Figures 10 and 13).

3.2. Relation to Prior Studies

Prior work by our lab has investigated calcium channel regulation through a host of various mechanisms. Groff and Smith investigated the influence of inactivation on calcium spark dynamics in a channel regulated by both calcium-activation and -inactivation [44]. Mazzag et al. demonstrated that residual calcium that accumulates in a subspace during channel openings can influence channel gating [27]. Perhaps more relevant to this study of how concentration fluctuations depend on the subspace volume and influence average rates of calcium binding, channel gating, and calcium influx, Smith and coworkers previously investigated how the number of subspace domains and the number of channels per subspace can influence cellular responses. Williams et al. demonstrated that a population of subspace domains can be represented by a probability density approach and can be utilized to simulate global calcium dynamics [45]. Hartman et al. utilized a model of a small number of coupled calcium activated channels to predict the global calcium release dynamics in response to pharmacological modification of single channel kinetics [46]. However, this study is the first to compare a model of calcium channel regulation accounting for the finite subspace volume (and using a discrete representation of the number of subspace calcium ions) with the corresponding ODE formulation that assumes a large system size (and uses a continuous representation of calcium concentration).

Only a few previous studies have utilized a discrete representation of calcium ions in the context of cardiac myocyte subspace dynamics. Winslow and colleagues simulated the spatial location of discrete diffusing calcium ions, as well as the spatial structure and geometry of the L-type calcium channel and ryanodine receptor in the cardiac dyad [47]. They demonstrate that stochastic fluctuations produce variability in the L-type calcium channel-ryanodine receptor signaling interactions (specifically excitation-contraction coupling gain), but their analysis does not distinguish between the influence of fluctuations due to channel gating, calcium diffusion, and small calcium ion number. Similar to this study, von Wegner and Fink presented a stochastic model of the L-type calcium channel, incorporating calcium diffusion, buffering, and channel gating and conductance, and they demonstrated how calcium concentration fluctuations could influence downstream signalling pathways [48]. Our results are novel in their focus on the influence of subspace volume and the kinetics of calcium-regulation of an L-type channel. Most importantly, we provide a thorough analysis of the deviation of the approximate deterministic description from the full stochastic model and clarify the conditions leading to large versus small deviations.

Previous studies have modeled biochemical reaction networks using master equations and compared results with deterministic ODE models. McQuarrie demonstrated in 1963 that for first-order reactions, the expected steady-state concentrations derived from the chemical master equation and deterministic ODEs agree [5]. In Section 2.1, calcium influx from an unregulated channel is modeled using zeroth- and first-order reactions and, consequently, the stochastic and deterministic descriptions must agree. Our observation that concentration fluctuations increased as the subspace volume became smaller is consistent with well-understood principles of statistical physics and should come as no surprise [23].

Darvey et al. demonstrated for several generic second-order reactions, the expected concentration computed from the chemical master equation may deviate from the corresponding ODE model [20]. The deviation is typically negative (i.e., Δ<0), with greatest suppression when concentration fluctuations are large. Other recent studies have demonstrated that the concentrations of species in stochastic biochemical networks can deviate from deterministic ODE descriptions. In agreement with our findings, the deviation is often negative [18, 21, 28], although positive deviation was observed in some biochemical systems [12, 49]. Our findings are consistent with Darvey et al., in that Δ had the greatest magnitude when either the subspace volume or calcium influx rate was small (Figure 10) (both result in larger concentration fluctuations, see (14)).

We found that the small system size deviation was particularly complex in cases where the deterministic ODE, that is, the model appropriate for the large system size limit, is bistable (Figure 8). Lestas et al. recently investigated bistability/bimodality in a network of gene regulation and demonstrated that bistability in the deterministic ODE model did not imply bimodality in the discrete system and, conversely, bimodality in the discrete system did not imply bistability in the corresponding ODEs [50]. We obtained similar results, as bimodality in distribution was not present in the bistable system for 𝑣=𝑣₀ (Figure 4(a)) but was present for 𝑣=8𝑣₀ (Figure 4(b)). Conversely, for the calmodulin-regulated channel, bimodality was present in the distribution for the monostable system (Figure 11). Interestingly, for the bistable system, E[C] computed from the discrete model need not be well approximated by either of the two stable equilibria in the deterministic model; rather, E[C] is given by an intermediate value and can have a complex dependence of subspace volume (Figure 8(C)). But it is important to note that the small system size deviation does not require a bistable deterministic model. The deviation can be quite pronounced even in a monostable deterministic model (Figures 8(B), 9, and 10).

3.3. Limitations

The two-state kinetic models of the calcium channel introduced in Section 2.1 is minimal and should be interpreted as phenomenological (as opposed to statistical) model of single channel kinetics, that is, the topology and parameters of this model were not obtained by fitting to patch clamp recordings [51]. On the other hand, the kinetic model for regulation of the calcium channel presented in Section 2.4 is state-of-the-art. Both minimal and physiologically realistic channel models are affected by the decision to account for (or neglect) fluctuations in calcium concentration that result from the small number of ions in the subspace.

The most significant limitation in the model formulation is our neglect of spatial dynamics of calcium diffusion within the dyadic subspace and the details of the spatial arrangement of the ryanodine receptors [47, 52]. However, for the purposes of the present study, that is, investigation of the influence of concentration fluctuations on the regulation of calcium influx, a nonspatial Markov chain model that includes subspace volume as a model parameter and accounts for the finite number of calcium ions in the domain is sufficient.

Another limitation of the present work is that we focus on stationary statistics, for example, the expected value of the steady-state subspace calcium concentration, in our analysis of the deviation of continuous ODE description from the discrete stochastic formulation. Future studies could address how transient dynamics, for example, the cellular response to a depolarizing voltage step, excitation-contraction coupling gain, and so forth, are affected by calcium concentration fluctuations resulting from small subspace volume.

4. Conclusions

Our findings demonstrate the physiological relevance of concentration fluctuations in both minimal and realistic models of a calcium-regulated channels associated with subspaces of small volume. The take home message is: concentration fluctuations do not “average out” in a manner that causes stochastic and deterministic descriptions of subspace dynamics to be equivalent. Future studies will investigate how subspace calcium concentration fluctuations may influence global calcium dynamics and plasma membrane electrical activity in physiological and pathophysiological conditions.

Conflict of Interests

The authors do not have any conflict of interests.

Acknowledgments

The work was supported in part by National Science Foundation Grant DMS-1121606 and the Biomathematics Initiative at The College of William & Mary. The authors acknowledge stimulating discussions with graduate students Xiao Wang and Mary Mohr. This study was inspired by comments made at the 2012 UC Davis Cardiovascular Symposium.