This Article Abstract Full Text (PDF) Supplement All Versions of this Article: biophysj.106.097469v1 93/3/709    most recent Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via Google Scholar Google Scholar Articles by Maurya, M. R. Articles by Subramaniam, S. Search for Related Content PubMed PubMed Citation Articles by Maurya, M. R. Articles by Subramaniam, S.

A Kinetic Model for Calcium Dynamics in RAW 264.7 Cells: 1. Mechanisms, Parameters, and Subpopulational Variability

Mano Ram Maurya ^* and Shankar Subramaniam ^*  

^* Department of Bioengineering; Department of Chemistry and Biochemistry; and Graduate Program in Bioinformatics, University of California, San Diego, California

Correspondence: Address reprint requests to Shankar Subramaniam, Dept. of Bioengineering, University of California, San Diego, 9500 Gilman Dr., La Jolla, CA 92093-0412. Tel.: 858-822-0986; Fax: 858-822-5722; E-mail: shankar{at}ucsd.edu.

	ABSTRACT

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
Calcium (Ca²⁺) is an important second messenger and has been the subject of numerous experimental measurements and mechanistic studies in intracellular signaling. Calcium profile can also serve as a useful cellular phenotype. Kinetic models of calcium dynamics provide quantitative insights into the calcium signaling networks. We report here the development of a complex kinetic model for calcium dynamics in RAW 264.7 cells stimulated by the C5a ligand. The model is developed using the vast number of measurements of in vivo calcium dynamics carried out in the Alliance for Cellular Signaling (AfCS) Laboratories. Ligand binding, phospholipase C-ß (PLC-ß) activation, inositol 1,4,5-trisphosphate (IP₃) receptor (IP₃R) dynamics, and calcium exchange with mitochondria and extracellular matrix have all been incorporated into the model. The experimental data include data from both native and knockdown cell lines. Subpopulational variability in measurements is addressed by allowing nonkinetic parameters to vary across datasets. The model predicts temporal response of Ca²⁺ concentration for various doses of C5a under different initial conditions. The optimized parameters for IP₃R dynamics are in agreement with the legacy data. Further, the half-maximal effect concentration of C5a and the predicted dose response are comparable to those seen in AfCS measurements. Sensitivity analysis shows that the model is robust to parametric perturbations.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
Cytosolic calcium is a second messenger and plays an important role in intracellular signaling (1). Cytosolic Ca²⁺ is involved in regulating numerous cellular functions by regulating the activity of proteins such as calmodulin (2), calreticulin (3–5), and calcineurin (6). Dynamic changes in intracellular calcium serve both as an important indicator of cellular events and a quantitative measure of cellular response to stimuli.

Basic mechanism of calcium dynamics
In a eukaryotic cell, the average cytosolic Ca²⁺ concentration ([Ca²⁺]_i) is maintained low (0.05–0.5 µM), whereas Ca²⁺ concentrations in extracellular space and endoplasmic (or sarcoplasmic) reticulum (ER/SR) ([Ca²⁺]_ER) are several thousands times higher (7). This large concentration gradient of Ca²⁺ between different cellular compartments is utilized by cells to generate rapid intracellular Ca²⁺ changes through receptor-mediated mechanisms. In nonexcitable cells, such as macrophages, ligand-induced release of calcium from the ER is the main initiator of calcium dynamics. In excitable cells, other sources that initiate dynamics include the calcium influx through voltage-gated channels on the plasma membrane (PM) (8–11). Both these eventually lead to increased activation and opening of the inositol 1,4,5-trisphosphate (IP₃) receptor (IP₃R) channels on the ER membrane through either increased hydrolysis of phosphatidylinositol 4,5-bisphosphate (PIP₂) into IP₃ (12–14) or local increase in [Ca²⁺]_i leading to calcium-induced calcium release (CICR) (15–17). Calcium release from the ER leads to increased Ca²⁺ in the cytosol. Most of the calcium released binds to various proteins, such as calmodulin (CaM). Calcium is also pumped back to the ER by the sarco(endo)plasmic reticulum calcium ATPase (SERCA) pump. Some calcium is also expelled to the extracellular space through an Na⁺/Ca²⁺ exchanger (NCX) and the plasma membrane calcium ATPase (PMCA) pump. Calcium exchange with the mitochondria also has been observed at elevated levels of [Ca²⁺]_i. If Ca²⁺ concentration in the ER becomes quite low, calcium can enter the cell as calcium-release-activated current (CRAC) through store-operated channels (SOC) on the plasma membrane (18–20). Finally, calcium response can propagate to nearby cells connected through gap-junctions resulting in calcium waves (20).

Variation across different cell types
Cytosolic Ca²⁺ can be measured relatively easily using fluorescence and dye-based methods (21,22) at high sample rates (one sample every few seconds) (23). Recent reviews by Van Den Brink et al. (21) and Takahashi et al. (22) deal with elaborate descriptions of various methods of measuring free [Ca²⁺] and calcium fluxes between various compartments, such as the cytosol, the ER, and the mitochondria. Calcium levels can vary substantially from one cell type to another. For example, in excitable cells such as neuronal cells and muscle cells, for which most calcium models and measurement techniques have been developed (13,21,24–28), basal [Ca²⁺]_i is in the 0.2–0.5 µM range, whereas in nonexcitable cells such as macrophages, including RAW 264.7 cells (29,30), it is between 0.03 and 0.1 µM. In macrophages, due to their nonexcitable nature, oscillations are generally not observed, as evidenced by the Alliance for Cellular Signaling (AfCS) calcium data (29) and discussed in (30–32). In fact, the maximal changes in [Ca²⁺]_i (spatially averaged) are also only 0.05–0.1 µM, except when cells are treated with strong activators such as the bacterial toxin lipopolysaccharide (30).

The need for a model for calcium dynamics in RAW 264.7 cells
Recently, AfCS has conducted a large array of experiments to measure the temporal response of cytosolic calcium in RAW 264.7 cells in response to various ligands. RAW 264.7 cells are a macrophage-like, Abelson leukemia virus transformed cell line derived from BALB/c mice (AfCS protocol PP00000226 (23,33)). Considering the importance of calcium signaling, substantial efforts have been devoted to develop mathematical models for calcium dynamics. However, much of the work has focused on excitable cells, such as neuronal cells (13,27) and muscle cells (10,11,21,24–26,28), in which calcium oscillations are prevalent. Recent efforts have also focused on nonexcitable cells (12,13,34). Although literature is replete with simple, as well as complex, models of calcium signaling, including those by Fink et al. (35), Mishra and Bhalla (13), and Lemon et al. (12), there are no models of calcium dynamics in RAW 264.7 cells in the published literature. The model presented by Mishra and Bhalla (13), which incorporates the mechanisms for G-protein signaling, IP₃R channels, IP₃ metabolism, and binding of calcium with many proteins such as calmodulin, is the most comprehensive detailed kinetic model for mammalian cells. Most of the parameters used for simulation of the model are taken from the literature and are robust for making qualitative comparisons for various cell types. However, these parameter values may not be able to capture the fine quantitative features of data for a given cell type. In fact, a comparison of the time course predicted by Lemon et al. (12) for a nonspecific cell type and by Fink et al. (35) for neuronal cells with the time course in AfCS experiments reveals that the effective time constant in AfCS data for the rise phase is several times larger than in Lemon et al. (12) and Fink et al. (35), whereas the decay phase is more gradual than in Fink et al. (35) and much faster than in Lemon et al. (12). Although precise quantitative comparisons/predictions of features such as peak height are not mandated from calculations, it is essential to capture overall temporal variations for a model to be accurate and predictive for describing calcium dynamics. Thus, there is a clear need for some cell-specific constraining of the parameters. One way to carry out constraining of the parameters is to tweak each parameter in a trial-and-error fashion. The alternative and arguably the best way is to use cell-specific in vivo data to develop a systems-level model. The known biochemistry and other constraints from legacy data must be included to maintain the qualitative generality of the results. RAW 264.7 cells are different from other cells such as muscle cells not only in size but also in the basal level of [Ca²⁺]_i, as stated earlier. Further, AfCS experiments have shown that calcium response in RAW cells is nonoscillatory and exhibits just one peak after applying a ligand. Thus, it is reasonable to expect differences in the values of some of the lumped/effective kinetic parameters besides differences in the species concentrations in various cell types. Cell-specific modeling of calcium signaling, including the estimation/measurement of parameters, is quite important for real-life applications such as drug-dose response, and applications of systems biology. Also, in systemic models that include many processes and parameters, many different parameter ranges can yield different qualitative responses. Hence, even if one is interested in only qualitative results, some constraining of the parameters using cell-specific data is desired.

Recent AfCS experiments have revealed a new and very challenging problem. Different cell populations (cloned from the same parent cell), when triggered by a stimulus of the same strength, result in quantitatively as well as qualitatively different responses (different peak heights, rise times, etc.). According to our current understanding, the underlying cause is the variation in the unmeasurable concentrations of hundreds of species (components) inside the cells. This variability could be due to intrinsic sources, such as noise in gene expression (36), as well as extrinsic noise, such as unsynchronized cell cycles in different experiments. Though some of these fluctuations across cells get averaged out in cell populations, a nontrivial variation is still observed from one population to another. This variation is observed even in the basal levels in unstimulated and stimulated control cell lines and in knockdown cell lines, in addition to the variability after stimulation with ligands. If the basal levels were similar across different subpopulations (datasets), a common set of parameter values could have been used and then a knockdown could be modeled by modifying the parameter related to the knockdown process, e.g., the concentration of the relevant protein. The substantial "subpopulational variability" in the basal level required us to employ a special modeling approach in which some parameters are allowed to vary across different datasets. The AfCS data also includes a large set of perturbation data on the knockdown of proteins affecting the temporal characteristics of calcium dynamics. This data has been used for both modeling and validation of our context-specific models. The ligand used in this study is the complement 5a (C5a), a potent mediator of inflammation (37–44). An ordinary differential equation-based model is developed to understand calcium dynamics in RAW264.7 cells.

The first part of this two-part article deals with the development of a kinetic model using the AfCS data exhibiting subpopulational variability and prediction of dose response of C5a. Part 2 deals with the prediction of knockdown response and long-term response using the model developed in this article. Our modeling approach allowed us to utilize in vivo data exhibiting variability across different subpopulations. The resulting model elucidates the central role of phospholipase C (PLC) ß (PLCß) and IP₃ for calcium dynamics. Using the model, the efficacy of C5a is characterized and the robustness of calcium response is analyzed using global sensitivity analysis.

	MATERIALS AND METHODS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
Experimental data
We have utilized the time-course data on [Ca²⁺]_i response in macrophage RAW 264.7 cell populations made available by the AfCS (45). Owing to limitations of space, we refer to the AfCS website for the necessary details pertaining to experimental data (45). In brief, free [Ca²⁺]_i is measured using the Ca²⁺-sensitive fluorescent dye fura-2 (AfCS protocol PP00000211 (23)). As described on the AfCS Web site (23) and in the literature (46,47), free [Ca²⁺]_i is computed from the ratio of the intensities of fluorescent light at two wavelengths (340 nm and 380 nm). The intensities are actually related to the amount of calcium-bound fura-2. However, using the rapid buffering approximation for fura-2, as well as other buffers present in the cytosol, it is possible to back out free [Ca²⁺]_i from the intensity measurements. In the Results section, we show that the rapid buffering approximation is valid, i.e., for a small perturbation in [Ca²⁺]_i, the new equilibrated free [Ca²⁺]_i is reached much earlier than the next measurement. We use data from native cells as well as from knockdown cell lines. The knockdowns were carried out using short hairpin RNA interference (shRNAi). Data from four time-courses are used for constraining/estimating the parameters. Additional description of the data used is given in Supplementary Material (Preprocessing of time-course data).

Mathematical model of calcium dynamics in RAW 264.7 cells
A simplified model for calcium signaling that incorporates most of the mechanisms included in the existing models has been developed (Fig. 1). This simple model, however, includes several mechanisms explicitly so that knockdown of important proteins, such as G-protein-coupled receptor (GPCR) kinase (GRK), Arrestin, G_ß, and G_,i, can be modeled quantitatively. We note that although these mechanisms have been modeled by several researchers to study G-protein signaling (48–51), they have not been included in most models for calcium signaling. Although the model by Lemon et al. (12) does include these details, the authors assumed a constant ligand concentration, which is appropriate only for saturating concentrations.

View larger version (35K): [in this window] [in a new window]   FIGURE 1  A simplified model for calcium signaling including calcium influx, ER, and mitochondrial exchange and storage, used in the conceptual-model-based computation. (A) Overall schematic model. The ligand C5a binds to its receptor C5aR on the plasma membrane, activating G-protein Gi. The free subunit Gß binds to and activates PLCß, which hydrolyzes PIP2 into IP3 and DAG. IP3 binds to its receptor on the ER membrane and the IP3R channels open to release calcium into the cytosol. Other calcium fluxes (e.g., with mitochondria and extracellular space) are also shown. (B) The mechanisms for the receptor module (box 1), the GTPase cycle module (box 2), and IP3-generation module (box 3), and the feedback effects (boxes 1 and 4). PIP2, phosphatidylinositol 4,5-bisphosphate; IP3, inositol 1,4,5-trisphosphate; IP3R, IP3 receptor; IP3,p, a lumped product of IP3 phosphorylation; Cai, cytosolic Ca2+; ER, endoplasmic reticulum; CICR, calcium-induced calcium release; SERCA, sarco(endo)plasmic reticulum calcium ATPase; PMCA, plasma membrane calcium ATPase; NCX, Na+/Ca2+ exchanger; mit (subscript), mitochondria; L, ligand C5a; R, receptor C5aR; GRK, G-protein-coupled receptor kinase; CaM, calmodulin; PLCß, phospholipase C-ß; GAP, GTPase activating protein; RGS, regulator of G-protein signaling; DAG, diacylglycerol; PKC, protein kinase C; Pi, phosphate.

A text-based prototype reaction interpreter has been developed that can read both simple and complex reactions. The simple reactions have law of mass action rate expression, whereas lumped reactions follow Michaelis-Menten (M-M) kinetics or Hill dynamics. More complex formats include modulation (law of mass action, M-M, or Hill dependence). Some examples of the reactions are presented in Table S1 (Supplementary Material). Lumped-reaction formats are quite useful for model simplification.

Various modules involved in ligand-induced release of calcium from ER are discussed below. In short, binding of the ligand (C5a, denoted by L) to its receptor (C5aR, denoted by R) leads to the activation of the receptor, which in turn leads to the activation of G-protein G_i. Free G_ß binds with PLCß and activates it, resulting in increased hydrolysis of PIP₂ into IP₃. Increase in the IP₃ level results in increased release of calcium from the ER.

Receptor module
Reactions 1–11 (Fig. 1 B, box 1) are as follows. Reaction 1 is activation of the receptor due to ligand binding. Reaction 2 represents binding of GRK to G_ß and is very fast. Reactions 3 (catalyzed by GRK) and 4 (catalyzed by GRK.G_ß) represent desensitization of the ligand-bound active receptor (L.R), due to its phosphorylation. The activity of GRK is enhanced by the PKC.DAG.Ca_i complex due to increased phosphorylation, resulting in a negative feedback on [L.R] (53,54). This effect on the rate of reaction 3 is modeled as an enzymatic activation by calcium to avoid explicit modeling of the complexes. Reaction 5 is the dissociation of ligand and phosphorylated receptor (12,49,50). Reactions 6 and 7 are for internalization of the ligand-bound phosphorylated receptor, L.Rp (55). Reaction 7 is catalyzed by Arrestin. Reactions 8–10 together constitute receptor recovery. Reaction 8 is the dissociation of the internalized ligand-bound phosphorylated receptor, L.Ri, into the free ligand and the phosphorylated receptor (internalized form), R_p,i (56). The internalized free ligand gets ubiquitinated and degraded. R_p,i then gets dephosphorylated and recycled back to the surface through reaction 9 (12,48). A small part of R_p,i also can get degraded (reaction 10). Reaction 11 is fresh receptor generation (48,51). To the best of our knowledge, in most existing models, GRK and arrestin (Arr) have not been explicitly included in calcium signaling studies, probably because knockdown data on these proteins have not been quantitatively analyzed with respect to kinetic modeling. The rate constant for reaction 4 is much higher (100-fold) than for reaction 3, so that most phosphorylation occurs due to GRK.G_ß when ligand is present and due to GRK during the basal state. Similarly, internalization (reactions 6 and 7) is enhanced substantially by arrestin (buffered). Reactions 3, 4, and 7 are lumped-enzymatic reactions, as described below:

GTPase cycle module
Reactions 12–16 (Fig. 1 B, box 2) depict GTPase cycle. Reactions 12 and 13 take place in the absence of active receptor and GTPase activating protein (GAP), respectively. Reaction 15 is similar to reaction 12 but is catalyzed by L.R, and reaction 16 is similar to reaction 13 but is catalyzed by GAP (A, RGS) (57). A more detailed description of the GTPase cycle, similar to the three-cube model in Fig. 1 A of Bornheimer et al. (58) but with simplification of reactions into (Fig. 1 B, reaction 12) and explicit modeling of G_ß and G_,i, was included in a detailed model, but its fit to one set of basic experimental data (stimulation with 250 nM C5a) was similar to the compact model of Fig. 1 B. Hence, to reduce computational complexity, only the simpler model (Fig. 1 B) is used for further analysis with knockdown data. Additional discussion on this simplification is presented in Supplementary Material (Approximations and the lumped/simplified mechanisms in the model). Reactions 15 and 16 are lumped-enzymatic reactions, where T is GTP, D is GDP, and A is the GAP RGS.

IP₃ module
In the basal state, most of the IP₃ is generated due to slow hydrolysis of PIP₂ (reaction 17), since free G_ß is present in very small amounts. Upon G-protein activation, dissociated G_ß binds to PLCß. Cytosolic Ca²⁺ (Fig. 1 B, box 3, Ca_i) can bind to both PLCß and PLCß.G_ß. Binding affinities for G_ß and calcium-bound forms of PLCß are 10 times higher than that of free PLCß (13). Each of PLCß.G_ß, PLCß.Ca_i and PLCß.G_ß.Ca_i catalyze hydrolysis of PIP₂, but PLCß.G_ß.Ca_i is the most potent. In our model, for simplification, a lumped-enzymatic reaction is used to model the enhancement due to PLCß.G_ß.Ca_i (see reaction 18 (Fig. 1 B, box 3)):

As explained in Supplementary Material (Approximations and the lumped/simplified mechanisms in the model), the rate expression of the above reaction is a close approximation to an expression derived based on the assumption that all the four reversible reactions related to PLCß are in equilibrium. Further, it is assumed that the enzymatic effect of PLCß.G_ß and PLCß.Ca_i is captured through suitable increase in the value of the rate constant. The amount of G_ß bound to PLCß complexes (estimated to be 10% of free G_ß) is ignored for G_ß balance (Supplementary Material, Details of explicitly modeled reactions) in the reduced model; if all other species are included, then little reduction is achieved, and therefore it would be appropriate to model all reactions explicitly. Since PLCß is not used in any other reaction in our model, PLCß denotes the total PLCß (buffered).

Reactions 19 and 20 (Fig. 1 A) are simplified and highly lumped representations of IP₃ metabolism, i.e., degradation/conversion to/from other inositol phosphates and back to PIP₂, with only one intermediate pseudospecies, namely, IP_3,p or IP₃ product (12). Since oscillation is not observed in the experimental data on C5a stimulation of RAW cells, this simple representation is considered sufficient. As Mishra and Bhalla (13) have noted, such a simple description may be insufficient to model oscillatory response, even though Marhl et al. (59) have shown that interactions with the mitochondria (included in our model) can result in oscillations, especially above cytosolic [Ca²⁺] of 0.5 µM.

It should be noted here that the free G_ß subunit is responsible for the activation of PLCß, as opposed to activation by the free subunit of a G-protein, e.g., G_qT (G-protein G_q) upon activation of P₂Y₂ receptors, as presented by Lemon et al. (12). For the G-protein G_i, the subunit, i.e., G_iT, does not bind to PLCß. G_ß also has been implicated in calcium oscillations during fertilization (60).

Feedback effects from calmodulin
Calmodulin binds with intracellular Ca²⁺ (Ca_i) (Fig. 1 B, box 4, reaction 21), and the resulting complex binds with GRK (reaction 22), reducing the effective amount of free GRK that can bind G_ß. The result is reduced phosphorylation of the active receptor, and thus, this constitutes a functional positive feedback.

Mathematical representation of the model
The state variables are described by a set of ordinary differential equations (61) involving the Ca²⁺ fluxes between different cellular compartments and other fluxes due to reactions. The state variables used to model the details of ligand-induced generation of IP₃ are [L], [R], [L.R], [G_ß], [GRK], [L.Rp], [Rp], [L.Ri], [R_p,i], [R_pool], [G_,iT], [G_,iD], [PIP₂], [IP₃], and [CaM] (15 of them). Calcium dynamics introduces four more state variables, [Ca²⁺]_i, [Ca²⁺]_ER, h, and [Ca²⁺]_mit. [Ca²⁺]_ER and [Ca²⁺]_mit denote the concentration of free Ca²⁺ in the ER and mitochondria, respectively; h is explained later. Thus, the model has 19 state variables. The quantities of all chemical species are in terms of their concentrations, normalized with respect to a unit volume of the cytosol. Since, during optimization, most unknowns, e.g., parameters and initial conditions, are specified in terms of a range rather than unique values, whenever needed, a cytosolic volume of 10 pL (picoliter) or cell diameter 27 µm, typical for large macrophages, has been used.

The reaction fluxes, the related parameters, and the differential equations for the 15 state variables related to the reactions are given in Supplementary Material, Details of explicitly modeled reactions. Moieties (conservation relations) have been used to reduce the number of state variables by computing the concentrations of some species. The expressions for the other four state variables and fluxes, related to the dynamics of calcium, are given below.

Equations for calcium dynamics

(1)

are the ratio of free calcium to total (free and bound) calcium in the cytosol, ER, and mitochondria, respectively, assuming fast buffering (equilibrium) with calcium-binding proteins. A fixed value of is used. is the ratio of the ER volume to the volume of the cytosol and is the ratio of the volume of the mitochondria to that of the cytosol. The expressions for ß_i and ß_ER are

(2)

where Pr_tot,e is the total concentration of endogenous buffer (Ca²⁺-binding proteins) in cytosol, K_m,e is the dissociation constant for binding of Ca²⁺ to endogenous buffer in the cytosol, Pr_tot,x and K_m,x are the corresponding quantities for binding of Ca²⁺ to the exogenous/mobile buffer (fura-2) in cytosol, Pr_tot,ER is the total buffer in the ER, and K_m,ER is the dissociation constant for binding of Ca²⁺ to the buffer in the ER.

   The flux of calcium from ER to cytosol through the IP₃R channel
De Young and Keizer (62) have developed a detailed model for activation/deactivation of IP₃R channels that captures the CICR effect and the biphasic behavior of the channel opening at low and extremely high calcium concentrations. Li and Rinzel (63) simplified the De Young and Keizer model (62) using time-scale analysis, which was also used by Fink et al. (35) to simulate spatial variations in neuroblastoma cells. The channel flux is given by

(3)

where P_open is the IP₃R channel instantaneous open fraction or probability and h is the fraction of IP₃R to which calcium is not bound to the inhibitory site. K_IP3 is the dissociation constant for binding of IP₃ on the IP₃R channel. K_act is the dissociation constant of Ca²⁺ binding to the activation site of the IP₃R channel. In the above expression, the second term arises due to the binding of calcium to the activation site on the IP₃R. Thus, at low calcium concentration, calcium has an activation effect on J_ch. However, at higher levels of calcium, the equilibrium value of h decreases with increasing [Ca²⁺]_i (see the differential equation for h). This is due to binding of calcium to the low-affinity inhibitory site on the IP₃R channel and results in biphasic response of J_ch with respect to [Ca²⁺]_i.

   Flux through the SERCA pump back to ER

(4)

where is the maximal rate of SERCA pump uptake and is the dissociation constant of Ca²⁺ binding to the SERCA pump (12,34,35,64, 65).

   Leakage flux from the ER

(5)

where (s^–1) is the rate constant of Ca²⁺ leak flux through the ER membrane (12,59).

   IP₃-dependent flux across the plasma membrane from the ECM into the cytosol
This is a combined effect of flux through various channels, including the SOC, which open in response to the depletion of calcium in the ER (20). Overall, this flux acts as a positive feedback mechanism to maintain cytosolic and ER calcium levels in anticipation of the emptying of the ER and the calcium efflux to the ECM (18–20,66,67). Since includes the SOC flux, a separate SOC flux is not used.

(6)

   Leakage flux from the ECM to cytosol (20)

   Calcium efflux from cytosol to ECM through the PMCA and NCX

(7)

is due to two pumps, a low-capacity (high-affinity) and a high-capacity (low-affinity) pump (34). and are the maximum capacities of the low and high capacity PMCA pump, respectively, and and are the Ca²⁺ concentrations for half-maximal flux through the two pumps. Similarly, is the maximum flux through the Na⁺/Ca²⁺ exchanger and is the corresponding Michaelis-Menten constant for the binding of Ca²⁺.

   Calcium uptake into and efflux from mitochondria
Based on the model of Marhl et al. (59), calcium flux into mitochondria (68) is given by

(8)

where is the maximum permeability for uptake of Ca²⁺ by the mitochondrial uniporter and is the half-maximal [Ca²⁺]_i. Marhl et al. (59) have used a Hill coefficient of 8, since this uptake rate is quite steep with respect to increased concentration of the cytosolic Ca²⁺. However, we have used a much smaller value of 4 (still quite steep) to be able to fit experimental data. The efflux term is given by

(9)

where is the maximal rate for the Na⁺/Ca²⁺ exchanger and permeability transition pores (PTPs), is the corresponding half-maximal [Ca²⁺]_i for efflux and is the rate constant for a nonspecific leak flux.

In the expression for dh/dt (Eq. 1), Q is an effective Michaelis constant capturing the IP₃R channel's calcium-dependent inhibition, and is given by (63)

(10)

where is the dissociation constant for Ca²⁺ binding to the inhibition site and is the dissociation constant of IP₃ binding when the Ca²⁺ inhibitory site is occupied. Finally, is the rate of reaction 21, the binding of two Ca²⁺ to calmodulin (Fig. 1 B, box 4). The values of the parameters used in the above expressions are listed in Table S2 (Supplementary Material). Initial conditions for the state variables and the values of other parameters, such as the concentrations of buffered species/enzymes, are listed in Supplementary Material, Conservation parameters, buffered species and the initial conditions.

Solution of the ordinary differential equations, and parameter estimation
The procedures for solving the ordinary differential equations and parameter estimation using a hybrid stochastic-search-based algorithm that combines genetic algorithm (GA), differential evolution (DE), and particle-swarm optimization (PSO) are summarized in Supplementary Material. The objective function and the constraints used for parameter estimation are also given in Supplementary Material.

Prediction of dose response
We predicted the dose response of C5a using the best set of parameter values for the master dataset. The C5a dose is varied from 1 nM to 500 nM. Between different runs, the only change is in the ligand strength i.e., [C5a] at t = t_LA (time of ligand addition). To perform the simulation, the system is solved to get the initial basal state (phase 1), and then the system evolves at the basal conditions till t = t_LA (phase 2, part 1). At t = t_LA, [C5a] is reset to the dose strength, e.g., 30 nM, corresponding to the experimental dataset 1. All other initial conditions remain at the basal steady-state level, and simulation is carried out (phase 2, part 2).

Sensitivity analysis
To consider both small and large changes, each parameter (one at a time) is perturbed by multiplying its base value by factors [, , , 1, 2, 4, 8]. After simulation, the shift in the basal level (baseline shift), the peak-height compared to the respective basal level, and their sum, i.e., the peak height from the basal level for the base set, are computed. Since our fit-error function included both the good fit to basal level and good fit to response after ligand addition, a weighted average of normalized baseline shift B (weight 0.5), and normalized variation in peak height (H) is used to rank the parameters. B is computed as the maximum absolute baseline shift divided by the maximum peak height across the seven perturbations. H is computed as the maximum absolute difference in peak height across the seven perturbations divided by the maximum peak height. Let b_i be the baseline shift and h_i be the peak height for the perturbations. Then, mathematically,

where h_i is always nonnegative. The parameters are sorted in the decreasing order of (0.5 x B + H), i.e., from the most to least sensitive. It should be noted that only the response of [Ca²⁺]_i is considered here for sensitivity.

Variation between different subpopulations of cells
For a chosen strength of the ligand C5a, when the experiments were repeated with multiple subpopulations of cells, though the qualitative response was similar, substantial quantitative variation is observed in key features such as basal level, peak height, etc. This variation is present across different control datasets, as well as in the basal levels of knockdown cell lines. This variability is an inevitable fact concerning biological systems (69). It cannot be explained by using a single set of values of the kinetic parameters and the initial conditions, since concentrations and fluxes vary across subpopulations. As mentioned earlier, if the basal levels were similar across different subpopulations (datasets) without knockdown, a common set of parameter values could have been used. However, in the data used, variation in the basal level even without knockdown is as much as 30% of the peak height in some datasets, e.g., compare datasets 1 and 2 (Fig. 2, panel 5, and Fig. 3, panel 2, in Results).

View larger version (22K): [in this window] [in a new window]   FIGURE 2  Time course of several state variables including [Ca2+]i for dataset 1 (control for knockdown of RGS with 30 nM C5a), as predicted by using the best parameter values. (1) Fraction (percentage as compared to total surface receptors, Rtot) of free surface receptors ([R], solid line), ligand-bound active receptors ([L.R], dashed line), internalized ligand-bound phosphorylated receptors ([L.Ri], dash-dotted line), and internalized free receptors ([Rp,i], dotted line). Substantial change (decrease in [R] and increase in [L.R]) is observed soon after ligand addition (tLA = 0 s). Upon phosphorylation, the active receptor, L.R, gets internalized. Thus, [L.Ri] starts to increase after a few seconds and reaches its peak at 120 s. As L.Ri dissociates into the ligand and Rp,i, [Rp,i] shows an increase after a delay of 20 s. (2) Percentage of Gß in the free state (solid line) and bound to GRK (dash-dotted line), and percentage of active G-protein ([G,iT], dashed line). All three lag behind [L.R] and exhibit qualitatively similar profiles. (3) Indicators of IP3-receptor channel activity: Q, the effective Michaelis constant (solid line), and h (dashed line). As indicated by the arrow, the y axes for Q and h are on the left- and right-hand sides, respectively. (4) Concentration of IP3 (solid line, the left y axis) and the fraction of unhydrolyzed PIP2 (dashed line, the right y axis). (5) Fit of the model prediction (dashed line) to experimental data (solid line). (6) Time course of [Ca2+]ER (the left y axis) and [Ca2+]mit (labeled as [Ca2+]m on the right y axis).

View larger version (11K): [in this window] [in a new window]   FIGURE 3  Time course of state variables for dataset 2 (control for knockdown of Gi with 100 nM C5a). (1) [IP3] (solid line, left y axis) and the fraction of unhydrolyzed PIP2 (dashed line, right y axis). (2) Fit of the model prediction (dashed line) to the experimental data (solid line). (3) Time course of [Ca2+]ER (solid line, left y axis) and [Ca2+]mit (dashed, right y axis).

ß

,i

This flexibility requires the generation of simulation/optimization code that is dependent on the specification of the data. This is the only effective approach to utilize such data for kinetic modeling. Thus, a total of 19 parameters are allowed to vary from one subpopulation to another. Several other combinations of fewer parameters were found to be unsatisfactory. The parameters are separated into two classes, nonvarying and varying. Corresponding to the chosen datasets, first the master list of parameters is related to one basic dataset (referred to as "Master dataset"). Here, the first dataset (the control data on RGS10 knockdown) is chosen as the master dataset (Supplementary Material, Table S3). The total number of parameters to be estimated that are associated with the master dataset is 65. This includes k_f,9, which was optimized between 7.0 x 10^–4 and 0.001 but later set to 0.001 to get nearly full recovery of the receptor in 30–35 min. Of these 65 parameters, 46 are common to all subpopulations, whereas the 19 varying parameters are instantiated for each dataset. Essentially, for each additional dataset, the list of varying parameters is appended to the list. The list of varying parameters is tracked so that during simulation stage, appropriate parameter assignment can be made for each dataset. The list of varying parameters cannot include a calculated parameter. Any calculated parameters are recalculated to maintain the constraints on the parameters.

With four datasets, the total number of unknown parameters becomes 46 + 4 x 19 = 122, an increase by a factor of 2. The total amount of experimental data has increased by a factor of 4. Hence, at least in principle, this scheme results in more constraints on the parameter space. In reality, the result could be somewhat different, because with nonlinear dependence on the parameters, even a small increase in the dimensionality of the search space could result in substantial effective freedom in the search space. Since the constraints described in Supplementary Material are applied to all datasets, this problem gets alleviated to some extent.

There is one potential problem in letting the variable parameters vary for different subpopulations across the entire range used for optimization. In reality, the values of a variable parameter across different subpopulations are within a certain factor of each other. For example, is optimized between 0.01 and 0.5. Its value for different subpopulations could be in the range 0.05–0.12, within a factor of 2.4 of each other. Thus, constraints must be imposed on how much such parameters could differ from each other across different populations. We generate the values of such parameters () for the master dataset and then require that the corresponding values for other datasets (populations), (for the k^th dataset), should be within a factor range of 0.5–1.5 (user-defined) as compared to ensuring that such parameters cannot deviate by more than a factor of 3 (= 1.5/0.5) between any two datasets. Thus, the additional parameters for other datasets are not the actual values of the parameters (), but instead are the corresponding ratios (). Using the ratio is easily computed as The resulting value is clipped at the lower or upper bound on the parameter value and is then ready for use in simulation. The calculated parameters are recomputed for each dataset with these corresponding values of the other parameters. A major advantage of this approach of handling the variability is that the additional ratio parameters to be estimated are all in a small range of (0.5–1.5) instead of being dependent on the optimization range of the actual parameter itself. This approach provides an elegant way of relating the variable parameters (intrinsic biochemical properties) of a knockdown dataset and the corresponding control dataset.

Modeling of knockdown data
In the shRNAi-based knockdown cell lines, cells reach a new initial basal state (AfCS protocol PP00000211 (23)). This new basal state is a manifestation of the knockdown of the protein of interest. Hence, to model the knockdown mathematically, once the values of the variable parameters are assigned through the approach described in the previous section, the knockdowns are emulated. Any calculated parameters that could be dependent on the value of the variable being altered due to knockdown are computed only then, and the initial basal steady state is computed. The value with knockdown is computed as where x is the base value of the affected variable (initial condition, concentration of a buffered protein, etc.) that would have been effective in the absence of knockdown, is the value with knockdown, and %knockdown is the level of knockdown expressed in percentage.

It should be noted that for the ratio constraint on G_ß,tot/G_,i,tot (Supplementary Material, The objective function and the constraints), the values before knockdown are used. After the knockdown, this ratio constraint need not be satisfied, especially with the knockdown of G_ß or G_,i (G_,i2 and G_,i3). Further, the G-proteins from other families such as G_,q are not modeled. It is assumed that the total amount of all G-proteins in the cell is around five times that of G_,i. Since the experimental data on the knockdown of G_,i is based upon the actual cell, with no effect on other G-proteins, for modeling purposes, the effective knockdown level is reduced by a factor of 5. Without this scaling, absurd values of the parameters and response were obtained.

Modeling and prediction of knockdown response requires computation of a new basal state for nontraining scenarios before simulating the model with stimulus as described in Supplementary Material, Computation of new basal state for knockdown scenarios.

	RESULTS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
Constrained parameter values and fit to experimental data
The detailed results of parameter estimation using the hybrid genetic algorithm-differential evolution-particle swarm optimization approach are presented in Supplementary Material. The essential results are presented here briefly. The best parameter values of the common parameters and the base values of the variable parameters are listed in Table 1. Table S6 (Supplementary Material) lists two other such parameter-value sets that fit the experimental data nearly to the same degree as the best set and satisfy all constraints. Some parameters, such as and vary by a factor of 1.5 or more across the three sets (Table S6), suggesting some redundancy and nonlinear relationships among the parameters.

Constrained parameter values versus values based on legacy data
As mentioned earlier, parameter-constraining using in vivo data was necessary, since the legacy values based on in vitro data in different cell types reported in the literature could not predict experimental data even within an order of magnitude. It is possible that the difference in the values of the individual parameters based on in vitro or legacy data and in vivo data may be small, but this difference is present in many parameters whose overall effect becomes very large. Here we compare these values for some of the parameters to justify the use of in vivo data for developing computational models.

Parameters for the receptor activation/desensitization
The bounds for these parameters, listed in Table S4 (Supplementary Material), were based on Hoffman et al. (48), Pumiglia et al. (70), Woolf and Linderman (50), Chen et al. (71,72), Daaka et al. (73), and Sklar and co-workers (74,75). None of the values were measured or designed for the RAW 264.7 cell system. To cite a few examples, Hoffman et al. (48) used k_f,1 = 84 ± 19 µM/s, whereas Sklar et al. (74) used 17–27 µM/s for the N-formyl peptide receptor system. We used the bounds 10–100 µM/s. The optimized/constrained value of k_f,1 is 58.9 µM/s for the best set (Table 1). For the two other sets, it is 60 µM/s (Table S6). The bounds for k_b,1, chosen similarly, were 0.1–1 s^–1. Hoffman et al. (48) used k_b,1 = 0.37 ± 0.1 s^–1 for the N-formyl peptide receptor system. The optimized value is 0.2 s^–1 (Tables 1 and S6), well within a factor of 2, even though they are different systems. We had also used the constraint that was based upon the values suggested by Linderman and colleagues (48–50). This constraint on the K_D is also within the value suggested specifically for C5a-C5aR, viz. 1–100 nM for neutrophils and eosinophils, by Leslie (37). Thus, a fair amount of legacy information was used to ensure that the constrained parameter values satisfy the biochemistry embedded in the literature values. The list goes on for the parameters related to reactions 2–11. For brevity, we refer to Tables 1, S4, and S6 and the constraints listed in Supplementary Material.

Parameters for G-protein module (reactions 12–16)
The bounds for most of the parameters of these reactions were chosen based upon the values suggested by Bornheimer et al. (58), Ross and Wilkie (76), and Biddlecome et al. (77). Since reactions 12 and 15 use the lumped scheme (Supplementary Material), the parameter values suggested in Bornheimer et al. (58), Ross and Wilkie (76), and Biddlecome et al. (77) were scaled accordingly to compute the bounds. For reaction 12 (basal activity), the bounds on k_f,12 were 1.0 x 10^–7 to 3.0 x 10^–6 µM/s (Table S4), whereas the optimized value is 3.67 x 10^–7 for the best set (Table 1). Due to lumping, and scaling of the parameter, this value cannot be directly compared, but the values of the raw parameters were based on a large body of literature (58). For reaction 13 (intrinsic GTPase activity), the bounds on k_f,13 were 0.05–0.07 s^–1 (Table S4), as suggested by Ross and Wilkie (76) and Sprang (78). The optimized value, 0.05 s^–1, actually matches exactly with the value for the G-protein G_,i1 (78,79). For RGS-catalyzed hydrolysis (reaction 16), Ross and Wilkie (76) report k_cat,16 = 15–25 s^–1 for G_,q-PLCß and G_,q-RGS4 systems. We allowed a range of 1–50 s^–1, since G-protein G_,i is from a different family (Table S4). The optimized value is 1.54 s^–1 (Table 1). This lower value reflects the observation of Ross and Wilkie (76) that RGS has a selectivity for G_,q over G_,i. Such corroborations increase our confidence that the parameter values are in the right range.

Parameters for the IP₃ module
Since the mechanisms proposed by Li and Rinzel (63) are directly used to model the release of Ca²⁺ from ER, the related parameters were constrained within the bounds gleaned from their work (63). The notation used by Fink et al. (35) is used here. Despite some relaxation to compensate for the simplification, the optimized values are quite close to those suggested by Li and Rinzel (63). Several examples are given here. The optimized value of the on-rate constant for Ca²⁺ binding to an inhibitory site, 0.104 µM/s (Table 1) is close to the value 0.2 µM/s suggested by De Young and Keizer (62) and used by Li and Rinzel (63), even though a wide bound of 0.1–4.0 was specified (Table S2), since Fink et al. (35) had used a value of 2.7 µM/s. The same is true for the K_D for Ca²⁺ binding to an inhibitory site. The bounds were 1–1.5 s^–1 (Table S2). The optimized value is 1.0 s^–1, whereas Li and Rinzel (63) have used 1.05 s^–1. Fink et al. (35) had used a much lower value of 0.2 s^–1, but we biased it to Li and Rinzel's value (63) by specifying a lower bound of 1. For (K_M for activation of IP₃R by IP₃), Hofer et al. (20) used 0.3 µM and Li and Rinzel (63) used 0.13 µM (Table S2). We optimized in the range 0.1–1, and the optimized value of 0.136 (Table 1) agrees well with data in Li and Rinzel (63). Similar agreement is found for the parameters and (Tables S2 and 1).

Parameters for calcium dynamics
There are 20 parameters related to various processes such as Ca²⁺ binding to cytosolic and ER proteins, and exchange with the ER across the plasma membrane and the mitochondria. The most important parameters in regulating the calcium dynamics are related to the exchange with ER ( and ) and binding to the calcium-binding proteins ( and ). The bounds for (maximum intrinsic channel flux per µM of [Ca²⁺]_ER) were based upon comparison of the corresponding values used by Marhl et al. (59), Lemon et al. (12), and Fink et al. (35) to get a peak flux of 1.5 µM/s with [Ca²⁺]_ER = 200–500 µM (Table S2). There is wide difference in the values used by Marhl et al. (59) (4100 s^–1) and Lemon et al. (12) (575 s^–1). Both these are very high, since the authors have used much lower [Ca²⁺]_ER. We get an optimized value of 0.189 s^–1. This low value in our model is justified, since we have much higher [Ca²⁺]_ER and much lower [Ca²⁺]_i. As an example, Fink et al. (35) use 3500 µM/s as maximum flux and peak height of the [Ca²⁺]_i time-course is as much as 1 µM. In our case, peak height is 0.06 µM (experimental data), lower by a factor of 17, and the equivalent maximum flux is 30–50 µM/s, lower by a factor of 70–100. After binding to the proteins, the remaining excess calcium is cleared primarily by the SERCA pump, which follows Hill dynamics with = 2. The ultrasensitivity and the saturation nature of the SERCA pump justify the need to get a 70–100 times lower maximum flux to get 17 times lower peak height of [Ca²⁺]_i. The optimized value of (half-saturation level for SERCA pump), 0.754 µM, is about double the value used by Lemon et al. (12), whereas (maximum flux through the SERCA pump), 114 µM/s, is 2.5 times higher than that used by Lemon et al. (12). Larger was necessary to get a flux of 1 µM/s, since [Ca²⁺]_i is about two to four times lower. Such differences in the parameter values, as compared to values reported previously for other specific or nonspecific cell types, justify the in vivo data-based constraining of the parameters. Although using the values reported in the literature for other cell types will result in similar gross qualitative shapes (e.g., first rise and then return to basal levels), the quantitative differences can be of more than an order of magnitude as explained above. Even the finer details of qualitative shape, such as the time constants for the rise and decay phases, can differ significantly, further strengthening the case for cell-specific constraining of the parameters. The bounds for the parameters and were based on values used by others (12,35,52,59). The optimized value of 191 µM, is about the same as used in Marhl et al. (59) and Lemon et al. (12) (120–150 µM). The value of 2.43 µM, is lower than reported in Lemon et al. (12) but much higher than 0.1 µM, the value used by Marhl et al. (59). Thus, there is huge variation in the literature values. Since Ca²⁺ binds to the endogenous buffer with a low affinity, a value >1 µM is justified. On the contrary, Ca²⁺ has high affinity for the exogenous buffer. For fura-2, is mostly reported. In our model, we relaxed it somewhat to partially compensate for uncertainty in and but still the optimized value of 0.139 µM for the best set is not much different. In fact, it is 0.2 µM for set 2. The lower values of and could be partially due to the lower [Ca²⁺]_i. Somewhat higher could offset this, but the fact that there are several pools of proteins with somewhat different K_D values points to the approximations made here.

Variation in parameters across datasets due to sub-populational variability
Due to subpopulational variability within the best set itself, the values of the variable parameters for the different datasets are substantially different, as is evident from the values of the ratio parameters listed in Table 2. In Table 2, ratio values are listed for datasets 2–4. For dataset 1, this value is 1 by default. As mentioned earlier (Materials and Methods), for datasets 2–4, the actual value of the variable is computed by multiplying the raw ratio () by the value of the parameter for dataset 1 () and then limiting the resulting value to the lower or upper bound, i.e., if and if Then the true value of the ratio parameter, reported in Table 2, is computed as Of all the 19 ratio parameters, the range of variation (MAX/MIN ratio across the four datasets) is <1.5 for eight parameters. For the parameters and this range is 2, which is not surprising since all three exert large control on calcium dynamics through the size factor, maximal flux through the IP₃R channel, and the expulsion of calcium to the extracellular space. For other parameters, the range of variation is between 1.5 and 2. Overall, these values suggest that substantial variation is exhibited by different subpopulations of cells. The specifics of how these differences help explain the observed variety of calcium response in different datasets is presented below in the context of fit to experimental data.