^{1,a)}, Mano Ram Maurya

^{2,b)}, Daniel M. Tartakovsky

^{1,c),d)}and Shankar Subramaniam

^{2,3,c),e)}

### Abstract

Deterministic models of biochemical processes at the subcellular level might become inadequate when a cascade of chemical reactions is induced by a few molecules. Inherent randomness of such phenomena calls for the use of stochastic simulations. However, being computationally intensive, such simulations become infeasible for large and complex reaction networks. To improve their computational efficiency in handling these networks, we present a hybrid approach, in which slow reactions and fluxes are handled through exact stochastic simulation and their fast counterparts are treated partially deterministically through chemical Langevin equation. The classification of reactions as fast or slow is accompanied by the assumption that in the time-scale of fast reactions, slow reactions do not occur and hence do not affect the probability of the state. Our new approach also handles reactions with complex rate expressions such as Michaelis–Menten kinetics. Fluxes which cannot be modeled explicitly through reactions, such as flux of from endoplasmic reticulum to the cytosol through inositol 1,4,5-trisphosphate receptor channels, are handled deterministically. The proposed hybrid algorithm is used to model the regulation of the dynamics of cytosolic calcium ions in mouse macrophage RAW 264.7 cells. At relatively large number of molecules, the response characteristics obtained with the stochastic and deterministic simulations coincide, which validates our approach in the limit of large numbers. At low doses, the response characteristics of some key chemical species, such as levels of cytosolic calcium, predicted with stochastic simulations, differ quantitatively from their deterministic counterparts. These observations are ubiquitous throughout dose response, sensitivity, and gene-knockdown response analyses. While the relative differences between the peak-heights of the cytosolic time-courses obtained from stochastic (mean of 16 realizations) and deterministic simulations are merely 1%–4% for most perturbations, it is specially sensitive to levels of (relative difference as large as 90% at very low ).

We would like to acknowledge the UCSD Triton Resource of San Diego Supercomputer Center (SDSC) used in this work. This research was supported by the National Heart, Lung and Blood Institute (NHLBI) Grant No. 5 R33 HL087375-02 (S.S.), National Science Foundation (NSF) Grant No. DBI-0641037 (S.S.), and the NSF collaborative Grant No. DBI-0835541 (S.S.). Conceived the modeling study: S.S., D.M.T., and M.R.M. Developed the algorithm for hybrid stochastic simulation: T.J.C. and M.R.M. Wrote the computer program and analyzed the results: T.J.C. and M.R.M. Wrote the paper: T.J.C., M.R.M., D.M.T., and S.S. Supervised the overall research: S.S., D.M.T., and M.R.M.

I. INTRODUCTION

II. DYNAMICS OF CYTOSOLIC CALCIUM

A. Biological mechanisms and pathways

B. Mathematical representations of calcium dynamics

III. MATERIALS AND METHODS

A. The mathematical model of cytosolic calcium dynamics

B. Comparison of computational efficiency of stochastic simulation algorithms

C. A multiscale hybrid approach

1. Multiscale approach

2. Deterministic modeling of nonreaction fluxes

3. Reactions with complex rate expressions

D. Application to cytosolic calcium dynamics in RAW cells

IV. RESULTS

A. Dose response

B. Convergence of stochastic simulations at low doses

C. Random variability of the response at low doses

D. Sensitivity analysis

E. Calcium response to protein knockdown

V. SUMMARY AND DISCUSSION

A. Methodological components

B. Sensitivity analysis

C. Knockdown (KD) analysis

D. Stochastic effects at low molecular numbers

E. Deriving statistics from stochastic simulation

VI. SUPPLEMENTARY MATERIAL

### Key Topics

- Calcium
- 69.0
- Biochemical reactions
- 17.0
- Proteins
- 17.0
- Reaction kinetics modeling
- 15.0
- Reaction rate constants
- 10.0

## Figures

A simplified model for calcium signaling including calcium influx, ER, and mitochondrial exchange and storage [diagram in panel B taken from Maurya and Subramaniam (Ref. 3) *with permission from Biophysical Journal*]. (a) Ligand Complement 5a (C5a) binds to its receptor on plasma membrane (PM) and activates G protein . The free binds to and increases its activity which accelerates the phosphorylation of into and DAG. binds to its receptor on the ER membrane. Thus, calcium from the ER is released into the cytosol. Other fluxes between cytosol and mitochondria or ECM are also shown. (b) Receptor module (box 1), GTPase cycle module (box 2), generation module (box 3), and feedback module (box 4); ECM, extracellular matrix; , phosphatidylinositol 4,5-bisphosphate; , inositol 1,4,5-trisphosphate; , a lumped product of phosphorylation; , cytosolic ; Pr, proteins; ER, endoplasmic reticulum; ATP, adenosine triphosphate; ADP, adenosine diphosphate; SERCA, sarco(endo) plasmic reticulum calcium ATPase; PMCA, plasma membrane calcium ATPase; NCX, exchanger; L, ligand C5a; R, receptor C5aR; GRK, G-protein-coupled receptor kinase; CaM, calmodulin; , phospholipase C-; GAP, GTPase activating protein; RGS, regulator of G-protein signaling; DAG, diacylglycerol; PKC, protein kinase C; Pi, phosphate.

A simplified model for calcium signaling including calcium influx, ER, and mitochondrial exchange and storage [diagram in panel B taken from Maurya and Subramaniam (Ref. 3) *with permission from Biophysical Journal*]. (a) Ligand Complement 5a (C5a) binds to its receptor on plasma membrane (PM) and activates G protein . The free binds to and increases its activity which accelerates the phosphorylation of into and DAG. binds to its receptor on the ER membrane. Thus, calcium from the ER is released into the cytosol. Other fluxes between cytosol and mitochondria or ECM are also shown. (b) Receptor module (box 1), GTPase cycle module (box 2), generation module (box 3), and feedback module (box 4); ECM, extracellular matrix; , phosphatidylinositol 4,5-bisphosphate; , inositol 1,4,5-trisphosphate; , a lumped product of phosphorylation; , cytosolic ; Pr, proteins; ER, endoplasmic reticulum; ATP, adenosine triphosphate; ADP, adenosine diphosphate; SERCA, sarco(endo) plasmic reticulum calcium ATPase; PMCA, plasma membrane calcium ATPase; NCX, exchanger; L, ligand C5a; R, receptor C5aR; GRK, G-protein-coupled receptor kinase; CaM, calmodulin; , phospholipase C-; GAP, GTPase activating protein; RGS, regulator of G-protein signaling; DAG, diacylglycerol; PKC, protein kinase C; Pi, phosphate.

Temporal evolution of the concentrations of substrate and product computed using the Gillespie, tau-leap, and CLE approaches. (a) shows time-course of one realization from each method. (b) and (c) show the time-course of mean and standard deviation from 1024 realizations, which show excellent agreement among the three different methods. (d)–(f) show histograms and probability distribution of the number of molecules of S sampled at . The shapes of the three histograms are very similar.

Temporal evolution of the concentrations of substrate and product computed using the Gillespie, tau-leap, and CLE approaches. (a) shows time-course of one realization from each method. (b) and (c) show the time-course of mean and standard deviation from 1024 realizations, which show excellent agreement among the three different methods. (d)–(f) show histograms and probability distribution of the number of molecules of S sampled at . The shapes of the three histograms are very similar.

Dose response. (a) Comparison between ensemble average of 16 realizations and deterministic simulation. For better contrast, the time-course from deterministic simulation is shifted by 100 s. (b) Comparison between ensemble average and individual realizations in stochastic simulation for 0.1% (of 30 nM) strength of the ligand C5a. (c) Comparison of the dose response (peak-heights): The difference is quite small as compared to the scale of peak-height. (d) At lower doses, the NRD is larger indicating the stochastic effects. The NRD decreases with increasing dose as the number of the molecules of C5a becomes several hundreds or more.

Dose response. (a) Comparison between ensemble average of 16 realizations and deterministic simulation. For better contrast, the time-course from deterministic simulation is shifted by 100 s. (b) Comparison between ensemble average and individual realizations in stochastic simulation for 0.1% (of 30 nM) strength of the ligand C5a. (c) Comparison of the dose response (peak-heights): The difference is quite small as compared to the scale of peak-height. (d) At lower doses, the NRD is larger indicating the stochastic effects. The NRD decreases with increasing dose as the number of the molecules of C5a becomes several hundreds or more.

Revelation of stochastic effects at low doses. (a)–(d) Distributions of peak-height for the 0.1% dose of C5a computed from 16, 64, 256, and 512 realizations, respectively. The dotted vertical line represents the mean value and the solid curves denote theoretical Gaussian distributions. As the number of realizations increases, the shape of the histogram approaches a Gaussian distribution. (e)–(h) The mean of 4, 8, 16, or 32 realizations was computed and 1024 such mean values were generated. All the four histograms are similar to a Gaussian distribution and the standard deviation from these distributions indeed decreased proportional to , being the number of realizations used to compute the mean. (i) The standard deviation computed from 16 realizations for several doses. Contrary to the expectation, higher doses result in larger absolute standard deviations. (j) The normalized standard deviation decreases as the dose is increased, signifying the effect of randomness at lower doses.

Revelation of stochastic effects at low doses. (a)–(d) Distributions of peak-height for the 0.1% dose of C5a computed from 16, 64, 256, and 512 realizations, respectively. The dotted vertical line represents the mean value and the solid curves denote theoretical Gaussian distributions. As the number of realizations increases, the shape of the histogram approaches a Gaussian distribution. (e)–(h) The mean of 4, 8, 16, or 32 realizations was computed and 1024 such mean values were generated. All the four histograms are similar to a Gaussian distribution and the standard deviation from these distributions indeed decreased proportional to , being the number of realizations used to compute the mean. (i) The standard deviation computed from 16 realizations for several doses. Contrary to the expectation, higher doses result in larger absolute standard deviations. (j) The normalized standard deviation decreases as the dose is increased, signifying the effect of randomness at lower doses.

Sensitivity analysis. (a)–(c) *Response of* *to changes in* . The decrease in the peak-height due to decrease in is much more pronounced than that caused by the same decrease of IC:[R]. (c) NRD is extremely high at very low , suggesting significant stochastic effects at low numbers of molecules of . (b) and (c) also show the effect of using different numbers of realizations for computing the mean. Such differences are small (see text) indicating that 16 realizations are sufficient for computing the mean.

Sensitivity analysis. (a)–(c) *Response of* *to changes in* . The decrease in the peak-height due to decrease in is much more pronounced than that caused by the same decrease of IC:[R]. (c) NRD is extremely high at very low , suggesting significant stochastic effects at low numbers of molecules of . (b) and (c) also show the effect of using different numbers of realizations for computing the mean. Such differences are small (see text) indicating that 16 realizations are sufficient for computing the mean.

Knockdown response of . [(a)–(b)] The response to the 50%, 80%, 90%, and 99% knockdown of for 0.1% and 10% levels of IC:[R], respectively. As the knockdown rate of increases, both the basal level and peak-height of decrease because the production decreases due to decrease in . (c) Peak-height of response corresponding to different combinations of the and IC:[R] levels. Peak-height increases with high amount of IC:[R] and . (d) NRD is insignificant and decreases as doses of R and increase.

Knockdown response of . [(a)–(b)] The response to the 50%, 80%, 90%, and 99% knockdown of for 0.1% and 10% levels of IC:[R], respectively. As the knockdown rate of increases, both the basal level and peak-height of decrease because the production decreases due to decrease in . (c) Peak-height of response corresponding to different combinations of the and IC:[R] levels. Peak-height increases with high amount of IC:[R] and . (d) NRD is insignificant and decreases as doses of R and increase.

Knockdown response of GRK. [(a) and (b)] The response to the 50%, 80%, 90%, and 99% knockdown of GRK for 0.1% and 10% levels of IC:[R], respectively. (c) Peak-height of response corresponding to different combinations of [GRK] and IC:[R] levels. (d) NRD is insignificant, reaching its maximum of about 1.5% at low IC:[R].

Knockdown response of GRK. [(a) and (b)] The response to the 50%, 80%, 90%, and 99% knockdown of GRK for 0.1% and 10% levels of IC:[R], respectively. (c) Peak-height of response corresponding to different combinations of [GRK] and IC:[R] levels. (d) NRD is insignificant, reaching its maximum of about 1.5% at low IC:[R].

The response to the simultaneous knockdown of GRK and gene/protein related to . Knockdown of GRK and reduction of have opposite effects on the response. The response is much more sensitive to knockdown of GRK than to decrease in .

The response to the simultaneous knockdown of GRK and gene/protein related to . Knockdown of GRK and reduction of have opposite effects on the response. The response is much more sensitive to knockdown of GRK than to decrease in .

## Tables

The run-time scalability of the Gillespie, tau-leap, and chemical Langevin equation algorithms as a function of the number of molecules.

The run-time scalability of the Gillespie, tau-leap, and chemical Langevin equation algorithms as a function of the number of molecules.

Criteria used to identify slow and fast reactions and corresponding numerical method. Columns 2 and 3 list the scale and simulation method in the scale (method) format.

Criteria used to identify slow and fast reactions and corresponding numerical method. Columns 2 and 3 list the scale and simulation method in the scale (method) format.

Summary of results of KD response. The change in the features of calcium response listed is for increase in KD-level (decrease in IC:[.] of the protein). Qualitative nature of the features is mostly independent of the level of [R].

Summary of results of KD response. The change in the features of calcium response listed is for increase in KD-level (decrease in IC:[.] of the protein). Qualitative nature of the features is mostly independent of the level of [R].

Article metrics loading...

Full text loading...

Commenting has been disabled for this content