Information transfer in signaling pathways: A study using coupled simulated and experimental data
 Jürgen Pahle^{1, 4}Email author,
 Anne K Green^{2},
 C Jane Dixon^{3} and
 Ursula Kummer^{1, 4}
DOI: 10.1186/147121059139
© Pahle et al; licensee BioMed Central Ltd. 2008
Received: 22 November 2007
Accepted: 04 March 2008
Published: 04 March 2008
Abstract
Background
The topology of signaling cascades has been studied in quite some detail. However, how information is processed exactly is still relatively unknown. Since quite diverse information has to be transported by one and the same signaling cascade (e.g. in case of different agonists), it is clear that the underlying mechanism is more complex than a simple binary switch which relies on the mere presence or absence of a particular species. Therefore, finding means to analyze the information transferred will help in deciphering how information is processed exactly in the cell. Using the informationtheoretic measure transfer entropy, we studied the properties of information transfer in an example case, namely calcium signaling under different cellular conditions. Transfer entropy is an asymmetric and dynamic measure of the dependence of two (nonlinear) stochastic processes. We used calcium signaling since it is a wellstudied example of complex cellular signaling. It has been suggested that specific information is encoded in the amplitude, frequency and waveform of the oscillatory Ca^{2+}signal.
Results
We set up a computational framework to study information transfer, e.g. for calcium signaling at different levels of activation and different particle numbers in the system. We stochastically coupled simulated and experimentally measured calcium signals to simulated target proteins and used kernel density methods to estimate the transfer entropy from these bivariate time series. We found that, most of the time, the transfer entropy increases with increasing particle numbers. In systems with only few particles, faithful information transfer is hampered by random fluctuations. The transfer entropy also seems to be slightly correlated to the complexity (spiking, bursting or irregular oscillations) of the signal. Finally, we discuss a number of peculiarities of our approach in detail.
Conclusion
This study presents the first application of transfer entropy to biochemical signaling pathways. We could quantify the information transferred from simulated/experimentally measured calcium signals to a target enzyme under different cellular conditions. Our approach, comprising stochastic coupling and using the informationtheoretic measure transfer entropy, could also be a valuable tool for the analysis of other signaling pathways.
Background
The simulation of complex biochemical networks has become very important to gain insight into the dynamic behavior of cellular processes [1–3]. Signaling pathways, in particular, often evade intuitive, and therefore rather static, explanations because of their highly nonlinear dynamics and many crosslinks. However, despite the emergence of sophisticated highthroughput and in vivo imaging techniques, there is still a lack of highquality singlecell multivariate data.
Such data would be very helpful in elucidating the nuts and bolts of many signaling mechanisms. In this study we use calcium signaling as an example. Calcium signaling represents one of the most versatile secondmessenger pathways and, in many cell types, Ca^{2+} (calcium) ions control a variety of cell functions from fertilization, secretion, enzyme activation and gene expression to cell death [4, 5].
An intriguing fact is that, even in nonexcitable cells like hepatocytes, the concentration of cytosolic calcium can display regular (spiking) or more complex (bursting) oscillations or prolonged elevated levels [6] after stimulation by an agonist and depending on the nature of this agonist. This oscillatory behavior is not only believed to save the cell from the toxic effects of sustained high cytosolic calcium levels and from desensitization, but has also been shown to increase the efficiency of calcium signaling [7]. In addition to these temporal patterns of calcium dynamics, interesting spatiotemporal patterns (e.g. calcium puffs and waves) have been described [4, 8]. However, we will concern ourselves in this study with temporal patterns only.
Due to its importance for the functioning of many cell types and its interesting dynamics [9], calcium signal transduction has attracted numerous theoretical studies. Many different models of calcium signaling have been proposed, ranging from simple onepool models [10] to more elaborate ones [11] incorporating many different processes. For a review on calcium models, see Schuster et al. 2002 [12].
Since a range of different agonists such as hormones (e.g. vasopressin) or nucleotides (e.g. ATP) trigger calcium responses and, on the other hand, a range of different targets (e.g. Ca^{2+} dependent proteins such as calmodulin, CaM kinase II, protein kinase C, phosphorylase kinase or transcription factors e.g. NFAT or NFκ B) exist in the cell [13], specific information is likely to be encoded in the calcium signal and decoded again later on. It has been proposed that information might be encoded in the amplitude, frequency, duration, waveform or timing of calcium oscillations and the search for this calcium code has attracted a number of experimental and theoretical studies (for a review, see [14]).
On the experimental side, mainly the frequency decoding of spiking calcium oscillations has been examined. De Koninck and Schulman 1998 [15] demonstrated the sensitivity of immobilized CaM kinase II to Ca^{2+} oscillation frequency by in vitro rapid superfusion. Li et al. 1998 [16] found that NFAT is activated optimally at a Ca^{2+} oscillation frequency of about 1/min and Dolmetsch et al. 1998 [7] studied the differential regulation of Tcell NFAT and NFκ B by Ca^{2+}oscillations of different frequencies. The interesting work of Oancea and Meyer 1998 [17] describes the activation of protein kinase C γ (PKC γ) by diacylglycerol (DAG) combined with highfrequency Ca^{2+} spikes, which points to a joint code of calcium and DAG in that case.
Most theoretical studies also limit themselves to the spiking mode of calcium oscillations. Dupont et al. 2003 [18] could successfully reproduce the findings of [15] in a model. Gall et al. 2000 [19] examined the activation of liver glycogen phosphorylase by modeling a de/phosphorylation cycle. Salazar et al. 2004 [20] studied the activation of target proteins by Ca^{2+} oscillations in terms of efficiency, speed and specificity. Marhl et al. [21, 22] investigated the decoding of timelimited calcium oscillations by downstream proteins. Recently, the bursting mode of Ca^{2+}oscillations has been investigated by Larsen et al. 2004 [11] and Schuster et al. 2005 [23]. Using a simple model of calcium oscillations [11] and artificially generated calcium bursts [23] respectively to drive protein activation, these studies showed that specific information can be encoded in the waveform of bursting oscillations and thus that different proteins can be activated differentially at the same time. Rozi and Jia [24] studied the activation of glycogen phosphorylase by spiking as well as bursting calcium oscillations.
Even though informationtheoretic measures [25] are in widespread use for physiological data [26] and neural information transfer [27], their application to biochemical systems is restricted to only relatively few studies. For instance, Prank et al. 1998 [28] and Kropp et al. 2005 [29] studied the encoding of hormonal signals in intracellular calcium signals using the socalled coding fraction and mutual information. The authors drive a deterministic model of calcium spiking with a specific form of generated noise and estimate the amount of information transferred. In [30] the same group couples a deterministic model of CaM kinase activation to experimentally measured data from HIT (hamster insulinsecreting tumor)cells, but here the results are not analyzed in an informationtheoretic manner.
We propose to use the informationtheoretic measure transfer entropy [31] to estimate the information transferred by spiking or bursting calcium oscillations under different conditions. Transfer entropy has advantages over conventional measures such as (timelagged) correlations, in that it detects all statistical dependencies (linear and nonlinear), it is asymmetric, i.e. it distinguishes between information source and target, and it considers shared information due to a common history of the source and target by using conditioned transition probabilities. Transfer entropy has been applied to physiological data [26, 32], financial time series [33], geological data [34] and others [35, 36], but, so far, not to biochemical data. We used both simulated and experimentally measured time series for the estimation of transfer entropy. The simulated data were generated by a stochastic version of the simple calcium oscillations model in [37], extended by a stochastically simulated activation of target protein. We set up a framework for stochastic simulation of the calcium system, stochastic coupling of the enzyme activation process and estimation of the transfer entropy using kernel density estimation methods. We used this framework to investigate calcium information transfer in systems with different levels of activation and particle numbers.
Since multivariate experimental data is scarce, we devised a method, inspired by hybrid deterministic/stochastic simulation techniques, which allows the stochastic coupling of the enzyme activation process to arbitrary univariate calcium time series. We took experimental data from singlecell measurements on rat hepatocytes and coupled the activation of the stochastically simulated enzyme to them in order to get bivariate data. Finally we used these semiexperimental data as input for the estimation of the information transfer.
Results
The concentration of the active form of a simulated Ca^{2+}dependent enzyme, which was stochastically coupled to the calcium data, is also shown. We implemented a stochastic coupling scheme to be able to couple the simulated enzyme to arbitrary, simulated or experimental, calcium time series. This method is described in detail in Methods.
We also tested our estimation process by using surrogate data (constrained realizations, [39]). We estimated the transfer entropy of time series in which the temporal order of the calcium signal was destroyed by shuffling the samples (data not shown). This removed all dependencies, while the marginal probability distributions were preserved. Indeed, here the estimated transfer entropy showed values near zero (~0.02 – 0.07).
In the case of under or overstimulation (k_{2} = 1 or k_{2} = 3.2), the system is in a (noisy) steady state and this results in low values for the transfer entropy. For k_{2} = 1 the calcium concentration is near its resting level, which is far below the K_{ M }value of the enzyme. No enzyme gets activated and no information can be transferred. For k_{2} = 3.2 the calcium steady state lies above the enzyme's K_{ M }value and the amount of active enzyme reaches its maximum. In contrast to understimulation, here the information transfer is not exactly zero, even though it takes low values of ~0.2. The reason for this is that now the noisy steady state is near the K_{ M }value of the enzyme and it can pick up some random fluctuations in calcium concentration. If the system is in an oscillatory mode, such as spiking (k_{2} = 2) or bursting (k_{2} = 2.5, 2.85), the transfer entropy increases with increasing volume until it seems to flatten out for volumes above 5 × 10^{10}, as shown above.
An interesting effect can be observed for k_{2} = 2.99 and k_{2} = 3, where the deterministic limits of the calcium dynamics are elevated oscillations and an elevated steady state, respectively. However, the stochastic system shows irregular behavior with small volumes. For high volumes, oscillations are observed even for k_{2} = 3. For both parameter values, the generally very high level of transfer entropy is due to the position of the center of their oscillations. It is near the K_{ M }value of the enzyme, so that the enzyme is responsive even to minute variations in the calcium level. Interestingly, for k_{2} = 3 the transfer entropy shows a maximum at the volumes 1 × 10^{9} and 5 × 10^{10}. An explanation for this effect is that for the higher volume 5 × 10^{9}, the system is already near the deterministic limit, which is just a rather uninteresting elevated steady state with relatively low information transfer. On the other hand, for smaller volumes, the information transfer gets degraded because of increasing stochastic fluctuations. Those fluctuations are especially pronounced in this parameter range, because the sensitivity of the system (measured by the divergence) is high (see [40] for details). Those two opposed trends lead to a maximum in a range where the system is still oscillatory, but not yet too noisy.
If we look at the estimates for a volume of 5 × 10^{9} only (the biggest systems considered in this study), there is a slight increase in estimated transfer entropy from spiking to increasingly complex bursting oscillations (see Table 1). The transfer entropy is very high for elevated oscillations near the enzyme's K_{ M }value and it drops to a very low value in the case of an elevated steady state, e.g. overstimulation. Intuitively, one would think that the information transfer should be correlated to the complexity (spiking, bursting or irregular oscillations) of the calcium oscillations, since more complex input signals can potentially carry more information. However, this can only be hinted at from our experiments. One should be wary not to overinterpret the absolute numbers, since we found them very much dependent on the estimation process used. Also, they are subject to statistical fluctuations. Furthermore, the enzyme is most sensitive for calcium levels near its K_{ M }value. For the input signal to generate a high information transfer, it is important to meet that range. The transfer entropy nicely detects this for the oscillatory regime with k_{2} = 2.99 and high volumes, where the oscillations exactly meet the K_{ M }value. Here the estimated transfer entropy is high, even though the dynamics is apparently less complex than in the bursting mode.
Different dynamical behavior and maximal transfer entropy values. Maximum values of the transfer entropy (TE) for different stimulation strengths k_{2} and their respective dynamic regime in the simulated system with volume 5 × 10^{9}.
k _{2}  Dynamic behavior  TE 

1  Understimulation  0.00 
2  Spiking  0.52 
2.5  Bursting  0.59 
2.85  Bursting  0.60 
2.99  Elevated oscillations  0.95 
3.2  Overstimulation  0.15 
The significantly higher transfer entropy values of the simulated system can partly be explained by the existence of two episodes in the experimental data without bursts (The calciummobilizing agonist was absent from the experimental medium for the duration of these two episodes). We removed these episodes and repeated the estimation which yielded a transfer entropy maximum of roughly 0.39 bit. An explanation for the remaining discrepancy is that the simulated bursts have a considerably longer duration than the bursts in real hepatocytes. Therefore, the calcium signal spends more time within the sensitive region of the enzyme (near the K_{ M }value) which clearly increases information transfer.
Discussion
In the following we will motivate the choice of several technical elements as well as discuss their strengths and limitations.
Stochastic coupling procedure
Stochastic fluctuations in cellular systems are not just random noise, but can even change the dynamics of the system [41] as was seen, for instance, in our simulations for parameter values near bifurcation points (k_{2} = 2.99 and small volumes). Therefore it is important to consider random effects (and the effects of the system size on those fluctuations) when modeling systems with relatively low particle numbers, e.g. signal transduction pathways.
It should be noted here that, even in those cases where stochastic effects do not change the dynamics significantly, deterministic coupling of a biochemical reaction system to experimental data [30] is not appropriate for our purposes. The estimation of transfer entropy diverges for coarsegrained continuous systems and increasing resolution if the coupling between the processes is deterministic [31]. Therefore our stochastic coupling scheme of the simulated enzyme to calcium time series is absolutely essential for this study.
Since the experimentally measured calcium concentration is only known at a discrete set of points in time and therefore we assumed it to be constant between two samplings, the coupling of a simulated enzyme to those time series can only be an approximation. However, it is apparent that, in the limit of a sampling frequency of the given time series near the frequency of reaction events in the system and a measurement resolution in the range of single particles, our method converges to the mathematically exact solution. For nearly every practical case, neither the number of samples nor the resolution will satisfy these theoretical conditions. To make sure that this fact did not compromise our results, we compared simulated data where the enzyme was only coupled to a calcium time series with data that was calculated by exact stochastic simulation of the whole system, i.e. calcium dynamics plus enzyme activation, and where no approximation was involved (data not shown). For the parameter values and sampling times we used, our results were not changed considerably by the approximate coupling.
One shortcoming of the stochastic coupling procedure described here is that it is a oneway process. Obviously, the input calcium time series is fixed and can not be changed during the process and so possible feedback of the target enzyme on the calcium system, e.g. calcium buffering by proteins or feedback via protein kinase C, has to be neglected.
Choice of model parameters
The volume of a hepatocyte is about 2 pL [42]. Assuming that the cytosol, where the free Ca^{2+} is located, takes up about half of the total volume of the cell and that, in the case of bursting, the calcium level peaks around 1 μ M, this results in a particle number of about 600 000. This particle number roughly corresponds to a volume of 10^{10} in the arbitrary units of the calcium model used. Therefore our results lie well in the range of physiologically meaningful parameters. Also the parameters of the simulated enzyme have been chosen to be, at least, biologically plausible. Most of the time calcium binding to enzymes occurs cooperatively, as e.g. with calmodulin. Calmodulin has four binding sites with high affinity (K_{ d }≈ 0.1 – 1 μ M) for Ca^{2+}. For this reason, we, like the authors of other numerical studies [11, 23, 30], employ a Hill term of 4th order. The K_{ M }value of the simulated enzyme lies between the calcium resting level and the amplitude of secondary peaks, in the case of bursting oscillations.
The reason for choosing this calcium model instead of a more detailed one was that, even though it is very simple in terms of size and kinetic functions, it can show both spiking and bursting behavior in addition to (elevated) steady states, thereby mimicking the dynamics of real cells after stimulation by different agonists (see [37] for details). Most other models cannot show bursting oscillations. It also was relatively easy to implement and fast to simulate stochastically. Nevertheless, the generation of some of the time series with high particle numbers required computation times in the range of several days. In fact, the purpose of this study was not to analyze this specific calcium model and therefore the approach presented here is not restricted to that model. It should also be mentioned that our framework can easily be applied to arbitrary enzyme regulation mechanisms, provided that they allow stochastic simulation of the Gillespie type, i.e. a propensity can be assigned to every possible event in the system.
One problem of the calcium model we used is that the amplitudes of the oscillations vary for different dynamic modes (see Fig. 1), whereas in real hepatocytes the amplitudes of calcium oscillations have been reported to be independent of the type of oscillations. Also the duration of bursts is longer than in experiments which, we believe, led to the discrepancy in transfer entropy between simulated and experimental data. To mitigate these issues we plan to use more realistic calcium models with more constant oscillation amplitudes, e.g. [11] in the future.
Estimation of transfer entropy
Often, (timelagged) correlations are used to quantify the coherence of two observables. However, correlations can only indicate linear relations, not nonlinear ones. Therefore mutual information has been developed which is sensitive to all statistical dependencies [43]. Unfortunately, this measure is still (like correlations) symmetric and cannot distinguish between information sources and targets.
The transfer entropy, on the other hand, is explicitly asymmetric because it uses conditioned transition probabilities. As stated by Schreiber [31, p. 461], "transfer entropy is able to distinguish effectively driving and responding elements and to detect asymmetry in the interaction of subsystems." In addition, the use of transition probabilities makes it a dynamic measure, meaning that it can account for the history of the processes. This, together with its ability to consider linear and nonlinear dependencies, renders it appropriate for use on nonlinear signal transduction systems.
We found that a major issue with this measure is that it is not trivial to estimate it from time series in a reliable way and that the estimation is quite dataintensive.
One crucial point is that the processes have to be ergodic to allow for the estimation of the probability densities from one time series alone. Also they must be Markovian. In other words, their histories of length k and l (see Methods), which are taken into account, must be longer than possible correlation times. This is very important, because the transfer entropy detects the deviation from the Markov property. One simple example where this condition would not be fulfilled is when we just reversed the direction and estimated the transfer entropy from the enzyme signal to the calcium signal (P_{act} → Ca^{2+}). We saw already that in our setting there can be no feedback from the enzyme to the calcium system and thus no information can be transferred this way. Because the transfer entropy is a directional measure and can distinguish between information transferred in one and the other direction, one would naively think that it should equal zero (plus statistical fluctuations) here. This, however, is not the case, because the calcium signal alone is not Markovian. In fact, in the model it is coupled to G_{ α }and PLC and their influence will lead to a transfer entropy which is not zero. There are two possible solutions to this issue: a) consider the whole system (Ca^{2+}, G_{ α }and PLC) or condition on all coupled subsystems, or b) take into account a long enough history for the processes in which all relevant information is already embedded.
In all practical applications of the transfer entropy, especially with purely experimental data, one has to fix the lengths of the two signal histories (k and l) with care. Since the characteristic timescale of autodependencies in measured data is not known a priori, they can not be regarded as stemming from an orderone Markov process. Therefore, one should estimate the transfer entropy using different values for the order of the underlying processes, and longer histories should be preferred. However, the often very limited amount of data renders this avenue infeasible in many cases, since kernel estimation would have to be applied to distribution functions in four and more dimensions. One possible resort here would be coarsegraining of the time series and the use of the discrete version of the transfer entropy. In the present study we restricted ourselves to the orderone case, the reason being that, in our case, the coupled protein is actually described by a Markov process of order one and is not dependent on previous values. Therefore, a history length of 1 (k = l = 1) suffices.
Kernel density estimation is known to be very dependent on the choice of a correct kernel bandwidth ε. Rules of thumb exist for the optimal bandwidth of (univariate) Gaussian kernels [44] which, however, are said to often lead to oversmoothing. Little has been done for multivariate kernels however. Instead of just using one bandwidth, we scanned the estimated transfer entropy over a range of different ε values and checked for the range of bandwidths where the estimates are independent of ε, e.g. a plateau is visible in the scans (Fig. 3). If there is a definite plateau, its values are simultaneously the maximal values of the scan. Due to this and because the estimated transfer entropy was observed to underestimate the true value [45], we chose to take the maxima of the scans as estimates of the transfer entropy.
The calcium signal and the enzyme signal have different ranges of values. Therefore we normalized the time series to have mean 0.0 and standard deviation 1.0 prior to the estimation, which allowed us to use the same ε in both spaces. This is justified, because the (continuous) transfer entropy is independent of coordinate transformations [45].
To improve our calculations, we used a Theiler window approach and excluded all estimates where only less than a required minimal number of neighbors could be found. This avoided spurious effects caused by temporal correlations and dampened statistical fluctuations, respectively. In this study we mainly used rectangular kernels. However, we also tried Gaussian kernels (data not shown), which did not change our results considerably.
Transfer entropy is an averaged measure, i.e. it describes the information transfer over the whole observation interval. We observed that periods in the experimental calcium time series without bursts (Fig. 6) decreased the overall transfer entropy. Therefore, if the processes under study are expected to show some kind of locking or unlocking episodes, which we would dub statistical locking, the measure would have to be calculated on smaller (disjoint or overlapping) windows in order to see possible changes over time. Care has to be taken, though, that the windows are big enough to get a sound statistical basis for the estimation.
We want to stress that the absolute values of our transfer entropy estimates are, of course, dependent on the parameters of the estimation procedure. In particular, the minimum number of neighbors needed for a sample to be considered plays a major role here. Setting this number to values greater than 1 helps to diminish statistical fluctuations, but can create a bias towards zero if there are not enough samples available. Therefore, one should be cautious when interpreting these values and should not mix results coming from different estimation procedures without justification. We only compared estimates where the estimation parameters, the type of kernel and the length of the input time series were the same. We attributed the discrepancy in estimated transfer entropy of simulated and experimentally measured data to lacking realism of the simple calcium oscillations model used. Hence, we note here that transfer entropy could very well be employed as a measure of realism for signaling pathway models. We envisage its use in biochemical modeling where models are optimized so as to have the same information transfer as observed in experiment.
Also, regarding the rates of information transfer we estimated in this study, one should be cautious. Even though they can provide a useful basis for hypotheses on the functioning of cellular signal transduction, it is not known what fraction of the information that can maximally be transferred is actually used by downstream cellular processes. Because it is not yet clear what features of the calcium signal really carry relevant information, we used an informationtheoretic approach. It potentially measures all the information from the calcium signal that can be found in the protein signal. In addition, this modelfree approach facilitates direct comparisons between simulated and experimentally measured data.
Nevertheless, specific information transferred from calcium to cellular processes could, in principle, be estimated by extending the simple model to include these processes under consideration and estimating the transfer entropy directly between calcium and the observables of these processes. This includes the detection, analysis and quantification of possible crosstalk between different signaling pathways.
A general framework
It should be mentioned that there are many potential variants and extensions of the estimation algorithm (simple or adaptive histograms, adaptive kernel density estimation, likelihood estimators and others [44]), which we could not cover here. However, regardless of the algorithm used, the basic strength of the informationtheoretic approach is that it is modelfree. This allows the direct comparison of simulated and experimental data.
Conclusion
In this study we combined methodologies from different fields in order to elucidate the cellular information transfer via Ca^{2+} signaling. The main ingredients we used are:

Modeling and simulation of calcium signal transduction, in particular stochastic approaches.

Stochastic coupling of a Ca^{2+}dependent protein to experimental and simulated data.

The socalled transfer entropy introduced by Schreiber [31] and its estimation using kernel density estimation techniques.
We developed and implemented a framework for the analysis of both simulated and experimentally measured Ca^{2+} time series with the informationtheoretic measure transfer entropy. This involved the stochastic coupling of a simulated enzyme to arbitrary calcium time series and the estimation of the transfer entropy of those bivariate data using kernel density estimation methods.
We investigated the information transfer from the calcium signal to the target enzyme under different conditions, namely different particle numbers (by varying the volume) and different calcium dynamics (corresponding to different stimuli). We found that, most of the time, information transfer increases with increasing particle numbers in the system. However, this increase is different in systems with different dynamic modes (spiking, bursting, etc.). More complex dynamic modes (bursting or irregular oscillations) tend to result in higher values of the transfer entropy. We observed that the input signal has to lie in the sensitive range, e.g. near the K_{ M }value of the enzyme, for the information transfer to be efficient. We also estimated the transfer entropy based on experimental data from hepatocytes. The values of these estimates are significantly lower than those from comparable simulated data. The major reason for this seems to be the unphysiologically long duration of simulated bursts. Further study is needed to investigate that in detail.
Even though the estimation of transfer entropy from time series is tricky and there are still some unsolved issues, it is a promising tool not only for the quantification of information transfer in biochemical networks, but also, for instance, to distinguish between different stochastic time series where a pure visual investigation is difficult. The direct comparison of two stochastic trajectories is difficult: Not the actual trajectory is important, but the features of it, that are essential for the correct functioning of the cell. In the case of calcium signaling, they are the ones that can be decoded by downstream elements.
Each dynamic state exhibits its own sensitivity to random fluctuations [40] and this should be reflected in the faster degradation of information transfer if the sensitivity is high. Therefore, one possible application of this approach could be the detection of the transition between stochastic and quasideterministic behavior, in cases where it is difficult to be identified by visual inspection alone. We saw one example of that already in the case of k_{2} = 2.99 (see Results), where the stochastic behavior is qualitatively different from the deterministic limit and where the transfer entropy could detect this transition. Another application could be information theorybased model fitting where models are optimized so as to have the same information transfer as observed experimentally.
It is worth mentioning that our framework is not at all limited to calcium signaling. Stochastic coupling and/or estimation of transfer entropy from biochemical data can be easily applied to other biochemical models and pathways.
Our approach can also be extended in a number of ways. On the technical side, for example, the estimation of transfer entropy from limited data sets should be improved. This could include the automatic determination of an optimal kernel bandwidth, the use of different kernels or alternative probability density estimation methods.
On the biological side, we plan to investigate the type and amount of information carried by the different properties of the calcium signal (amplitude, frequency, duration, shape, timing), because it is not yet clear which of those are really used in cells. For instance, enzymes which cannot decode a specific signal should lead to a transfer entropy value of almost zero. On the contrary, the transfer entropy is expected to show significant higher values for enzymes which are sensitive to the input signal. Thus we hope that the transfer entropy can give valuable hints for further theoretical and experimental studies. Furthermore, we want to use our framework to study different enzyme regulation mechanisms and to analyze other signaling pathways including their possible crosstalks.
Methods
Model
Model of calcium oscillations. Simple model of calcium oscillations [37]. Parameters: k_{1} = 0.212, k_{3} = 1.52, K_{4} = 0.19, k_{5} = 4.88, K_{6} = 1.18, k_{7} = 1.24, k_{8} = 32.24, K_{9} = 29.09, k_{10} = 13.58, k_{11} = 153, K_{12} = 0.16. k_{2} is bifurcation parameter and is set to different values in simulations depending on the desired behavior.
Reactions  Kinetics  

→  G_{ α }  v_{1} = k_{1}  
$\stackrel{{\text{G}}_{\alpha}}{\to}$  G_{ α }  v_{2} = k_{2} × [G_{ α }]  
G_{ α }  $\stackrel{\text{PLC}}{\to}$  ${\nu}_{3}=\frac{{k}_{3}\times [{\text{G}}_{\alpha}]\times [\text{PLC}]}{{K}_{4}+[{\text{G}}_{\alpha}]}$  
G_{ α }  $\stackrel{{\text{Ca}}^{\text{2+}}}{\to}$  ${\nu}_{4}=\frac{{k}_{5}\times [{\text{G}}_{\alpha}]\times [{\text{Ca}}^{2+}]}{{K}_{6}+[{\text{G}}_{\alpha}]}$  
$\stackrel{{\text{G}}_{\alpha}}{\to}$  PLC  v_{5} = k_{7} × [G_{ α }]  
PLC  →  ${\nu}_{6}=\frac{{k}_{8}\times [\text{PLC}]}{{K}_{9}+[\text{PLC}]}$  
$\stackrel{{\text{G}}_{\alpha}}{\to}$  Ca^{2+}  v_{7} = k_{10} × [G_{ α }]  
Ca^{2+}  →  ${\nu}_{8}=\frac{{k}_{11}\times [{\text{Ca}}^{\text{2+}}]}{{K}_{12}+[{\text{Ca}}^{\text{2+}}]}$ 
The parameter k_{2} represents the stimulation strength and serves as bifurcation parameter to vary the dynamic behavior of the model. In addition, we changed the numbers of particles present by varying the volume of the system.
Activation of the inactive form of protein P_{inact} to its active form P_{act} by Ca^{2+} is modeled by a Hill term of order four while deactivation follows mass action kinetics (Eq. (1)). The parameters were set to k_{act} = 0.025, k_{inact} = 0.005, K_{ M }= 1.0, P_{tot} = 5.0 and p = 4.
Stochastic simulation and coupling
Since the seminal work of Gillespie 1976 [38, 46] several algorithms have been proposed to calculate trajectories governed by the Chemical Master Equation [47, 48]. These trajectories are instances of the underlying stochastic process and consider the random fluctuations correctly. The different methods are mathematically equivalent and differ only in their algorithmic implementation. There are a number of software packages supporting stochastic simulation [49, 50].
This function describes the probability P(τ, μx, t) that starting in state x at timepoint t the next reaction in the system will occur after time τ and will be of type R_{ μ }. One can see that the reaction times are exponentially distributed and thus the process is a homogeneous Poisson process. By iteratively drawing pseudorandom numbers according to this reaction probability density function and updating the system state, trajectories can be calculated in a Monte Carlo scheme.
Since we not only want to analyze simulated calcium dynamics, but also intend to couple measured calcium time series to our enzyme activation model, we have to take that influence into account. The coupled calcium system exerts an influence on the reaction propensities а_{ μ }in the protein model and thus they can no longer be considered constant between two reaction events. Mathematically, this is equivalent to changing the homogeneous Poisson process into an inhomogeneous one and therefore the pure stochastic simulation methods cannot be used in this case.
One can sample this inhomogeneous Poisson process by integrating the differential reaction probabilities over time. Whenever a stopping criterion for one of the reactions is reached, the integration is interrupted and the corresponding reaction event is instantiated. This method has been used in hybrid stochastic/deterministic simulation methods [52, 53], where the set of reactions is partitioned into a stochastically simulated and a deterministically simulated subset. During the simulation, the influence of the (fast) deterministic subset on the stochastic subset has to be considered.
When we couple a time series to a stochastically simulated system, we do not know the states of the system which produced the time series between two samples. Therefore a reasonable presumption is to assume piecewise constant particle numbers between two sample times. In this special case we can use, for instance, Gillespie's Direct Method, with the modification that recalculation of all the a_{ μ }in the system, which are dependent on the coupled time series, is needed whenever a new sample is observed in the time series. This approximation is discussed in section Discussion.
We implemented this simple method in C++code dynamically linked to Octave (Version 2.9.9 on Linux) [54]. Our implementation accepts time series with arbitrary sampling times, both evenly and unevenly sampled.
Transfer Entropy
The transfer entropy has KullbackLeibler divergence form and measures the deviation of process I from its Markov process behavior of order k due to the interaction with process J. In this study, we set k = l = 1, which is justified in section Discussion. One should keep in mind though, that in the general case longer history lengths might be required. Setting those parameters correctly is crucial for a reliable estimation. If this is the case, probability densities in spaces of dimension >3 must be estimated.
For the estimation of the transfer entropy, usually either the time series is coarsegrained by histogrambased methods and the transfer entropy is estimated on the symbolic time series or kernel density estimation [44] methods are used.
with Θ the Heaviside function.
To avoid spurious effects caused by temporal correlations, we employed a Theiler window approach which excluded all neighbors that were too close in time. In addition, in order to dampen statistical fluctuations, only samples that had a userdefined minimal number of neighbors were considered. The kernel density estimation procedure was implemented using a twodimensional boxassisted neighborsearching algorithm [55], which resulted in a five to six fold speedup compared to the naive implementation.
We scanned the transfer entropy of the simulated calcium model and the coupled enzyme for different volumes (between 1 × 10^{13} and 5 × 10^{9} [arbitrary units]), corresponding to different particle numbers in the system (roughly between 600 and 30 000 000 during primary peaks), and for different dynamics, e.g. different values of the bifurcation parameter k_{2} (1 understimulation, 2 spiking, 2.5, 2.85 bursting, 2.99 irregular, 3, 3.2 overstimulation). The same kernel bandwidth were used in the calcium and the protein concentration spaces, but the data was normalized to have mean 0.0 and standard deviation 1.0 prior to the estimation.
Experiments
Single hepatocytes were isolated from fed male Wistarstrain rats (150–250 g) by collagenase perfusion as described previously [56]. The cells were harvested and incubated at 37°C at low density (10^{3} cells/ml) in 2% type IX agarose in William's medium E (WME). Single hepatocytes were prepared for microinjection with the bioluminescent Ca^{2+} indicator aequorin as described previously [57]. The injected cell was transferred to a perfusable cup held at 37°C, positioned under a cooled, lownoise photomultiplier, and continuously superfused with WME, to which agonists were added. Photon counts were sampled every 50 ms by computer. At the end of an experiment, the total aequorin content of each cell was determined by discharging the aequorin by lysing the cell. The data were normalized retrospectively by computer, by calculating the photon counts per second divided by the total counts remaining. The computed fractional rate of aequorin consumption could then be plotted as [Ca^{2+}]_{ i }using in vitro calibration data and exponential smoothing with time constants: for resting [Ca^{2+]}_{ i }, 12 s; for transients, 1 s.
In addition, we transformed these data to be roughly in the same range of values ([0, 10]) as that of the simulated data.
Materials
Aequorin was provided by Prof. O. Shimomura (Marine Biological Laboratory, Woods Hole, MA, U.S.A). Collagenase was obtained from Roche Diagnostics (Lewes, U.K.) and WME from Invitrogen (Paisley, U.K.). Agarose and agonists were purchased from SigmaAldrich (Poole, U.K.).
Declarations
Acknowledgements
JP and UK would like to thank the Klaus Tschira Foundation and the BMBF for financial support. Also we would like to thank Sven Sahle for fruitful discussions and Katrin Hübner, Markus Müller and Jocelyn Faberman for helpful comments on the manuscript.
Authors’ Affiliations
References
 Endy D, Brent R: Modelling cellular behavior. Nature 2001, 409: 391–395. 10.1038/35053181View ArticlePubMedGoogle Scholar
 Érdi P, Tóth J: Mathematical models of chemical reactions – Theory and applications of deterministic and stochastic models. Princeton, NJ: Princeton University Press; 1989.Google Scholar
 Heinrich R, Schuster S: The Regulation of Cellular Systems. New York, NY: Chapman and Hall; 1996.View ArticleGoogle Scholar
 Berridge M, Bootman M, Lipp P: Calcium – a life and death signal. Nature 1998, 395: 645–648. 10.1038/27094View ArticlePubMedGoogle Scholar
 Carafoli E: Calcium signaling: A tale for all seasons. PNAS 2002, 99(3):1115–1122. 10.1073/pnas.032427999PubMed CentralView ArticlePubMedGoogle Scholar
 Woods N, Cuthbertson K, Cobbold P: Repetitive transient rises in cytoplasmic free calcium in hormonestimulated hepatocytes. Nature 1986, 319: 600–602. 10.1038/319600a0View ArticlePubMedGoogle Scholar
 Dolmetsch R, Xu K, Lewis R: Calcium oscillations increase the efficiency and specificity of gene expression. Nature 1998, 392: 933–936. 10.1038/31960View ArticlePubMedGoogle Scholar
 Falcke M: Deterministic and stochastic models of intracellular Ca^{2+}waves. New Journal of Physics 2003, 5: 96.1–96.28. 10.1088/13672630/5/1/396View ArticleGoogle Scholar
 Keener J, Sneyd J: Mathematical Physiology. Springer 2001 chap. Calcium Dynamics;Google Scholar
 Somogyi R, Stucki J: Hormoneinduced Calcium Oscillations in Liver Cells Can Be Explained by a Simple One Pool Model. J Biol Chem 1991, 266(17):11068–11077.PubMedGoogle Scholar
 Larsen A, Olsen L, Kummer U: On the encoding and decoding of calcium signals in hepatocytes. Biophys Chem 2004, 107: 83–99. 10.1016/j.bpc.2003.08.010View ArticlePubMedGoogle Scholar
 Schuster S, Marhl M, Höfer T: Modelling of simple and complex calcium oscillations – From singlecell responses to intercellular signalling. Eur J Biochem 2002, 269: 1333–1355. 10.1046/j.00142956.2001.02720.xView ArticlePubMedGoogle Scholar
 Celio M, (Ed): Guidebook to the CalciumBinding Proteins. Oxford: Oxford University Press; 1996.Google Scholar
 Larsen A, Kummer U: Information Processing in Calcium Signal Transduction. In Understanding Calcium Dynamics. Volume 623. Edited by: Falcke M, Malchow D. Berlin: Springer Verlag; 2003:153–178.View ArticleGoogle Scholar
 De Koninck P, Schulman H: Sensitivity of CaM Kinase II to the Frequency of Ca^{2+}Oscillations. Science 1998, 279: 227–230. 10.1126/science.279.5348.227View ArticlePubMedGoogle Scholar
 Li WH, Llopis J, Whitney M, Zlokarnik G, Tsien R: Cellpermeant caged InsP_{3}ester shows that Ca^{2+}spike frequency can optimize gene expression. Nature 1998, 392: 936–941. 10.1038/31965View ArticlePubMedGoogle Scholar
 Oancea E, Meyer T: Protein Kinase C as a Molecular Machine for Decoding Calcium and Diacylglycerol Signals. Cell 1998, 95(3):307–318. 10.1016/S00928674(00)817638View ArticlePubMedGoogle Scholar
 Dupont G, Houart G, Koninck PD: Sensitivity of CaM kinase II to the frequency of Ca^{2+}oscillations: a simple model. Cell Calcium 2003, 34: 485–497. 10.1016/S01434160(03)001520View ArticlePubMedGoogle Scholar
 Gall D, Baus E, Dupont G: Activation of the Liver Glycogen Phosphorylase by Ca^{2+}Oscillations: a Theoretical Study. J Theor Biol 2000, 207: 445–454. 10.1006/jtbi.2000.2139View ArticlePubMedGoogle Scholar
 Salazar C, Politi A, Höfer T: Decoding of calcium oscillations by phosphorylation cycles. In Proceedings of Fourth International Workshop on Bioinformatics and Systems Biology Edited by: Mamitsuka H, Smith T, Holzhütter H, Kanehisa M, DeLisi C, Heinrich R, Miyano S, Kyoto. 2004, 50–51.Google Scholar
 Marhl M, Perc M, Schuster S: Selective regulation of cellular processes via protein cascades acting as bandpass filters for timelimited oscillations. FEBS Letters 2005, 579: 5461–5465.View ArticlePubMedGoogle Scholar
 Marhl M, Perc M, Schuster S: A minimal model for decoding of timelimited Ca^{2+}oscillations. Biophys Chem 2006, 120: 161–167. 10.1016/j.bpc.2005.11.005View ArticlePubMedGoogle Scholar
 Schuster S, Knoke B, Marhl M: Differential regulation of proteins by bursting calcium oscillations – a theoretical study. BioSystems 2005, 81: 49–63. 10.1016/j.biosystems.2005.02.004View ArticlePubMedGoogle Scholar
 Rozi A, Jia Y: A theoretical study of effects of cytosolic Ca^{2+}oscillations on activation of glycogen phosphorylase. Biophys Chem 2003, 106: 193–202. 10.1016/S03014622(03)001923View ArticlePubMedGoogle Scholar
 Weaver W, Shannon C: The Mathematical Theory of Communication. Urbana and Chicago, IL: University of Illinois Press; 1949.Google Scholar
 Imas O, Ropella K, Ward B, Wood J, Hudetz A: Volatile anesthetics disrupt frontalposterior recurrent information transfer at gamma frequencies in rat. Neuroscience Lett 2005, 387(3):145–50. 10.1016/j.neulet.2005.06.018View ArticleGoogle Scholar
 Borst A, Theunissen F: Information theory and neural coding. Nature Neuroscience 1999, 2(11):947–957. 10.1038/14731View ArticlePubMedGoogle Scholar
 Prank K, Schöfl C, Läer L, Wagner M, von zur Mühlen A, Brabant G, Gabbiani F: Coding of timevarying hormonal signals in intracellular calcium spike trains. Pac Symp Biocomput 1998, 633–44.Google Scholar
 Kropp M, Gabbiani F, Prank K: Differential coding of humoral stimuli by timing and amplitude of intracellular calcium spike trains. IEE ProcSyst Biol 2005, 152(4):263–268. 10.1049/ipsyb:20050040View ArticleGoogle Scholar
 Prank K, Läer L, von zur Mühlen A, Brabant G, Schöfl C: Decoding of intracellular calcium spike trains. Europhys Lett 1998, 42(2):143–147. 10.1209/epl/i1998002202View ArticleGoogle Scholar
 Schreiber T: Measuring information transfer. Phys Rev Lett 2000, 85(2):461–4. 10.1103/PhysRevLett.85.461View ArticlePubMedGoogle Scholar
 Katura T, Tanaka N, Obata A, Sato H, Maki A: Quantitative evaluation of interrelations between spontaneous lowfrequency oscillations in cerebral hemodynamics and systemic cardiovascular dynamics. Neuroimage 2006, 31(4):1592–600. 10.1016/j.neuroimage.2006.02.010View ArticlePubMedGoogle Scholar
 Marschinski R, Kantz H: Analysing the information flow between financial time series – An improved estimator for transfer entropy. Europ Phys J B 2002, 30(2):275–281. 10.1140/epjb/e2002003792View ArticleGoogle Scholar
 Materassi M, Wernik A, Yordanova E: Determining the verse of magnetic turbulent cascades in the Earth's magnetospheric cusp via transfer entropy analysis: preliminary results. Nonlin Processes Geophys 2007, 14: 153–161.View ArticleGoogle Scholar
 Nichols J, Seaver M, Trickey S, Salvino L, Pecora D: Detecting impact damage in experimental composite structures: an informationtheoretic approach. Smart Mater Struct 2006, 15(2):424–434. 10.1088/09641726/15/2/023View ArticleGoogle Scholar
 Lungarella M, Sporns O: Mapping Information Flow in Sensorimotor Networks. PLOS Comp Biol 2006, 2(10):e144. 10.1371/journal.pcbi.0020144View ArticleGoogle Scholar
 Kummer U, Olsen L, Dixon C, Green A, BornbergBauer E, Baier G: Switching from Simple to Complex Oscillations in Calcium Signaling. Biophys J 2000, 79: 1188–1195.PubMed CentralView ArticlePubMedGoogle Scholar
 Gillespie D: A General Method for Numerically Simulating the Stochastic Time Evolution of Coupled Chemical Reactions. J Comp Phys 1976, 22: 403–434. 10.1016/00219991(76)900413View ArticleGoogle Scholar
 Schreiber T, Schmitz A: Surrogate time series. Physica D 2000, 142(3–4):346–382. 10.1016/S01672789(00)000439View ArticleGoogle Scholar
 Kummer U, Krajnc B, Pahle J, Green A, Dixon C, Marhl M: Transition from Stochastic to Deterministic Behavior in Calcium Oscillations. Biophys J 2005, 89(3):1603–1611. 10.1529/biophysj.104.057216PubMed CentralView ArticlePubMedGoogle Scholar
 Rao C, Wolf D, Arkin A: Control, exploitation and tolerance of intracellular noise. Nature 2002, 420: 231–237. 10.1038/nature01258View ArticlePubMedGoogle Scholar
 Savage V, Allen A, Brown J, Gillooly J, Herman A, Woodruff W, West G: Scaling of number, size, and metabolic rate of cells with body size in mammals. PNAS 2007, 104(11):4718–4723. 10.1073/pnas.0611235104PubMed CentralView ArticlePubMedGoogle Scholar
 Steuer R, Kurths J, Daub C, Weise J, Selbig J: The mutual information: Detecting and evaluating dependencies between variables. Bioinformatics 2002, 18(Suppl 2):S231S240.View ArticlePubMedGoogle Scholar
 Silverman B: Density estimation for statistics and data analysis. London: Chapman & Hall/CRC; 1986.View ArticleGoogle Scholar
 Kaiser A, Schreiber T: Information transfer in continuous processes. Physica D 2002, 166(1–2):43–62. 10.1016/S01672789(02)004323View ArticleGoogle Scholar
 Gillespie D: Exact Stochastic Simulation of Coupled Chemical Reactions. J Phys Chem 1977, 81(25):2340–2361. 10.1021/j100540a008View ArticleGoogle Scholar
 Gibson M, Bruck J: Efficient Exact Stochastic Simulation of Chemical Systems with Many Species and Many Channels. J Phys Chem A 2000, 104(9):1876–1889. 10.1021/jp993732qView ArticleGoogle Scholar
 Cao Y, Li H, Petzold L: Efficient formulation of the stochastic simulation algorithm for chemically reacting systems. J Chem Phys 2004, 121(9):4059–4067. 10.1063/1.1778376View ArticlePubMedGoogle Scholar
 Copasi[http://www.copasi.org]
 Dizzy[http://magnet.systemsbiology.net/software/Dizzy]
 Gillespie D: Marcov Processes – An Introduction for Physical Scientists. San Diego, CA: Academic Press, Inc; 1992.Google Scholar
 Salis H, Kaznessis Y: Accurate hybrid stochastic simulation of a system of coupled chemical or biochemical reactions. J Chem Phys 2005, 122(5):054103. 10.1063/1.1835951View ArticleGoogle Scholar
 Alfonsi A, Cancès E, Turinici G, Di Ventura B, Huisinga W: Adaptive Simulation of Hybrid Stochastic and Deterministic Models for Biochemical Systems. ESAIM Proceedings 2005, 14: 1–13.Google Scholar
 Octave[http://www.octave.org]
 Kantz H, Schreiber T: Nonlinear Time Series Analysis. Cambridge, UK: Cambridge University Press; 1997.Google Scholar
 Dixon C, Cobbold P, Green A: Actions of ADP, but not ATP, on cytosolic free Ca^{2+}in single rat hepatocytes mimicked by 2methylthioATP. Br J Pharmacol 1995, 116: 1979–1984.PubMed CentralView ArticlePubMedGoogle Scholar
 Cobbold P, Lee J: Aequorin measurements of cytoplasmic free calcium. In Cellular Calcium: A Practical Approach. Edited by: McCormack J, Cobbold P. Oxford: I.R.L. Press; 1991:55–81.Google Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.