Computational and Mathematical Methods in Medicine

Volume 2012 (2012), Article ID 897371, 17 pages

doi:10.1155/2012/897371

Research Article

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

Seth H. Weinberg and Gregory D. Smith

Department of Applied Science, The College of William and Mary, Williamsburg, VA 23187, USA

Received 29 June 2012; Accepted 17 September 2012

Academic Editor: Ling Xia

Copyright © 2012 Seth H. Weinberg and Gregory D. Smith. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

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

1. Introduction

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

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

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

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

2. Methods and Results

2.1. Calcium Influx and Subspace Calcium Concentration Fluctuations

We begin with the case of a single calcium channel that is associated with a spatially restricted subspace but not regulated by subspace calcium (Figure 1 and Section 2.1). The description of the model in the absence of calcium regulation simplifies the initial presentation of the model and allows us to illustrate general properties of subspace calcium concentration fluctuations. Subsequently, we present a more complete model formulation that includes calcium-regulated calcium influx (Figure 3 and Section 2.2). For simplicity, we neglect the presence of endogenous calcium binding proteins and assume a constant flux of calcium, denoted by , into the subspace. The subspace calcium concentration is passively coupled with relaxation rate  ms⁻¹ to the constant bulk concentration of  μ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 is a continuous real-valued quantity.

Figure 1: Diagram of the components and fluxes in a minimal subspace model. Calcium influx (in units of μM/s) leads to increased calcium concentration (units of μM) in a diadic subspace of volume (liters). Subspace calcium moves to the bulk passively via diffusion at rate (given by 0.01 ms⁻¹). Bulk calcium at the concentration  μM returns to the subspace at the same rate. The equilibration time of subspace calcium is ms [27].

2.1.1. Stochastic Model

In the corresponding stochastic description of calcium influx into a diadic subspace, the state variable is the number of calcium ions in the subspace (a discrete quantity that we will denote by , where the caret (hat) indicates that 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 ) depends on both and the subspace volume (), that is,

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

where the index that labels states, , ranges over all possible numbers of calcium ions in the subspace and the constant is proportional to the subspace volume, that is,

2.1.2. Master Equation and Steady-State Probability Distribution

If we write , 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

Note that the correspondence between the rate constants in the deterministic (1) and stochastic (3)–(5) models is established by substituting in (1) to find the rate of change of the number of calcium ions in the deterministic model, that is,

This equation indicates that increases at rate (due to influx and diffusion from the bulk), a value that is independent of the number of calcium ions in the subspace. At the same time, decreases at rate , a value that is proportional to because each ion has an opportunity to diffuse into the bulk. Consequently, the transition rates leading out of state in the stochastic model are given by for the to transitions and for the to transitions.

To find the steady-state probability distribution of , we set the left hand sides of (5) to zero to obtain

from which it follows that is a Poisson distribution with parameter , that is,

2.1.3. Analysis of Concentration Fluctuations

To see how the subspace calcium concentration fluctuations predicted by this minimal model depend on the parameters , , , and , recall that the mean and variance of the Poisson distribution (8) is equal to the parameter and, consequently, the steady-state expected number of calcium ions in the subspace is given by

where the last equality defines as follows:

Using (9) and the fact that implies , we can identify as the expected subspace calcium concentration:

Similarly, the steady-state variance of the number of calcium ions in the subspace is

and implies that the variance of the subspace calcium concentration is

Note that the coefficient of variation of and are identical and inversely proportional to subspace volume, that is, 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 ).

Figure 2 illustrates fluctuation amplitudes in the minimal subspace model by plotting the steady-state probability distribution of and (left and right columns, resp.). In the first row, using subspace volume of liters and influx rate of  μM/ms, the expected calcium concentration is  μM, and the expected number of subspace calcium ions is . In both cases the coefficient of variation is (the spread of the distributions as illustrated is due to the different -axis scales). The following rows of Figure 2 show that in a subspace three or ten times larger ( or ), the coefficient of variation drops to 0.11 and 0.058, respectively, when the calcium influx rate is scaled to result in the same expected calcium concentration ( fixed, see (14)). As might be expected, concentration fluctuations in the stochastic model are more pronounced for small volumes and become negligible for large volumes, because as for fixed .

Figure 2: Steady-state probability distribution of the number of calcium ions (, left column) and subspace calcium concentration (, right column) for subspace volume of liters and subspaces that are 3 and 10 times larger. Parameters:  μM/ms,  ms⁻¹,  μM; the steady-state expected subspace calcium concentration is  μM.

Figure 3: Diagram of the components and fluxes in a subspace model that includes calcium-regulated calcium influx. A single calcium channel (with two calcium binding sites) is associated with a subspace of volume . The calcium influx rate is and when calcium is unbound and bound, respectively, and the transition rates between these states are and , where is the subspace calcium concentration. Subspace calcium is passively coupled at rate to the bulk cytosol with constant concentration .

Most importantly, the deterministic and stochastic descriptions of this minimal subspace model agree in the following sense: the expected value of the fluctuating calcium concentration in the stochastic model 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

By conservation of probability, the zeroth moment and the first moment is the expected number of calcium ions in the subspace (9),

The second moment 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 () and, furthermore,

where we have used . Setting the left hand sides of these equations to zero, we see that steady-state first and second moment are and , consistent with (11) and (17).

2.2. Stochastic Subspace Model with Calcium Regulation

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

2.2.1. Stochastic Model

Let us denote the states of the stochastic system by and , where 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,

where . The rate of calcium binding to the channel in the stochastic model is inversely proportional to the square of the volume, because of the concentration dependence of the association reaction (). The downward transitions between states and include the combinatorial coefficient, , double the number of ways that two indistinguishable calcium ions can be chosen from the ions in the subspace. This factor of two is required so that the microscopic propensity agrees with the macroscopic rate for large and with fixed, that is,

an expression that approaches as [28].

2.2.2. Master Equation

Let us write to indicate the probability that at time the channel is calcium-free and . Similarly, is the probability that 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:

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 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 (). 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 are related to the conditional variances via

2.2.3. Moment Calculation

By differentiating (22) with respect to time and substituting for the time derivatives using (21), it can be shown that the time-derivatives of the zeroth-moments, and —that is, the probability of the channel being in the calcium free or bound state—are given by

where we note that and . In the same way, the equations for the first moments, and , 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 .

Note that as the volume increases (), becomes negligible compared to , while as the conditional variance goes to zero (). 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, 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

and the probability that a calcium-activated channel is open, , may depend on the subspace volume. In order to analyze the effect of small system size and concentration fluctuations at steady-state, we integrated (21) and determined the probability distributions for various model parameters.

Figures 4(a) and 4(b) show the probability distribution for and 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,

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 (). Note that the eight-fold increase in system size leads to a significant increase in the channel open probability, that is, = 0.23 and 0.78 for and , respectively. Thus, the open probability of the channel is significantly influenced by the subspace volume, in spite of the fact that the calcium influx rate is scaled so that in the absence of calcium-regulation there is no effect of volume ( constant, as in Section 2.1). Comparison of Figures 4(a) and 4(b) also shows a qualitative change in the probability distribution of the subspace calcium concentration (unimodal when , bimodal when ).

(a)

(b)

(c)

Figure 4: Subspace volume-dependence of calcium fluctuations and open probability of a calcium-activated channel. Steady-state probability distribution for = (a) and (b) for the calcium-activated channel ( μM,  μM⁻² ms⁻¹) [29, 30]. (c) Steady-state and for integer multiples of the unitary volume and different rate constants for calcium binding (0.005 to 0.15 μM⁻² ms⁻¹) with fixed.

Figure 4(c) shows the expected calcium concentration (30) and open probability () for the calcium-activated channel as a function of subspace volume and different rate constants for calcium binding (fixed dissociation constant ). Both and increase with subspace volume , that is, the restricted volume of a physiological subspace leads to an open probability and expected calcium concentration that is less than predicted in the corresponding (approximate) continuous description.

Both the open probability and expected calcium concentration asymptotically approach values in a range that are easily precalculated. For example, and must be greater than the values obtained by assuming channel gating is extremely slow, in which case and , because the transition from the calcium-free to -bound channel usually occurs with a subspace that is equilibrated with bulk calcium. In addition, and are always less than the values obtained under the assumption of rapid channel binding, values given by simultaneous solution of and . 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

where is calculated using a system volume of and is the same quantity calculated in the large system size limit (, numerically estimated using ). Figure 5 shows the small system deviation as a function of unitary subspace volume and influx parameter for a calcium-activated channel. In all cases, was negative, meaning that for was suppressed below the large system size limit and increased with volume (cf. Figure 4(c)). For small , was near zero. For intermediate values of (5–10 M), the suppression was quite large (). As increased, becomes less negative and approaches zero. In general, as increased above this range, the suppression ultimately becomes negligible.

(a)

(b)

Figure 5: Percentage small system deviation (, (32)) as a function of unitary subspace volume and influx parameter for a single calcium-activated channel ( μM,  μM⁻² ms⁻¹) and calcium-inactivated channel ( μM,  μM⁻² ms⁻¹).

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

2.3. Calcium Regulation of Multiple Channels

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

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

2.3.1. Deterministic Model

Assuming as before that two free calcium ions bind to channel to form the complex , 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 and . Because the total (calcium-free plus-bound) concentration of channels, , is a constant determined by initial conditions, we have eliminated the equation for . At steady-state the channels will be in equilibrium with subspace calcium, that is, .Thus, in the case of a calcium-activated channel , the steady-state calcium concentration satisfies

while in the case of a calcium-inactivated channel (, ),

Figure 7 shows bifurcation diagrams for the steady-state calcium concentration in both cases. For the calcium-activated channel there is a range of that leads to bistability (Figure 7(a)), while no bistable regime exists for a calcium-inactivated channel (Figure 7(b)).

(a)

(b)

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

2.3.2. Stochastic Model

Following the notation developed in the previous section, we write for and and, where is the total number of channels (for integer , when ). The state-transition diagram for the Markov process (not shown) is analogous to (19) but with rows as opposed to two. The master equation for the dynamics of the calcium channel and subspace calcium concentration is

where and

In (38), it is understood that terms in the master equation involving negative indices (i.e., or ) evaluate to zero.

2.3.3. Concentration Fluctuations

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

Figure 8: Subspace volume-dependence of concentration fluctuation and channel open probability for multiple calcium-activated channels. (A) Steady-state probability distribution for = (a), (b), and (c) ( μM,  μM⁻² ms⁻¹). For each panel, the dashed black line denotes the conditional expected concentration (). The steady-state probability distribution is shown for each possible number of closed () and open () channels. (B) (a) Steady-state and for the monostable system as a function of for rate constants of calcium binding ( to  μM⁻² ms⁻¹). The fast/large and slow system limits are shown in red and blue, respectively. (b) Steady-state and for the bistable system ( μM) as a function of using = 0.005 to 0.015 μM⁻² ms⁻¹. In the bistable system, the larger of the two stable equilibrium (large system limit) is shown in red. The smaller equilibrium is approximately equal to the slow system limit (shown in blue).

Figure 8(B)a shows and for subspace volumes given by different discrete multiples of the unitary volume . Using parameters that lead to a monostable deterministic ODE system, we find, similar to the case of the single channel (Figure 4), a significant deviation between the expected calcium concentration and open probability for a small subspace as compared to the large system limit (35). As expected, both and approached the fast/large system limit as increased. This also occurs for fixed with increasing , that is, the rate constant for calcium binding. For fixed , smaller values of can cause to approach , 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 and 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, and were suppressed below the fast/large system limit (Figure 9(B)). is often negative, but became negligible as increased. Similarly, as increased, both and approached the fast/large system limit.

Figure 9: Subspace volume-dependence of concentration fluctuation and channel open probability for multiple calcium-inactivated channels. (A) Steady-state probability distribution for = (a) , (b) , and (c) ( μM,  μM⁻² ms⁻¹). (B) Steady-state and as a function of ( to 0.15 μM⁻² ms⁻¹).

Figure 10 summarizes the dependence of the small system deviation (, (32)) on the unitary subspace volume () and calcium influx parameter () for the scaling that involves multiple calcium-activated and -inactivated channels. In all cases, was negative, meaning that 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 is small (i.e., small for the calcium-activated channel and large for the calcium-inactivated channel). In general, as the unitary volume is increased, there is less suppression compared to the large system size limit.

(a)

(b)

Figure 10: Percentage small system deviation () as a function of unitary subspace volume and influx parameter for multiple (a) calcium-activated channels ( μM,  μM⁻² ms⁻¹) and (b) calcium-inactivated channels ( μM,  μM⁻² ms⁻¹).

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

where states and are calcium-free, states and are calcium-bound, and state determines the fraction of channels activated (or inactivated) by calmodulin. When it is assumed that a single calmodulin molecule is colocalized with the calcium channel (as in Section 2.2), the master equation takes the following form:

where for a calmodulin-activated channel and are given by (31).

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

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

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

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

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

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

3. Discussion

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

3.1. Physiological Implications

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

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

3.2. Relation to Prior Studies

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

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

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

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

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

3.3. Limitations

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

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

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

4. Conclusions

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

Conflict of Interests

The authors do not have any conflict of interests.

Acknowledgments

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

References

T. J. Kamp and J. W. Hell, “Regulation of cardiac L-type calcium channels by protein kinase A and protein kinase C,” Circulation Research, vol. 87, no. 12, pp. 1095–1102, 2000.

W. C. Rose, C. W. Balke, W. G. Wier, and E. Marban, “Macroscopic and unitary properties of physiological ion flux through L-type Ca2+ channels in guinea-pig heart cells,” Journal of Physiology, vol. 456, pp. 267–284, 1992.

J. L. Greenstein and R. L. Winslow, “An integrative model of the cardiac ventricular myocyte incorporating local control of Ca2+ release,” Biophysical Journal, vol. 83, no. 6, pp. 2918–2945, 2002.

V. Iyer, R. Mazhari, and R. L. Winslow, “A computational model of the human left-ventricular epicardial myocyte,” Biophysical Journal, vol. 87, no. 3, pp. 1507–1525, 2004.

D. A. McQuarrie, “Kinetics of small systems. I,” The Journal of Chemical Physics, vol. 38, pp. 433–436, 1963.

J. Liang and H. Qian, “Computational cellular dynamics based on the chemical master equation: a challenge for understanding complexity,” Journal of Computer Science and Technology, vol. 25, no. 1, pp. 154–168, 2009.

H. Qian, “Nonlinear stochastic dynamics of mesoscopic homogeneous biochemical reaction systems—an analytical theory,” Nonlinearity, vol. 24, no. 6, pp. R19–R49, 2011.

H. Qian and L. M. Bishop, “The chemical master equation approach to nonequilibrium Steady-state of open biochemical systems: linear single-molecule enzyme kinetics and nonlinear biochemical reaction networks,” International Journal of Molecular Sciences, vol. 11, no. 9, pp. 3472–3500, 2010.

C. J. Gadgil, “Size-independent differences between the mean of discrete stochastic systems and the corresponding continuous deterministic systems,” Bulletin of Mathematical Biology, vol. 71, no. 7, pp. 1599–1611, 2009.

E. Arslan and I. J. Laurenzi, “Kinetics of autocatalysis in small systems,” Journal of Chemical Physics, vol. 128, no. 1, Article ID 015101, 2008.

D. J. Higham, “Modeling and simulating chemical reactions,” SIAM Review, vol. 50, no. 2, pp. 347–368, 2008.

C. A. Gómez-Uribe and G. C. Verghese, “Mass fluctuation kinetics: capturing stochastic effects in systems of chemical reactions through coupled mean-variance computations,” Journal of Chemical Physics, vol. 126, no. 2, Article ID 024109, 2007.

M. Calder, A. Duguid, S. Gilmore, and J. Hillston, “Stronger computational modelling of signalling pathways using both continuous and discrete-state methods,” in Computational Methods in Systems Biology, vol. 4210 of Lecture Notes in Computer Science, pp. 63–77, 2006.

O. Wolkenhauer, M. Ullah, W. Kolch, and K. H. Cho, “Modeling and simulation of intracellular dynamics: choosing an appropriate framework,” IEEE Transactions on Nanobioscience, vol. 3, no. 3, pp. 200–207, 2004.

T. E. Turner, S. Schnell, and K. Burrage, “Stochastic approaches for modelling in vivo reactions,” Computational Biology and Chemistry, vol. 28, no. 3, pp. 165–178, 2004.

K. Vasudeva and U. S. Bhalla, “Adaptive stochastic-deterministic chemical kinetic simulations,” Bioinformatics, vol. 20, no. 1, pp. 78–84, 2004.

H. Qian, S. Saffarian, and E. L. Elson, “Concentration fluctuations in a mesoscopic oscillating chemical reaction system,” Proceedings of the National Academy of Sciences of the United States of America, vol. 99, no. 16, pp. 10376–10381, 2002.

I. J. Laurenzi, “An analytical solution of the stochastic master equation for reversible bimolecular reaction kinetics,” Journal of Chemical Physics, vol. 113, no. 8, pp. 3315–3322, 2000.

H. Wang and Q. Li, “Master equation analysis of deterministic chemical chaos,” Journal of Chemical Physics, vol. 108, no. 18, pp. 7555–7559, 1998.

I. G. Darvey, B. W. Ninham, and P. J. Staff, “Stochastic models for second-order chemical reaction kinetics. The equilibrium state,” The Journal of Chemical Physics, vol. 45, no. 6, pp. 2145–2155, 1966.

Q. Zheng and J. Ross, “Comparison of deterministic and stochastic kinetics for nonlinear systems,” The Journal of Chemical Physics, vol. 94, no. 5, pp. 3644–3648, 1991.

J. Goutsias, “Classical versus stochastic kinetics modeling of biochemical reaction systems,” Biophysical Journal, vol. 92, no. 7, pp. 2350–2365, 2007.

J. E. Keizer, Statistical Thermodynamics of Nonequilibrium Processes, Springer, Berlin, Germany, 1987.

S. Wray and T. Burdyga, “Sarcoplasmic reticulum function in smooth muscle,” Physiological Reviews, vol. 90, no. 1, pp. 113–178, 2010.

R. L. Winslow and J. L. Greenstein, “Cardiac myocytes and local signaling in nano-domains,” Progress in Biophysics and Molecular Biology, vol. 107, no. 1, pp. 48–59, 2011.

M. B. Cannell and C. H. T. Kong, “Local control in cardiac E-C coupling,” Journal of Molecular and Cellular Cardiology, vol. 52, pp. 298–303, 2011.

B. Mazzag, C. J. Tignanelli, and G. D. Smith, “The effect of residual Ca2+ on the stochastic gating of Ca 2+-regulated Ca2+ channel models,” Journal of Theoretical Biology, vol. 235, no. 1, pp. 121–150, 2005.

J. Wu, B. Vidakovic, and E. O. Voit, “Constructing stochastic models from deterministic process equations by propensity adjustment,” BMC Systems Biology, vol. 5, no. 187, pp. 1–21, 2011.

I. Györke and S. Györke, “Regulation of the cardiac ryanodine receptor channel by luminal Ca2+ involves luminal Ca2+ sensing sites,” Biophysical Journal, vol. 75, no. 6, pp. 2801–2810, 1998.

M. A. Huertas and G. D. Smith, “The dynamics of luminal depletion and the stochastic gating of Ca2 +-activated Ca2 + channels and release sites,” Journal of Theoretical Biology, vol. 246, no. 2, pp. 332–354, 2007.

M. R. Tadross, I. E. Dick, and D. T. Yue, “Mechanism of local and global Ca2+ sensing by calmodulin in complex with a Ca2+ channel,” Cell, vol. 133, no. 7, pp. 1228–1240, 2008.

H. Liang, C. D. DeMaria, M. G. Erickson, M. X. Mori, B. A. Alseikhan, and D. T. Yue, “Unified mechanisms of Ca2+ regulation across the Ca 2+ channel family,” Neuron, vol. 39, no. 6, pp. 951–960, 2003.

Y. L. Hashambhoy, J. L. Greenstein, and R. L. Winslow, “Role of CaMKII in RyR leak, EC coupling and action potential duration: a computational model,” Journal of Molecular and Cellular Cardiology, vol. 49, no. 4, pp. 617–624, 2010.

J. Heijman, P. G. A. Volders, R. L. Westra, and Y. Rudy, “Local control of β-adrenergic stimulation: effects on ventricular myocyte electrophysiology and Ca2+-transient,” Journal of Molecular and Cellular Cardiology, vol. 50, no. 5, pp. 863–871, 2011.

S. Pendyam, A. Mohan, P. W. Kalivas, and S. S. Nair, “Computational model of extracellular glutamate in the nucleus accumbens incorporates neuroadaptations by chronic cocaine,” Neuroscience, vol. 158, no. 4, pp. 1266–1276, 2009.

D. A. Rusakov, M. Y. Min, G. G. Skibo, L. P. Savchenko, M. G. Stewart, and D. M. Kullmann, “Role of the synaptic microenvironment in functional modification of synaptic transmission,” Neurophysiology, vol. 31, pp. 99–101, 1999.

D. M. Bers, “Cardiac excitation-contraction coupling,” Nature, vol. 415, no. 6868, pp. 198–205, 2002.

W. Chen, G. Aistrup, J. Andrew Wasserstrom, and Y. Shiferaw, “A mathematical model of spontaneous calcium release in cardiac myocytes,” American Journal of Physiology, vol. 300, no. 5, pp. H1794–H1805, 2011.

M. X. Mori, M. G. Erickson, and D. T. Yue, “Functional stoichiometry and local enrichment of calmodulin interacting with Ca2+ channels,” Science, vol. 304, no. 5669, pp. 432–435, 2004.

D. L. Stenoien, T. V. Knyushko, M. P. Londono et al., “Cellular trafficking of phospholamban and formation of functional sarcoplasmic reticulum during myocyte differentiation,” American Journal of Physiology, vol. 292, no. 6, pp. C2084–C2094, 2007.

J. Krupnick and J. Benovic, “The role of receptor kinases and arrestins in G protein-coupled receptor regulation,” Annual Review of Pharmacology and Toxicology, vol. 39, Article ID 2890319, 1998.

M. Narita, S. Shimizu, T. Ito et al., “Bax interacts with the permeability transition pore to induce permeability transition and cytochrome c release in isolated mitochondria,” Proceedings of the National Academy of Sciences of the United States of America, vol. 95, no. 25, pp. 14681–14686, 1998.

L. S. Song, E. A. Sobie, S. McCulle, W. J. Lederer, C. W. Balke, and H. Cheng, “Orphaned ryanodine receptors in the failing heart,” Proceedings of the National Academy of Sciences of the United States of America, vol. 103, no. 11, pp. 4305–4310, 2006.

J. R. Groff and G. D. Smith, “Calcium-dependent inactivation and the dynamics of calcium puffs and sparks,” Journal of Theoretical Biology, vol. 253, no. 3, pp. 483–499, 2008.

G. S. B. Williams, M. A. Huertas, E. A. Sobie, M. S. Jafri, and G. D. Smith, “A probability density approach to modeling local control of calcium-induced calcium release in cardiac myocytes,” Biophysical Journal, vol. 92, no. 7, pp. 2311–2328, 2007.

J. M. Hartman, E. A. Sobie, and G. D. Smith, “Spontaneous Ca2+ sparks and Ca2+ homeostasis in a minimal model of permeabilized ventricular myocytes,” American Journal of Physiology, vol. 299, no. 6, pp. H1996–H2008, 2010.

R. L. Winslow, A. Tanskanen, M. Chen, and J. L. Greenstein, “Multiscale modeling of calcium signaling in the cardiac dyad,” Annals of the New York Academy of Sciences, vol. 1080, pp. 362–375, 2006.

F. Von Wegner and R. H. A. Fink, “Stochastic simulation of calcium microdomains in the vicinity of an L-type calcium channel,” European Biophysics Journal, vol. 39, no. 7, pp. 1079–1088, 2010.

R. Grima, “Study of the accuracy of moment-closure approximations for stochastic chemical kinetics,” The Journal of Chemical Physics, vol. 136, Article ID 154105, 2012.

I. Lestas, J. Paulsson, N. E. Ross, and G. Vinnicombe, “Noise in gene regulatory networks,” IEEE Transactions on Automatic Control, vol. 53, pp. 189–200, 2008.

B. Sakmann and E. Neher, Single-Channel Recording, Springer, 2nd edition, 1995.

J. R. Groff and G. D. Smith, “Ryanodine receptor allosteric coupling and the dynamics of calcium sparks,” Biophysical Journal, vol. 95, no. 1, pp. 135–154, 2008.