## Abstract

Ca^{2+} is a ubiquitous intracellular messenger that regulates diverse cellular activities. Extracellular stimuli often evoke sequences of intracellular Ca^{2+} spikes, and spike frequency may encode stimulus intensity. However, the timing of spikes within a cell is random because each interspike interval has a large stochastic component. In human embryonic kidney (HEK) 293 cells and rat primary hepatocytes, we found that the average interspike interval also varied between individual cells. To evaluate how individual cells reliably encoded stimuli when Ca^{2+} spikes exhibited such unpredictability, we combined Ca^{2+} imaging of single cells with mathematical analyses of the Ca^{2+} spikes evoked by receptors that stimulate formation of inositol 1,4,5-trisphosphate (IP_{3}). This analysis revealed that signal-to-noise ratios were improved by slow recovery from feedback inhibition of Ca^{2+} spiking operating at the whole-cell level and that they were robust against perturbations of the signaling pathway. Despite variability in the frequency of Ca^{2+} spikes between cells, steps in stimulus intensity caused the stochastic period of the interspike interval to change by the same factor in all cells. These fold changes reliably encoded changes in stimulus intensity, and they resulted in an exponential dependence of average interspike interval on stimulation strength. We conclude that Ca^{2+} spikes enable reliable signaling in a cell population despite randomness and cell-to-cell variability, because global feedback reduces noise, and changes in stimulus intensity are represented by fold changes in the stochastic period of the interspike interval.

## INTRODUCTION

Signaling from receptors in the plasma membrane requires a strong correlation between the extracellular stimulus and downstream intracellular events if information is not to be lost. The mechanisms by which changes in extracellular stimulus intensity are reliably converted into graded changes in cellular activity have not been fully resolved. It is unclear, for example, whether individual cells reliably transmit changes in the intensity of an extracellular stimulus to a graded change in cellular activity or whether the correlation between stimulus intensity and cellular response is realized largely by the average behavior of many cells (*1*). The complex spatiotemporal organization of the changes in cytosolic free Ca^{2+} concentration ([Ca^{2+}]_{i}) evoked by receptors that stimulate formation of inositol 1,4,5-trisphosphate (IP_{3}) enables Ca^{2+} to regulate many cellular events (*2*–*8*). IP_{3} receptors (IP_{3}Rs) are intracellular Ca^{2+} channels located on the endoplasmic reticulum (ER). Opening of IP_{3}Rs is stimulated by IP_{3} and Ca^{2+} (*9*), allowing them to propagate Ca^{2+} signals regeneratively. As IP_{3} concentrations increase, openings of single IP_{3}Rs lead to coordinated opening of clustered IP_{3}Rs, and then to cytosolic Ca^{2+} waves (*10*). Repetitive initiation of these waves generates sequences of Ca^{2+} spikes, the frequency of which often increases with stimulus intensity (*2*, *11*, *12*). This repetitive spiking behavior is not limited to Ca^{2+} signaling because there are examples of other signaling pathways in which sustained stimuli evoke pulsatile intracellular signals or responses (*13*, *14*). The location, amplitude, duration, and frequency of Ca^{2+} signals are likely to convey information (*2*–*8*, *12*). However, most supporting evidence comes from biochemical analyses (*6*), experimentally imposed Ca^{2+} spikes (*4*, *6*), or analysis of cell populations (*4*). Few examples directly demonstrate that individual cells can mount graded responses to changes in extracellular stimulus intensity by decoding Ca^{2+} spikes (*5*, *7*, *8*).

In various cells, extracellular stimuli evoke trains of Ca^{2+} spikes for which the interspike intervals (ISIs) are not predictable, because each ISI has a large random component (*15*). Thus, we investigated whether there are features of the relationship between stimulus and ISI that enabled reliable encoding of stimulus intensity in the spike frequency. We analyzed Ca^{2+} spikes by imaging single cells exposed to ligands that stimulate IP_{3}-evoked Ca^{2+} signals. We then performed mathematical analysis to identify properties of the sequences of spikes that correlated with stimulus intensity, and tested those properties for robustness against variability between cells.

## RESULTS

### Ca^{2+} spike sequences exhibit temporal randomness and large cell-to-cell variability

We performed Ca^{2+} imaging of individual cells exposed to ligands that activate phospholipase C (PLC) through G protein (heterotrimeric guanine nucleotide–binding protein)–coupled receptors (GPCRs), thereby stimulating production of IP_{3} and Ca^{2+} release from the ER (Fig. 1A). We imaged human embryonic kidney (HEK) 293 cells exposed to carbachol (CCh), an agonist of muscarinic acetylcholine receptors, or rat primary hepatocytes exposed to either vasopressin or phenylephrine, an agonist of α_{1}-adrenoceptors (Fig. 1B, upper). We measured the ISI for sequences of spikes occurring in individual cells (Fig. 1B, middle), and plotted the average ISI (*T*_{av}) and its SD (σ) for each condition (Fig. 1B, bottom) (abbreviations and symbols are listed in table S1).

Repetitive Ca^{2+} spiking in stimulated HEK293 cells, hepatocytes, and many other cells (*2*, *11*, *12*, *16*) creates a false impression that Ca^{2+} spikes are predictable (*15*, *17*). All biological processes include some variability (*18*), but we found that in stimulated HEK293 cells and hepatocytes, a distinguishing feature of the Ca^{2+} spikes was the correlation between the variability of the ISI within each cell (σ) and *T*_{av}. This relationship between σ and *T*_{av} has been described previously for Ca^{2+} spikes in HEK293 cells, astrocytes, microglia, processed lipoaspirate (PLA) cells (*15*), and endothelial progenitor cells (*17*). The correlation is captured in the *T*_{av}-σ relation (Fig. 1B, bottom):*T*_{min} is the sum of the spike duration and refractory period, and σ_{min} is the SD of ISI sequences with *T*_{av} = *T*_{min}. The ISI has three components (Fig. 1C): spike duration, refractory period, and an interval before the next spike initiates that terminates stochastically (*15*). For Ca^{2+} spikes, the relationship between *T*_{av} and σ reveals the contribution of such a stochastic process to the ISI. This stochastic period shortens as the stimulus intensity increases until it becomes so brief that a spike initiates almost as soon as a cell emerges from the refractory period (*15*). With high stimulus intensities, spike sequences would have an average ISI equal to *T*_{min}, and the ISIs would become almost uniform, with a small σ (σ_{min}). Under all other conditions, the stochastic component makes a major contribution to the ISI.

We observed that there is also variability in *T*_{av} between cells treated with the same stimulus (Fig. 1D). Under various conditions, this variability between cells often exceeded *T*_{av} for the cell population, indicating that there is no consistent relationship between stimulus intensity and *T*_{av} that applies to all cells. These data indicated that there are at least two potential impediments to transferring information reliably from extracellular signals through repetitive Ca^{2+} spikes: individual cells differ in their sensitivities to stimuli, and within cells, Ca^{2+} spike sequences have a stochastic or random component.

### The signal-to-noise ratio is increased by negative feedback and robust against cell-to-cell variability

The ISI comprises fixed components (*T*_{min}: the spike duration and the refractory period) and a stochastic period (Fig. 1C, bottom). By analyzing the temporal randomness of the ISI, we can determine how the probability of a spike occurring changes with time after the refractory period. Therefore, we examined many sequences of stimulus-evoked Ca^{2+} spikes to determine the variability of the ISI. If the probability of a spike occurring stepped immediately to its final value, the occurrence of each spike would be unaffected by the timing of any preceding Ca^{2+} spike. The slope of the *T*_{av}-σ relation (α) would then be 1 (*15*), that is, the ISIs would be described by a Poisson distribution. Once the refractory period had passed, the timing of the next Ca^{2+} spike would then behave like the random decay of radioactive atoms. Under this condition, the signal-to-noise ratio (α^{−1}), and thus the information content of the signal, would have their minimal values (*19*, *20*). There is, however, an interval after the refractory period when the cell recovers from negative feedback inhibition toward the maximal probability of firing another Ca^{2+} spike (eq. S4, fig. S1, and text S1: Mathematical description of ISI distributions) (*15*, *17*, *19*). This lingering inhibition after the refractory period delays the time at which the information content of the signal becomes minimal. The time scale of recovery from feedback inhibition is, therefore, critical in determining the reliability of signaling: slow recovery from feedback inhibition improves reliability by increasing signal-to-noise ratios (reducing α) (eqs. S3 and S6, fig. S2, and text S1: Mathematical description of ISI distributions) (*21*, *22*).

For hepatocytes and HEK293 cells, α is less than 1 (Fig. 1B), indicating the involvement of feedback inhibition with a recovery time after the refractory period that is at least as long as the ISI (*23*). Because the ISI can range from tens of seconds to several minutes, the recovery is too slow to result from known properties of clustered IP_{3}Rs: the intervals between Ca^{2+} signals evoked by clusters of IP_{3}Rs, often described as Ca^{2+} puffs, are no more than a few seconds (*23*). The feedback inhibition must, therefore, be “global,” that is, a consequence of the Ca^{2+} spike invading the entire cell, rather than the result of local Ca^{2+} signaling from clusters of IP_{3}Rs. Such global mechanisms could include regulation of IP_{3} metabolism (*12*), changes in the sensitivity of IP_{3}Rs to IP_{3} or Ca^{2+} (*9*), or store depletion. The value of α is different for HEK293 cells and hepatocytes, and different for hepatocytes stimulated with different stimuli (Fig. 1B). That different stimuli produced different values of α suggested that the slowly reversing feedback mechanism must regulate signaling events specifically associated with GPCRs.

Because a high signal-to-noise ratio is advantageous, we examined how the ratio was affected by perturbations of steps within the Ca^{2+} signaling pathway. Theoretical analyses predict that α should be unaffected by such perturbations unless they affect the global feedback mechanism (*21*, *24*), the molecular targets of which have not yet been identified (*12*). Thus, changes in PLC activity, for example, would not be expected to affect α unless PLC contributed to the global feedback. Conversely, changing the Ca^{2+} sensitivity of a sensor mediating the feedback by, for example, altering the activity of IP_{3} 3-kinase, a kinase that inactivates IP_{3} (*25*), or influencing coupling of the GPCR to G proteins and downstream signaling events would affect α. Because the theoretical predictions have not been investigated experimentally, we perturbed various steps in the GPCR-mediated Ca^{2+} signaling pathway to determine their effects on ISI. We then calculated α to determine the consequences of the perturbations on the signal-to-noise ratio.

In the HEK293 cells, CCh activates muscarinic acetylcholine receptors coupled to Gα_{q}, which stimulate PLC and thus IP_{3}-mediated Ca^{2+} release from the ER. In the same cells, a stably expressed receptor for human parathyroid hormone (PTH) couples to Gα_{s}, which stimulates adenylyl cyclase to promote accumulation of adenosine 3′,5′-monophosphate (cAMP). cAMP directly increases the sensitivity of IP_{3}Rs to IP_{3} (*26*). We analyzed the ISI of spikes triggered by CCh in the presence or absence of U73122 to inhibit PLC (*27*), PTH to increase IP_{3}R sensitivity (*26*), or cyclopiazonic acid (CPA) to inhibit the sarcoplasmic/endoplasmic reticulum Ca^{2+}–ATPase (SERCA) (*28*) that pumps Ca^{2+} into the ER (Fig. 2A). Ca^{2+} imaging analysis showed that these manipulations changed the frequency of CCh-evoked Ca^{2+} spikes: PTH increased their frequency (shortened *T*_{av}), whereas U73122 and CPA decreased Ca^{2+} spike frequency (increased *T*_{av}) (Fig. 2, B and C, and fig. S3A). Despite substantial changes in *T*_{av}, none of these perturbations changed α (Fig. 2, D and E). Varying the concentration of CCh also changed *T*_{av} without affecting α (fig. S3B). Thus, cells can produce Ca^{2+} spikes with a constant signal-to-noise ratio despite perturbations to the signaling pathway that alter the frequency of the spikes. We found that the signal-to-noise ratio is specific for each cell type (Figs. 1B and 2E) and that within hepatocytes, it varies with the stimulus (Fig. 2E). Thus, the downstream mechanisms that decode Ca^{2+} spikes receive signals in which the signal-to-noise ratio is defined, and it is robust to variations within the pathway that triggers the spikes.

### Stimulation steps are encoded by fold changes in the average stochastic period of the ISI

We found that Ca^{2+} spikes in two different cell types responding to three different extracellular stimuli exhibited a stochastic component in the ISI (Figs. 1B and 2D), and *T*_{av} for cells exposed to identical stimuli varied widely between individual cells (Fig. 1D). With this amount of variability, how might Ca^{2+} spikes encode extracellular stimulus intensities?

In sensory perception, Weber’s law proposes that the ability to detect a small change in stimulus is a constant fraction of the initial stimulus (*29*). Therefore, we evaluated whether fold changes in the stochastic period of the ISI might encode stimulus intensities in a way that was insensitive to the variability between cells. Such fold changes in ISI would more reliably encode stimulus intensities than absolute values only if the fold changes were less variable than the absolute values. As a starting point, we define the fold change (β) as the change of the stochastic period relative to its initial value:*T*_{av1} and *T*_{av2} are average ISIs with different stimulus intensities, and *T*_{min} is the same as in Eq. 1 (see Fig. 1C and table S1). Equation 2 implies a linear relationship between Δ*T*_{av} = *T*_{av1} − *T*_{av2} and *T*_{av1} − *T*_{min}. The linear relationship is also predicted by an independently derived stochastic model, which is illustrated in Fig. 3A, and includes only basic properties of IP_{3}R clusters (*21*) (see Materials and Methods). We refer to this relationship between the change of the stochastic period (Δ*T*_{av}) relative to its initial value as the “encoding relation,” because as we show later, β reliably encodes changes in stimulus intensity:

We tested whether this relation (Eq. 3) applies to Ca^{2+} spiking using a paired stimulation protocol in which we stimulated HEK293 cells with one concentration of CCh and then switched the medium to one with a higher concentration (Fig. 3B). The amplitude of the Ca^{2+} spikes remained similar when the concentration of CCh was increased: for the largest step (30 to 200 μM CCh), the peak amplitude was 154 ± 21 nM (*n* = 23 cells, 292 spikes) for the first stimulation, and 155 ± 18 nM (632 spikes) for the second. Changes in stimulus intensity are not, therefore, encoded by spike amplitudes. Instead, *T*_{av} changed with stimulus intensity in accord with Eq. 3, and the stimulation step affected β, where β is the slope of the lines (Fig. 3C and table S2). The smallest step in CCh concentration (from 50 μM to 100 μM) caused only small changes in *T*_{av} (Fig. 3C right, red line). The practicable duration of experiments limits the reliability with which we can measure such small Δ*T*_{av}, because the length of the spike sequences sets the precision with which we can determine *T*_{av}. Thus, the linearity of the relationship between Δ*T*_{av} and *T*_{av1} − *T*_{min} was less clear for this small stimulation step than for larger steps. Nevertheless, even with such small changes in stimulus intensity, β increased with the size of the concentration step (Fig. 3C right, table S2), and the value of β obtained fits well in the context of the additional concentration-response data (Figs. 3H and 4A). Analysis of hepatocytes sequentially stimulated with 0.6 μM and then 1 μM phenylephrine showed a similar linear relationship between Δ*T*_{av} and *T*_{av1} − *T*_{min} (Fig. 3D). We assessed the linearity of the relationship between Δ*T*_{av} and *T*_{av1} − *T*_{min} using Pearson’s correlation coefficient (ρ) and analysis of explained uncertainty (*u*_{ex}) (see Materials and Methods) (*30*). A value of 1 for ρ or *u*_{ex} establishes a perfect linear correlation. These analyses (Fig. 3E) confirmed that the relation between *T*_{av1} − *T*_{min} and Δ*T*_{av} is linear and that β increases with the stimulus step in agreement with our theoretical predictions (Fig. 3A).

We assessed whether stimulation steps were more reliably encoded by fold changes (β) in the average stochastic period of the ISI (*T*_{av} – *T*_{min}) or by absolute responses (*T*_{av}) by quantifying the variability of each (Fig. 3F), using the data shown in Fig. 3 (C and D). To illustrate the methods used, we describe this analysis for HEK293 cells stimulated first with 30 μM and then with 150 μM CCh as an example (Fig. 3G). We plotted *T*_{av1} and *T*_{av2} from each single cell. If absolute responses more reliably encoded the stimulation steps, all values of *T*_{av2} should be similar to the population average (that is, *T*_{av2} = *T*_{pop2}, dashed line in Fig. 3G). If fold changes more reliably encoded the stimulation steps, the data should obey Eq. 2, which we rewrite for the purpose of this analysis as *T*_{av2} = (1 − β)*T*_{av1} + β*T*_{min} (Fig. 3G, solid line). To determine whether stimulation steps were more reliably encoded by absolute responses or fold changes, we compared the root mean square distances of the data points from these lines. We divided the distance by *T*_{pop2} to obtain the coefficient of variation (CV), which can then be compared across the different experiments (a to g) in Fig. 3F (see Materials and Methods, Eqs. 9 and 10). On the basis of this analysis, we found that the relative deviation of *T*_{av2} from its population average *T*_{pop2} [CV(*T*_{pop2}), Eq. 9] was consistently larger than its relative deviation from the encoding relation [CV(β), Eq. 10] (Fig. 3F). Thus, we concluded that β represents individual cell behavior better than does the population average (*T*_{pop}), and β encodes stimulation more reliably than does the average ISI.

Because spike amplitudes and durations were unaffected by stimulus intensity, the integral ratio (IR), that is, the ratio of the area beneath the Ca^{2+} spikes occurring during the stationary phases of responses to the first stimulus relative to that for the second stimulus, is given by Eq. 11 (see Materials and Methods). To determine whether absolute values or fold changes in the integrated Ca^{2+} signals more reliably encoded differences in stimulation intensity, we compared the CV(*T*_{pop2}) with the CV of the integral ratio CV(IR) for HEK293 cells and hepatocytes stimulated with paired steps in CCh or phenylephrine concentration (Fig. 3F). This comparison showed that fold changes of the integrated Ca^{2+} signal are less variable than is *T*_{av2}. Thus, this analysis indicated that fold changes in the average stochastic period of the ISI and fold changes in the integrated Ca^{2+} signal more reliably encode stimulus changes than does the average ISI (*T*_{av}).

### Fold changes reliably encode stimulus intensity through an exponential relationship between the stimulus concentration and response

For HEK293 cells exposed to paired steps in CCh concentration, we observed that β depended only on the step size (Δ[CCh]) and not the initial CCh concentration. This is illustrated in Fig. 3H, where the relationship between Δ[CCh] and β was the same whether the first challenge was with 30 μM (red symbols) or 50 μM CCh (blue). We showed in eq. S7 (see text S2: Mathematical derivation of the concentration-response relation) that this observation and Eq. 2 result in the differential equation:

The solution to Eq. 4 is the concentration-response relation, which predicts an exponential dependence of *T*_{av} on CCh concentration:

*T*_{av,ref} is the value of the average ISI measured at a reference CCh concentration ([CCh]_{ref}), and γ describes the sensitivity of the stochastic period of *T*_{av} to CCh. Cell-to-cell variability appears in Eq. 5 in the variability of *T*_{av,ref}, which captures differences between individual cells in the response of the average stochastic period to CCh. Because Eq. 5 applies to average ISIs, it does not conflict with the randomness of individual ISIs. Inserting Eq. 5 into Eq. 2 shows that these results entail an exponential dependence of β on stimulation step Δ[CCh] = [CCh] − [CCh]_{ref}:

Equation 6 describes the measured data well: The relationship between Δ[CCh] and the experimentally determined fold changes (β) fitted to the exponential function using the fit parameter γ confirmed the reliability with which β describes cell behavior (Fig. 3H).

Analysis of the effects of different CCh concentrations on the population average (*T*_{pop}) of *T*_{av} provided additional support for our suggestion that Eq. 5 appropriately describes the concentration-response relationship. If Eq. 5 correctly describes single-cell behavior, all cells contributing to the population average T_{pop} must obey the same exponential dependence. Consequently, *T*_{pop} is not the sum of exponentials; rather, it obeys a single exponential function:

Therefore, we analyzed the dependence of *T*_{pop} derived from analysis of the Ca^{2+} spikes evoked by different concentrations of CCh in HEK293 cells (Fig. 4A) and found that the relationship was well described by the single exponential function (Eq. 7) with values for γ and *T*_{min} obtained from Fig. 3H and Fig. 3C, respectively. This fit of the experimental data to Eq. 7 is inconsistent with β and γ varying substantially between individual cells.

Equation 7 follows from the encoding relation (Eq. 3) and the independence of β from the initial stimulus intensity (Fig. 3H, and see text S2: Mathematical derivation of the concentration-response relation). Conversely, if we had started with this exponential concentration-response relationship (Eq. 7) and found that it accurately described the responses of cell populations, it follows that Eqs. 2, 5, and 6 would apply, and that β would be independent of the initial stimulus intensity. To test if Eq. 7 accurately described experimentally determined concentration-response relations, we used our measurements of *T*_{pop} from CCh-stimulated HEK293 cells (Fig. 4A) and published data for hepatocytes (*31*) and insect salivary glands (*32*). Data from Rooney *et al*. (*31*) provided measurements of the frequency of the Ca^{2+} spikes evoked by different concentrations of phenylephrine or vasopressin in hepatocytes (Fig. 4, B and C). Data from Rapp and Berridge (*32*) with blowfly salivary glands provided measurements of the frequency of the Ca^{2+}-regulated changes in transepithelial membrane potential evoked by different concentrations of 5-hydroxytryptamine (5-HT) (Fig. 4D). For the published data, we derived *T*_{pop} from the frequency of the spiking responses, and then used least-square nonlinear regression to fit the relationships between stimulus concentration and *T*_{pop} (Fig. 4, B to D), with both *T*_{min} and γ as fitting parameters (the values are provided in the legend of Fig. 4). In each case, the effects of stimulus intensity on *T*_{pop} were well described by Eq. 7: The concentration-response relation was exponential, the response to maximal stimulation was the sum of spike duration and refractory period (*T*_{min}, Fig. 1C), and stimulation controlled the average stochastic period of the ISI (*T*_{av} − *T*_{min}).

## DISCUSSION

All biological systems are variable (*18*, *33*–*35*). Heterogeneity between cells can have beneficial effects, such as contributing to adaptability (*36*), increasing the range of sensitivities to extracellular stimuli (*1*), and providing robustness (*37*). However, heterogeneity may also distort relationships between stimulus and response. We found that the frequency of Ca^{2+} spikes in individual cells varied with stimulus intensity, but that no consistent relationship applied to all cells in the population (Fig. 1D). We identified two general features of receptor-regulated Ca^{2+} spiking that enable effective encoding of stimulus intensity.

First, Ca^{2+} spikes occurred with a reliable signal-to-noise ratio. Our observation that the slope (α) of the *T*_{av}-σ relation was constant between cells exposed to the same type of stimulus formed the basis for this conclusion (Fig. 2E and fig. S3B). Furthermore, α determines the quality of signal transmission (*19*): A small α enables more information to be transmitted by spike sequences. We showed that the value of α was unaffected by stimulus intensity (fig. S3B), as predicted earlier by mathematical modeling (*21*), nor was α affected by perturbations of the signaling pathway, like inhibition of PLC, sensitization of IP_{3}Rs, or inhibition of Ca^{2+} uptake by the ER (Fig. 2), that do not affect recovery from global negative feedback. The value of α, and thus the signal-to-noise ratio, is predicted to depend on the time scale of recovery from global negative feedback (*19*, *21*), which our results suggested is likely to operate close to GPCRs, because α was different for different extracellular stimuli (Fig. 1E).

Second, changes in extracellular stimulus intensity (ligand concentration) were encoded by fold changes in the average stochastic period of the ISI (*T*_{av} − *T*_{min}), rather than by *T*_{av} itself. The fold change (β) was similar for each cell, but it varied with stimulus step size, and β correlated better with stimulus intensity than did *T*_{pop} (Fig. 3F). Many effectors, including transcription factors, integrate Ca^{2+} signals (*1*, *3*, *17*, *38*). Fold changes in the stochastic period of the ISI lead to integrated Ca^{2+} signals that also occur as approximate fold changes (IR in Fig. 3F). Thus, downstream effectors that respond to integrated Ca^{2+} signals will also benefit from the reliability of β (*3*, *4*, *38*). Our observation that stimulus intensities are encoded by fold changes in Ca^{2+} signals may provide the explanation for downstream cellular activities responding with fold changes (*3*, *4*, *38*).

Our findings resulted in a fundamental equation (Eq. 4), which can be applied to all signals that generate a fold change in response. The solution to this equation for Ca^{2+} signaling (Eqs. 5 and 7) defines an exponential dependence of *T*_{pop} on stimulus intensity, with individual cells having similar exponents (γ). Our analysis of three different cell types responding to four different ligands suggested that this relationship is generally applicable to GPCR-stimulated PLC signaling pathways (Fig. 4). Other signaling pathways also exhibit fold changes in responses after stimulation (*14*): epidermal growth factor–mediated regulation of extracellular signal–regulated kinase 2 (ERK2) has been reported to produce fold changes in the peak of nuclear ERK2 (*39*), Wnt signaling produces fold changes in the abundance of the transcriptional regulator β-catenin (*40*), and tumor necrosis factor stimulates fold changes in the nuclear abundance of the transcription factor nuclear factor κB (NF-κB) (*41*). We suggest, therefore, that measuring the fold changes according to Eq. 2 for nuclear ERK2, β-catenin, or nuclear NF-κB and solving the corresponding Eq. 4 would provide the concentration-response relations for these systems (eq. S8).

We determined that, for Ca^{2+} spiking, temporal randomness, cell-to-cell variability, and fold changes are unified by an exponential concentration-response relationship that applies to individual cells and the population average (Eqs. 5 and 7). The same features—temporal randomness (*42*), cell variability (*39*, *42*), and fold changes (*39*)—have been reported for ERK2 activity, but without a mathematical analysis leading to the concentration-response relation. These analogies between very different signaling systems suggest that we have identified a signaling concept that extends beyond Ca^{2+} spiking.

We conclude that despite variability between cells and in the ISIs within cells, there are generic features of Ca^{2+} spikes that enable repetitive spikes to encode stimulus intensities. Slow recovery from global feedback provides spike sequences with reliable signal-to-noise ratios. Changes in stimulus intensity are encoded by fold changes (β) in the average stochastic period of the ISI, which generate exponential concentration-response relationships defined by the stimulus sensitivity (γ).

## MATERIALS AND METHODS

### Materials

Cell culture media, fura-2 AM, and fluo-4 AM were from Invitrogen. Human PTH (residues 1 to 34) was from Bachem. Other chemicals, including Arg^{8}-vasopressin, were from Sigma-Aldrich.

### Single-cell imaging of [Ca^{2+}]_{i} in HEK293 cells

HEK293 cells stably transfected with human type 1 PTH receptor (*26*) were cultured in Dulbecco’s modified Eagle’s medium/Ham’s F-12 supplemented with l-glutamine (2 mM), fetal calf serum (10%), and G-418 (800 μg/ml). For imaging experiments, the cells were plated onto 22-mm round glass coverslips coated with 0.01% poly-l-lysine. Cells were loaded with fura-2 AM (2 μM) for 45 min at 20°C in Hepes-buffered saline (HBS), washed, and imaged after a further 45 min. HBS had the following composition: 135 mM NaCl, 5.9 mM KCl, 1.2 mM MgCl_{2}, 1.5 mM CaCl_{2}, 11.6 mM Hepes, 11.5 mM glucose, pH 7.3. Single-cell fluorescence measurements were performed at 20°C in HBS as previously described (*26*). Fura-2 fluorescence ratios were collected at 5-s intervals and calibrated to [Ca^{2+}]_{i} after correction for background fluorescence determined after addition of MnCl_{2} (1 mM) and ionomycin (1 μM) to quench fura-2 fluorescence.

### Single-cell imaging of [Ca^{2+}]_{i} in hepatocytes

Hepatocytes were isolated by collagenase perfusion of livers from adult male Sprague-Dawley rats and used on the same day for Ca^{2+} imaging measurements (*5*, *31*). Cells were plated on poly-l-lysine–coated coverslips for 30 min in HBS solution (HBSS). HBSS comprised 121 mM NaCl, 5 mM NaHCO_{3}, 4.7 mM KCl, 1.2 mM KH_{2}PO_{4}, 2 mM CaCl_{2}, 1.2 mM MgSO_{4}, 10 mM glucose, 100 μM sulfobromophthalein, 0.25% (w/v) fatty acid–free bovine serum albumin, 25 mM Hepes, pH 7.4, at 37°C. Cells were loaded with fura-2 AM (8 μM) for 30 min in the presence of 0.02% (v/v) Pluronic F-127 and immediately transferred to an imaging chamber at 37°C. Pairs of fura-2 fluorescence images using excitation at 340 and 380 nm were collected at 3-s intervals with a cooled charge-coupled device camera and corrected for background fluorescence by release of cytosolic fura-2 using digitonin (25 μg/ml) before calculation of fluorescence ratios for individual cells.

### Analysis of Ca^{2+} spikes

In some sequences of Ca^{2+} spikes, the stationary phase was preceded by an initial transient in the ISI (see Fig. 3B for an example). Only the stationary components of spike sequences were used for analyses of ISI. In some recordings, there was a linear trend through the entire ISI sequence (see top panels of Fig. 1B). Such a trend, which may reflect slow changes in cellular sensitivity, is not the result of a stochastic process, but it would exaggerate the randomness of ISIs (increase σ) if it were not removed. Therefore, when calculating σ, these trends were removed by subtracting a linear fit to the trend. All other analyses (for example, calculation of *T*_{av}) were performed without this correction.

After normalization using a running average, spikes were identified as an increase in *F*_{340}/*F*_{380} fluorescence of at least 20% from the basal fluorescence ratio. Because long stationary sequences are required to determine σ reliably, only cells in which at least 12 spikes were recorded were used to determine *T*_{av}-σ relations. Within the text, *T*_{av} refers to the average ISI recorded from a single cell (SD, σ), and T_{pop} is the average of the *T*_{av}s from a population of identically stimulated cells (see table S1 for definitions and abbreviations). Where only averages of ISIs were analyzed, stationary trains with fewer Ca^{2+} spikes (>7) were accepted. Table S2 summarizes the properties of the data used.

Paired Student’s *t* tests were used to compare *T*_{av} values from cells exposed to sequential stimuli. *F* tests were used to compare the slopes of *T*_{av}-σ relations. Statistical tests used Prism (version 5, GraphPad Software); *P* < 0.05 was considered significant.

### Calculation of the encoding relation

The building blocks of intracellular Ca^{2+} signals are Ca^{2+} puffs. These are brief, local increases in [Ca^{2+}]_{i} that result from the coordinated opening of a few clustered IP_{3}Rs (*43*). Ca^{2+} spikes occur when the probability of initiating a Ca^{2+} puff increases, allowing regenerative propagation of a Ca^{2+} wave across the cell. In previous publications, we developed a stochastic model for IP_{3}-evoked Ca^{2+} signals in terms of the distributions of interpuff intervals and puff durations (*20*, *21*), which we use here to calculate the *T*_{av} values entering the encoding relations in Fig. 3A in the same way as described in (*20*, *21*). The distributions of both the interpuff intervals and puff durations are generalized exponential functions (see eq. S1 with *T*_{min} = 0 for this type of distribution). In the model, each cell has identical clusters of IP_{3}Rs, the parameters that determine the distribution of puff durations are fixed, and the values of the parameters of the interpuff interval distributions specify individual cells. The dependence of the interpuff interval distribution parameters on IP_{3} concentration was calculated from the DeYoung-Keizer model (*21*). The parameter ν (see eq. S1) consists of a constant and an additive part that depends on IP_{3} concentration. Both have values typically in the range from 1.0 to ~4.0. We captured cell-to-cell variability by varying the constant part in a way that provided ν values from 1 to 3.85 at the first stimulation level. The range of values of *T*_{av1} in Fig. 3A arises from this cell-to-cell variability. Stimulation steps were modeled by increasing the IP_{3} concentration from 0.5 μM by the steps specified in Fig. 3A. Average ISIs were calculated as the average first passage time for the transition from 0 to 4 open clusters. All other parameter values are given in (*21*).

### Statistical analysis of the encoding relation

The encoding relation (Eq. 3) predicts a linear relationship between individual *T*_{av1} and *T*_{av2} values. The linear fit was assessed in two ways. Pearson’s correlation coefficient (ρ) measures how close data points are to a straight line, whereas the explained uncertainty (*u*_{ex}) evaluates how much of the variance in the experimental data can be explained by fitting to a model like Eq. 3. In a regression model, it is generally true that *u*_{ex} is defined as the fraction of the total variance that remains if the regression error is subtracted (*31*):

The explained uncertainty thereby reports the extent to which a model adequately describes the experimental variation. For ρ and *u*_{ex}, a value of 1 indicates a perfect linear correlation between values. We used both ρ and *u*_{ex} to assess the linearity of the *T*_{av1}-Δ*T*_{av} relations for HEK293 cells and hepatocytes exposed to steps in stimulus concentration using our paired stimulus paradigm.

Our analysis of whether absolute responses (*T*_{av}) or fold changes (β) provided the most reliable encoding of stimuli was performed for the fold change β = Δ*T*_{av}/(*T*_{av1} − *T*_{min}) (Eq. 3). The analysis is equally valid for the fold change of stochastic periods (*T*_{av2} − *T*_{min})/(*T*_{av1} − *T*_{min}), which is equal to 1 − β (see Eq. 2). We compared the coefficients of variation CV(*T*_{pop2}) and CV(β):

*j*) responding to a second stimulus. *T*_{pop2} is the average *T*_{av2} from a population of *N*_{cell} cells.

*j*) responding to a first and second stimulus, respectively. *T*_{min} is defined in Eq. 1. β is determined by Δ*T*_{av}/(*T*_{av1} − *T*_{min}) (see Eqs. 2 and 3). The expression in the numerator of Eq. 9 is the root mean square orthogonal distance of the data points from the line *T*_{av2} = *T*_{pop2}, and the numerator in Eq. 10 is the root mean square orthogonal distance of the data points from the line *T*_{av2} = (1 − β)*T*_{av1} + β*T*_{min}.

In general, the integrated Ca^{2+} signal comprises the integral over the initial transient and over the stationary part of the spike sequence. Because transients vary between cells and with stimulation conditions, no general statement about their integral is possible. For paired stimuli, the ratio of the integrated stationary signal approximately equals the ratio of Ca^{2+} spike frequencies, because spike amplitudes were unaffected by stimulus intensity, and we did not detect any obvious change in spike duration. The IR of the stationary parts of the spike trains is therefore

Mathematical descriptions of ISI distributions are provided in Supplementary Materials.

## SUPPLEMENTARY MATERIALS

www.sciencesignaling.org/cgi/content/full/7/331/ra59/DC1

Text S1. Mathematical description of ISI distributions.

Text S2. Mathematical derivation of the concentration-response relation.

Fig. S1. ISI distributions for hepatocytes.

Fig. S2. *T*_{av}-σ relation.

Fig. S3. Responses of individual HEK293 cells to stimulation.

Table S1. Symbols and abbreviations.

Table S2. Selection of data for paired stimulation experiments.

## REFERENCES AND NOTES

**Funding:**Supported by grants from the Wellcome Trust (101844) and Biotechnology and Biological Sciences Research Council (BB/H009736) to C.W.T., DFG (GRK 1772) to G.M., EMBO ASTF (Application Short-Term Fellowship) 398-09 to K.T., and Infusino Endowment to A.P.T.

**Author contributions:**M.F., C.W.T., and A.S. designed the study. K.T., S.C.T., V.L.P., and A.M. collected experimental data and analyzed the data with C.W.T. and A.P.T. Mathematical analyses were performed by K.T., G.M., and A.S. and were supervised by M.F. The paper was written by C.W.T. and M.F. with input from all other authors.

**Competing interests:**The authors declare that they have no competing interests.