ARTICLE IN PRESS
Neurocomputing 69 (2006) 1066–1069 www.elsevier.com/locate/neucom
A novel approach to model neuronal signal transduction using stochastic differential equations Tiina Manninena,b,, Marja-Leena Linneb, Keijo Ruohonena a
Institute of Mathematics, Tampere University of Technology, P.O. Box 553, FI-33101 Tampere, Finland Institute of Signal Processing, Tampere University of Technology, P.O. Box 553, FI-33101 Tampere, Finland
b
Available online 7 February 2006
Abstract We introduce a new approach to model the behavior of neuronal signal transduction networks using stochastic differential equations. We present first a mathematical formulation for a stochastic model of protein kinase C pathway. Different kinds of numerical integration methods, including the explicit and implicit Euler–Maruyama methods, are used to solve the Ito^ form of the stochastic model. Stochastic models may provide more realistic representations for the chemical species in signal transduction networks compared to deterministic models. Our approach has the advantage of being computationally less demanding in the context of large-scale stochastic simulations than other approaches where individual chemical interactions are simulated stochastically. r 2006 Elsevier B.V. All rights reserved. Keywords: Stochastic differential equation; Itoˆ form; Euler–Maruyama method; Protein kinase C; Signal transduction
1. Introduction In recent years, stochastic differential equation (SDE) models have been increasingly used in many research fields, such as economics, finance, chemistry, and biology. It is known that all molecular reactions in living organisms are essentially random processes and stochastic models provide therefore a more detailed description of reactions than deterministic models [12]. To our knowledge, SDE models have not been extensively utilized in neuroscience to model signal transduction networks. In this work, we model the behavior of neuronal signal transduction networks using the mass action law (see also [8]). To obtain SDE models, we incorporate stochasticity into deterministic differential equation models describing macroscopic behavior. Another distinct way of incorporating stochasticity is to model individual chemical interactions as Markov processes, using for example the Gillespie algorithm, to describe microscopic behavior. In this paper, we introduce a new mathematical formulation for incorporating stochasticity into reaction rates and assume constant-intensity Gaussian white noise. Corresponding author.
E-mail address: tiina.manninen@tut.fi (T. Manninen). 0925-2312/$ - see front matter r 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.neucom.2005.12.047
Previously, we have studied a deterministic model of protein kinase C (PKC) signal transduction pathway [1,8,10]. The pathway is interesting because it has been found to be important in memory formation in neurons [1,2]. The model consists of 10 reversible reactions and 15 different interacting chemical species, types of molecules and ions. Three of the model species correspond to second messengers: calcium ion ðCa2þ Þ, arachidonic acid (AA), and diacylglycerol (DAG). The model has six computational intermediates which enable the study of concentrations before they are summed to the end product, active PKC (PKCa). We have also implemented and simulated a stochastic model for a larger network which is composed of four kinases (PKC, mitogen-activated protein kinase (MAPK), protein kinase A (PKA), and calmodulin-dependent protein kinase II (CaMKII)) and numerous regulatory pathways involved in synaptic signaling [1,10]. However, in this paper we concentrate on the PKC pathway model as an example to illustrate our stochastic approach. 2. Methods We investigate the numerical solutions for the Ito^ [6,9] form of SDEs in the context of signal transduction
ARTICLE IN PRESS T. Manninen et al. / Neurocomputing 69 (2006) 1066–1069
networks. The Ito^ SDE is given by
0.35
dXðtÞ ¼ SvðK; XðtÞ; YðtÞÞ dt þ b dWðtÞ,
0.30
Xð0Þ ¼ X0 ,
0.25
where the deterministic equation is modified by adding a stochastic component b dW. In (1), X: Rþ ! Rm describes m variables (concentrations of chemical species) and X0 describes the initial concentration values. Symbol Rþ denotes non-negative real numbers. Function Y: Rþ ! Rk describes k model inputs (concentrations of second messengers). Matrix S 2 Zmn is the stoichiometric matrix k and n is the number of reactions. Function v: Rm þ Rþ ! n R describes the deterministic reaction rates. The reaction rates depend on rate constants K, variables X, and model inputs Y according to the mass action law. Herepffiffib is a standard deviation parameter in units of M= s and WðtÞN pffiffi m ð0; tIÞ is m-dimensional Brownian motion in units of s. The solution for (1) at time t is given by Z t Z t XðtÞ ¼ X0 þ SvðK; XðsÞ; YðsÞÞ ds þ b dWðsÞ: ð2Þ
0.20
0.10
-0.05 -0.10
The behavior of the deterministic PKC model has been simulated using different kinds of inputs for Ca2þ , AA, and DAG [8]. It is known that concentrations of many intracellular molecules and ions, such as Ca2þ ions, behave in a sinusoidal manner [4,7,11]. Several studies indicate that intracellular Ca2þ concentration varies between 107 and 105 M [3,4,7,11]. Furthermore, calcium oscillations with period of 10 s to 2 min have been observed. We assume a similar behavior for AA and DAG. Therefore, artificially generated sine waves are used here as inputs, for example, Yi ¼ ðsin 2pt=pi þ 1:1Þ 106 M where period pi is 50 s for Ca, 10 s for AA, and 100 s for DAG. The simulated deterministic responses of PKCa are used here as a reference for the stochastic realizations of PKCa. The simulation conditions, including initial concentrations and rate constants, are those used in [8]. In this study, all tested integration methods are found to produce identical results with same Dt, as expected. Here the implicit Euler–Maruyama method is used. The first issue to address is to determine an appropriate value for b. Fig. 1 shows simulation results for the concentration of
0
10
20
30
40
50 t (s)
60
70
80
90
100
Fig. 1. Varying value of parameter b in stochastic PKC model. Dt is 214 s pffiffi and b has values of ½1 5 10 109 M= s.
0.30 Sample mean of 10000 realizations 3 individual realizations
0.25 0.20 0.15 0.10 0.05 0
3. Results
β is 1⋅10-9 β is 5⋅10-9 β is 10⋅10-9
0
0
The simplest integration method for solving (2) is the Euler–Maruyama method [5,6,12]. The efficiency of different integration methods depends on the SDE model, for example, is the model stiff or nonstiff. A model is stiff if it has multiple widely differing time scales. Implicit methods are more suitable for stiff SDEs. We use different kinds of integration methods, including the explicit and implicit Euler–Maruyama methods [5,6], to solve our stochastic model. The SDE model and numerical integration methods are implemented and simulated in MATLABs programming environment.
0.15
0.05
[PKCa] (µM)
0
[PKCa] (µM)
ð1Þ
1067
0
10
20
30
40
50 t (s)
60
70
80
90
100
Fig. 2. Sample mean of 10,000 realizations and three individual realizations of PKCa in stochastic PKC model. Dt is 214 s and b is pffiffi 109 M= s.
PKCa ([PKCa]) using different values of parameter b. We conclude that the larger the value of b, the larger the variation of the stochastic realizations of PKCa around the deterministic one. Occasionally, large values of b are found to give negative concentration values. Compared to the deterministic model, also high concentration values can occur. Therefore, the values of b for Dt ¼ p 214 ffiffi s are 9 predicted to be lower than or equal to 10 M= s. In the future, experimentally obtained time-series data will be used as a basis of determining suitable values of b and Dt. Next we simulate the stochastic PKC model 10,000 times in order to validate that different realizations of the stochastic PKC model are bounded and exhibit similar type of behavior than the pffiffideterministic response. We use 214 s as Dt and 109 M= s as b. Fig. 2 shows the sample
ARTICLE IN PRESS T. Manninen et al. / Neurocomputing 69 (2006) 1066–1069
1068
unknown values of rate constants based on obtained time-series data. Later, stochastic signal transduction models will be integrated with excitability models to understand information processing in neurons.
0.25
[PKCa] (µM)
0.20
Acknowledgements
0.15
This study was financially supported by Tampere University of Technology Graduate School, the Academy of Finland, Grant 104508, as well as the Ulla Tuominen Foundation and the Foundation of Technology for T.M., and the Academy of Finland, Grants 104508, 106030, and 107694 for M.-L.L.
0.10
0.05
0
0
50
100 150 200 250 300 350 400 450 500 t (s)
Fig. 3. Sample mean and sample standard deviation of 10,000 realizations pffiffi of PKCa in stochastic PKC model. Dt is 214 s and b is 109 M= s.
^ of the 10,000 realizations and three individual mean ðlÞ realizations. Furthermore, Fig. 3 shows the sample mean ^ The sample standard and sample standard deviation ðrÞ. deviation of the 10,000 realizations is calculated in every timepoint, but in Fig. 3 the sample standard deviation is ^ m^ þ s] ^ in every 20 s. We verify that plotted with bars [m^ s, the deterministic solution and the sample mean of the stochastic solutions for PKCa are identical. Fig. 3 also shows that the sample standard deviation increases when simulation proceeds. 4. Discussion and conclusions In this paper, we introduce a mathematical formulation for incorporating stochasticity into deterministic signal transduction models. Based on our studies, SDE models may provide more realistic representations for the chemical species in signal transduction networks when compared to deterministic models. The SDE approach, with appropriate choices of b and Dt, may prove useful in modeling systems where the number of interacting chemical species is low. Specifically, the SDE approach introduced in this work is computationally more efficient compared to other approaches where individual chemical interactions are simulated stochastically. In the future, experimentally obtained time-series data will assist in determining where the stochasticity should be incorporated and what values of timestep and parameter b should be used. If biological noise has larger variance than model realizations, the model needs to be improved, perhaps by using a varying b, to guarantee that negative concentration values are not possible. Furthermore, possible instability points need to be prevented mathematically taking into account experimentally obtained data. The development of stochastic parameter estimation methods will help in automatically determining the
References [1] U.S. Bhalla, R. Iyengar, Emergent properties of networks of biological signaling pathways, Science 283 (1999) 381–387 hhttp:// doqcs.ncbs.res.in/i. [2] T.V.P. Bliss, G.L. Collingridge, A synaptic model of memory: long-term potentiation in the hippocampus, Nature 361 (1993) 31–39. [3] E. De Schutter, P. Smolen, Calcium dynamics in large neuronal models, in: C. Koch, I. Segev (Eds.), Methods in Neuronal Modeling: From Ions to Networks, second ed., The MIT Press, Cambridge, 2001, pp. 211–250. [4] G.W. De Young, J. Keizer, A single-pool inositol 1,4,5,-triphosphate-receptor-based model for agonist-stimulated oscillations in Ca2þ concentration, Proc. Natl. Acad. Sci. USA 89 (1992) 9895–9899. [5] D.J. Higham, An algorithmic introduction to numerical simulation of stochastic differential equations, SIAM Rev. 43 (3) (2001) 525–546. [6] P. Kloeden, E. Platen, Numerical Solution of Stochastic Differential Equations, third ed., Springer, Berlin, 1999. [7] A.M. Lawrie, M.E. Graham, P. Thorn, D.V. Gallacher, R.D. Burgoyne, Synchronous calcium oscillations in cerebellar granule cells in culture mediated by NMDA receptors, NeuroReport 4 (1993) 539–542. [8] T. Manninen, A. Saarinen, M.-L. Linne, Simulation study of differential equation model for protein kinase C signaling, Report 3, Institute of Signal Processing and Institute of Mathematics, Tampere University of Technology, Finland, 2004. hwww.cs.tut.fi/ manninet/Report2004.pdfi. [9] B. Øksendal, Stochastic Differential Equations, fifth ed., Springer, Berlin, 1998. [10] S. Sivakumaran, S. Hariharaputran, J. Mishra, U.S. Bhalla, The Database of Quantitative Cellular Signaling: management and analysis of chemical kinetic models of signaling networks, Bioinformatics 19 (2003) 408–415. [11] Y. Tang, H.G. Othmer, Frequency encoding in excitable systems with applications to calcium oscillations, Proc. Natl. Acad. Sci. USA 92 (1995) 7869–7873. [12] T.E. Turner, S. Schnell, K. Burrage, Stochastic approaches for modelling in vivo reactions, Comput. Biol. Chem. 28 (2004) 165–178. Tiina Manninen received her MS in mathematics from Tampere University of Technology, Finland, in 2003. She is presently doing her Ph.D. studies in the Institute of Mathematics and the Institute of Signal Processing at Tampere University of Technology. She is specializing in modeling and simulation of signal transduction in neurons with an emphasis on automated model parameter estimation.
ARTICLE IN PRESS T. Manninen et al. / Neurocomputing 69 (2006) 1066–1069
Marja-Leena Linne received her MS and Ph.D. from Tampere University of Technology, Finland, in 1993 and 2001, respectively. She is presently Research Fellow of the Academy of Finland concentrating on computational neuroscience and systems biology research in the Institute of Signal Processing at Tampere University of Technology.
1069
Keijo Ruohonen is professor of applied mathematics in the Institute of Mathematics at Tampere University of Technology, Finland. He received his BS, MS, and Ph.D. in mathematics from the University of Turku, Finland. His research interests include computability and algorithms, and statistics.