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

Authors: Seth H. Weinberg; Gregory D. Smith
Journal: Computational and Mathematical Methods in Medicine
 (2012)

---

## 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−17 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.


---

## Body

## 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−17 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). Section2.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 (Figure1 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−1 to the constant bulk concentration of c∞=0.1 μM [27]. These assumptions lead to the following deterministic model of subspace calcium dynamics:
(1)dcdt=α-β(c-c∞),
where the influx rate α has units of concentration per time (e.g., μM/ms) and the calcium concentration c 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 c (units of μM) in a diadic subspace of volume v (liters). Subspace calcium moves to the bulk passively via diffusion at rate β (given by 0.01 ms−1). Bulk calcium at the concentration c∞=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 byC^∈{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 (v), that is,
(2)C=C^v.
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
(3)0⇌βα-1⇌2βα-2⋯n-1⇌nβα-n⇌(n+1)βα-n+1⋯,
where the index that labels states, n∈{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)α-=v(α+βc∞).


#### 2.1.2. Master Equation and Steady-State Probability Distribution

If we writepn(t)=Pr{C^(t)=n}, 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
(5)dp0dt=-α-p0+βp1,dpndt=-(α-+nβ)pn+α-pn-1+(n+1)βpn+1,n=1,2,….
Note that the correspondence between the rate constants in the deterministic (1) and stochastic (3)–(5) models is established by substituting c=c^/v in (1) to find the rate of change of the number of calcium ions in the deterministic model, that is,
(6)dc^dt=α--βc^.
This equation indicates that c^ 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, c^ decreases at rate βc^, a value that is proportional to c^ because each ion has an opportunity to diffuse into the bulk. Consequently, the transition rates leading out of state C^=n in the stochastic model are given by α- for the C^=n to n+1 transitions and βn for the C^=n to n-1 transitions.To find the steady-state probability distribution ofC^, we set the left hand sides of (5) to zero to obtain
(7)nβpn=α-pn-1,n=1,2,…
from which it follows that {pn} is a Poisson distribution with parameter λ=α-/β, that is,
(8)pn=e-λλnn!.


#### 2.1.3. Analysis of Concentration Fluctuations

To see how the subspace calcium concentration fluctuations predicted by this minimal model depend on the parametersα, β, c∞, and v, 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
(9)E[C^]=∑n=0∞‍npn=λ=α-β=v(αβ+c∞)=vc*,
where the last equality defines c* as follows:
(10)c*=αβ+c∞.
Using (9) and the fact that C=C^/v implies E[C]=E[C^]/v, we can identify c* as the expected subspace calcium concentration:
(11)E[C]=c*.
Similarly, the steady-state variance of the number of calcium ions in the subspace is
(12)Var[C^]=∑n=0∞‍(n-E[C])2pn=vc*,
and Var[C]=Var[C^]/v2 implies that the variance of the subspace calcium concentration is
(13)Var[C]=c*v.
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/vc* and, similarly,
(14)CV[C]=Var[C]E[C]=Var[C^]/v2E[C^]/v=1vc*.
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 v).Figure2 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 v=v0=10-17 liters and influx rate of α=0.049 μM/ms, the expected calcium concentration is E[C]=c*=α/β+c∞=5 μM, and the expected number of subspace calcium ions is E[C^]=v0c*=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 x-axis scales). The following rows of Figure 2 show that in a subspace three or ten times larger (v=3v0 or 10v0), 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 (c* 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/vc*→0 as v→∞ for fixed c*.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 v0=10-17 liters and subspaces that are 3 and 10 times larger. Parameters: α=0.049 μM/ms, β=0.01 ms−1, c∞=0.1 μM; the steady-state expected subspace calcium concentration is  E[C]=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 volumev. The calcium influx rate is α0 and α1 when calcium is unbound and bound, respectively, and the transition rates between these states are k+c2 and k-, where c is the subspace calcium concentration. Subspace calcium is passively coupled at rate β to the bulk cytosol with constant concentration c∞.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 modelE[C]=c*=α/β+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(15)μq=∑n=0∞‍nqpn.
By conservation of probability, the zeroth moment μ0=1 and the first moment is the expected number of calcium ions in the subspace (9),
(16)μ1=E[C^].
The second moment μ2 is related to the variance of the number of calcium ions via
(17)Var[C^]=μ2-(μ1)2.
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 (dμ0/dt=0) and, furthermore,
(18)dμ1dt=α--βμ1,dμ2dt=α-+(2α-+β)μ1-2βμ2,
where we have used μ0=1. Setting the left hand sides of these equations to zero, we see that steady-state first and second moment are μ1=α-/β and μ2=α-/β+(α-/β)2, 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 Figure3). We assume that calcium binding instantaneously modifies the conductance of the channel, that is, the rate of calcium influx into the domain is α0 when the channel is calcium-free and α1 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 k+c2 and k-, respectively, (Figure 3, curved arrows). Note that the rate constant k- has units of ms−1, k+ has units of  μM−2 ms−1, and the dissociation constant for calcium binding, denoted by κ, has units of μM and is given by κ2=k-/k+.


#### 2.2.1. Stochastic Model

Let us denote the states of the stochastic system by(n,0) and (n,1), where n∈{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,(19)(0,0)⇌βα-0(1,0)⇌2βα-0(2,0)⇌3βα-0(3,0)⇌4βα-0(4,0)⋯2k-+⥯k-6k-+⥯k-12k-+⥯k-(0,1)⇌βα-1(1,1)⇌2βα-1(2,1)⋯,wherek-+=k+/v2. 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 (k+c2=k+c^2/v2). The downward transitions between states (n,0) and (n-2,1) include the combinatorial coefficient, n(n-1), double the number of ways that two indistinguishable calcium ions can be chosen from the n ions in the subspace. This factor of two is required so that the microscopic propensity k-+ agrees with the macroscopic rate k+c2 for large n and v with c=n/v fixed, that is,
(20)n(n-1)k-+=n(n-1)k+v2=k+(c2-cv),
an expression that approaches k+c2 as v→∞ [28].


#### 2.2.2. Master Equation

Let us writepn0(t) to indicate the probability that at time t the channel is calcium-free and C^=n. Similarly, pn1(t) is the probability that C^(t)=n 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:
(21)dpn0dt=-[α-0+nβ+n(n-1)k-+]pn0+α-0pn-10+(n+1)βpn+10+k-pn-21,dpn1dt=-[α-1+nβ+k-]pn1+α-1pn-11+(n+1)βpn+11+(n+2)(n+1)k-+pn+20.
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:
(22)μq0/1=∑n=0∞‍nqpn0/1,
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 (μ00+μ01=1). The expected number of calcium ions in the subspace conditioned on the channel being calcium free or bound, respectively, is given by
(23)E0/1[C^]=∑n=0∞npn0/1∑n=0∞pn0/1=μ10/1μ00/1.
Similarly the second moments μ20/1 are related to the conditional variances via
(24)Var0/1[C^]=μ20/1μ00/1-(μ10/1μ00/1)2.


#### 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, μ00 and μ01—that is, the probability of the channel being in the calcium free or bound state—are given by
(25)dμ00dt=-k-+μ20+k-+μ10+k-μ01,(26)dμ01dt=k-+μ20-k-+μ10-k-μ01,
where we note that dμ00/dt+dμ01/dt=0 and μ00+μ01=1. In the same way, the equations for the first moments, μ10 and μ11, are found to be
(27)dμ10dt=α-0μ00-βμ10-k-+μ30+k-+μ20+k-μ11+2k-μ01,dμ11dt=α-1μ01-βμ11+k-+μ30-3k-+μ20+2k-+μ10-k-μ11.
Setting the left hand side of (25) to zero, we find that the steady-state probability of a calcium-bound channel is
(28)μ01=μ20-μ10κ2v2=E0[C2]-E0[C]/vκ2+E0[C2]-E0[C]/v,
where in the second equality we have used μ20=μ00E0[C^2]=μ00v2E0[C2].Note that as the volume increases (v→∞), E0[C]/v becomes negligible compared to E0[C2], while E0[C2]→(E0[C])2  as the conditional variance goes to zero (Var0[C]→0). Thus, in the large system limit, the probability that the channel is in the calcium-bound state is given by
(29)μ01=(E0[C])2κ2+(E0[C])2.
In the case of a calcium-activated channel, μ01 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(30)E[C]=μ00E0[C]+μ01E1[C]
and the probability that a calcium-activated channel is open, popen=μ01, 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 {pn0/1} for various model parameters.Figures4(a) and 4(b) show the probability distribution for v=v0 and 8v0 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)α-0=vβc∞<vβc*=v(α+βc∞)=α-1.
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 (E1[C]>E0[C]). Note that the eight-fold increase in system size leads to a significant increase in the channel open probability, that is, popen=μ01 = 0.23 and 0.78 for v=v0 and 8v0, 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 v=v0, bimodal when v=8v0).Subspace volume-dependence of calcium fluctuations and open probability of a calcium-activated channel. Steady-state probability distribution forv = v0 (a) and 8v0 (b) for the calcium-activated channel (κ=2 μM, k+=0.05 μM−2 ms−1) [29, 30]. (c) Steady-state E[C] and popen=μ01 for integer multiples of the unitary volume v0 and different rate constants for calcium binding k+ (0.005 to 0.15 μM−2 ms−1) with κ fixed.
(a)
(b)
(c)Figure4(c) shows the expected calcium concentration E[C] (30) and open probability (popen=μ01) for the calcium-activated channel as a function of subspace volume v and different rate constants for calcium binding k+ (fixed dissociation constant κ). Both E[C] and popen increase with subspace volume v, 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 popen must be greater than the values obtained by assuming channel gating is extremely slow, in which case E[C]≈c∞ and popen≈c∞2/(κ2+c∞2), 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 popen are always less than the values obtained under the assumption of rapid channel binding, values given by simultaneous solution of popen=c2/(κ2+c2) and c=popenα/β+c∞. 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
(32)Δ=E[C]-E[C]∞E[C]∞,
where E[C] is calculated using a system volume of v=v0 and E[C]∞ is the same quantity calculated in the large system size limit (v→∞, numerically estimated using v⩾10v0). Figure 5 shows the small system deviation as a function of unitary subspace volume v0 and influx parameter c* for a calcium-activated channel. In all cases, Δ was negative, meaning that E[C] for v=v0 was suppressed below the large system size limit and increased with volume (cf. Figure 4(c)). For small c*, Δ was near zero. For intermediate values of c* (5–10 μM), the suppression was quite large (Δ≈-60%). As v0 increased, Δ becomes less negative and approaches zero. In general, as c* increased above this range, the suppression ultimately becomes negligible.Percentage small system deviation (Δ, (32)) as a function of unitary subspace volume v0 and influx parameter c* for a single calcium-activated channel (κ=2 μM, k+=0.005 μM−2 ms−1) and calcium-inactivated channel (κ=0.63 μM, k+=0.05 μM−2 ms−1).
(a)
(b)For comparison, Figure5(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 c* (1–3 μM), the magnitude of Δ increased with c*, while above this range Δ was essentially independent of c*. As with the calcium-activated channel, the magnitude of Δ decreased and approached zero as v0 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 Figure6). 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 α0 when all channels are calcium-free and α1 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 v, 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 v, 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 ionsC bind to channel B to form the complex C2B, we can write the following kinetic scheme:
(33)2C+B⇌k-k+C2B.
The deterministic ODE system that applies in the case of a large subspace volume is
(34)dcdt=α0bbt+α1bt-bbt-β(c-c∞)-k+c2b+k-(bt-b),dbdt=-k+c2b+k-(bt-b),
where we have written c=[C] and b=[B]. Because the total (calcium-free plus-bound) concentration of channels, bt=[B]+[C2B], is a constant determined by initial conditions, we have eliminated the equation for [C2B]. At steady-state the channels will be in equilibrium with subspace calcium, that is, b/bt=κ2/κ2+c2.Thus, in the case of a calcium-activated channel (α0=0,α1>0), the steady-state calcium concentration satisfies
(35)0=α1c2κ2+c2-β(c-c∞),
while in the case of a calcium-inactivated channel (α0>0, α1=0),
(36)0=α0κ2κ2+c2-β(c-c∞).
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)).Bifurcation diagram showing the steady-state calcium concentrationc 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.
(a)
(b)


#### 2.3.2. Stochastic Model

Following the notation developed in the previous section, we writepnm=P{C^=n,C2B^=m,t}=P{C^=n,B^=b^t-m} for n∈{0,1,…,∞} and m∈{0,1,…,b^t} and, where b^t is the total number of channels (for integer ℓ, b^t=ℓ when v=ℓv0). The state-transition diagram for the Markov process (not shown) is analogous to (19) but with b^t+1 rows as opposed to two. The master equation for the dynamics of the calcium channel and subspace calcium concentration is
(37)dpnmdt=-[α-m+nβ+mk-+n(n-1)(b^t-m)k-+]pnm+α-mpn-1m+(n+1)βpn+1m+(m+1)k-pn-2m+1+(n+2)(n+1)(b^t-m+1)k-+pn+2m-1,
where α-m=v(αm+βc∞) and
(38)αm=α0b^t-mb^t+α1mb^t.
In (38), it is understood that terms in the master equation involving negative indices (i.e., n or m<0) evaluate to zero.


#### 2.3.3. Concentration Fluctuations

Figure8(A) shows the steady-state probability distribution for v=v0, 2v0 and 4v0 for a calcium-activated channel with dissociation constant chosen so that the deterministic system is monostable (κ=0.45 μM). For v=v0, 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 c∞, while for the open channel, the conditional mean is near c*. For v=2v0 and 4v0, 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 c∞ and c*, 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 forv = v0 (a), 2v0 (b), and 4v0 (c) (κ=0.45 μM, k+=5×10-4 μM−2 ms−1). For each panel, the dashed black line denotes the conditional expected concentration (Em[C]). The steady-state probability distribution is shown for each possible number of closed (NC) and open (NO) channels. (B) (a) Steady-state E[C] and popen for the monostable system as a function of v for rate constants of calcium binding (k+=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 popen for the bistable system (κ=2 μM) as a function of v using k+= 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).Figure8(B)a shows E[C] and popen for subspace volumes v given by different discrete multiples of the unitary volume v0. 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 popen approached the fast/large system limit as v increased. This also occurs for fixed v with increasing k+, that is, the rate constant for calcium binding. For fixed κ, smaller values of k+ 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 popen is more complex (Figure 8(B)b). Interestingly, the small system deviation in this case is often a biphasic function of system volume.Figure9 shows analogous results for calcium-inactivated calcium influx. As with the calcium-activated channel, E[C] and popen were suppressed below the fast/large system limit (Figure 9(B)). Δ is often negative, but became negligible as k+ increased. Similarly, as v increased, both E[C] and popen 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 forv = (a) v0, (b) 2v0, and (c) 4v0 (κ=0.63 μM, k+=0.005 μM−2 ms−1). (B) Steady-state E[C] and popen as a function of v (k+=0.0015 to 0.15 μM−2 ms−1).Figure10 summarizes the dependence of the small system deviation (Δ, (32)) on the unitary subspace volume (v0) and calcium influx parameter (c*) 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 popen is small (i.e., small c* for the calcium-activated channel and large c* for the calcium-inactivated channel). In general, as the unitary volume v0 is increased, there is less suppression compared to the large system size limit.Percentage small system deviation (Δ) as a function of unitary subspace volume v0 and influx parameter c* for multiple (a) calcium-activated channels (κ=0.45 μM, k+=0.005 μM−2 ms−1) and (b) calcium-inactivated channels (κ=0.63 μM, k+=0.05 μM−2 ms−1).
(a)
(b)


### 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(39)S1⇌δ-γ-S2⇌k-k+c2S3⇌δ+γ+S4,
where states S1 and S2 are calcium-free, states S3 and S4 are calcium-bound, and state S4 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:
(40)dpn1dt=-(α-0+nβ+γ-)pn1+α-0pn-11+(n+1)βpn+11+δ-pn2,dpn2dt=-[α-0+nβ+δ-+n(n-1)k-+]pn2+α-0pn-12+(n+1)βpn+12+γ-pn1+k-pn-23,dpn3dt=-(α-0+nβ+γ+k-)pn3+α-0pn-13+(n+1)βpn+13+δ+pn4+(n+2)(n+1)k-+pn+22,dpn4dt=-(α-1+nβ+δ+)pn4+α-1pn-14+(n+1)βpn+14+γ+pn3,
where for a calmodulin-activated channel α-0 and α-1 are given by (31).Figure11(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 (S4) had the greatest steady-state probability. Because the calmodulin binding was slow, each conditional distribution had their respective largest peaks near the slow limit (c∞ for states S1 and S2, c* for S3 and S4) (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 S4 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 (S1 and S2) were bimodel, with peaks near c∞ and c*.Figure12 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 v0 and c* was similar to our simplified model (cf. Figure 5(a)), decreasing in magnitude and approaching 0 as v0 or c* 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 (Figure5(b)). For both the slow CaM and SQS parameters, Δ decreased as c* or v0 increased. Additionally, for both the calmodulin-activated and -inactivated channels, Δ had greater dependence on v0 for the SQS parameters, which is consistent with calmodulin being more sensitive to the bulk concentration (since increasing v0 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 c* and v0 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 c* increased. However, in contrast with the parameter space using SQS parameters, Δ was fairly insensitive to v0 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.


## 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 (Figure1 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−1 to the constant bulk concentration of c∞=0.1 μM [27]. These assumptions lead to the following deterministic model of subspace calcium dynamics:
(1)dcdt=α-β(c-c∞),
where the influx rate α has units of concentration per time (e.g., μM/ms) and the calcium concentration c 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 c (units of μM) in a diadic subspace of volume v (liters). Subspace calcium moves to the bulk passively via diffusion at rate β (given by 0.01 ms−1). Bulk calcium at the concentration c∞=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 byC^∈{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 (v), that is,
(2)C=C^v.
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
(3)0⇌βα-1⇌2βα-2⋯n-1⇌nβα-n⇌(n+1)βα-n+1⋯,
where the index that labels states, n∈{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)α-=v(α+βc∞).


### 2.1.2. Master Equation and Steady-State Probability Distribution

If we writepn(t)=Pr{C^(t)=n}, 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
(5)dp0dt=-α-p0+βp1,dpndt=-(α-+nβ)pn+α-pn-1+(n+1)βpn+1,n=1,2,….
Note that the correspondence between the rate constants in the deterministic (1) and stochastic (3)–(5) models is established by substituting c=c^/v in (1) to find the rate of change of the number of calcium ions in the deterministic model, that is,
(6)dc^dt=α--βc^.
This equation indicates that c^ 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, c^ decreases at rate βc^, a value that is proportional to c^ because each ion has an opportunity to diffuse into the bulk. Consequently, the transition rates leading out of state C^=n in the stochastic model are given by α- for the C^=n to n+1 transitions and βn for the C^=n to n-1 transitions.To find the steady-state probability distribution ofC^, we set the left hand sides of (5) to zero to obtain
(7)nβpn=α-pn-1,n=1,2,…
from which it follows that {pn} is a Poisson distribution with parameter λ=α-/β, that is,
(8)pn=e-λλnn!.


### 2.1.3. Analysis of Concentration Fluctuations

To see how the subspace calcium concentration fluctuations predicted by this minimal model depend on the parametersα, β, c∞, and v, 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
(9)E[C^]=∑n=0∞‍npn=λ=α-β=v(αβ+c∞)=vc*,
where the last equality defines c* as follows:
(10)c*=αβ+c∞.
Using (9) and the fact that C=C^/v implies E[C]=E[C^]/v, we can identify c* as the expected subspace calcium concentration:
(11)E[C]=c*.
Similarly, the steady-state variance of the number of calcium ions in the subspace is
(12)Var[C^]=∑n=0∞‍(n-E[C])2pn=vc*,
and Var[C]=Var[C^]/v2 implies that the variance of the subspace calcium concentration is
(13)Var[C]=c*v.
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/vc* and, similarly,
(14)CV[C]=Var[C]E[C]=Var[C^]/v2E[C^]/v=1vc*.
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 v).Figure2 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 v=v0=10-17 liters and influx rate of α=0.049 μM/ms, the expected calcium concentration is E[C]=c*=α/β+c∞=5 μM, and the expected number of subspace calcium ions is E[C^]=v0c*=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 x-axis scales). The following rows of Figure 2 show that in a subspace three or ten times larger (v=3v0 or 10v0), 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 (c* 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/vc*→0 as v→∞ for fixed c*.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 v0=10-17 liters and subspaces that are 3 and 10 times larger. Parameters: α=0.049 μM/ms, β=0.01 ms−1, c∞=0.1 μM; the steady-state expected subspace calcium concentration is  E[C]=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 volumev. The calcium influx rate is α0 and α1 when calcium is unbound and bound, respectively, and the transition rates between these states are k+c2 and k-, where c is the subspace calcium concentration. Subspace calcium is passively coupled at rate β to the bulk cytosol with constant concentration c∞.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 modelE[C]=c*=α/β+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(15)μq=∑n=0∞‍nqpn.
By conservation of probability, the zeroth moment μ0=1 and the first moment is the expected number of calcium ions in the subspace (9),
(16)μ1=E[C^].
The second moment μ2 is related to the variance of the number of calcium ions via
(17)Var[C^]=μ2-(μ1)2.
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 (dμ0/dt=0) and, furthermore,
(18)dμ1dt=α--βμ1,dμ2dt=α-+(2α-+β)μ1-2βμ2,
where we have used μ0=1. Setting the left hand sides of these equations to zero, we see that steady-state first and second moment are μ1=α-/β and μ2=α-/β+(α-/β)2, consistent with (11) and (17).


## 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 byC^∈{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 (v), that is,
(2)C=C^v.
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
(3)0⇌βα-1⇌2βα-2⋯n-1⇌nβα-n⇌(n+1)βα-n+1⋯,
where the index that labels states, n∈{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)α-=v(α+βc∞).


## 2.1.2. Master Equation and Steady-State Probability Distribution

If we writepn(t)=Pr{C^(t)=n}, 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
(5)dp0dt=-α-p0+βp1,dpndt=-(α-+nβ)pn+α-pn-1+(n+1)βpn+1,n=1,2,….
Note that the correspondence between the rate constants in the deterministic (1) and stochastic (3)–(5) models is established by substituting c=c^/v in (1) to find the rate of change of the number of calcium ions in the deterministic model, that is,
(6)dc^dt=α--βc^.
This equation indicates that c^ 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, c^ decreases at rate βc^, a value that is proportional to c^ because each ion has an opportunity to diffuse into the bulk. Consequently, the transition rates leading out of state C^=n in the stochastic model are given by α- for the C^=n to n+1 transitions and βn for the C^=n to n-1 transitions.To find the steady-state probability distribution ofC^, we set the left hand sides of (5) to zero to obtain
(7)nβpn=α-pn-1,n=1,2,…
from which it follows that {pn} is a Poisson distribution with parameter λ=α-/β, that is,
(8)pn=e-λλnn!.


## 2.1.3. Analysis of Concentration Fluctuations

To see how the subspace calcium concentration fluctuations predicted by this minimal model depend on the parametersα, β, c∞, and v, 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
(9)E[C^]=∑n=0∞‍npn=λ=α-β=v(αβ+c∞)=vc*,
where the last equality defines c* as follows:
(10)c*=αβ+c∞.
Using (9) and the fact that C=C^/v implies E[C]=E[C^]/v, we can identify c* as the expected subspace calcium concentration:
(11)E[C]=c*.
Similarly, the steady-state variance of the number of calcium ions in the subspace is
(12)Var[C^]=∑n=0∞‍(n-E[C])2pn=vc*,
and Var[C]=Var[C^]/v2 implies that the variance of the subspace calcium concentration is
(13)Var[C]=c*v.
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/vc* and, similarly,
(14)CV[C]=Var[C]E[C]=Var[C^]/v2E[C^]/v=1vc*.
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 v).Figure2 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 v=v0=10-17 liters and influx rate of α=0.049 μM/ms, the expected calcium concentration is E[C]=c*=α/β+c∞=5 μM, and the expected number of subspace calcium ions is E[C^]=v0c*=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 x-axis scales). The following rows of Figure 2 show that in a subspace three or ten times larger (v=3v0 or 10v0), 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 (c* 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/vc*→0 as v→∞ for fixed c*.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 v0=10-17 liters and subspaces that are 3 and 10 times larger. Parameters: α=0.049 μM/ms, β=0.01 ms−1, c∞=0.1 μM; the steady-state expected subspace calcium concentration is  E[C]=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 volumev. The calcium influx rate is α0 and α1 when calcium is unbound and bound, respectively, and the transition rates between these states are k+c2 and k-, where c is the subspace calcium concentration. Subspace calcium is passively coupled at rate β to the bulk cytosol with constant concentration c∞.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 modelE[C]=c*=α/β+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(15)μq=∑n=0∞‍nqpn.
By conservation of probability, the zeroth moment μ0=1 and the first moment is the expected number of calcium ions in the subspace (9),
(16)μ1=E[C^].
The second moment μ2 is related to the variance of the number of calcium ions via
(17)Var[C^]=μ2-(μ1)2.
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 (dμ0/dt=0) and, furthermore,
(18)dμ1dt=α--βμ1,dμ2dt=α-+(2α-+β)μ1-2βμ2,
where we have used μ0=1. Setting the left hand sides of these equations to zero, we see that steady-state first and second moment are μ1=α-/β and μ2=α-/β+(α-/β)2, 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 Figure3). We assume that calcium binding instantaneously modifies the conductance of the channel, that is, the rate of calcium influx into the domain is α0 when the channel is calcium-free and α1 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 k+c2 and k-, respectively, (Figure 3, curved arrows). Note that the rate constant k- has units of ms−1, k+ has units of  μM−2 ms−1, and the dissociation constant for calcium binding, denoted by κ, has units of μM and is given by κ2=k-/k+.


### 2.2.1. Stochastic Model

Let us denote the states of the stochastic system by(n,0) and (n,1), where n∈{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,(19)(0,0)⇌βα-0(1,0)⇌2βα-0(2,0)⇌3βα-0(3,0)⇌4βα-0(4,0)⋯2k-+⥯k-6k-+⥯k-12k-+⥯k-(0,1)⇌βα-1(1,1)⇌2βα-1(2,1)⋯,wherek-+=k+/v2. 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 (k+c2=k+c^2/v2). The downward transitions between states (n,0) and (n-2,1) include the combinatorial coefficient, n(n-1), double the number of ways that two indistinguishable calcium ions can be chosen from the n ions in the subspace. This factor of two is required so that the microscopic propensity k-+ agrees with the macroscopic rate k+c2 for large n and v with c=n/v fixed, that is,
(20)n(n-1)k-+=n(n-1)k+v2=k+(c2-cv),
an expression that approaches k+c2 as v→∞ [28].


### 2.2.2. Master Equation

Let us writepn0(t) to indicate the probability that at time t the channel is calcium-free and C^=n. Similarly, pn1(t) is the probability that C^(t)=n 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:
(21)dpn0dt=-[α-0+nβ+n(n-1)k-+]pn0+α-0pn-10+(n+1)βpn+10+k-pn-21,dpn1dt=-[α-1+nβ+k-]pn1+α-1pn-11+(n+1)βpn+11+(n+2)(n+1)k-+pn+20.
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:
(22)μq0/1=∑n=0∞‍nqpn0/1,
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 (μ00+μ01=1). The expected number of calcium ions in the subspace conditioned on the channel being calcium free or bound, respectively, is given by
(23)E0/1[C^]=∑n=0∞npn0/1∑n=0∞pn0/1=μ10/1μ00/1.
Similarly the second moments μ20/1 are related to the conditional variances via
(24)Var0/1[C^]=μ20/1μ00/1-(μ10/1μ00/1)2.


### 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, μ00 and μ01—that is, the probability of the channel being in the calcium free or bound state—are given by
(25)dμ00dt=-k-+μ20+k-+μ10+k-μ01,(26)dμ01dt=k-+μ20-k-+μ10-k-μ01,
where we note that dμ00/dt+dμ01/dt=0 and μ00+μ01=1. In the same way, the equations for the first moments, μ10 and μ11, are found to be
(27)dμ10dt=α-0μ00-βμ10-k-+μ30+k-+μ20+k-μ11+2k-μ01,dμ11dt=α-1μ01-βμ11+k-+μ30-3k-+μ20+2k-+μ10-k-μ11.
Setting the left hand side of (25) to zero, we find that the steady-state probability of a calcium-bound channel is
(28)μ01=μ20-μ10κ2v2=E0[C2]-E0[C]/vκ2+E0[C2]-E0[C]/v,
where in the second equality we have used μ20=μ00E0[C^2]=μ00v2E0[C2].Note that as the volume increases (v→∞), E0[C]/v becomes negligible compared to E0[C2], while E0[C2]→(E0[C])2  as the conditional variance goes to zero (Var0[C]→0). Thus, in the large system limit, the probability that the channel is in the calcium-bound state is given by
(29)μ01=(E0[C])2κ2+(E0[C])2.
In the case of a calcium-activated channel, μ01 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(30)E[C]=μ00E0[C]+μ01E1[C]
and the probability that a calcium-activated channel is open, popen=μ01, 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 {pn0/1} for various model parameters.Figures4(a) and 4(b) show the probability distribution for v=v0 and 8v0 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)α-0=vβc∞<vβc*=v(α+βc∞)=α-1.
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 (E1[C]>E0[C]). Note that the eight-fold increase in system size leads to a significant increase in the channel open probability, that is, popen=μ01 = 0.23 and 0.78 for v=v0 and 8v0, 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 v=v0, bimodal when v=8v0).Subspace volume-dependence of calcium fluctuations and open probability of a calcium-activated channel. Steady-state probability distribution forv = v0 (a) and 8v0 (b) for the calcium-activated channel (κ=2 μM, k+=0.05 μM−2 ms−1) [29, 30]. (c) Steady-state E[C] and popen=μ01 for integer multiples of the unitary volume v0 and different rate constants for calcium binding k+ (0.005 to 0.15 μM−2 ms−1) with κ fixed.
(a)
(b)
(c)Figure4(c) shows the expected calcium concentration E[C] (30) and open probability (popen=μ01) for the calcium-activated channel as a function of subspace volume v and different rate constants for calcium binding k+ (fixed dissociation constant κ). Both E[C] and popen increase with subspace volume v, 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 popen must be greater than the values obtained by assuming channel gating is extremely slow, in which case E[C]≈c∞ and popen≈c∞2/(κ2+c∞2), 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 popen are always less than the values obtained under the assumption of rapid channel binding, values given by simultaneous solution of popen=c2/(κ2+c2) and c=popenα/β+c∞. 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
(32)Δ=E[C]-E[C]∞E[C]∞,
where E[C] is calculated using a system volume of v=v0 and E[C]∞ is the same quantity calculated in the large system size limit (v→∞, numerically estimated using v⩾10v0). Figure 5 shows the small system deviation as a function of unitary subspace volume v0 and influx parameter c* for a calcium-activated channel. In all cases, Δ was negative, meaning that E[C] for v=v0 was suppressed below the large system size limit and increased with volume (cf. Figure 4(c)). For small c*, Δ was near zero. For intermediate values of c* (5–10 μM), the suppression was quite large (Δ≈-60%). As v0 increased, Δ becomes less negative and approaches zero. In general, as c* increased above this range, the suppression ultimately becomes negligible.Percentage small system deviation (Δ, (32)) as a function of unitary subspace volume v0 and influx parameter c* for a single calcium-activated channel (κ=2 μM, k+=0.005 μM−2 ms−1) and calcium-inactivated channel (κ=0.63 μM, k+=0.05 μM−2 ms−1).
(a)
(b)For comparison, Figure5(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 c* (1–3 μM), the magnitude of Δ increased with c*, while above this range Δ was essentially independent of c*. As with the calcium-activated channel, the magnitude of Δ decreased and approached zero as v0 increased.


## 2.2.1. Stochastic Model

Let us denote the states of the stochastic system by(n,0) and (n,1), where n∈{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,(19)(0,0)⇌βα-0(1,0)⇌2βα-0(2,0)⇌3βα-0(3,0)⇌4βα-0(4,0)⋯2k-+⥯k-6k-+⥯k-12k-+⥯k-(0,1)⇌βα-1(1,1)⇌2βα-1(2,1)⋯,wherek-+=k+/v2. 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 (k+c2=k+c^2/v2). The downward transitions between states (n,0) and (n-2,1) include the combinatorial coefficient, n(n-1), double the number of ways that two indistinguishable calcium ions can be chosen from the n ions in the subspace. This factor of two is required so that the microscopic propensity k-+ agrees with the macroscopic rate k+c2 for large n and v with c=n/v fixed, that is,
(20)n(n-1)k-+=n(n-1)k+v2=k+(c2-cv),
an expression that approaches k+c2 as v→∞ [28].


## 2.2.2. Master Equation

Let us writepn0(t) to indicate the probability that at time t the channel is calcium-free and C^=n. Similarly, pn1(t) is the probability that C^(t)=n 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:
(21)dpn0dt=-[α-0+nβ+n(n-1)k-+]pn0+α-0pn-10+(n+1)βpn+10+k-pn-21,dpn1dt=-[α-1+nβ+k-]pn1+α-1pn-11+(n+1)βpn+11+(n+2)(n+1)k-+pn+20.
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:
(22)μq0/1=∑n=0∞‍nqpn0/1,
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 (μ00+μ01=1). The expected number of calcium ions in the subspace conditioned on the channel being calcium free or bound, respectively, is given by
(23)E0/1[C^]=∑n=0∞npn0/1∑n=0∞pn0/1=μ10/1μ00/1.
Similarly the second moments μ20/1 are related to the conditional variances via
(24)Var0/1[C^]=μ20/1μ00/1-(μ10/1μ00/1)2.


## 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, μ00 and μ01—that is, the probability of the channel being in the calcium free or bound state—are given by
(25)dμ00dt=-k-+μ20+k-+μ10+k-μ01,(26)dμ01dt=k-+μ20-k-+μ10-k-μ01,
where we note that dμ00/dt+dμ01/dt=0 and μ00+μ01=1. In the same way, the equations for the first moments, μ10 and μ11, are found to be
(27)dμ10dt=α-0μ00-βμ10-k-+μ30+k-+μ20+k-μ11+2k-μ01,dμ11dt=α-1μ01-βμ11+k-+μ30-3k-+μ20+2k-+μ10-k-μ11.
Setting the left hand side of (25) to zero, we find that the steady-state probability of a calcium-bound channel is
(28)μ01=μ20-μ10κ2v2=E0[C2]-E0[C]/vκ2+E0[C2]-E0[C]/v,
where in the second equality we have used μ20=μ00E0[C^2]=μ00v2E0[C2].Note that as the volume increases (v→∞), E0[C]/v becomes negligible compared to E0[C2], while E0[C2]→(E0[C])2  as the conditional variance goes to zero (Var0[C]→0). Thus, in the large system limit, the probability that the channel is in the calcium-bound state is given by
(29)μ01=(E0[C])2κ2+(E0[C])2.
In the case of a calcium-activated channel, μ01 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(30)E[C]=μ00E0[C]+μ01E1[C]
and the probability that a calcium-activated channel is open, popen=μ01, 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 {pn0/1} for various model parameters.Figures4(a) and 4(b) show the probability distribution for v=v0 and 8v0 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)α-0=vβc∞<vβc*=v(α+βc∞)=α-1.
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 (E1[C]>E0[C]). Note that the eight-fold increase in system size leads to a significant increase in the channel open probability, that is, popen=μ01 = 0.23 and 0.78 for v=v0 and 8v0, 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 v=v0, bimodal when v=8v0).Subspace volume-dependence of calcium fluctuations and open probability of a calcium-activated channel. Steady-state probability distribution forv = v0 (a) and 8v0 (b) for the calcium-activated channel (κ=2 μM, k+=0.05 μM−2 ms−1) [29, 30]. (c) Steady-state E[C] and popen=μ01 for integer multiples of the unitary volume v0 and different rate constants for calcium binding k+ (0.005 to 0.15 μM−2 ms−1) with κ fixed.
(a)
(b)
(c)Figure4(c) shows the expected calcium concentration E[C] (30) and open probability (popen=μ01) for the calcium-activated channel as a function of subspace volume v and different rate constants for calcium binding k+ (fixed dissociation constant κ). Both E[C] and popen increase with subspace volume v, 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 popen must be greater than the values obtained by assuming channel gating is extremely slow, in which case E[C]≈c∞ and popen≈c∞2/(κ2+c∞2), 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 popen are always less than the values obtained under the assumption of rapid channel binding, values given by simultaneous solution of popen=c2/(κ2+c2) and c=popenα/β+c∞. 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
(32)Δ=E[C]-E[C]∞E[C]∞,
where E[C] is calculated using a system volume of v=v0 and E[C]∞ is the same quantity calculated in the large system size limit (v→∞, numerically estimated using v⩾10v0). Figure 5 shows the small system deviation as a function of unitary subspace volume v0 and influx parameter c* for a calcium-activated channel. In all cases, Δ was negative, meaning that E[C] for v=v0 was suppressed below the large system size limit and increased with volume (cf. Figure 4(c)). For small c*, Δ was near zero. For intermediate values of c* (5–10 μM), the suppression was quite large (Δ≈-60%). As v0 increased, Δ becomes less negative and approaches zero. In general, as c* increased above this range, the suppression ultimately becomes negligible.Percentage small system deviation (Δ, (32)) as a function of unitary subspace volume v0 and influx parameter c* for a single calcium-activated channel (κ=2 μM, k+=0.005 μM−2 ms−1) and calcium-inactivated channel (κ=0.63 μM, k+=0.05 μM−2 ms−1).
(a)
(b)For comparison, Figure5(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 c* (1–3 μM), the magnitude of Δ increased with c*, while above this range Δ was essentially independent of c*. As with the calcium-activated channel, the magnitude of Δ decreased and approached zero as v0 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 Figure6). 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 α0 when all channels are calcium-free and α1 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 v, 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 v, 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 ionsC bind to channel B to form the complex C2B, we can write the following kinetic scheme:
(33)2C+B⇌k-k+C2B.
The deterministic ODE system that applies in the case of a large subspace volume is
(34)dcdt=α0bbt+α1bt-bbt-β(c-c∞)-k+c2b+k-(bt-b),dbdt=-k+c2b+k-(bt-b),
where we have written c=[C] and b=[B]. Because the total (calcium-free plus-bound) concentration of channels, bt=[B]+[C2B], is a constant determined by initial conditions, we have eliminated the equation for [C2B]. At steady-state the channels will be in equilibrium with subspace calcium, that is, b/bt=κ2/κ2+c2.Thus, in the case of a calcium-activated channel (α0=0,α1>0), the steady-state calcium concentration satisfies
(35)0=α1c2κ2+c2-β(c-c∞),
while in the case of a calcium-inactivated channel (α0>0, α1=0),
(36)0=α0κ2κ2+c2-β(c-c∞).
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)).Bifurcation diagram showing the steady-state calcium concentrationc 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.
(a)
(b)


### 2.3.2. Stochastic Model

Following the notation developed in the previous section, we writepnm=P{C^=n,C2B^=m,t}=P{C^=n,B^=b^t-m} for n∈{0,1,…,∞} and m∈{0,1,…,b^t} and, where b^t is the total number of channels (for integer ℓ, b^t=ℓ when v=ℓv0). The state-transition diagram for the Markov process (not shown) is analogous to (19) but with b^t+1 rows as opposed to two. The master equation for the dynamics of the calcium channel and subspace calcium concentration is
(37)dpnmdt=-[α-m+nβ+mk-+n(n-1)(b^t-m)k-+]pnm+α-mpn-1m+(n+1)βpn+1m+(m+1)k-pn-2m+1+(n+2)(n+1)(b^t-m+1)k-+pn+2m-1,
where α-m=v(αm+βc∞) and
(38)αm=α0b^t-mb^t+α1mb^t.
In (38), it is understood that terms in the master equation involving negative indices (i.e., n or m<0) evaluate to zero.


### 2.3.3. Concentration Fluctuations

Figure8(A) shows the steady-state probability distribution for v=v0, 2v0 and 4v0 for a calcium-activated channel with dissociation constant chosen so that the deterministic system is monostable (κ=0.45 μM). For v=v0, 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 c∞, while for the open channel, the conditional mean is near c*. For v=2v0 and 4v0, 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 c∞ and c*, 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 forv = v0 (a), 2v0 (b), and 4v0 (c) (κ=0.45 μM, k+=5×10-4 μM−2 ms−1). For each panel, the dashed black line denotes the conditional expected concentration (Em[C]). The steady-state probability distribution is shown for each possible number of closed (NC) and open (NO) channels. (B) (a) Steady-state E[C] and popen for the monostable system as a function of v for rate constants of calcium binding (k+=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 popen for the bistable system (κ=2 μM) as a function of v using k+= 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).Figure8(B)a shows E[C] and popen for subspace volumes v given by different discrete multiples of the unitary volume v0. 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 popen approached the fast/large system limit as v increased. This also occurs for fixed v with increasing k+, that is, the rate constant for calcium binding. For fixed κ, smaller values of k+ 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 popen is more complex (Figure 8(B)b). Interestingly, the small system deviation in this case is often a biphasic function of system volume.Figure9 shows analogous results for calcium-inactivated calcium influx. As with the calcium-activated channel, E[C] and popen were suppressed below the fast/large system limit (Figure 9(B)). Δ is often negative, but became negligible as k+ increased. Similarly, as v increased, both E[C] and popen 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 forv = (a) v0, (b) 2v0, and (c) 4v0 (κ=0.63 μM, k+=0.005 μM−2 ms−1). (B) Steady-state E[C] and popen as a function of v (k+=0.0015 to 0.15 μM−2 ms−1).Figure10 summarizes the dependence of the small system deviation (Δ, (32)) on the unitary subspace volume (v0) and calcium influx parameter (c*) 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 popen is small (i.e., small c* for the calcium-activated channel and large c* for the calcium-inactivated channel). In general, as the unitary volume v0 is increased, there is less suppression compared to the large system size limit.Percentage small system deviation (Δ) as a function of unitary subspace volume v0 and influx parameter c* for multiple (a) calcium-activated channels (κ=0.45 μM, k+=0.005 μM−2 ms−1) and (b) calcium-inactivated channels (κ=0.63 μM, k+=0.05 μM−2 ms−1).
(a)
(b)


## 2.3.1. Deterministic Model

Assuming as before that two free calcium ionsC bind to channel B to form the complex C2B, we can write the following kinetic scheme:
(33)2C+B⇌k-k+C2B.
The deterministic ODE system that applies in the case of a large subspace volume is
(34)dcdt=α0bbt+α1bt-bbt-β(c-c∞)-k+c2b+k-(bt-b),dbdt=-k+c2b+k-(bt-b),
where we have written c=[C] and b=[B]. Because the total (calcium-free plus-bound) concentration of channels, bt=[B]+[C2B], is a constant determined by initial conditions, we have eliminated the equation for [C2B]. At steady-state the channels will be in equilibrium with subspace calcium, that is, b/bt=κ2/κ2+c2.Thus, in the case of a calcium-activated channel (α0=0,α1>0), the steady-state calcium concentration satisfies
(35)0=α1c2κ2+c2-β(c-c∞),
while in the case of a calcium-inactivated channel (α0>0, α1=0),
(36)0=α0κ2κ2+c2-β(c-c∞).
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)).Bifurcation diagram showing the steady-state calcium concentrationc 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.
(a)
(b)


## 2.3.2. Stochastic Model

Following the notation developed in the previous section, we writepnm=P{C^=n,C2B^=m,t}=P{C^=n,B^=b^t-m} for n∈{0,1,…,∞} and m∈{0,1,…,b^t} and, where b^t is the total number of channels (for integer ℓ, b^t=ℓ when v=ℓv0). The state-transition diagram for the Markov process (not shown) is analogous to (19) but with b^t+1 rows as opposed to two. The master equation for the dynamics of the calcium channel and subspace calcium concentration is
(37)dpnmdt=-[α-m+nβ+mk-+n(n-1)(b^t-m)k-+]pnm+α-mpn-1m+(n+1)βpn+1m+(m+1)k-pn-2m+1+(n+2)(n+1)(b^t-m+1)k-+pn+2m-1,
where α-m=v(αm+βc∞) and
(38)αm=α0b^t-mb^t+α1mb^t.
In (38), it is understood that terms in the master equation involving negative indices (i.e., n or m<0) evaluate to zero.


## 2.3.3. Concentration Fluctuations

Figure8(A) shows the steady-state probability distribution for v=v0, 2v0 and 4v0 for a calcium-activated channel with dissociation constant chosen so that the deterministic system is monostable (κ=0.45 μM). For v=v0, 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 c∞, while for the open channel, the conditional mean is near c*. For v=2v0 and 4v0, 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 c∞ and c*, 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 forv = v0 (a), 2v0 (b), and 4v0 (c) (κ=0.45 μM, k+=5×10-4 μM−2 ms−1). For each panel, the dashed black line denotes the conditional expected concentration (Em[C]). The steady-state probability distribution is shown for each possible number of closed (NC) and open (NO) channels. (B) (a) Steady-state E[C] and popen for the monostable system as a function of v for rate constants of calcium binding (k+=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 popen for the bistable system (κ=2 μM) as a function of v using k+= 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).Figure8(B)a shows E[C] and popen for subspace volumes v given by different discrete multiples of the unitary volume v0. 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 popen approached the fast/large system limit as v increased. This also occurs for fixed v with increasing k+, that is, the rate constant for calcium binding. For fixed κ, smaller values of k+ 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 popen is more complex (Figure 8(B)b). Interestingly, the small system deviation in this case is often a biphasic function of system volume.Figure9 shows analogous results for calcium-inactivated calcium influx. As with the calcium-activated channel, E[C] and popen were suppressed below the fast/large system limit (Figure 9(B)). Δ is often negative, but became negligible as k+ increased. Similarly, as v increased, both E[C] and popen 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 forv = (a) v0, (b) 2v0, and (c) 4v0 (κ=0.63 μM, k+=0.005 μM−2 ms−1). (B) Steady-state E[C] and popen as a function of v (k+=0.0015 to 0.15 μM−2 ms−1).Figure10 summarizes the dependence of the small system deviation (Δ, (32)) on the unitary subspace volume (v0) and calcium influx parameter (c*) 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 popen is small (i.e., small c* for the calcium-activated channel and large c* for the calcium-inactivated channel). In general, as the unitary volume v0 is increased, there is less suppression compared to the large system size limit.Percentage small system deviation (Δ) as a function of unitary subspace volume v0 and influx parameter c* for multiple (a) calcium-activated channels (κ=0.45 μM, k+=0.005 μM−2 ms−1) and (b) calcium-inactivated channels (κ=0.63 μM, k+=0.05 μM−2 ms−1).
(a)
(b)


## 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(39)S1⇌δ-γ-S2⇌k-k+c2S3⇌δ+γ+S4,
where states S1 and S2 are calcium-free, states S3 and S4 are calcium-bound, and state S4 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:
(40)dpn1dt=-(α-0+nβ+γ-)pn1+α-0pn-11+(n+1)βpn+11+δ-pn2,dpn2dt=-[α-0+nβ+δ-+n(n-1)k-+]pn2+α-0pn-12+(n+1)βpn+12+γ-pn1+k-pn-23,dpn3dt=-(α-0+nβ+γ+k-)pn3+α-0pn-13+(n+1)βpn+13+δ+pn4+(n+2)(n+1)k-+pn+22,dpn4dt=-(α-1+nβ+δ+)pn4+α-1pn-14+(n+1)βpn+14+γ+pn3,
where for a calmodulin-activated channel α-0 and α-1 are given by (31).Figure11(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 (S4) had the greatest steady-state probability. Because the calmodulin binding was slow, each conditional distribution had their respective largest peaks near the slow limit (c∞ for states S1 and S2, c* for S3 and S4) (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 S4 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 (S1 and S2) were bimodel, with peaks near c∞ and c*.Figure12 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 v0 and c* was similar to our simplified model (cf. Figure 5(a)), decreasing in magnitude and approaching 0 as v0 or c* 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 (Figure5(b)). For both the slow CaM and SQS parameters, Δ decreased as c* or v0 increased. Additionally, for both the calmodulin-activated and -inactivated channels, Δ had greater dependence on v0 for the SQS parameters, which is consistent with calmodulin being more sensitive to the bulk concentration (since increasing v0 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 c* and v0 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 c* increased. However, in contrast with the parameter space using SQS parameters, Δ was fairly insensitive to v0 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 v0 and c*), 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 (Figure8). 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 v=v0 (Figure 4(a)) but was present for v=8v0 (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 Section2.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 onstationary 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.


## 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 v0 and c*), 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 (Figure8). 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 v=v0 (Figure 4(a)) but was present for v=8v0 (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 Section2.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 onstationary 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.



---

*Source: 897371-2012-12-29.xml*