Calcium responses have been observed as spikes of the whole-cell calcium concentration in numerous cell types and are essential for translating extracellular stimuli into cellular responses. While there are several suggestions for how this encoding is achieved, we still lack a comprehensive theory. To achieve this goal it is necessary to reliably predict the temporal evolution of calcium spike sequences for a given stimulus. Here, we propose a modelling framework that allows us to quantitatively describe the timing of calcium spikes. Using a Bayesian approach, we show that Gaussian processes model calcium spike rates with high fidelity and perform better than standard tools such as peri-stimulus time histograms and kernel smoothing. We employ our modelling concept to analyse calcium spike sequences from dynamically-stimulated HEK293T cells. Under these conditions, different cells often experience diverse stimulus time courses, which is a situation likely to occur in vivo. This single cell variability and the concomitant small number of calcium spikes per cell pose a significant modelling challenge, but we demonstrate that Gaussian processes can successfully describe calcium spike rates in these circumstances. Our results therefore pave the way towards a statistical description of heterogeneous calcium oscillations in a dynamic environment.
Transient changes in the intracellular calcium (Ca2+) concentration have long been associated with the activation of plasma membrane receptors [1]. Since the seminal work by Woods et al. [2] linking the frequency of cytosolic Ca2+ oscillations in hepatocytes to the concentration of various hormones, both experimental and theoretical studies have provided compelling evidence for encoding extracellular stimuli into intracellular Ca2+ oscillations [3–14]. In whole-cell recordings, Ca2+ oscillations are usually observed as sequences of spikes of the intracellular Ca2+ concentration. A prominent feature of Ca2+ spike sequences is that they are random. Ca2+ spikes only occur with some probability that generally changes over time. For example, there are distributions of inter-spike intervals (ISIs) for agonist induced Ca2+ oscillations in HEK293 cells and spontaneous Ca2+ oscillations in astrocytes, microglia and PLA cells instead of a single value [15]. When astrocytes are transiently stimulated with ATP three times, with a recovery period between stimuli, the observed Ca2+ spikes display any number of response patterns, from no spikes to three [16]. To elucidate the principles that govern the translation of extracellular cues into changes of the intracellular Ca2+ concentration therefore requires faithfully capturing the stochasticity in Ca2+ spike generation. To date, the main modelling approach to investigate stochastic Ca2+ spike generation has been based on numerical solutions of differential equations, both ordinary and partial [11,12,15,17–30]. In these studies, the randomness of Ca2+ spikes results from the stochastic behaviour of Ca2+ releasing channels, such as the inositol-1,4, 5-trisphosphate (InsP3) receptor (InsP3R), and their interactions. The random dynamics of the InsP3R is then either described by coupling a Markov chain for the InsP3R to the differential equations, or by assuming a Langevin-type equation. All of these approaches require detailed models of the InsP3R (with often a considerable number of rate constants), and other assumptions such as the number of InsP3Rs per cluster and the spatial distribution of InsP3Rs, Ca2+ pumps and Ca2+ buffers in the case of partial differential equations. However, such mechanistic detail, which has been instrumental in advancing our understanding of Ca2+ spikes, often comes at considerable computational costs. It might be therefore desirable to change the perspective from the mechanistic bottom-up approach to a top-down view, in which cellular Ca2+ spikes are described directly. The mathematical concept that has proven particularly useful for this endeavour is the theory of point processes [31]. Indeed, in [15] a time-dependent conditional rate for the generation of a Ca2+ spike was introduced, and its two parameters (a time scale and an amplitude) were determined from experiments on four different cell types. Subsequent work [32,33] demonstrated in more detail that for constant stimulation, the time scale was cell type specific, while different cells of the same type could be distinguished by their amplitude. These modelling approaches have not confronted the issue of the dynamic nature of cell stimulation that occurs under physiological conditions. Cells in vivo experience a complex and dynamically-changing environment, where signals frequently arrive in a time-varying manner—such as transient release of neurotransmitters or oscillations in the concentration of circulating hormones. Diffusion of messengers through tissue (such as away from a blood vessel) also introduce spatial variation in signal strength, meaning that cell populations encounter a complex spatiotemporal pattern of stimulation. An efficient mechanism for modelling Ca2+ responses in such a heterogeneous population has not yet been devised. Our approach offers a solution to this issue and is based on combining point processes with Bayesian inference. Bayesian concepts have been used with great success across various disciplines (see e. g. [34] for an overview). One advantage of utilising Bayesian ideas is that model parameters can be effectively constrained by observed data and can be estimated in a controlled fashion. For example, each set of parameter values comes with its own probability that informs us about how likely this set represents the observed data. In the past, the combination of Bayesian inference and point processes has been successfully applied to action potential spike trains in neurons [35,36], but to our knowledge, this is the first time that Ca2+ spike sequences have been analysed in this way. While we can draw on these previous results, the substantial differences between action potential spike trains and Ca2+ spike sequences (e. g. the number of spikes per train or the time scales of spikes) have required significant attention. A particular characteristic of Ca2+ spikes is that their generation depends on the cellular Ca2+ spike history and hence is non-Markovian. This results from both the stochastic nature of Ca2+ spike formation as well as the dynamic variation in the cellular signalling micro-environment. To model such history dependence, we follow a concept introduced in [37], which effectively turns a non-Markovian Ca2+ spike sequence into a Markovian one. This is achieved with the help of a so-called intensity function x (t), for which we provide a definition and more details in the Materials and Methods section. Since the intensity function is directly inferred from individual Ca2+ spike sequences, it is specific to each cell. This results from two facts. Firstly, Ca2+ spikes are shaped by the cellular composition of the Ca2+ signalling toolkit, i. e. the expression levels and spatial arrangement of Ca2+ channels, pumps, transporters and buffers. Secondly, each cell experiences a different signalling micro-environment as illustrated above. Autocrine and paracrine signalling modify the original signal further. In addition, any time-dependence of x (t) that originates from a time varying signal is compounded by dynamic changes of the Ca2+ signalling apparatus, such as Ca2+ store refilling, adaptation or desensitisation. In the present paper we demonstrate how to estimate x (t) in the presence of such single cell variability. This will allow us to make progress in two main directions. Firstly, one of the most highly discussed questions in computational biology is concerned with scaling dynamics from single cells to tissue. Taking into account our current understanding of intracellular Ca2+ signalling, this will require the simulation of spatially extended single cells driven by fluctuating Ca2+ releasing channels. The computational costs for such studies are extraordinary, and how to sample the associated high dimensional parameter space remains an open challenge. On the other hand, generating Ca2+ spike sequences from an intensity function is computationally cheap and hence puts researchers in an advantageous position to obtain high quality statistical insights into tissue level dynamics. Secondly, while our approach is statistical, the Ca2+ ISI distribution that we determine as part of estimating x (t) has direct mechanistic interpretations. This may further our understanding of how exactly Ca2+ spikes emerge at the single cell level from the orchestrated action of the Ca2+ signalling toolkit. For instance, one line of argument suggests that Ca2+ spikes result from the co-ordinated interplay of a certain number of Ca2+ puffs. Our statistical analysis may provide quantitative estimates for this.
A main motivation for our study is the need to model the heterogeneity of Ca2+ responses observed in cell populations in a computationally efficient manner. To illustrate the problem, the results in Fig 1 show HEK293T cells challenged with a solution containing 100 μM carbachol. Fig 1B illustrates the large variability of Ca2+ responses observed between different cells even when the stimulus is stationary. While some cells only detect the onset of the stimulus (third trace), other cells exhibit Ca2+ spikes during the entire stimulation period. However, the spike characteristics vary significantly. Some cells show irregular Ca2+ spikes (first trace), while other cells settle into an almost regular pattern (fourth trace). It is also worth noting that the frequency of Ca2+ spikes spans a considerable range (cf. first and second trace) and that all cells display a decrease in peak amplitude, which is indicative of adaptation. Fig 1C and S1 Video provide further evidence for the large heterogeneity in the timing of Ca2+ spikes for a constant stimulus. The first step in developing a top-down model that reproduces this heterogeneity is to determine the probability distribution that most accurately describes recorded ISIs. Throughout this study, Ca2+ spikes are treated as all-or-nothing events and any information on the width or amplitude of a Ca2+ spike is excluded. The reason why we can introduce an ISI distribution and hence treat successive ISIs as independent—instead of trying to fit a time-dependent ISI distribution—comes from the introduction of the intensity function x (t). We tested three possible candidate ISI statistics (see Materials and methods for details): an inhomogeneous Poisson (IP), an inhomogeneous Gamma (IG) and an inhomogeneous inverse Gaussian (IIG) distribution. The IP distribution often serves as starting point for analysing spiking behaviour as it is the most basic statistical distribution. While the parsimony of the IP distribution has undoubtedly helped in establishing a large body of mathematical results, real world data often exhibit more complex statistics. The IG distribution provides a natural extension of the IP distribution, in that it contains the IP distribution as a special case: putting γ = 1 in (5) recovers (10). The shape parameter γ endows the IG distribution with more flexibility, which has proven fruitful in numerous applications. Conceptually, spikes in general and Ca2+ spikes in particular have been described as first passage events [38]. One of the most fundamental models, which contains positive drift and random motion only, gives rise to the IIG distribution. To ascertain which ISI distribution describes the results in Fig 1 best, we transformed the recorded ISIs using the time rescaling theorem and analysed the results in a Kolmogorov-Smirnov plot (see Material and methods). The Kolmogorov-Smirnov plot allows the visual inspection of whether two probability distributions are the same by plotting their respective cumulative distribution functions against one another. If the two distributions are identical, the cumulative distributions functions coincide and plotting one against the other results in a straight line with slope 1. Fig 2A shows results for the Kolmogorov-Smirnov plot for the data taken from Fig 1. Here, u corresponds to the cumulative distribution of the transformed data, which depends on the chosen ISI distribution, while s stems from a theoretical prediction based on the time rescaling theorem (see Materials and methods). If we have identified the correct ISI distribution, data points should cluster around a straight line with a 45° slope. Note that in our analysis, each cell was treated individually and no data amalgamation took place. Results for the IP and IIG clearly deviate from a line with slope 1, while the data for the IG exhibit much less deviation. This suggests that an IG, but not an IP or IIG, better describes the ISI statistics of Ca2+ spikes. The box plots in Fig 2C provide further quantitative evidence. They demonstrate that the slopes obtained from individual cells are concentrated close to 1 and exhibit a significantly smaller variability for the IG compared to the IP and IIG. To further corroborate these findings, we analysed data from cells exposed to a lower concentration of carbachol (10μM). Fig 2B and 2D illustrate that again the IG distribution captures the ISI statistics most closely. In addition to Kolmogorov-Smirnov plots, we also tested our data with a quantile-quantile plot (see S1 Fig and Materials and methods). As with the Kolmogorov-Smirnov plot the correct ISI distribution leads to data points that accumulate around the 45° line. Panels A and B in S1 Fig show that this is the case for the IG, but not for the IP and IIG, hence confirming our results from the Kolmogorov-Smirnov plot. The box plots in panels C and D of S1 Fig show that again the slopes obtained from single cells for the IG model are much closer to one and exhibit less variability than those for the IP and IIG. Taken together, these results indicate that the patterns of Ca2+ spikes observed across the population of cells can most accurately be reproduced with an IG distribution. This approach will therefore be used as the basis of our analysis. It is worth noting that when testing the different ISI distributions, we also obtain an estimate for the intensity function x (t). As illustrated by Eq (11), the probability for a specific ISI depends on x (t). Therefore, all ISI distributions in this study have to be understood as being conditioned on x (t). In the next section, we show how to estimate x (t) from measured Ca2+ spike sequences. In addition to the ISI distribution, we also need to know the intensity function x (t) to fully describe Ca2+ spike sequences. In the following, we will use three different approaches to estimate x (t): peri-stimulus time histograms (PSTHs), kernel smoothing (KS), and Gaussian processes (GPs) combined with Bayesian inference. An illustration of a GP is shown in Fig 3. In contrast to a deterministic curve, a GP generates infinitely many curves (3 possible candidates are shown), and the statistics of these curves is the organising principle. At each time point, the values of a GP are Gaussian distributed, and the mean and standard deviation can change over time. To guarantee a controlled comparison between PSTHs, KS and GPs, we will fix an intensity function, generate surrogate Ca2+ spike sequences from it and then estimate how well the above methods recapture the original intensity function (see Materials and methods for details). Before proceeding it is worth noting that the intensity function that we need to estimate is identical to the Ca2+ spike rate that we find from using either PSTHs or KS [36]. In other words, if we can obtain a high quality estimate for the Ca2+ spike rate from either PSTHs or KS, we have a very good estimate for x (t). However, this usually requires a large number of Ca2+ spike sequences, which is an issue that we will address below. Since we can identify the Ca2+ spike rate with the intensity function, we will use both terms interchangeably. Note, however, that this Ca2+ spike rate is different from the conditional intensity function defined in Eq (15), which is often used in generating spike trains. A major objective of this model is to reproduce Ca2+ signalling patterns during complex stimulation conditions. For example, physiological patterns of hormone or neurotransmitter release are time-varying, rather than the step-changes used in typical experiments (such as Fig 1). We therefore tested the performance of candidates for x (t) for reproducing dynamically-changing signals—specifically sinusoidal oscillations. As a first choice, we considered a regularly oscillating intensity function xdet (t) = 0. 5 cos (t) + 0. 5 cos (0. 5t) + 1. Fig 4A shows Ca2+ spike sequences generated from xdet, and Fig 4C reveals that both PSTH and KS capture the original intensity function very well. By using a specific functional form for x (t) as in xdet (t), we make strong assumptions about the intensity function. A more flexible and versatile approach is based on GPs. In Fig 4B we plot Ca2+ spike sequences generated from one xGP candidate, while Fig 4D depicts the estimation of xGP from a PSTH and KS. As with xdet we find very good agreement between the original and estimated Ca2+ spike rate. This strategy supposes that all Ca2+ spike sequences can be combined into a single large dataset, but in a physiological context this will not always be true. To be a more useful tool, the model should be able to simulate the diversity of responses expected in a more complex environment, where the stimulus varies in both space and time. Under these circumstances, the number of cells that receive an equivalent stimulus would be more limited, and so it is important to assess how the number of spike sequences available for parameter estimation affects the accuracy of predicting the Ca2+ spike rate. Consequently, we randomly picked groups of 1,2, 4 and 7 Ca2+ spike sequences and computed the Ca2+ spike rate based on a GP and KS. In Fig 5 we show results for xGP. Since the Ca2+ spike rate estimation obtained from KS depends on the bandwidths σ (see Eq (17) ), we employed different σ values. For a single Ca2+ spike sequence (Fig 5A) the estimated Ca2+ spike rates differ visibly from the theoretical one, and the smallest bandwidth leads to spurious oscillations. As we increase the number of Ca2+ spike sequences the estimated Ca2+ spike rates capture the true Ca2+ spike rate more faithfully. The light blue area in each panel delineates the 95% confidence interval, which we obtain as a by-product from the GP optimisation. Overall, all estimates lie within this confidence interval except the one for σ ^ in Fig 5A (see Materials and methods for the definition of σ ^). We obtain similar results for xdet (t) as illustrated with S2 Fig. We quantified the accuracy of predicting the Ca2+ spike rate by computing the normalised L2 norm of the difference between the known and estimated Ca2+ spike rate (see Materials and methods). Fig 6 shows that for a given method, the L2 norm decreases as we increase the number of Ca2+ spike sequences, which corresponds to better predictions. When we fix the number of Ca2+ spike sequences, GPs yield a better estimate. The improvement is particularly evident when comparing the Ca2+ spike rate estimation based on σ ^ at small numbers of Ca2+ spike sequences. We further tested that our results did not depend on the particular choice of Ca2+ spike sequences nor on the details of the surrogate generator. For the latter, we compared three different approaches: inverse sampling, a Bernoulli process and time rescaling. We generated a number of Ca2+ spike sequences with each method and then estimated the Ca2+ spike rate using the same methods as in Fig 6, i. e. KS with different bandwidths and GPs. S3 Fig shows box plots of the normalised L2 norm between the estimated and the true Ca2+ spike rate xdet. We find that for all three methods, the normalised L2 norm is generally smallest for GP estimates. To test the statistical significance of this result, we computed the corresponding p-values as shown in S1 Table using the non-parametric Mann-Whitney test. Based on the common assumption that a finding is statistically significant if p < 0. 05, estimates using GPs perform statistically better than KS since the largest p value was 0. 0375. We repeated the analysis for xGP and report the normalised L2 norm in S4 Fig. While the GP performs clearly better than KS with bandwidths of σ ^ and 35, the distributional results for bandwidths of 52 and 70 look similar to those of the GP. This is also confirmed by the p-values shown in S1 Table, where some exceed the threshold of 0. 05. This indicates that KS can approach the performance of GPs. However, since there are no a priori estimates for this optimal bandwidths for a given scenario, GPs provide the more robust estimation method. The results so far provide strong evidence that an intensity function derived from a GP allows accurate prediction of Ca2+ spike patterns even when the estimate is based on small numbers of Ca2+ spike sequences. We next applied our Bayesian approach to a more complex experimental system, designed to reproduce some of the stimulus heterogeneity expected in vivo. We used a microfluidics chamber to deliver sinusoidal changes in carbachol concentration to HEK293T cells (see Materials and methods). The concentration of carbachol varied in both space and time. Fig 7A shows a complex concentration surface throughout the chamber at a fixed time point and illustrates how the sharp interface between high and low agonist concentration on the left side of the chamber widens as the flow progresses through the chamber. In Fig 7B we plot agonist concentration time courses sampled at four positions along the transverse direction of the microfluidics chamber, which demonstrates the stimulation heterogeneity that cells experience depending on their position within the chamber. We also include the corresponding Ca2+ spike sequences, which again display significant variability. The goal of this experiment was to generate an environment in which a population of cells is exposed to agonist in a manner that varies with the cells’ distance from the stimulus source and with a dynamic mechanism of delivery (by analogy to a circulating hormone diffusing from a blood vessel to underlying cells, for example). This scenario is the context in which a Bayesian framework is most useful, as it can predict spiking Ca2+ responses to a complex but physiologically meaningful stimulus profile. In Fig 7C–7E, we show clustered stimulus curves, the corresponding Ca2+ spike sequences and the estimated intensity functions, respectively. We grouped cells that experienced similar agonist concentration profiles to allow for a meaningful comparison of the resultant intensity functions. To determine how similar stimulus traces are, we computed the weights of the three leading principal components (see Materials and methods). Results are shown in Fig 8, where we employed k-means [39] to detect possible clusters. Data points that belong to the same group are plotted in the same colour, and these colours correspond to those used in Fig 7C–7E. Overall, 4 distinctive groups represented the data best, containing 10,13,13 and 3 cells, respectively. It is worth noting that while clustering was performed on the stimuli time courses the Ca2+ spike sequences show a consistent pattern in that they are generally more similar within a given group than between groups. The black lines in Fig 7C and 7E denote the mean stimulus and mean intensity function, respectively. We observed that the intensity functions broadly mirror the global behaviour of the stimulus. The mean intensity function for all responses showed a single peak with a similar time course to the stimulus. The amplitude of stimulus and intensity function are also well matched. However, there are also observable differences. For example, individual cells showed intensity functions with a more complex time course than the mean (such as multiple peaks). Furthermore, while the mean stimulus was symmetrical, the mean intensity functions could exhibit notable asymmetries. For example, for weak stimuli, the rising phase of the intensity function may be markedly slower than the falling phase (Fig 7C and 7E; yellow and green traces). This most likely reflects the excitable character of intracellular Ca2+ signalling [30,40]. For weaker stimulation, it takes longer to reach the threshold for generating a Ca2+ spike, hence the intensity function grows more slowly. The quicker decrease results from the stimulus dropping below the Ca2+ spike generating threshold quickly after reaching its maximum, hence prohibiting further Ca2+ spikes. The faster increase of the intensity function for stronger stimuli (grey traces) lends further support for this interpretation, as the Ca2+ spike generating threshold is reached more quickly. To illustrate the value of the mean intensity functions shown in Fig 7E, we generated surrogate Ca2+ spike sequences from them using the IG distribution and plotted them in Fig 9C. For ease of comparison, we also show the measured Ca2+ spike sequences in Fig 9B. We first note that the simulated Ca2+ spike sequences resemble the measured ones. This is also confirmed by the histrograms in Fig 9D, which exhibit large overlaps between the experimental and theoretical Ca2+ spike sequences. To quantify how similar the two histograms are, we computed the histogram distance given by H (R, S) = ∑ i min (R i, S i) max (∑ i R i, ∑ i S i), (1) where Ri and Si denote the histogram count in the ith bin of the recorded and simulated data, respectively. The histogram distance is bounded between 0 and 1, and the closer it is to 1 the more similar the histograms are. The histograms coincide if H = 1. We found that H = 0. 86,0. 94,0. 91 and 0. 89 (from top to bottom), which confirms our visual inspection that the histograms vary little between recorded and simulated Ca2+ spike sequences. Moreover, Ca2+ spike sequences generated from one intensity function exhibit a certain degree of heterogeneity, which is consistent with our experimental findings. Taken together, these results show that without explicitly including any information about the stimulus into the estimation of the intensity functions, our approach yields intensity functions that reflect the characteristics of the stimulus and that are consistent with experimentally recorded Ca2+ spike sequences.
In this work we have developed a mathematical framework to quantitatively describe the heterogeneous timing of Ca2+ spikes in a cell population subject to time-varying stimulation. At the heart of this new approach is the use of Bayesian inference to determine the most likely intensity function and hence the most likely Ca2+ spike rate for a given stimulus. As part of this estimation process, we found that the statistics of Ca2+ ISIs are best captured by an IG distribution. Importantly, knowledge of the intensity function and the ISI statistics suffices to completely describe Ca2+ spiking. Since generating Ca2+ spike sequences from an ISI distribution and intensity function is computationally significantly cheaper than solving partial differential equations for cellular Ca2+ transport, this approach is ideally suited for numerically studying large numbers of cells. The estimation of inhomogeneous single cell behaviour also puts us in an ideal position to ascertain whether or not there is signal processing at the cell population level. Indeed, numerous examples exist where the average population behaviour is not shared by any cell (see e. g. [55]). These incongruous dynamics also warrant investigations into population invariances, where cell populations respond consistently in the same manner, albeit with completely heterogeneous single cell behaviour [56,57]. By reliably estimating single cell Ca2+ dynamics, the present study provides a stepping stone towards answering these questions for intracellular Ca2+ signalling.
We here follow the exposition in [37] for the definition of the intensity function. Assume that Ca2+ spikes occur at times y1 < y2 < … < yN. Let p (v) denote the probability density for a general renewal process on v ∈ (0, ∞), i. e. p (v) dv is the probability for an event in [v, v + dv], and subsequent events are independent. For ya > 0, let y correspond to a time variable on (ya, ∞) and X be a one-to-one mapping X (y) = v of (ya, ∞) to (0, ∞). Conservation of probability then entails that p (y) = | d v d y | p (v) = | X ′ (y) | p (X (y) ). (2) In other words, the probability density for a Ca2+ spike at yi can be computed from the renewal probability density p if we know the mapping X. A convenient form of X is X (y) = ∫ y a y x (u) d u, (3) which satisfies the conditions above and where x is called the intensity function, which is the object that we need to estimate. Eq (3) can be interpreted as rescaling the original time y such that Ca2+ ISIs become independent and identically distributed in the new time [58]. Given two subsequent Ca2+ spikes times yi−1 and yi in the original time, the ISI in the new time is X (y i - 1, y i) = ∫ y i - 1 y i x (u) d u. (4) Since it is only through the introduction of the intensity function x (t) that Ca2+ ISIs become Markov, we introduce the notation p (yi, yi−1|x), which corresponds to the ISI probability density given x (t). Note that formally the conditional ISI probability density is defined as the joint conditional probability density for spikes at yi and yi−1 (and hence no spike in [yi−1, yi]) given an intensity function x (t). We will employ three different choices for the ISI probability density: an inhomogeneous Gamma distribution p (y i, y i − 1 | x) = γ x (y i) Γ (γ) [ γ X (y i − 1, y i) ] γ − 1 e − γ X (y i − 1, y i), (5) where γ > 0 denotes the shape parameter and Γ is the Gamma function; an inhomogeneous inverse Gaussian distribution p (y i, y i − 1 | x) = x (y i) 2 π X 3 (y i − 1, y i) exp { − (X (y i − 1, y i) − α) 2 2 α 2 X (y i − 1, y i) }, (6) where α > 0 is the location parameter; and an IP distribution p (y i, y i - 1 | x) = x (y i) e - X (y i - 1, y i). (7) The time-dependent intensity function x (t) is modelled as a Gaussian Process (GP) [34,59]. A GP is uniquely defined by its mean μ (t) and covariance function Σ (t1, t2). While there are many possible choices for Σ [34,60], we employ the widely used squared exponential (SE) kernel Σ (t 1, t 2) = σ f 2 e - κ (t 1 - t 2) 2 2 + δ (t 1 - t 2) σ v 2, (8) where κ measures the smoothness of the GP and σf controls its variance. The last term allows us to model additional noise sources. We originally included σ v 2 as a hyperparameter in the optimisation. However, we consistently found small values for σ v 2 and hence decided to fix it at a presentative value of σ v 2 = 10 - 4. We collect the spike times in a sequence of N Ca2+ spikes in a vector y = {y1, …, yN}. For consistency, we set y0 = 0. Through the introduction of an intensity function x (t), the joint probability density for a spike sequence y given x (t) factorises and reads as [36] p (y | x) = p 1 (y 1 | x) p T (T, y N | x) ∏ i = 2 N p (y i, y i - 1 | x). (9) Here, p1 (y1|x) represents the conditional probability density of finding the first spike at time y1. We also take into account that the observation time T usually exceeds the last spike time through the term pT (T, yN|x), which denotes the conditional probability that no spike occurs after yN. The statistics for p1 and pT are often based on an inhomogeneous Poisson (IP) process, i. e. p 1 (y 1 | x) = x (y 1) e - X (0, y 1), p T (T, y N | x) = e - X (y N, T), (10) where X is given by Eq (3). For practical purposes, we discretise time with a time step Δ such that T = nΔ [36]. When working with experimental spike trains, we set Δ equal to the inverse of the recording frame rate. A spiking time yi can then be expressed as yi = li Δ for an appropriate l i ∈ N. By setting xi = x (iΔ) and using Eq (9) with e. g Eq (5), we obtain the probability density for a spike sequence for the inhomogeneous Gamma distribution as p (y | x, θ) = x l 1 e − X ^ 0,1 e − X ^ N, n ∏ i = 2 N γ x l i Γ (γ) [ γ X ^ i − 1, i ] γ − 1 e − γ X ^ i − 1, i, (11) where X ^ i, j = Δ ∑ k = l i l j x k and l0 = 0, ln = n. By introducing θ on the left hand side, we make explicit the dependence of the probability density on the hyperparameters θ, which in this case are θ = {γ, κ, σf}. The most probable intensity function x* (t) given a spike train y is determined by x* = argmaxx≥0 p (x|y). Under the assumption that the nodal value x* is close to its mean, we have x * ≈ ∫ x θ * p (θ | y) d θ = 1 Z ∫ x θ * F (y, x θ *, θ) d θ, (12) where x θ * = argmax x ≥ 0 p (x | y, θ) = argmax x ≥ 0 p (y | x, θ) p (x | θ). (13) To evaluate the first integral in Eq (12) we note that p (θ | y) = p (θ) p (y) ∫ p (y | x, θ) p (x | θ) d x = F (y, x θ *, θ) p (y), (14) with F (y, x θ *, θ) = p (θ) p (y | x θ *, θ) p (x θ * | θ) / | Λ * + Σ - 1 | and Λ * = - L x 2 log p (y | x θ *, θ), where we used Laplace’s approximation for the integral as shown in S1 Appendix. We further introduced the notation L x 2 f to denote the Hessian of f with respect to x (t) and Z = ∫ F (y, x θ *, θ) d θ = p (y) / (2 π) n / 2. Let q (t|yk, x), t > yk denote the conditional intensity function, i. e. q (t|yk, x) dt is the probability for a spike in [t, t + dt] given an intensity function x (t) and the last spike at yk. We can express q (t|yk, x) in terms of the ISI probability density as [31] q (t | y k, x) = p (t, y k | x) 1 - ∫ y k t p (s, y k | x) d s. (15) The time rescaling theorem then states that the rescaled ISIs [31,47,61,62] τ k = ∫ y k - 1 y k q (s | y k - 1, x) d s, (16) are independent and identically distributed exponential random variables with mean one if y is a realisation from a point process with conditional intensity function q (t|yk, x). Suppose there are K rescaled ISIs. For a quantile-quantile plot [46], we order the τk in ascending order giving rise to the new ISIs τ ˜ n. We then plot the quantiles of the distribution of the τ ˜ n against the quantiles of an exponential distribution with unit rate, which are given by τ ^ n = - ln (1 - s n) with sn = (n − 0. 5) /K. For the Kolmogorov-Smirnov, plot [62], we define the random variable u k = 1 - e - τ k and then plot the ordered set of the uk against the cumulative distribution function of the uniform distribution, i. e. F (x) = x for 0 ≤ x ≤ 1, sampled at sn. The Ca2+ spike rate is estimated from m spike sequences via kernel smoothing (KS) through [49,63] r = 1 m ∑ j = 1 m ∑ i = 1 N j f (t - y i j, σ), (17) where y i j denotes the ith spike time in the jth Ca2+ spike sequence yj, and Nj is the total number of spikes in yj. The function f represents the kernel, and we chose a Gaussian of the form f (t, σ) = 1 2 π σ 2 exp (- t 2 2 σ 2). (18) The parameter σ is referred to as the bandwidth of the kernel. In case we work with a large number of independent Ca2+ spike sequences yj, we can use an optimal bandwith [49,50]. To evaluate how well a given method (e. g. Bayesian inference or KS) approximates the true Ca2+ spike rate used to generate surrogate data, we evaluated the normalised L2 norm as L 2 = [ ∫ 0 t (r ^ (t) − r ˜ (t) ) 2 d t ] 1 / 2 [ ∫ 0 t r ˜ (t) d t ] − 1, (19) where r ˜ and r ^ denote the known and estimated Ca2+ spike rate, respectively. We arrange the stimuli experienced by individual cells in a matrix X such that each row corresponds to a single stimulus time course. We then compute the singular value decomposition of X, i. e. X = UΣVt, where t denotes transposition. The columns of V correspond to the eigenvectors of XtX, and Σ is a diagonal matrix that holds the singular values of X. The weights of the principal components of the stimuli time courses are the rows of XV = UΣ. The k-means algorithm requires the number k of clusters as input and then determines the members of each cluster by minimising the error function [64] E = ∑ i = 1 k ∑ x ∈ C i ‖ x - μ i ‖ 2. (20) Here, x are the data points, C1, …, Ck are the k disjoint clusters and μi is the centroid of the ith cluster. We varied k and visually inspected the clustering. For consistency, we also clustered the data using other algorithms such as mean shift, spectral clustering and density-based spatial clustering of applications with noise. While there were minor differences between the suggested clusters, the overall clustering structure remained the same. To see which is the most likely ISI statistics, we apply the following protocol to every single cell from the experiment shown in Fig 1: To test the performance of Ca2+ spike rate estimation, we generated surrogate data from an IG for the two different intensity functions xdet (t) and xGP (t) using inverse sampling, a Bernoulli process based on the conditional intensity function in Eq (15) and time rescaling [47,65,66]. A key factor in estimating Ca2+ spike rates from PSTHs and KS is the choice of a bin width and bandwidth, respectively. For a large number of Ca2+ spike sequences, optimal estimates exist [49–51], and we use them for Fig 4. In case of only a few Ca2+ spike sequences with a small number of spikes per sequence, as in Fig 5, no estimates for a bin width or bandwidth exist. We therefore employed a bandwidth that was approximately equal to the optimal bandwidth determined in Fig 4 as well as bandwidths 1. 5 and 2 times larger than this. In addition, we used the same formal expression as for the optimal value, which resulted in the bandwidth σ ^. Note that σ ^ differs from the optimal bandwidth in Fig 4, since it explicitly depends on the number of Ca2+ spike sequences. A bespoke perfusion system connected to a 3-port microfluidics device [67] was used to expose cultured HEK293T cells to varying concentrations of the muscarinic receptor agonist, carbachol. The HEK293T cell line was a gift from Dr N. Holliday, University of Nottingham, that had been frozen after passage 28 of the original stock. After thawing, cells were used for up to a further ten passages. Cells were seeded at a density of 105 cells/ml in the central micro-channel of the microfluidic devices, in DMEM D6429 growth media (Invitrogen, Paisley, UK) containing 10% fetal calf serum. Cells were loaded inside the microchannels with 1 μM of the Ca2+ indicator Fluo5F-AM for 30 min, followed by washout with imaging buffer (135 mM NaCl, 3 mM KCl, 10 mM HEPES, 15 mM D-glucose, 2 mM MgSO4 and 2 mM CaCl2) for at least a further 30 min. To stimulate the cells, the flow rates of two inlet channels into the microchannel were varied, allowing the interface between the two solutions to be shifted laterally across the chamber. One inlet stream contained the agonist (100 μM carbachol) and Alexa Fluor 594 (AF594,2 nM; to allow monitoring of agonist concentration in proportion to AF594 fluorescence). The second inlet contained buffer alone. The interface formed between the two solutions due to laminar flow was shifted across the width of the microchannel by controlled changes in the fractional flow rates for each stream, with total flow being constant. In combination with the shifting interface position, the concentration gradient formed by diffusional collapse of the interface as the co-flow progresses through the channel length results in a spatiotemporal gradient in agonist concentration throughout the channel. This method enables the exposure of cells to pre-defined, time-varying changes in agonist concentration, from simple step-changes to complex waveforms. During dynamic stimulation with agonist, AF594 and Fluo-5F AM indicators were excited sequentially (100 ms exposure, 1 Hz frame rate) using a pE2 LED system (excitation peaks 470 nm and 565 nm; CoolLED, Andover UK). Emission was detected at 535 ± 50 nm and 565 ± 20 nm with an ORCA-R2 camera (Hamamatsu, Welwyn Garden City, UK). A time-series analyser plugin to ImageJ (Wayne Rasband, National Institutes of Health, Bethesda, MD, available at http: //rsb. info. nih. gov/ij) was used to manually define circular regions of interest (ROI) centred on each cell. Mean Fluo-5F emission intensity of pixels falling within each ROI was quantified and expressed as the ratio of fluorescence at time t divided by mean intensity from a 25 s window prior to the first increase in stimulus concentration (F/F0). The baseline window is selected as the window with minimum standard deviation from sliding 25 s windows taken from 0 to 120 s (before increase in stimulus concentration). Fluorescence of AF594 was quantified as the mean fluorescence intensity of pixels falling within each ROI being quantified; therefore each cell has a Ca2+ response measure and an associated stimulation profile.
Upon stimulation a large number of cell types respond with transient increases of the intracellular calcium concentration, which often take the form of repetitive spikes. It is therefore believed that calcium spikes play a central role in cellular signal transduction. A critical feature of these calcium spikes is that they occur randomly, which raises the question of how we can predict the timing of calcium spikes. We here show that by using Bayesian ideas and concepts from stochastic processes, we can quantitatively compute the calcium spike rate for a given stimulus. Our analysis also demonstrates that traditional methods for spike rate estimation perform less favourably compared to a Bayesian approach when small numbers of cells are investigated. To test our methodology under conditions that closely mimic those experienced in vivo we challenged cells with agonist concentrations that vary both in space and time. We find that cells that experience similar stimulus profiles are described by similar calcium spike rates. This suggests that calcium spike rates may constitute a quantitative description of whole-cell calcium spiking that reflects both the randomness and the spatiotemporal organisation of the calcium signalling machinery.
carbachol
action potentials
medicine and health sciences
engineering and technology
paracrine signaling
drugs
membrane potential
electrophysiology
neuroscience
endocrine physiology
probability distribution
mathematics
pharmacology
statistical distributions
fluidics
microfluidics
autocrine signaling
probability density
endocrinology
probability theory
signal transduction
cell biology
physiology
biology and life sciences
physical sciences
cell signaling
neurophysiology
A Bayesian approach to modelling heterogeneous calcium responses in cell populations