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

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.


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 2 Computational and Mathematical Methods in Medicine dynamics [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][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 calciumregulated calcium influx and sarcoplasmic reticulum calcium release in cardiac myocytes [24][25][26], we sought to determine precisely when the conventional deterministic formulation of these processes are a valid approximation. When is it appropriate to model the dynamics of subspace calcium using deterministic ODEs? When does one require a stochastic description that accounts for the fundamentally discrete nature of calcium-regulated calcium influx?
To answer this question, we constructed and analyzed a minimal Markov model of a calcium-regulated calcium channel and associated subspace. We compared the expected steady-state subspace calcium concentration in this stochastic model (a result that accounts for the small subspace volume) with the result obtained using the corresponding deterministic ODE model (an approximation that assumes large system size). Section 2.1 introduces our model formulation and shows the agreement between deterministic and stochastic descriptions when subspace calcium does not regulate calcium influx. However, when calcium binding regulates channel activity (through either activation or inactivation), the deterministic and stochastic descriptions often disagree (Sections 2.2 and 2.3). In general, the effect of concentration fluctuations in a spatially restricted calcium domain with a calcium-regulated calcium influx pathway (e.g., a stochastically gating L-type calcium channel) is only well-approximated by the deterministic description when the subspace volume is sufficiently (unphysiologically) large or the kinetics of calcium binding to the calciumregulated 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).

Calcium Influx and Subspace Calcium Concentration
Fluctuations. We begin with the case of a single calcium channel that is associated with a spatially restricted subspace but not regulated by subspace calcium (Figure 1 and Section 2.1). The description of the model in the absence  of calcium regulation simplifies the initial presentation of the model and allows us to illustrate general properties of subspace calcium concentration fluctuations. Subsequently, we present a more complete model formulation that includes calcium-regulated calcium influx ( Figure 3 and Section 2.2). For simplicity, we neglect the presence of endogenous calcium binding proteins and assume a constant flux of calcium, denoted by α, into the subspace. The subspace calcium concentration is passively coupled with relaxation rate β 0.01 ms 1 to the constant bulk concentration of c 0.1 μM [27]. These assumptions lead to the following deterministic model of subspace calcium dynamics: 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.

Stochastic Model.
In the corresponding stochastic description of calcium influx into a diadic subspace, the state variable is the number of calcium ions in the subspace (a discrete quantity that we will denote by þ C è0, 1, . . . , í, where the caret (hat) indicates that þ C is a dimensionless number of molecules rather than concentration, and the capitalization indicates a random variable. The fluctuating subspace calcium concentration (also a random variable, denoted by C) depends on both þ C and the subspace volume (v), that is, Using this relationship, it is straightforward to derive the transition rates between the discrete states of the stochastic Computational and Mathematical Methods in Medicine 3 model that are consistent with (1). The resulting state-transition diagram for the stochastic model is 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,

Master Equation and
Steady-State Probability Distribution. If we write p n 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 dp 0 dt αp 0 βp 1 , dp n dt α nβ p n αp n 1 n 1 βp n 1 , n 1, 2, . . . . (5) 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, 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 of þ C, we set the left hand sides of (5) to zero to obtain nβp n αp n 1 , n 1, 2, . . . (7) from which it follows that èp n í is a Poisson distribution with parameter λ α β, that is, p n e λ λ n n! .
where the last equality defines c as follows: 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: Similarly, the steady-state variance of the number of calcium ions in the subspace is and Var C¥ Var þ C¥ v 2 implies that the variance of the subspace calcium concentration is 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, 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). Figure 2 illustrates fluctuation amplitudes in the minimal subspace model by plotting the steady-state probability distribution of þ C and C (left and right columns, resp.). In the first row, using subspace volume of v v 0 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¥ v 0 c 30. In both cases the coefficient of variation is 1 ö 30 0.18 (the spread of the distributions as illustrated is due to the different x-axis scales). The following rows of Figure 2 show that in a subspace three or ten times larger (v 3v 0 or 10v 0 ), 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  following sense: the expected value of the fluctuating calcium concentration in the stochastic model E 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).

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 By conservation of probability, the zeroth moment μ 0 1 and the first moment is the expected number of calcium ions in the subspace (9), The second moment μ 2 is related to the variance of the number of calcium ions via 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, 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).

Stochastic Subspace Model with Calcium Regulation.
This section augments the subspace model presented above to include calcium regulation of a calcium channel (see Figure 3). We assume that calcium binding instantaneously modifies the conductance of the channel, that is, the rate of calcium influx into the domain is α 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 c 2 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 Figure 3: Diagram of the components and fluxes in a subspace model that includes calcium-regulated calcium influx. A single calcium channel (with two calcium binding sites) is associated with a subspace of volume v. The calcium influx rate is α0 and α1 when calcium is unbound and bound, respectively, and the transition rates between these states are k c 2 and k , where c is the subspace calcium concentration. Subspace calcium is passively coupled at rate β to the bulk cytosol with constant concentration c .

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 statetransition diagram for the stochastic subspace model with calcium influx, where k k v 2 . 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 c 2 k þ c 2 v 2 ). 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 c 2 for large n and v with c n v fixed, that is, an expression that approaches k c 2 as v [28].

Master Equation.
Let us write p 0 n t to indicate the probability that at time t the channel is calcium-free and þ C n. Similarly, p 1 n 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: dp 0 n dt ¢α 0 nβ n n 1 k §p 0 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: where the superscript 0 1 indicates either index occurring on both the left and right hand sides of the equality. Note ). The expected number of calcium ions in the subspace conditioned on the channel being calcium free or bound, respectively, is given by Similarly the second moments μ 0 1 2 are related to the conditional variances via (22) with respect to time and substituting for the time derivatives using (21), it can be shown that the time-derivatives of the zerothmoments, μ 0 0 and μ 1 0 -that is, the probability of the channel being in the calcium free or bound state-are given by

Moment Calculation. By differentiating
where we note that dμ 0 0 dt dμ 1 0 dt 0 and μ 0 0 μ 1 0 1. In the same way, the equations for the first moments, μ 0 1 and μ 1 1 , are found to be Setting the left hand side of (25) to zero, we find that the steady-state probability of a calcium-bound channel is where in the second equality we have used μ 0 Note that as the volume increases (v ), E 0 C¥ v becomes negligible compared to E 0 C 2 ¥, while E 0 C 2 ¥ E 0 C¥ 2 as the conditional variance goes to zero (Var 0 C¥ 0). Thus, in the large system limit, the probability that the channel is in the calcium-bound state is given by In the case of a calcium-activated channel, μ 1 0 is the open probability.

Analysis of Concentration Fluctuations.
The moment analysis in the previous section suggests that the expected calcium concentration in the subspace given by (30) and the probability that a calcium-activated channel is open, p open μ 1 0 , 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 èp show the probability distribution for v v 0 and 8v 0 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,  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 where E C¥ is calculated using a system volume of v v 0 and E C¥ is the same quantity calculated in the large system size limit (v , numerically estimated using v 10v 0 ).  Figure 5 shows the small system deviation as a function of unitary subspace volume v 0 and influx parameter c for a calcium-activated channel. In all cases, Δ was negative, meaning that E C¥ for v v 0 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 v 0 increased, Δ becomes less negative and approaches zero. In general, as c increased above this range, the suppression ultimately becomes negligible. For comparison, Figure 5(b) shows the small system deviation for a calcium-inactivated channel. In general, the magnitude of Δ for the calcium-inactivated channel was smaller than the calcium-activated channel. For small 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 v 0 increased.

Calcium Regulation of Multiple Channels.
The previous section analyzed the effect of subspace volume when the influx pathway involves calcium regulation of a single channel. In this section, we assume that the total number of channels increases with subspace volume (see Figure 6). As before, we assume that calcium binding instantaneously  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. 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.

Deterministic Model.
Assuming as before that two free calcium ions C bind to channel B to form the complex C 2 B, we can write the following kinetic scheme: The deterministic ODE system that applies in the case of a large subspace volume is where we have written c C¥ and b B¥. Because the total (calcium-free plus-bound) concentration of channels, b t B¥ C 2 B¥, is a constant determined by initial conditions, we have eliminated the equation for C 2 B¥. At steady-state the channels will be in equilibrium with subspace calcium, that is, b b t κ 2 κ 2 c 2 .Thus, in the case of a calciumactivated channel α 0 0, α 1 0 , the steady-state calcium concentration satisfies while in the case of a calcium-inactivated channel (α 0 0, (36) Figure 7 shows bifurcation diagrams for the steady-state calcium concentration in both cases. For the calciumactivated channel there is a range of κ that leads to bistability (Figure 7(a)), while no bistable regime exists for a calciuminactivated channel (Figure 7(b)).

Stochastic Model.
Following the notation developed in the previous section, we write p m n Pè þ C n, C 2 B 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 v 0 ). 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 where α m v α m βc and In (38), it is understood that terms in the master equation involving negative indices (i.e., n or m 0) evaluate to zero. Figure 8(A) shows the steady-state probability distribution for v v 0 , 2v 0 and 4v 0 for a calcium-activated channel with dissociation constant chosen so that the deterministic system is monostable (κ 0.45 μM). For v v 0 , 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 2v 0 and 4v 0 , 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(B)a shows E C¥ and p open for subspace volumes v given by different discrete multiples of the unitary volume v 0 . 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 Computational and Mathematical Methods in Medicine  expected calcium concentration and open probability for a small subspace as compared to the large system limit (35). As expected, both E C¥ and p open 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 p open is more complex (Figure 8(B)b). Interestingly, the small system deviation in this case is often a biphasic function of system volume. Figure 9 shows analogous results for calcium-inactivated calcium influx. As with the calcium-activated channel, E C¥ and p open 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 p open approached the fast/large system limit. Figure 10 summarizes the dependence of the small system deviation (Δ, (32)) on the unitary subspace volume (v 0 ) 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 p open is small (i.e., small c for the calciumactivated channel and large c for the calcium-inactivated channel). In general, as the unitary volume v 0 is increased, there is less suppression compared to the large system size limit.

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 calciumactivated 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 calmodulinmediated 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.
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 where states S 1 and S 2 are calcium-free, states S 3 and S 4 are calcium-bound, and state S 4 determines the fraction of channels activated (or inactivated) by calmodulin. When it is assumed that a single calmodulin molecule is colocalized . 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).  with the calcium channel (as in Section 2.2), the master equation takes the following form: where for a calmodulin-activated channel α 0 and α 1 are given by (31). Figure 11(B) shows the steady-state probability distribution numerically calculated from these equations. Using a parameter set referred to as "slow CaM," Tadross et al. showed that calmodulin was primarily sensitive to the local subspace calcium level (representing the C-lobe) when calcium binding to calmodulin was slow [31]. With "slow CaM" parameters, we found that calmodulin bound to the effector site (S 4 ) 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 S 1 and S 2 , c for S 3 and S 4 ) (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 S 4 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 (S 1 and S 2 ) were bimodel, with peaks near c and c . Figure 12 shows the small system size suppression Δ for both slow CaM and SQS parameter sets assuming a single calmodulin-regulated channel. As in our simplified model ( Figure 5), Δ was quite large in magnitude for some conditions (up to 30% suppression). For the calmodulin-activated channel, the dependence of Δ on v 0 and c was similar to our simplified model (cf. Figure 5(a)), decreasing in magnitude and approaching 0 as v 0 or c increased. The parameter space for the calmodulin-inactivated channel differed somewhat from our simplified model ( Figure 5(b)). For both the slow CaM and SQS parameters, Δ decreased as c or v 0 increased. Additionally, for both the calmodulin-activated and -inactivated channels, Δ had greater dependence on v 0 for the SQS parameters, which is consistent with calmodulin being more sensitive to the bulk concentration (since increasing v 0 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 v 0 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 v 0 using slow CaM parameters, which is consistent with calmodulin being, in this case, less sensitive to the bulk calcium concentration.

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 calmodulinmediated calcium channel regulation.

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  Figure 13: Small system deviation for multiple calmodulin-activated and -inactivated channels, using "slow CaM" and "SQS" parameters.
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 v 0 and c ), subspace volume can greatly influence the steady-state properties of stochastically gating channels (Figures 10 and 13).

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 zerothand 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 secondorder reactions, the expected concentration computed from the chemical master equation may deviate from the corresponding ODE model [20]. The deviation is typically negative (i.e., Δ 0), with greatest suppression when concentration fluctuations are large. Other recent studies have demonstrated that the concentrations of species in stochastic biochemical networks can deviate from deterministic ODE descriptions. In agreement with our findings, the deviation is often negative [18,21,28], although positive deviation was observed in some biochemical systems [12,49]. Our findings are consistent with Darvey et al., in that Δ had the greatest magnitude when either the subspace volume or calcium influx rate was small ( Figure 10) (both result in larger concentration fluctuations, see (14)).
We found that the small system size deviation was particularly complex in cases where the deterministic ODE, that is, the model appropriate for the large system size limit, is bistable ( Figure 8). Lestas et al. recently investigated bistability/bimodality in a network of gene regulation and demonstrated that bistability in the deterministic ODE model did not imply bimodality in the discrete system and, conversely, bimodality in the discrete system did not imply bistability in the corresponding ODEs [50]. We obtained similar results, as bimodality in distribution was not present in the bistable system for v v 0 (Figure 4(a)) but was present for v 8v 0 (Figure 4(b)). Conversely, for the calmodulinregulated 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).

Limitations.
The two-state kinetic models of the calcium channel introduced in Section 2.1 is minimal and should be interpreted as phenomenological (as opposed to statistical) model of single channel kinetics, that is, the topology and parameters of this model were not obtained by fitting to patch clamp recordings [51]. On the other hand, the kinetic model for regulation of the calcium channel presented in Section 2.4 is state-of-the-art. Both minimal and physiologically realistic channel models are affected by the decision to account for (or neglect) fluctuations in calcium concentration that result from the small number of ions in the subspace.
The most significant limitation in the model formulation is our neglect of spatial dynamics of calcium diffusion within the dyadic subspace and the details of the spatial arrangement of the ryanodine receptors [47,52]. However, for the purposes of the present study, that is, investigation of the influence of concentration fluctuations on the regulation of calcium influx, a nonspatial Markov chain model that includes subspace volume as a model parameter and accounts for the finite number of calcium ions in the domain is sufficient.
Another limitation of the present work is that we focus on stationary statistics, for example, the expected value of the steady-state subspace calcium concentration, in our analysis of the deviation of continuous ODE description from the discrete stochastic formulation. Future studies could address how transient dynamics, for example, the cellular response to a depolarizing voltage step, excitation-contraction coupling gain, and so forth, are affected by calcium concentration fluctuations resulting from small subspace volume.

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.