# Systematic Characterization of Dynamic Parameters of Intracellular Calcium Signals

^{1}Department of Physiology, McGill University, Montreal, QC, Canada^{2}Faculty of Dentistry, McGill University, Montreal, QC, Canada^{3}Shriners Hospital for Children-Canada, Montreal, QC, Canada

Dynamic processes, such as intracellular calcium signaling, are hallmark of cellular biology. As real-time imaging modalities become widespread, a need for analytical tools to reliably characterize time-series data without prior knowledge of the nature of the recordings becomes more pressing. The goal of this study is to develop a signal-processing algorithm for MATLAB that autonomously computes the parameters characterizing prominent single transient responses (TR) and/or multi-peaks responses (MPR). The algorithm corrects for signal contamination and decomposes experimental recordings into contributions from drift, TRs, and MPRs. It subsequently provides numerical estimates for the following parameters: time of onset after stimulus application, activation time (time for signal to increase from 10 to 90% of peak), and amplitude of response. It also provides characterization of the (i) TRs by quantifying their area under the curve (AUC), response duration (time between 1/2 amplitude on ascent and descent of the transient), and decay constant of the exponential decay region of the deactivation phase of the response, and (ii) MPRs by quantifying the number of peaks, mean peak magnitude, mean periodicity, standard deviation of periodicity, oscillatory persistence (time between first and last discernable peak), and duty cycle (fraction of period during which system is active) for all the peaks in the signal, as well as coherent oscillations (i.e., deterministic spikes). We demonstrate that the signal detection performance of this algorithm is in agreement with user-mediated detection and that parameter estimates obtained manually and algorithmically are correlated. We then apply this algorithm to study how metabolic acidosis affects purinergic (P2) receptor-mediated calcium signaling in osteoclast precursor cells. Our results reveal that acidosis significantly attenuates the amplitude and AUC calcium responses at high ATP concentrations. Collectively, our data validated this algorithm as a general framework for comprehensively analyzing dynamic time-series.

## Introduction

Cellular biology is vastly populated with dynamic processes, which can be altered dramatically or subtly by pathological causes. Calcium signals, characterized by fast and transient increases in cytosolic free calcium concentration ([Ca^{2+}]_{i}), which vary in amplitude and duration and can exhibit oscillatory dynamics with frequency-dependent downstream effects (Clapham, 2007), represent a prominent example of such dynamic processes (Figure 1A). To fully understand the data of such dynamic complexity, a robust methodology to analyse, and characterize these responses is necessary. Numerous studies have investigated [Ca^{2+}]_{i} dynamics, but the analysis have in many cases been limited to qualitative assessments (Cao et al., 1997; Frame and de Feijter, 1997; Jorgensen et al., 1997; Jørgensen et al., 2003; Isakson et al., 2001; Romanello and D'Andrea, 2001). Studies that have pursued quantitative analysis of calcium time-series reported a number of different, often non-overlapping characteristics of the response (Table 1). In cases where experiments were conducted on a smaller-scale, manual analysis was achievable. However, to achieve larger-scale analyses for experiments with hundreds of individual recordings, open-source signal-processing algorithms are required and becoming increasingly relied on to overcome these bottlenecks in productivity. None of 11 published algorithms we examined provided a comprehensive analysis of the entire response observed within a recording (Table 2). As a direct consequence of the lack of a standardized methodology to quantify such data-sets, findings from various studies are challenging to compare, relate, and generalize. Hence, the motivation of this study was to achieve faster analysis while standardizing the methodology involved, thereby minimizing user-bias, and ensuring consistency in the analysis of complex physiological signals. While such a tool may or may not change the conclusions of individual studies, it would improve comparability between different studies, and potentially enable meta-analysis of different experiments. We have therefore developed an algorithm that addresses these concerns and focused on the dynamic signals generated by purinergic (P2) receptors to demonstrate its utility.

**Figure 1. Characterization of dynamic calcium signals**. **(A)** Examples of heterogeneity of [Ca^{2+}]_{i} responses observed in various live-cell recordings; note many signature forms that TRs and MPRs can exhibit. **(B)** Analysis of parameters for single-peak TRs. Inset: Representative single-peak TRs. **(C)** Analysis of parameters for MPRs. Inset: Representative MPRs. Time of onset, t_{onset}; area under curve, AUC; full-width half-max, FWHM; activation time, t_{10–90%}; decay constant, τ_{decay}; periodicity, T; number of peaks, N_{osc}; oscillatory peak magnitude, E; width of oscillatory peaks, ξ^{peak}; point of inflection in deactivation phase of TR, *F*″(ρ). Arrows/text in red illustrate how parameters of interest are obtained.

Purinergic receptors that evoke intracellular responses upon extracellular stimulation with nucleotides, such as ATP and ADP, are known to induce complex [Ca^{2+}]_{i} signals. P2 receptors are subdivided into two families, P2X and P2Y receptors, which are omnipresent in virtually all mammalian tissue (Burnstock and Verkhratsky, 2009). The mammalian P2X receptor family, consisting of seven subtypes (P2X_{1–7}), are ionotropic ligand-gated cation channels that can permit the influx of extracellular calcium upon stimulation (Kaczmarek-Hajek et al., 2012). The mammalian P2Y receptor family, consisting of eight subtypes (P2Y_{1–2,4,6,11–14}), are metabotropic G-protein coupled receptors that can indirectly modulate the release of calcium from intracellular calcium stores through inositol triphosphate (von Kugelgen and Hoffmann, 2016). P2 receptors have been demonstrated to play an important role on bone physiology (Lenertz et al., 2015). Since individual bone cells commonly express multiple active P2 receptors (Gallagher and Buckley, 2002), responses to purinergic stimulation result in complex, concentration-dependent [Ca^{2+}]_{i} transients (Xing et al., 2016). While it remains difficult to experimentally isolate the contribution of individual receptors, a number of studies have demonstrated that various P2 receptor subtypes have distinct calcium response kinetics and signatures. For instance, various P2X receptors desensitize at distinct rates under sustained agonist stimulation (Koshimizu et al., 1999). P2X_{7}-mediated responses in particular are biphasic (Yan et al., 2010) and characterized by sustained [Ca^{2+}]_{i} elevation (Nobile et al., 2003). It is becoming increasingly clear, however, that P2 receptors cannot be studied and manipulated as individual components, but rather must be regarded as building blocks of a far more “dynamic architecture” that permits diverse functionality and flexibility (Volonte et al., 2006).

The goal of this study is to develop a universal signal-processing algorithm for MATLAB (MathWorks, Natick, MA) that would facilitate and standardize the parameter characterization of time series calcium imaging recordings containing prominent single transient responses (TR) and/or multi-peaked responses (MPR). All signals, no matter their complexity, can be reduced to a set of defined characteristics that describe the magnitude and kinetics of a given response. Based on our expertise and literature review (Tables 1, 2), we have selected the following parameters: time of onset after stimulus application (t_{onset}), activation time (time for signal to increases from 10 to 90% of peak; t_{10–90%}), and amplitude of response. Additionally, TRs are specifically described by their area under the curve (AUC), response duration (time between 1/2 amplitude on ascent and descent of the transient; FWHM), and decay constant of the exponential decay region of the deactivation phase of the response (τ_{decay}, Figure 1B) while MPRs are described by their number of peaks (N_{osc}), mean peak magnitude (E), mean periodicity (T), standard deviation of periodicity (σ_{T}), oscillatory persistence (time between first and last discernable peak; *l*_{osc}), and duty cycle (fraction of period during which system is active; ξ^{peak}/T, where ξ^{peak} is the width of the oscillatory peaks, Figure 1C). Since MPRs can be either stochastic or deterministic (Skupin et al., 2008; Dupont and Combettes, 2009; Dupont et al., 2011), the algorithm reports two sets of MPR parameters. The first describes MPR parameters for all the peaks present, while the second set reports the MPR parameters describing the subset of coherent oscillations, to omit the influence of stochastic processes, and to focus on the deterministic properties of the signal.

Live cell recordings will inevitably contain signal contaminations arising from experimental conditions and instrumentation, including (a) photochemical effects induced by the measurement process and (b) unrelated biological processes. While these imperfections are inherent to the experimental process, dynamic processes of interest can still be extracted from these recordings. This process in itself can be complicated and highly subjective depending on the extent to which the raw data are corrupted by noise and drift. Therefore, to reliably evaluate the magnitude and kinetics of these signals, we have developed a systematic way of first identifying unwanted signal contaminations, and then removing their effects when determining the parameters of interest. Following algorithm validation, we have also investigated the effect of acidosis on ATP-mediated [Ca^{2+}]_{i} responses in bone-marrow derived osteoclast precursors to demonstrate the efficacy of this algorithm in characterizing real-time cellular dynamics.

To implement the algorithm in MATLAB, the user is required to store the set of discrete sample points of the measured signal *F*(*t*), {*F*_{i}}, *i* = 1, 2, ⋯ , *N*, along with the corresponding discrete time points {*t*_{i}} in an Excel file (filename.xlsx). The algorithm can then be run using the MATLAB command

The MATLAB code required to execute this command along with examples of data and scripts are provided in Supplmentary Materials (Data Sheet 1).

## Materials and Methods

### Cell Culture

All procedures were approved by McGill University's Animal Care Committee and complied with the ethical guidelines of the Canadian Council on Animal Care. Bone marrow precursor cells were isolated from the femur and tibia of 6 week old FVB mice (Charles River), plated at a density of 7.5 × 10^{3} on 48-well glass-bottom plates (No. 1.5 Coverslip, 6 mm glass diameter, uncoated, MatTek Corp.) and cultured for 3 days in αMEM (12,000-022, GIBCO) supplemented with 10% FBS (080152, Wisent), 1% sodium pyruvate (600-110-UL, Wisent), 1% penicillin streptomycin (450-201-EL, Wisent), 50 ng/mL MCSF (300-25, Peprotech), and 50 ng/mL RANKL according to the protocol previously described (Boraschi-Diaz and Komarova, 2016).

### Intracellular Calcium Measurements

After 3 days of culture, osteoclast precursors were loaded with fura2-AM, a ratiometric fluorescent calcium dye (F1221, Invitrogen), incubated at room temperature for 30 min and washed twice with physiological solution (130 mM NaCl; 5 mM KCl; 1 mM MgCl_{2}; 1 mM CaCl_{2}; 10 mM glucose; 20 mM HEPES, pH 7.6). The final volume of 270 μL of physiological solution at pH 7.6 or pH 7.0 was added and cells were acclimatized for 10 min to reduce the effects of mechanical agitation that resulted from fura2-AM loading and washing. 10X ATP (Sigma) solutions were prepared in physiological solution at corresponding pH and 30 μL was added after 10 s of baseline [Ca^{2+}]_{i} recording to obtain a 1X dilution (i.e., final concentrations ranging from 1 μM to 10 mM ATP). [Ca^{2+}]_{i} was imaged for an additional 110 s at a sampling rate of 2 images per second using a fluorescent inverted microscope (T2000, Nikon). The excitation wavelength was alternated between 340 and 380 nm using an ultra-high-speed wavelength switching illumination system (Lambda DG-4, Quorum Technologies). Regions of interest (ROI) were manually defined and the ratio of the fluorescent emission at 510 nm, following 340 and 380 nm excitation (f340/f380), was calculated and exported using the imaging software (Velocity, Improvision). All data were imported into an excel spreadsheet for subsequent analysis.

### Validation and Statistical Analysis

Algorithm performance was evaluated using the algorithm generated figures for 450 individual signal fitting that enabled retrospective visual examination of both response-detection and quality of parameter fitting. Manual and automated estimates were compared using a correlation plot and Bland Altman analysis (Bland and Altman, 1986) to assess the degree of correlation and agreement, respectively. For correlation analysis, the line of exact linear correlation (i.e., *y* = *x*) is plotted as a reference to assess deviation of the linear regression curve from the desired 1:1 relationship between the manual and automated estimates. For the Bland Altman analysis, we compared the automated (*a*) and manual (*m*) parameter estimates of the i^{th} recording to obtain a Z-score, given by

where *x*_{i} = (*m*_{i} + *a*_{i})/2 is the average value of the estimated parameter, $\stackrel{\u0304}{x}$ is the mean value of *x* averaged over all recordings, and σ_{x} is the standard deviation of *x* overall recordings. Furthermore, the percent difference, Ξ_{i}, of the i^{th} recording is defined by

Ξ_{i} vs. *Z*_{i} were plotted to illustrate systematic biases. Negative values of Ξ_{i} were interpreted as manual estimates being greater than automated estimates, and vice versa. A quantitative estimate of the interval of agreement, within which 95% of differences lie, is defined by

where $\stackrel{\u0304}{\Xi}$ is the mean percent difference over all recordings and ${\sigma}_{\stackrel{\u0304}{\Xi}}$ is the standard deviation of $\stackrel{\u0304}{\Xi}$.

Experimental data were expressed as means ± S.E.M. Effect of ATP treatment under control conditions was evaluated using one-way ANOVA followed with Bonferroni *post hoc* test. The effect of acidosis was evaluated using two-way ANOVA with Bonferroni *post hoc* test. Results were accepted as significant at *p* < 0.05. Statistical analysis was performed in MATLAB.

## Results and Discussion

Although the notation used throughout the text implies that the recorded signals are fluorescence, the methodology remains the same for any other type of signals. For fluorescent recordings, the measured fluorescence *F* may consist of multiple parts: the drift, TR including the activation and deactivation phases, and the superimposed MPR. It can be expressed as the sum of the actual signal, *F*_{true}, and the normally distributed noise with a standard deviation σ.

To characterize parameters that reliably reflect *F*_{true}(*t*), *F*(*t*) is first preprocessed to remove the effects of noise (Section Noise Characterization) and to estimate the contributions of drift to *F*_{true}(*t*) (Section Baseline Drift). Next, the activation phase of TR is fit while simultaneously refining the estimated contribution of the drift. This approach allows us to determine if the recording is consistent with the expected model of a TR superimposed on a drifting baseline (i.e., whether activation phase is followed by a deactivation phase, Section Activation Fitting). If a TR is detected, we proceed by fitting the full set of TR model parameters simultaneously with the drift parameters [Section Transient Response (TR) Model]. In the case where there remain multiple significant deviations in the data from the TR model, we investigate and characterize the presence of oscillatory MPRs (Section Multi-Peaked Responses). The fitting of the TR is refined to remove the effects of the multiple peaks on the initial fit (Section Identifying Coherent Oscillations), in order to provide the best estimate of the baseline around which the MPRs oscillate. The deviations resulting from this secondary fitting of the TR are then analyzed to determine those resulting from coherent oscillatory processes (Section Characterizing Oscillatory Parameters). At each step throughout the fitting procedure, an updated estimate of the optimal set of parameters (e.g., the drift parameters) is obtained. These parameters are then used in a feed-forward manner, where the optimal set of parameters of the preceding fit is used as an initial guess for the subsequent step. This ensures that the algorithm produces high-fidelity fittings. Finally, the algorithm performance and utility is demonstrated with a new data set describing the effect of acidosis on ATP-induced calcium signaling in osteoclast precursors (Section Application to Pathophysiology).

### Noise Characterization

The first step in the processing of data is to evaluate four values that will be used for the remainder of the text: the derivative (*û*) of the noisy signal, the standard deviation of noise (σ_{−}), indices at which this noise is not prevalent (*j*), and the noise-to-signal ratio (ϕ). *û* will be used in Section Drift Delimitation to separate the TR from the underlying drift. Data points excluded by *j* (i.e., noise) will be omitted. The methodology detailed throughout the following section is an iterative procedure. In the instances where those quantities are used, we are referring to the value determined by the final iteration of the procedure.

#### Euler-Lagrange Formalism

The presence of noise in a recording renders naïve methods of derivative estimation inadequate (Chartrand, 2011). This is particularly exacerbated by the intermittent presence of large amplitude noise (spikes) related to the use of high gain settings on instrumentation. To reliably delimit (i.e., define the boundaries of) the drift in a recording of a noisy transient signal (Section Drift Delimitation), we have adapted the total-variational (TV) technique commonly used to estimate the first derivative of a signal contaminated with various types of noise (Chartrand and Staneva, 2008; Chartrand and Wohlberg, 2010; Chartrand, 2011). This technique performs better than the low-pass filter in distinguishing the drift from the transient response, as it does not indiscriminately remove high frequency components that affect the overall trend of the signal.

Our TV-based methodology seeks a function, $\widehat{\mathrm{\text{u}}}\left(\mathrm{\text{t}}\right)$, which represents the derivative of F_{true}(t), that solves the optimization problem.

The first term in Equation (2) is a regularization term which penalizes sudden changes in the derivative (to make the fitting smooth), the second term is an *L*^{2} data fidelity term, where *A* is the anti-differentiation operator (*Au* ≈ *F*_{true}), and α is a parameter dictating the balance between the two terms. In order to solve this minimization problem, we have to find the stationary solution to the following equation

derived from the Euler-Lagrange equation associated with Equation (2), where ${A}^{T}v\left(x\right)=\underset{x}{\overset{L}{\int}}vdx$ is the *L*^{2} -adjoint of *A*. Within the context of ratiometric fluorescent dyes (such as Fura2 AM used for [Ca^{2+}]_{i} recordings in this study) the recorded signal is the ratio of two Poisson random variables. The variance and the mean of such a signal follow a complex, and seemingly non-linear, function of the photon count rates at each wavelength. Since these rates are assumed to be unknown *a posteriori* in a recording, we cannot accurately determine how the noise is distributed. However, we will assume that instrumentation and experimental settings contribute to a noise distribution that is approximately Gaussian. Moreover, in the specific case of ratiometric dyes (Section Intracellular Calcium Measurements), we have found that true noise distribution is a complex function of time, but can be represented reliably using a time-dependent Gaussian noise model. Therefore, to reduce data-fidelity and conversely increase regularity in the regions of highest amplitude noise while accurately reproducing data in regions of lowest noise, we find instead a stationary solution to the following equation

where ψ(x, u) is an iteratively determined weighting function (as described below), and ε and η, are parameters introduced to avoid dividing by zero.

At the *n*^{th} iteration of the algorithm, we solve for *u*^{n+1} by setting the left-hand-side of Equation (4) to zero, and linearizing the problem through substituting every *u* appearing in the denominators by the value of *u*^{n} obtained from the previous iterate. For a more detailed descritption of the means used to solve this type of problem (see Vogel, 2002; Chartrand and Staneva, 2008; Chartrand and Wohlberg, 2010). It is known that an appropriate choice of the denominator offsets, ε and η, is necessary to produce acceptable minimizations (Chartrand, 2007), yet this choice is rarely considered beyond their status as parameters that must be tweaked to obtain acceptable results (Li et al., 2007; Chambolle et al., 2009; Oh et al., 2013). In what follows, we detail a methodology on how to determine the parameters α, ε, η, and the function ψ(*t*), based on the data *F*(*t*), and the derivative estimate *u*(*t*).

#### Dynamic Determination of Total-Variational Parameters

Given a set of fluorescence recordings {*F*_{i}} of length *N*, at the *n*^{th} iterate of the regularization algorithm, we identify the set of indices $j=\left\{\text{i}=1,2,\cdots \phantom{\rule{0.3em}{0ex}},\text{N}|{\psi}^{\text{n}}\left({t}_{\text{i}}\right)\ne \infty \right\}$ (as explained in Section Removal of Noise Spikes) whose data are not likely dominated by noise and thus should contribute to the fidelity term of Equation (4). Letting Δ_{i} = *F*_{i} − *F*_{i − 1}, we can define the weighting sequence ${g}_{\text{i}}=1-\left|{\Delta}_{\text{i}}\right|/\underset{\text{i}}{\text{max}}\text{\hspace{0.17em}}\left|{\Delta}_{\text{i}}\right|$ in order to provide an upper bound on the noise of the signal, given by

The weights *g*_{i} will tend to zero as Δ_{i} approach their maximum, and converge to a positive number (< 1) as |Δ_{i}| approach their minimum. We can thus conclude that the weighted average of |Δ_{j}| will identify the smallest differences as being the most informative of the magnitude of noise. Large discrete differenecs, whether they result from transient increases in the noise level or from the fact that the signal is non-stationary, contribute only modestly to the estimate ${\sigma}_{+}^{n}$. On the other hand, we can also estimate a lower bound on the noise using

where ${\zeta}_{\text{i}}^{n}={F}_{\text{i}}-{\left(A{u}^{n}\right)}_{\text{i}}$, i.e., by taking the difference between the data and the cumulative integral of *u*. Because we use the discrete differences Δ_{i} as our initial solution: ${u}_{\text{i}}^{0}={\Delta}_{\text{i}}$, the integral of *u*^{n} (*n* > 0) will likely diverge away from the data with each succesive iteration of the algortithm due to the action of the regularization term in Equation (4). This tends to result in ${\sigma}_{-}^{n}$ being smaller than ${\sigma}_{+}^{n}$, although this is not always strictly true. The use of the two different estimators for the noise allows for a more robust performance of the methodology, as both estimators are prone to becoming inaccurate in different scenarios. With an estimate of the noise, we can also estimate the noise-to-signal ratio by

The value of ϕ is a critical parameter in our algorithm as it discreminates between small and large values of various quantities. For example, it is employed to calculate an appropriate value of ε, defined in Equation (4), based on the scale of variations of small values of *u*′, as follows

The calculation of η requires defining another weighting sequence ${h}_{i}=\left|{u}^{\prime}{}_{i}\right|/\underset{\text{i}}{\text{max}}|{u}^{\prime}{}_{i}|$, which tends to zero when *u* is the most regular, as well as an estimate for the upper bound on the total error, χ, between the integral of *u* and the data in the least regular regions of the solution, given by

When χ is small, η must be made large enough to improve the smoothness of the fitting (at the expense of data-fidelity). This can be acheived by making η a decreasing function of χ. However, if a recording does not contain any rapid jumps or noise spikes (but is nonetheless noisy), such a relation between χ and η will not be sufficient to infer a proper choice of η given χ. Thus, we must include another term independent of χ which will produce modest data-fidelity for signals dominated by drift. We therefore define η to be

which includes the smallest non-zero scale of the weighted differences ${\delta}_{\text{i}}^{n}=\left(1-{h}_{\text{i}}^{n}\right){\Delta}_{\text{i}}^{n}$.

TV methods tend to smooth the fit when there are large amplitude jumps in the data, or where *u* is large. If a single large jump dominates the derivative, this can lead to excessive local smoothing, which can be resolved by having enhanced data-fidelity at that point. On the other hand, in the presence of large jumps in the data, small noise-driven fluctuations may be under-regularized. In this case, data-fidelity at these points must be reduced. The function ψ allows for local enhancement or reduction of data-fidelity. Unlike ε, η does not depend on ϕ, but data-fidelity must through the weighting function ψ (by making ψ proportional to ϕ). Therefore, we define ψ to be

where

${\tau}_{max}^{n}=\left(mean\left(\left|{\zeta}_{j}^{n}\right|\right)+std\left(\left|{\zeta}_{j}^{n}\right|\right)\right)/{\sigma}_{+}^{n}$ is a iterative error scale parameter, ${\upsilon}_{\text{i}}^{\text{n}}={\text{u}}_{\text{i}+1}^{\text{n}-1}+{(1-{\text{h}}_{\text{i}+1})}^{2}{\text{u}}_{\text{i}-1}^{\text{n}-1}$ is an estimator of *u*_{i} based on the adjacent values of *u*, and ω_{i} is a three-point moving average of δ_{i}. The exponent *J*_{1} emphasizes data-fidelity (regularity) when the derivative of the following point is large (small), whereas the exponent *J*_{2} emphasizes data-fidelity when Δ_{i} are small or when the derivative is near its maximum. When Δ_{i} is large and *u*_{i−1}, *u*_{i} are regular, on the other hand, *J*_{2} is small allowing regularity to propagate forward into regions of signal possessing large ampltidue noise (i.e., where the data is not informative). The balance between the two effects of *J*_{2} along with the χ -dependent terms of Equation (6) produce an acceptable compromise between the regularity of the fit and data-fidelity for recordings across a wide range of signal-to-noise ratio and extremely varied dynamics. Finally, once the maximum relative change between two iterations of the procedure becomes less than $\sqrt{\varphi}$, we consider the solution to have achieved quasi-stationarity and terminate the procedure.

#### Removal of Noise Spikes

When a recording exhibits intermittent periods of high amplitude noise (noise spikes), the data contaminated by these noise spikes is minimally informative. Detecting them allows for determining the indices *j* (Figure 2A). Within our regularization algorithm, this is done by (a) setting ψ to infinity at those time points in such a way that Equation (4) only penalizes irregularity at these time points, and (b) determining the fit at these points based on the surrounding (reliable) data. After each iteration, *n*, of the regularization algorithm, a smoother fit, *Au*^{n}, of the data is obtained. We also obtain a criterion that determines whether or not each point represents a noise spike based on a comparison between the residual differences, ζ^{n}, and a chosen threshold value. This threshold is specified using the parameter ${\tau}_{rm}^{n}$, given by

where $\xi =exp\left(-\left(\left({\sigma}_{+}^{n}-{\sigma}_{-}^{n}\right)/max\left\{{\sigma}_{+}^{n},{\sigma}_{-}^{n}\right\}\right)-{\left({\chi}^{n}\right)}^{2}/N\right)$ is a convergence parameter for the noise rejection method. Positive ${\zeta}_{\text{i}}^{n}$ are rejected if they are greater than ${\tau}_{rm}^{n}{\sigma}_{+}^{n}$, while negative ${\zeta}_{\text{i}}^{n}$ are rejected if they are less than ${\left({\tau}_{rm}^{n}\right)}^{2}{\sigma}_{+}^{n}$. The use of two thresholds is due to the asymmetry of the Poisson statistics underlying data collection using photodetectors. Rejection is achieved by setting ψ_{i} = ∞, which serves as the basis for defining the set of indices *j* (i.e., $j=\left\{\text{i}=1,2,\cdots \phantom{\rule{0.3em}{0ex}},\text{N}|{\psi}^{\text{n}}\left({t}_{\text{i}}\right)\ne \infty \right\}$ as stated before).

**Figure 2. Multitude of features that must be considered during parameter characterization. (A)** Experimental recordings are often contaminated with noise; red crosses represent noise detected and corrected. **(B)** Activation-coupled drift; Trajectory of baseline drift, d_{1}, can shift to secondary drift, d_{2} at time t_{2}. **(C)** Determination of activation in the presence of gradual activation. Region of activation is restricted to $\left[{t}_{{\text{i}}_{+}^{0}},{t}_{{\text{i}}_{local}^{0}}\right]$, within which the most likely point of activation, ${t}_{{\text{i}}_{act}^{0}}$, is statistically resolved. **(D)** Non-monotonic activation. Spurts of negligible activity, represented by the local maxima at ${t}_{{\text{i}}_{local}^{0}}$ and ${t}_{{\text{i}}_{local}^{1}}$, are disregarded as points of activation, in favor of ${t}_{{\text{i}}_{local}^{2}}$ which is characterized by the max peak, ${t}_{{\text{i}}_{\text{max}}^{0}}$. **(E)** Drift discrimination in cases where MPRs are superimposed on TRs. Underlying TR serves as a non-stationary baseline around which the MPR will oscillate within $\left[{t}_{2},{t}_{{\text{i}}_{end}^{0}}\right]$ until its contribution becomes negligible and the baseline drift, d_{2}, dominates the MPR baseline. **(F)** Complex drift; Change in the baseline drift, t_{shift}, prior to onset of TR is ignored in favor of d_{1} after the reweighting procedure outlined in Section Drift Fitting.

### Baseline Drift

When fitting data to specific functional forms, it is important to take into account temporal drifting of the baseline in a signal. The specific nature of the processes underlying this drift are not of particular interest here. Rather, we consider them as nuisance trends and aim to remove their effects from the data.

#### Drift Model

Some drifts are quite slow compared to the timeframe of the experimental recording, and thus can be fairly well represented by linear functions, whereas others are fast and better represented by an exponential decay function. Since, in a given recording, a number of processes will result in the observed drift, we postulate that, the drift throughout the signal can be well fit to a combination of linear and exponential functions

where, *a*_{ι}, τ_{ι}, and *m*_{ι} are the exponential amplitude, exponential time constant and the slope of the linear component of the ι^{th} drift within the signal, respectively (examples of fitted drifts shown in Figures 2B,E,F). With an appropriate choice of data, these are determined using the least squares fitting as discussed in Section Drift Fitting. We have found numerous cases where either the linear or exponential components were not justifed. However, this knowledge is unavailable to us prior to conducting manual or automated analysis of the data, and we cannot assume *a priori* a less general form than Equation (7). Thus, we have to rely on the fitting to optimize the contribution of the two compoenets in a data-dependent manner.

It is possible for the signal to exhibit (multiple) drifts with different trends separated perhaps by a TR. In order to capture this effect in a signal, we use a global drift model that combines the intial and secondary drifts (Figure 2B) in a semi-piecewise manner in which the initial drift, d_{1}, continues to contribute to the overall observed drift and succesfully captures the global behavior. This can be written as

where *t*_{2} is the time at which the secondary drift, d_{2}, begins and *z* is the offset at *t* = 0. Rather than assuming that the drift is similar for all recordings and attempting to construct a standard curve, we assume that a few points in each recording are highly informative of the drift in the baseline.

#### Drift Delimitation

To fit the drift model to the corresponding portion of the recorded signal, it is necessary to delimit the boundaries of TR by identifying the start of activation and end of deactivation. Firstly, we will aim to estimate the point in time at which the TR of a recording begins. Experimental TRs rarely activate abruptly and simultaneously for a field of imaged cells. Many factors play a role in the heterogeneity of responses observed such as variable diffusion fronts of applied agonist or heterogenous receptor expression among cells. These effects can manifest as very gradual rises or additional small amplitudes prior to a certain activation threshold being surpassed and a rapid activation phase being observed (Figures 2C,D). Therefore, analyzing the activation times of all components of the biological unit manually can be subjective. While a simple threshold value in the signal can be effective in detecting activation, the choice of the threshold requires some knowledge of the amplitude of the noise σ and is complicated by the presence of drift in the signal. Instead, we determine the end points of the time intervals dominated by the drift through statistical analysis of an estimate of the first derivative of the signal *û*(*t*). In other words, to distinguish the drift from TR, the derivative of the latter must change in a way that is more statisically significant than that of the former. The methodology used to obtain the estimate for the derivative is detailed in Section Noise Characterization.

Using the estimated time derivative, we aim to determine (i) the earliest possible time of activation ${t}_{{\text{i}}_{\text{act}}^{0}}$ (defined as the last time point exhibitng a significant increase in the first derivative before it reaches its maximum value), (ii) the time at which the activation reaches its maximum value *t*_{imax}, and (iii) the time at which the deactivation ends and the signal is once again dominated by the drift ${t}_{{\text{i}}_{\text{end}}^{0}}$.

i Assuming that the estimated first derivative during activation reaches a local maximum, we can find the most significant local maxima of *û* at the location

where we denote the first significant local maximum by ${i}_{local}^{0}\text{}=\text{}min\left\{{\text{i}}_{local}\right\}$. By focusing on the portion of the signal containing the first drift and the activation phase of the TR, we restrict our attention to the time interval in which the first derivative *û*(*t*) is non-negative to disentangle the effects of activation and drift in the data. In other words, we restrict our analysis to $\left[{t}_{{\text{i}}_{+}^{0}},{t}_{{\text{i}}_{\text{local}}^{0}}\right]$, where the derivative is non-negative and

Without prior knowledge about the sign of the derivative of the baseline drift, we cannot conclude that ${i}_{+}^{0}$ corresponds to the beginning of activation. To resolve this issue, we employ a statistical test to determine the likely time at which the activation occurs, located at the index

This methodology, based on the properties of the derivative around its first significant local maximum, generally picks out the first visually unambiguous activation (Figure 2C). As a result, it may be necessary to trim recordings where there are (large amplitude) artifacts prior to the activation of interest.

ii. TRs may be produced by the action of multiple active units (e.g., different receptor species) within the biological system under consideration, each having distinct properties and activation times. This leads to multiple delayed activations taking place over a broad range of time (Figure 2D). Due to the superposition of the drift in the baseline with these responses, it is entirely possible for recordings to be contaminated by strongly decreasing drift and for the expected maximum TR to not coincide with the actual maximum of the data. Thus, we have developed a method to search for the visually most likely point at which the TR reaches its maximum in the presence of a drift. For each one of the *N*_{local} significant local maxima of the first derivative along the activation phase, we find the location of the previous local minimum of the TV estimate of the data using

as well as the location of the next local maximum at

From the positions of the local extrema of the data, we can estimate the value of the baseline drift from each local minimum to the next local maximum using the linear extrapolation

where *A* is the operator of antidifferentiation with *Au* ≈ *F*_{true} (see Equation 1). Based on this, we then estimate the average rate of activation for each significant local maximum of the derivative according to

and select the first local maximum at the location

which has an average rate and magnitude within a statistically acceptable range (that excludes small outliers). ${t}_{{\text{i}}_{\text{max}}^{0}}$ is the time point at which the derivative reaches a local maximum. It may differ from the one that corresponds to the local maximum immediately following the first activation time point ${t}_{{\text{i}}_{\text{act}}^{0}}$ (Figure 2D). This is because ${i}_{act}^{0}$ depends solely on the derivative around its first local maximum, while ${i}_{max}^{0}$ takes into account an approximation to the average rate of change around all local maxima; the local maximum after ${i}_{act}^{0}$ should only differ from ${i}_{max}^{0}$ when there is a succession of activations and the first does not have the largest rate of activation.

Starting from this local maximum of the derivative, we seek the location of the first point in time ${t}_{{\text{i}}_{\text{max}}^{\text{end}}}$ where the change in the signal drops below the estimated noise level σ_{−} (see Section Dynamic Determination of Total-Variational Parameters)

and thus isolate the time interval in which the most significant activation occurs. Having isolated the most significant activation, we finally arrive at the location of the first estimate of the time where the response reaches its maximum

iii. Having identified the time at which activation is likely to begin ${t}_{{\text{i}}_{\text{act}}^{0}}$, we can now assume that the signal prior to this point is the drift. If the recording is of a sufficiently long duration, the response will return to baseline and the end of the recording should once again be dominated by the (secondary) drift (Figure 2A). To account for this drift, we need to estimate the time duration of deactivation. This is done in a manner nearly identical to how we determined the first time of activation ${t}_{{\text{i}}_{\text{act}}^{0}}$. However, due to the possibility of having MPR after the initiation of TR (Figure 2E), we cannot restrict ourselves to time intervals in which the first derivative is non-positive. To solve this issue, we define a set of time points after the presumed maximum of TR, located at i_{decay} = {i_{max}, …, *N*}, and construct a measure

to quantify how far the derivative deviates from its mean during the decay. This measure allows us to robustly detect the location of the time point where TR is negligible

The time point ${t}_{{\text{i}}_{\text{end}}^{0}}$ is then used to determine the start of the new drift.

#### Drift Fitting

To isolate the TR from the drift, it is necessary to generate an accurate fit for the drift. This is achieved by employing a succession of least square fits that progressively incorporates more data and models that account for additional components of the signal (including activation and deactivation phases of TR). The first step in this successive least-square-fitting method is to obtain preliminary estimates of the parameters θ_{drift} = [*a*_{1}, τ_{1}, *m*_{1}, *a*_{2}, τ_{2}, *m*_{2}] of the drift model. The MATLAB implementation of the non-linear least squares method (Marquardt, 1963; Moré, 1978) is used. More specifically, we initially minimize the error function between the drift model and the data

where the set of indices *k* is defined by $k=\left\{\left\{1,\dots ,{i}_{act}^{0}\right\}\bigcup \left\{{i}_{end}^{0},\dots ,N\right\}\right\}\bigcap j$ (${i}_{act}^{0}$ and ${i}_{end}^{0}$ are as defined in Section Drift Delimitation while *j* is defined in Section Removal of Noise Spikes). We denote the first estimate of parameters obtained from the minimizaiton of Equation (11) by ${\theta}_{drift}^{0}$. There is no guarantee, however, that the drift model described by Equation (8) can accurately represent the actual drift in the baseline. The emergence of drifting trends which are not related to the onset of TR (Figure 2F) may require the inclusion of more than two functions of the type described in Equation (7), yet the decision to include more drift terms, or alternatively to truncate the signal, would require manual intervention. To circumvent this limitation, we perform a second fit where the individual terms in Equation (11) are weighed according to a weight function *w*, and the first derivative of the data is taken into account. Although the set of parameters ${\theta}_{drift}^{0}$ will be able to produce the “general” trends in the data before and after TR, the presence of multiple drifting trends necessitates the use of the derivative estimate, *û*, to identify the segment(s) of the signal that actually follow the trend described in Equation (8) and remove the effects of others. The inclusion of *û*, however, in the sum of squared errors can lead to erroneous fittings when the drift is not well represented by Equation (8). The weight function *w* alleviates this problem by preventing the linear and exponential trends in Equation (7) from growing unjustifiably large.

To determine *w*, we require that it approaches zero when the match between *û* and $\u1e0a\left(t;{\theta}_{drift}^{0}\right)$ is minimal. To achieve this, we choose *w* to depend on $y=|\xfb-\u1e0a\left(t;{\theta}_{drift}^{0}\right)|$ as follows

We then apply the non-linear least squares fitting procedure to minimize the error function

and use ${\theta}_{drift}^{0}$ as an initial condition for the fitting procedure.

### Activation Fitting

For the activation phase of the response, we use the model

where *A*_{act} is the maximum of the Hill function, *n*_{act} is the Hill coefficient, β is the time at which the the Hill function reaches its half maximum, and *m*_{act} is the slope of the quasi-linear function that accounts for the trend which dominates at the end of activation. The Hill function allows for a rapid rise and switch in convexity due to many biological units being activated at once, whereas the linear trend allows for a delayed and slower rise induced by more units being progressively recruited into the generation of the signal. The values of *A*_{act} and *m*_{act} determine the magnitude of these two trends, whereas *n*_{act} and β affect primarily the timescales of switching between the two. The time at which the activation phase begins is denoted as *t*_{on}, and its estimate is confined to the time interval $\left[{t}_{{\text{i}}_{\text{act}}^{0}},{t}_{{\text{i}}_{\mathrm{max}}}\right]$. In order to obtain a preliminary estimate of the values of these parameters and those of the drift function, θ_{act} = [*t*_{on}, *A*_{act}, β, *n*_{act}, *m*_{act}, θ_{drift}], we minimize the sum of square errors between the activation data and the model along with their derivatives, given by the function

where ϕ is the noise-to-signal ratio as defined in Section Dynamic Determination of Total-Variational Parameters, and

is our activation data model that takes into account the effect of the drift *D* (defined in Equation 8) on the perceived activation. To abrogate the influence of noise spikes on measured parameters, the sum of squares in Equation (12) is evaluated at the set of indices *j* (defined in Section Dynamic Determination of Total-Variational Parameters) which are not dominated by noise. Moreover, we use ${\widehat{\theta}}_{drift}^{act}$ to denote the set of drift parameters obtained from the minimization of Equation (12).

#### Signal Detection

In order for the algorithm to resolve whether there is a discernable TR present in *F*_{true}(*t*), three conditions must be satisfied: (i) i_{local} must be a non-empty set, (ii) the initial drift estimate

must be below the TV data estimate by a detection threshold of 4σ_{−} for at least six data points, and (iii) the difference between the TV data estimate and the initial drift estimate must not be a strictly increasing function of time after ${t}_{max}^{0}$. These three criteria allow for the detection of the TR and further analysis of its characteristics [Section Transient Response (TR) Model].

To evaluate the signal detection performance of the algorithm, 450 individual traces of ATP-induced calcium responses were used as a validation set. Manual results were then compared to automated detection of TRs to assess extent of agreement between the two methods. Classical signal detection nomenclature (i.e., true positive or negative and false positive or negative) was intentionally avoided due to lack of certainty in determining the true presence of TRs in more ambiguous cases. We found that the automated and manual methods agreed in detecting a TR in 88.1% of cases, and they disagreed in 11.9% of cases (Figure 3A). Further dissection of these results showed that 64.7% of disagreement arose from the algorithm reporting an absence of TR while visual evaluation suggested otherwise (Figure 3A), indicating that the algorithm has a tendency to be more conservative than user-mediated assessments.

**Figure 3. Signal detection performance of algorithm. (A)** Signal detection analysis demonstrates relative selectivity and sensitivity of the algorithm as compared to user-mediated manual detection. “Agreement” refers to instances where the detection of a response was consistent between algorithm and manual methods. “Disagreement” refers to the contrary. “Present” indicates cases where a response was detected using the indicated method, while “absent” refers to cases where a response was not detected. **(B)** Representative time-series for pronounced (top) and obscure (bottom) signals. **(C)** Analysis of disagreement cases for pronounced and obscure signals. Automated detection, auto; manual detection, manu; sample size, n.

To determine whether there were particular types of recordings that contributed to these disagreements, time-series traces were qualitatively divided into two groups: Clean signals with clearly defined responses were classified as “pronounced” (Figure 3B, top), and signals containing ambiguous signals with low signal to noise ratio or strong drift were classified as “obscure” (Figure 3B, bottom). The total frequency of disagreement was 3.6 times greater for obscure signals compared to those classified as pronounced (17.4 vs. 4.8%; Figure 3C). Regardless of the group, the algorithm signal detection remained more conservative compared to the manual method.

### Transient Response (TR) Model

Transient cellular responses are generally complex with multiple time scales and amplitudes. They may, in fact, exhibit prolonged MPRs superimposed on a more acute response (see Figure 2E). For these reasons, a complete characterization of all possible TRs is unlikely to be attainable. In order to remain as general as possible, we propose modeling TR as a continuously differentiable piece-wise defined function that first increases during the activation phase and subsequently decreases during the deactivation phase. Due to the large number of parameters required for such a description and the automated nature of our fitting procedure, we decompose the fitting of the whole TR into a sequential fitting of the activation phase alone followed by a fit of both phases simultaeously. This yields significantly more reliable results with faster convergence rates over a wide gamut of input data, because it allows for information obtained during preliminary simple fits to be used in a progressively more complex manner. In order to capture the complex fluorescence response generated by the spatially separated units in a live cell, we use a combination of Hill functions and quasi-linear functions generated by the integral of the corresponding Hill functions.

#### Response Fitting

Following the least squares fitting of the activation data, we now seek to fit the entire recording (with the drift and TR) using a continuously differentiable function. A decreasing Hill Function is used to describe deactivation phase of the signal in a manner similar to Equation (11), as follows

where *A*_{de} is the amplitude of the Hill function, *n*_{de} is the Hill coefficient, γ is the time at which the Hill function reaches its half maximum, and *m*_{de} is the slope of the quasi-linear function that accounts for the trend dominating at the beginning of deactivation. The time when the response switches to the deactivation function is denoted by *t*_{de}. The parameter *m*_{de} is chosen such that the response function returns to zero by the end of the recording and is given by

where ℓ_{act} = *t*_{de}−*t*_{on} is the time duration of the activation phase of the response. If differentiability is not enforced at the point *t*_{de}, where the two functions *g*_{act} and *g*_{de} meet, then the fitting may contain sharp edges indicative of unconverged solution. To solve this issue, the continuity of the first derivative of these two functions, particularly at *t*_{de}, is achieved through a third-order Hermite spline (Traub, 1964) on the time interval [*t*_{de}−ς_{act}, *t*_{de}+ς_{de}], where ς_{act} and ς_{de} are two parameters that must satisfy ς_{act} < ℓ_{act} and ς_{de} < 2γ (see Data Sheet 2 in Supplementary Materials). The overall response model is thus given by

where θ_{resp} = [θ_{act}, *t*_{de}, ς_{act}, ς_{de}, *A*_{de}, γ, *n*_{de}θ_{drift}]. Given the response model, we define the global data model as

and minimize the error function

to obtain the fitting, where κ is a parameter quantifying the apparent coherence between the drift and reponse models (*D* and *g*_{resp}), given by

λ is a parameter defined by

and the values for the parameters *A*_{act}, *n*_{act}, β, *m*_{act}, and *t*_{on} are taken from the activation fitting. The first two terms in the sum of squares in Equation (14) are analgous to those in Equation (12), whereas the third term minimizes the AUC for the second drift function *D*_{2}. By including the coefficient κ in this third term, however, allows *D*_{2} to become more significant when there is a large mismatch in the value of the baseline between ${t}_{{\text{i}}_{\text{act}}^{0}}$ and *t*_{N} (that cannot be explained by *d*_{1}).

#### Activation Parameter Validation: t_{onset}, t_{10–90%}, and Amplitude

The three parameters, t_{onset}, t_{10–90%}, and amplitude, are considered together because they describe what happens at the activation phase of the TR, with no regard for the deactivation phase or MPR. The time at which a TR is discernable from baseline is defined as t_{onset} and is estimated directly as the parameter *t*_{on} of the response function *g*_{resp}. There was a strong linear agreement between the manual and automated-estimates, with a linear slope of 1.06 and a correlation coefficient (*r*^{2}) of 0.94 (Figure 4A, left). On average, automated estimates of t_{onset} were 4.6% greater than manual estimates, with limits of agreement ranging between −10 and 19% difference (Figure 4A, right). The t_{10–90%} is also estimated numerically from the response function *g*_{resp}. Using the response function allows to overcome issues arising from the subsampling of rapid dynamics by numerically evaluating on a time grid 10 times finer than the input times. Manual estimation of this parameter is contingent upon accurate estimation of the baseline and peak occurrences, both of which present potential sources of error, particularly for noisy signals with drift. The relationship between manual and automated estimates had a linear slope of 0.77 and a correlation coefficient of 0.77 (Figure 4B, left). The higher degree of scatter away from the line of exact linear correlation is reflected by the wider Bland–Altman interval of agreement, ranging from −89 to 54% difference between manual and automated estimates. Overall, there was a −17% difference between all paired estimates of t_{10–90%}, revealing that t_{10–90%} was manually overestimated compared to the automated estimates (Figure 4B, right). Amplitude estimates obtained by the manual and automated methods had a strong linear relationship with a slope of 1.02 and a correlation coefficient of 0.84 (Figure 4C, left). The limits of agreement, ranging from −26 to 25%, were narrow with a mean percent difference of −0.3% between all paired estimates of amplitude (Figure 4C, right).

**Figure 4. Validation of activation parameter measurements**. Parameters were estimated manually and then compared with automated-estimated values using correlation analysis (left: solid line, linear regression line; broken line, line of equality) and Bland-Altman analysis (right, solid line: mean percent difference, $\stackrel{\u0304}{\Xi}$; broken line: limits of agreement, $\stackrel{\u0304}{\Xi}$ ±1.96 ${\sigma}_{\stackrel{\u0304}{\Xi}}$). **(A)** Time of onset, t_{onset}. **(B)** Activation time, t_{10–90%}. **(C)** Amplitude. Arbitrary units, A.U. Insets: Visual representations of parameters measured.

#### TR Parameter Validation: AUC, FWHM, and τ_{decay}

Due to the inherent differences between TRs and MPRs and the approach taken with this algorithm, the AUC, FWHM, and t_{decay} are limited to describing TRs. Nevertheless, these parameters will be also reported in the presence of MPRs where they should be interpreted with the following considerations. (i) if a MPR is superimposed on a TR, the reported parameters describe the underlying TR, not the superimposed MPR. (ii) if TR presence is not detectable and MPR demonstrates a purely oscillatory response, the reported parameters characterize the first peak only. With these considerations in mind, manual evaluation of TR parameters was performed with a variety of signals, including MPRs. AUC estimates are manually determined using a geometric estimation of the area of a triangle whose vertices are at the start, peak, and end of the TR. Algorithmically, AUC was evaluated from the area under *g*_{resp}, using the trapezoidal rule (Rice, 1973) implemented in MATLAB using “trapz().” Comparing the manual and automated estimates demonstrated a linear relationship with a slope of 0.85 and a correlation coefficient of 0.79 (Figure 5A, left). On average, automated estimates were 4.6% larger than manual estimates with an interval of agreement ranging between −38 and 47% difference (Figure 5A, right). Considering the geometric-approach used for manual-estimation of AUC values, it is reasonable to assume that the error arose from manual limitations.

**Figure 5. Validation of TR parameter measurements**. Parameters were estimated manually and then compared with automated-estimated values using correlation analysis (left: solid line, linear regression line; broken line, line of equality) and Bland-Altman analysis (right: solid line, mean percent difference, $\stackrel{\u0304}{\Xi}$; broken line, limits of agreement, $\stackrel{\u0304}{\Xi}$ ±1.96 ${\sigma}_{\stackrel{\u0304}{\Xi}}$). **(A)** Area under curve, AUC. **(B)** Full-width half-max, FWHM. **(C)** Decay constant, τ_{decay}. Arbitrary units, A.U. Insets: Visual representations of parameters measured.

Due to the difficulty in determining the precise time at which the signal returns to its former baseline, it would be challenging to manually describe the duration and decay characteristics of the response. Incomplete recordings and background drift are largely responsible for generating such behavior. The full-width half-max (FWHM) is defined as the time elapsed between the two half-max coordinates of a peak. Our analysis of FWHM revealed that the linear relationship between manual and automated-estimates was strong, with a correlation coefficient of 0.87 and slope of 0.95 (Figure 5B, left). The Bland Altman analysis demonstrated that the agreement interval ranged from −20 to 22% difference with a mean percent difference of 0.85% between all paired estimates (Figure 5B, right).

Finally, to manually estimate the decay constant, the time of the inflection point (ρ) is visually estimated and the general trend of the data following ρ is represented by a mono-exponential decay. The decay constant is then determined by the time it takes for the signal to reduce to approximately 37% of its initial value (1/e). Algorithmically, the time ρ is determined by solving for the inflection point in the deactivation function, given by

The data following ρ is then fit to a mono-exponential decay function using least squares, to determine the time constant of decay. The slope of the linear agreement between manual and automated-estimates was 0.81 and the correlation coefficient was 0.75 (Figure 5C, left). The Bland Altman analysis revealed a strong systematic bias of −14% difference, with an agreement range of −63 to 34% difference, signifying that manual efforts to estimate the decay constant consistently overshot the values reported by the algorithm (Figure 5C, right).

### Multi-Peaked Responses

To isolate the characteristic parameters of MPRs that are frequently superimposed on TRs (Figure 2D), the TR model *g*_{resp} must be refined to serve as a non-stationary baseline around which the MPR will oscillate. This refinement is necessary because it is often the case that the TR model *g*_{resp} will produce sub-optimal fittings where the data deviates significantly from the TR fitting. Therefore, to accurately characterize the TR and quantify the properties of truly oscillatory MPRs, it is first necessary to adapt the least squares fitting procedure of Section Response Fitting to remove the effects of data points not well represented by *G*_{resp}. This is done by first identifying large deviations representing MPRs from the *G*_{resp} -fitting (obtained by minimizing Equation 14), and then reweighing the sum of squares in Equation (14) to remove the influence of those deviations from the fits. We subsequently perform a secondary fitting of *G*_{resp} to determine more accurately the baseline, delineating the TR, where the MPR-associated deviations originate from. Finally, we analyze the MPRs by determining whether they represent oscillations and, if so, quantify their properties.

To identify these MPR-associated deviations from the newly defined baseline, we first employ the MATLAB “findpeaks()” function. This finds the peaks and troughs of the significant deviations in the TV data estimate, *Aû*(*t*), from the TR estimate, ${G}_{resp}\left(t;{\widehat{\theta}}_{resp}\right)$, truncated at its half-maximums, where ${\widehat{\theta}}_{resp}$ is the optimal parameter set obtained from minimizing Equation (14). This truncation permits for the possibility that the onset of TR coincides with the first peak of the MPR-associated deviations. The implementation of peak-finding algorithm, on the other hand, allows for the specification of minimum heights and timing between peaks, set to be 6σ_{−} and 5 s, respectively. The algorithm yields the heights (*E*^{peak}, *E*^{trough}), the FWHM (ξ^{peak}, ξ^{trough}), and the times (${\stackrel{~}{t}}^{peak}$, ${\stackrel{~}{t}}^{trough}$) of significant peaks and troughs of the deviations, respectively (see Figure 2D). In total, there are Ñ = *N*_{peak}+*N*_{trough} of these deviations, including *N*_{peak} peaks and *N*_{trough} troughs. If Ñ ≤ 2, then the only deviation in the signal is the TR and the algorithm can terminate. Without prior knowledge of the nature of the MPR-associated deviations, it is very difficult to determine whether they result from trends which are above, below, or symmetric to the TR. To resolve this issue, we assume that the estimated baseline from Section Response Fitting underlies the signal in the absence of deviations. To incorporate this assumption algorithmically, we define two bias parameters based on the relative heights of the first peak and trough, as follows

These quantities are then used to calculate weighting functions for the data based on the properties of peaks

and troughs

The weighting functions ${\pi}_{\text{i}}^{peak},{\pi}_{\text{i}}^{trough}$ quantify the relative reliability of the data around each deviation based on how close it is to the fitting function *G*_{resp} and on its duration. We can also assess the reliability of each time point of the recording (including the TR, the drift, and any MPR-associated deviations present) by how well its derivative matches ${\u0120}_{resp}\left(t;{\widehat{\theta}}_{resp}\right)$. This is done using another weighting function, defined by

We combine these weighting functions using the criterion that for a data point to be reliable, it must have either a large value of π^{peak} or π^{trough}, and a large value of Γ. It is implemented in the weighting function Ω, as follows

With Ω_{i}, we can fit the TR reliably in the presence of significant deviations from the model of Equation (13). This is done by minimizing the error function

#### Identifying Coherent Oscillations

Not all MPRs correspond to periodic oscillations (Thurley et al., 2014). To address this, the algorithm reports two sets of MPR parameters, the first to describe all the peaks detected within a MPR, and the second to describe the subset of coherent peaks present within the same MPR. This section focuses on how the subset of coherent oscillatory peaks is identified. We use a clustering algorithm which is an unsupervised learning technique that enables for the identification of natural groupings or patterns with a defined data set. By minimizing Equation (15), we obtain the most reliable estimate of the TR (specified by the model *G*_{resp} and its optimal parameter set ${\widehat{\theta}}_{resp}^{\Omega}$), which we take to be the baseline of the MPR-associated deviations. Given the estimate ${G}_{resp}\left(t;{\theta}_{resp}^{\Omega}\right)$, we repeat the peak finding steps detailed in Section Multi-Peaked Responses. To determine whether or not the detected deviations represent oscillations, we use a Gaussian mixture model clustering technique. It groups together (in clusters) peaks and troughs with comparable periods, *T* (determined by the difference between two consecutive peak or trough times; i.e., ${\stackrel{~}{t}}^{peak}$ or ${\stackrel{~}{t}}^{trough}$) and FWHM, ξ. Two adjacent deviations are deemed to be coherent oscillations if they are grouped in the same cluster. In situations where the period or FWHM are modulated throughout time, Gaussian clustering technique may not be able to cluster all coherent oscillations adequately. We therefore process clusters by defining period- and FWHM-trends for all coherent oscillations. If this trend can accurately predict the period and FWHM of the first deviation of an adjacent cluster, then both clusters are deemed to form a set of coherent oscillations. This is repeated for all pairs of adjacent clusters, progressively updating the set of coherent oscillations with those previously deemed incoherent at prior steps.

This procedure use the Expectation Maximization (EM) clustering algorithm (McLachlan and Peel, 2005) to cluster the period and FWHM of a potential oscillatory MPRs, and to determine the optimal number of clusters using gap statistic (Tibshirani et al., 2001). It allows for a reliable separation of oscillatory data from recording artifacts or non-oscillatory MPRs with visually different properties. It also yields a set of ${\widehat{N}}_{peak}$ peaks (${\widehat{N}}_{trough}$ troughs) occurring at times ${\widehat{t}}^{peak}$ (${\widehat{t}}^{trough}$), which together represent coherent oscillations. Having determined the properties of the individual features making up the oscillations, they can be used to quantify the properties of the oscillations.

#### Characterizing Oscillatory Parameters

In order for the algorithm to report the oscillatory properties of a signal, the MPR-associated deviations must satisfy Ñ > 2, for the set of MPR parameters describing all detected peaks. If these deviations form a coherent set of oscillations, a second set of MPR parameters characterizing this coherent oscillatory behavior is also reported. In both instances, the following parameters will be reported: the number of oscillations (*N*_{osc}), the average magnitude of the oscillations (defined as *E*^{peak} + *E*^{trough}), the average period of oscillation (*T*), the standard deviation of the periodicity (σ_{T}), the total time for which the oscillations persist (defined as ${\ell}_{osc}={\widehat{t}}_{{\widehat{N}}_{peak}}^{peak}-{\widehat{t}}_{1}^{peak}$), and the mean duty cycle parameter (given by ξ^{peak}/*T*) (Smedler and Uhlén, 2014).

#### MPR Parameter Validation: N_{osc}, E, T, *l*_{osc}, and ξ^{peak}/T

Manual estimates of N_{osc} is determined by counting the number of discernable peaks within the signal. The slope of correlation was 0.85 with an *r*^{2}-score of 0.78 (Figure 6A, left). The mean difference between manual and automated estimates was 15% with an interval of agreement ranging from −39 to 69% (Figure 6A, right). The peak magnitude of the oscillations, E, is manually estimated by the mean change between peak maxima and their subsequent trough minima, after correcting for a non-stationary baseline that is often a consequence of a concurrent TR. For most signals the non-stationary baseline can be manually estimated to be linear. However, there are a few cases where an estimate of an exponential baseline is required. The correlation between manual and automated estimates was relatively strong, with an *r*^{2}-value of 0.92 and a slope of 1.08 (Figure 6B, left). The agreement analysis on the other hand revealed relatively no bias, with a mean difference of −2% and limit of agreement ranging from −48 to 44% difference (Figure 6B, right). The periodicity is manually estimated by the average time between adjacent peaks. The linear relationship was slightly weaker with a *r*^{2}-value of 0.55 and slope of 0.65 (Figure 6C, left). The mean difference between manual and automated estimates was negligible, at only 0.3%, indicating an absence systematic bias, and the limits of agreement spanned from −54 to 54% difference (Figure 6C, right). The standard deviation of periodicity was obtained from the same set of periods used to estimate the mean period. The linear slope was 0.77 and the *r*^{2}-value was 0.77 (Figure 6D, left). Similar to periodicity, the mean difference for the standard deviation of periodicity was a negligible −0.4%, with limits of agreement ranging from −61 to 62% difference (Figure 6D, right). Oscillatory persistence is chosen to describe how long oscillations are sustained within a given recording, and is estimated as the elapsed time between the first and last discernable peaks in the MPR. The correlation between manual and automated estimates was supported by a *r*^{2}-value of 0.79 and slope of 0.85 (Figure 6E, left). The mean difference between paired estimates was only −5.3% with a limit of agreement between −47 and 37% difference (Figure 6E, right). Finally, the duty cycle is manually estimated by the ratio between ξ^{peak} and T. ξ^{peak} is manually determined by the mean FWHM of individual oscillatory peaks and the same T value obtained above is used to calculate x^{peak} /T. The linear relationship between manual and automated estimated of ξ^{peak} /T was decidedly weak with a slope of 0.51 and *r*^{2} of 0.26 (Figure 6F, left). The Bland Altman analysis, however, suggests that there was a systematic bias that could explain the poorer correlation results. The mean difference between manual and automated estimated was −19% with a limit of agreement between −69 and 31% (Figure 6F, right).

**Figure 6. Validation of MPR parameter measurements**. Parameters were estimated manually and then compared with automated-estimated values using correlation analysis (left: solid line, linear regression line; broken line, line of equality) and Bland-Altman analysis (right: solid line, mean percent difference, $\stackrel{\u0304}{\Xi}$; broken line, limits of agreement, $\stackrel{\u0304}{\Xi}$ ±1.96 ${\sigma}_{\stackrel{\u0304}{\Xi}}$). **(A)** Number of oscillatory periods N_{osc}. **(B)** Magnitude of oscillatory peaks, E. **(C)** Periodicity, T. **(D)** Standard Deviation of Period, σ_{T}. **(E)** Oscillatory persistence, *l*_{osc}. **(F)** Duty cycle, ξ^{peak}/T. Arbitrary Units, A.U. Insets: Visual representations of parameters measured.

The MPR parameter validation described above focuses on the set of parameter estimates describing all the peaks in the MPR, rather than the subset of coherent oscillations. This is because manual detection of each peak in the MPR is less subjective than detecting only the coherent peaks in the MPR. Since the algorithm sub-selects the coherent oscillatory peaks from the initial set of identified deviations, the performance reported for the characterization of all peaks extends to the subset of coherent oscillations. Furthermore, as expected, the standard deviation of the periodicity is consistently lower for coherent oscillations when compare to the σ_{T} reported for all peaks in the same MPR (i.e., more regular periodicities result in lower standard deviations).

To ensure confidence in the reported MPR parameters, users of this algorithm are urged to visually verify the quality of the signal fittings to determine whether the algorithm is characterizing their peaks of interest, as these may not always coincide with the most prevalent oscillatory component of the signal (Thurley et al., 2014). Furthermore, N_{osc} reported for all peaks and coherent peaks can be compared to be aware of how many peaks were omitted during the clustering step. Collectively, the information reported for MPRs is sufficient for the informed analysis of a diverse selection of MPRs, including those that exhibit coherent oscillations and those that do not.

### Application to Pathophysiology

In the context of bone physiology, the deleterious consequences of disrupting extracellular nucleotide-mediated cross talk have been highlighted by the emergence of P2 receptor knockout mouse models (Lenertz et al., 2015). P2 receptors are particularly sensitive to changes in the extracellular milieu. Consequently, P2 receptor pathophysiology is often coupled to events that influence the extracellular composition, thereby compromising processes regulated by the P2 receptor network. In particular, changes in extracellular pH alter P2 receptor function (King et al., 1997; Gerevich et al., 2007; Wildman, 2009; Langfelder et al., 2015). Such conditions arise from pathological acidosis that is commonly caused by systemic acid-base disturbances, such as metabolic or respiratory acidosis (Krieger et al., 2004; Miller, 2012; Berend et al., 2014). More localized acidifications can also be associated with tumors (Martin and Jain, 1994; Kato et al., 2013). Since the skeleton is a common metastatic site for cancer, and participates in systemic buffering of protons, the effect of acidosis on the skeletal system is of particular interest. On the cellular level, acidosis promotes the activation of osteoclasts, resulting in elevated bone resorption which manifests itself in osteoporotic phenotype (Bushinsky and Frick, 2000; Krieger et al., 2003; Ahn et al., 2012; Gasser et al., 2014). However, it remains unclear whether the P2 receptor network plays a direct role in this cascade of events. The most immediate influence of acidosis on the P2-receptor network can be studied at the level of the [Ca^{2+}]_{i} response evoked immediately upon application of a purinergic agonist.

We investigated the effect of acidosis on ATP-mediated [Ca^{2+}]_{i} responses in bone-marrow derived osteoclast precursors to demonstrate the applicability of the developed algorithm. The application of ATP (100 nM to 10 mM) to the fura2-loaded osteoclast precursors evoked a [Ca^{2+}]_{i} TR in a dose-dependent manner in control and acidosis conditions (Figures 7A,B). The response amplitudes under acidic conditions were virtually indistinguishable from the control for ATP concentrations up to 10 μM. However, above this threshold concentration, the amplitude of the control responses continued to increase with rising concentrations of ATP, while [Ca^{2+}]_{i} responses under acidic conditions plateaued at 10 μM (Figure 7C). With respect to the AUC of the [Ca^{2+}]_{i} responses, the observed differences between the two conditions were more gradual with a diverging trend beginning as low as 1 μM ATP and becoming more prominent at high ATP (Figure 7D). Finally, acidosis was found to have no significant effect on the periodicity of the oscillatory responses (Figure 7E).

**Figure 7. Algorithm application in characterization of pathological states**. ATP (100 nM–10 mM) was applied to Fura-2 loaded osteoclast precursors, under control (pH 7.6) and acidosis (pH 7.0) conditions, and [Ca^{2+}]_{i} responses were recorded. Algorithm was used to obtain estimates for amplitude, AUC and periodicity. **(A)** Representative [Ca^{2+}]_{i} response traces for 100 μM and 1 mM ATP under control conditions. **(B)** Representative [Ca^{2+}]_{i} response traces for 100 μM and 1 mM ATP under acidosis conditions. **(C)** Amplitude dose-response curves. **(D)** AUC dose-response curves. **(E)** Period dose-response curves. For **(C–E)**, data are mean ± S.E.M. The effect of ATP under control conditions was examined using one-way ANOVA. The effect of acidosis on ATP-mediated responses was examined using two-way ANOVA. The Bonferroni test was used for post hoc multi-comparison analysis; ***p* < 0.01; ****p* < 0.001 indicate significant difference compared to the response to the lowest ATP concentration; ^{#}*p* < 0.05 indicates significant difference between responses to the same [ATP] in control and acidosis conditions.

These findings support that acidosis, while having no effect on ATP-mediated [Ca^{2+}]_{i} responses at lower ATP concentration, significantly attenuates the magnitude of [Ca^{2+}]_{i} transients responding to higher ATP concentrations (>10 μM ATP). Within the limited scope of this study that is focused on the development of a data analysis algorithm, we can only hypothesize on the mechanism by which these differences arise. One possibility is that the rise in extracellular [H^{+}] has a significant influence on the electro-chemical gradient across the cellular membrane, which may consequently alter the extent of calcium flux across certain ionotropic P2X receptors. Since the oscillations are commonly driven by inositol triphosphate-mediated release of calcium from internal calcium stores (i.e., isolated from extracellular [H^{+}]), it may explain why the oscillatory behavior is not affected by acidosis. Alternatively, there may exist a subset of P2 receptors that are sensitive to fluctuations in extracellular [H^{+}], while P2-receptors involved in oscillatory behavior and/or responses to lower ATP concentrations (≤10 μM) are resilient to such changes. Regardless of the underlying mechanism, these results highlight that the P2 receptor network can be differentially modulated by extracellular pH.

## Conclusions

This paper presents an autonomous signal-processing algorithm capable of robustly removing signal-contaminating noise and delineating the various components seen in a calcium response, including non-stationary drift, TRs, and MPRs (possibly caused by flickers, puffs, and sparks) sampled with at least twice the Nyquist frequency. By fitting piece-wise defined model functions to data, the algorithm also extracts estimates for the parameters that are relevant to the characterization of cellular transient dynamics. Any time-series recordings can be used as an input for the algorithm, provided that they resemble a single or multi-peak transient response. As demonstrated in the validation process, manual estimation of certain parameters has an inherent degree of subjectivity and measurement error associated with it. In particular, the manual evaluations of AUC values, decay constant and the time of onset, as well as most of the MPR-parameter values, were found to rely on subjective estimates and thus lacked true accuracy and consistency. Because of such limitation in the validation method, manual-estimates are to be recognized as representative estimates, rather than accurate values for these parameters. Consequently, validation method applied here should be considered as a comparison against imperfect estimates.

Nevertheless, our analysis of the automated method has verified that the algorithm performs within acceptable margins of agreement when compared to manual analysis. Regarding the response detection capabilities, the algorithm behaves conservatively compared to manual assessments, especially when presented with low-magnitude TRs or ambiguous response signals. Most importantly, our algorithm has been validated against experimental [Ca^{2+}]_{i} recording data, rather than simulated data, ensuring that the method is capable of handling variations in drift and noise that realistically reflect signal contaminations of experimental data acquisition. We have demonstrated that this automated methodology is effective in analyzing empirical data, providing quantitative insights about them and identifying differences between them.

A particularly unique feature of this algorithm is its capacity to characterize the magnitude and temporal characteristics of MPRs exhibiting stochastic and deterministic behavior. It is well established that a diverse amount of biochemical processes can be amplitude- and/or frequency-modulated (Adachi et al., 1999; Micali et al., 2015). To analyse such oscillatory data, the fast Fourier transform (FFT) is commonly used, which allows for the conversion of a signal from its time domain, into the frequency domain. Unfortunately, the variance in the frequency domain is proportional to the number of repetitive components in the time-domain. Therefore, if the oscillatory signals present few repetitive components then reliable resolution of the true periodicity of the signal is unachievable. To circumvent the limitations inherent to FFT, we apply the MATLAB “findpeaks()” function to identify peaks of interest. To isolate underlying coherent oscillations that are often present, we applied a clustering method. This is based on the principle of clustering deviations from baseline according to their temporal offset and respective FWHM. The advantage of this approach is that it allows for the reliable detection of periodic peaks, even in the presence of stochastic discharges, as is often the case in experimental recordings. Secondly, comparison of the set of MPR parameters for all peaks and subset of coherent peaks allows users to quantify the extent of stochastic activity within MPRs. Alternatively, the relationship between mean and standard deviation of periodicity in a MPR has been previously used to reveal the contribution of stochastic processes to the periodicity (Thurley et al., 2014). We anticipate this methodology will contribute to the comprehensive analysis of diverse MPRs.

Calcium signaling is by no means unique to the P2-receptor network, but rather represents the most ubiquitous and versatile messenger found in biological systems. All kinds of extracellular signals exploit calcium as a secondary messenger, including P2 agonists (i.e., ATP, ADP, UTP, and UDP), endothelin-1, oxotremorine-M, norephinephrine, thrombin, PDGF, bombensin (Balla et al., 1991; Palmer, 1994; Burnstock, 2004). The universal involvement of calcium ranges from basic physiological processes such as muscle contraction, neuronal discharge and pancreatic secretion, to early development events including mammalian egg fertilization and embryonic pattern formation (Berridge et al., 2000). Calcium signaling is also known to be impaired in various pathological states, as suggested for metabolic acidosis in this study, chronic renal failure (Massry et al., 1995), Alzheimer's (Brawek et al., 2014), Diabetes (Chen et al., 2015), and zinc deficiency (O'Dell and Browning, 2013). However, despite all that we know about calcium's role in biological processes, there remains ongoing debate on how calcium signals robustly encode information while still exhibiting a large degree of heterogeneity within and between various cellular populations. Many theories have been proposed to establish how information can be encoded. Some of these involve encoding information on the basis of calcium binding cooperativity (Larsen et al., 2004), amplitude and frequency modulation (De Pitta et al., 2009), changes in spike time variation (Thurley et al., 2014), and signal integration (Hannanta-anan and Chow, 2016). In order to reconcile these theories and establish a universal syntax for calcium-encoded information, tools such as this algorithm will aid in the large-scale analysis of experimental data sets required for the validation of mathematical models.

The consideration of signaling nuances that are specifically found in physiological signals, but may or may not be present in non-biological signals, was a critical step in the development of this algorithm. As demonstrated in this study, physiological signals were decomposed into their elementary components and mathematically generalized to enable for the computational reconstruction of a diverse range of signature forms. In doing so, we were able to provide a foundation for further modeling of the nonlinear multi-parametric physiological signals. This study demonstrates that the accurate description of complex physiological signals is non-trivial, but rather an extensive mathematical undertaking. Therefore, we believe that, beyond serving the purpose of a signal-processing tool, this algorithm will also contribute to future efforts to modeling physiological signals.

In summary, we have detailed an open-source MATLAB algorithm intended to facilitate the analysis of time-series recordings. With minimal user-input required, this tool dramatically decreases analysis time and ensures consistency in parameter characterization of complex physiological signals. This algorithm is capable of handling noise and drift and robustly characterizes the magnitude and kinetics of dynamic processes, outputting the amplitude, time of onset (t_{onset}), activation time (t_{10–90%}), full-width half-max (FWHM), AUC, and decay constant (τ_{decay}). In the presence of MPR, six additional parameters are characterized which include number of oscillations (N_{osc}), magnitude of oscillatory peaks (E), periodicity (T), standard deviation of periodicity (σ_{T}), oscillatory persistence (l_{osc}), and the duty cycle (ξ^{peak}/T). This algorithm is not limited to any specific data-type, but [Ca^{2+}]_{i} recordings represent an obvious application. In addition to calcium imaging, other imaging modalities such as adapted fluorescence resonance energy transfer (FRET) biosensors, real-time bioluminescence and voltage and current measurements can generate time-series data for which characterization of signal magnitude and kinetics can provide valuable information. As data acquisition becomes more efficient and data sets become increasingly complex, automated analysis will serve as an essential tool for conducting basic research and clinical screening.

## Author Contributions

Study conception and design: LM, NM, SK, AK. Algorithm development: LM, AK. Acquisition of data: NM. Analysis and interpretation of data: LM, NM. Drafting of Manuscript: LM, NM, SK, AK. All authors contributed to the critical revision of manuscript and approved the final version to be published.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Acknowledgments

This research work was supported by the Natural Sciences and Engineering Research Council (NSERC) grants to AK and SK. The authors would like to acknowledge the partial financial support provided to LM by Dr. Michael Mackey through his NSERC fund.

## Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/article/10.3389/fphys.2016.00525/full#supplementary-material

## Abbreviations

ADP, Adenosine diphosphate; ATP, Adenosine triphosphate; AUC, Area under curve; [Ca^{2+}]_{i}, Cytosolic free calcium concentration; D(t), Global drift model; E, Mean peak magnitude of MPR; FWHM, Full-width half-max; F(t), Measured fluorescent signal; *l*_{osc}, Oscillatory persistence; MPR, Multi-peak response; N_{osc}, Number of peaks in MPR; ξ^{peak}, Width of oscillatory peak; P2, Purinergic; σ_{T}, Standard deviation of periodicity; T, Period of oscillation; t_{10–90%}, Activation time; τ_{decay}, Decay constant; t_{onset}, Time of response onset; TR, Transient response; TV, Total-variational.

## References

Abu Khamidakh, A. E., Juuti-Uusitalo, K., Larsson, K., Skottman, H., and Hyttinen, J. (2013). Intercellular Ca(2+) wave propagation in human retinal pigment epithelium cells induced by mechanical stimulation. *Exp. Eye Res.* 108, 129–139. doi: 10.1016/j.exer.2013.01.009

Adachi, Y., Kindzelskii, A. L., Ohno, N., Yadomae, T., and Petty, H. R. (1999). Amplitude and frequency modulation of metabolic signals in leukocytes: synergistic role of IFN-γ in IL-6- and IL-2-mediated cell activation. *J. Immunol.* 163, 4367–4374.

Ahn, H., Kim, J. M., Lee, K., Kim, H., and Jeong, D. (2012). Extracellular acidosis accelerates bone resorption by enhancing osteoclast survival, adhesion, and migration. *Biochem. Biophys. Res. Commun.* 418, 144–148. doi: 10.1016/j.bbrc.2011.12.149

Appleby, P. A., Shabir, S., Southgate, J., and Walker, D. (2015). Sources of variability in cytosolic calcium transients triggered by stimulation of homogeneous uro-epithelial cell monolayers. *J. R. Soc. Interface* 12:20141403. doi: 10.1098/rsif.2014.1403

Balla, T., Sim, S. S., Iida, T., Choi, K. Y., Catt, K. J., and Rhee, S. G. (1991). Agonist-induced calcium signaling is impaired in fibroblasts overproducing inositol 1,3,4,5-tetrakisphosphate. *J. Biol. Chem.* 266, 24719–24726.

Berend, K., de Vries, A. P., and Gans, R. O. (2014). Physiological approach to assessment of acid-base disturbances. *N. Engl. J. Med.* 371, 1434–1445. doi: 10.1056/NEJMra1003327

Berridge, M. J., Lipp, P., and Bootman, M. D. (2000). The versatility and universality of calcium signalling. *Nat. Rev. Mol. Cell Biol.* 1, 11–21. doi: 10.1038/35036035

Bland, J. M., and Altman, D. G. (1986). Statistical methods for assessing agreement between two methods of clinical measurement. *Lancet* 1, 307–310. doi: 10.1016/S0140-6736(86)90837-8

Boraschi-Diaz, I., and Komarova, S. V. (2016). The protocol for the isolation and cryopreservation of osteoclast precursors from mouse bone marrow and spleen. *Cytotechnology* 68, 105–114. doi: 10.1007/s10616-014-9759-3

Brawek, B., Schwendele, B., Riester, K., Kohsaka, S., Lerdkrai, C., Liang, Y., et al. (2014). Impairment of *in vivo* calcium signaling in amyloid plaque-associated microglia. *Acta Neuropathol.* 127, 495–505. doi: 10.1007/s00401-013-1242-2

Bray, M. A., Geisse, N. A., and Parker, K. K. (2007). Multidimensional detection and analysis of Ca2+ sparks in cardiac myocytes. *Biophys. J.* 92, 4433–4443. doi: 10.1529/biophysj.106.089359

Burnstock, G. (2004). Introduction: P2 receptors. *Curr. Top. Med. Chem.* 4, 793–803. doi: 10.2174/1568026043451014

Burnstock, G., and Verkhratsky, A. (2009). Evolutionary origins of the purinergic signalling system. *Acta Physiol. (Oxf.)* 195, 415–447. doi: 10.1111/j.1748-1716.2009.01957.x

Bushinsky, D. A., and Frick, K. K. (2000). The effects of acid on bone. *Curr. Opin. Nephrol. Hypertens.* 9, 369–379. doi: 10.1097/00041552-200007000-00008

Cao, D., Lin, G., Westphale, E. M., Beyer, E. C., and Steinberg, T. H. (1997). Mechanisms for the coordination of intercellular calcium signaling in insulin-secreting cells. *J. Cell Sci.* 110(Pt 4), 497–504.

Chambolle, A., Caselles, V., Novaga, M., Cremers, D., and Pock, T. (2009). *An Introduction to Total Variation for Image Analysis*. Available online at: https://www.semanticscholar.org/paper/An-Introduction-to-Total-Variation-for-Image-Chambolle-Caselles/19d5136714aa8aa84129ca2b5f0fefca1ea6ce8d#citingPapers

Chartrand, R. (2007). “Nonconvex regularization for shape preservation,” in *Paper Presented at the 2007 IEEE International Conference on Image Processing* (San Antonio, TX). doi: 10.1109/ICIP.2007.4378949

Chartrand, R. (2011). Numerical differentiation of noisy, nonsmooth data. *ISRN Appl. Math.* 2011, 1–11. doi: 10.5402/2011/164564

Chartrand, R., and Staneva, V. (2008). A quasi-newton method for total variation regulatrization of images corrupted by non-gaussian noise. *IET Image Process.* 2, 1–6. doi: 10.1049/iet-ipr:20080017

Chartrand, R., and Wohlberg, B. (2010). “Total-variation regularization with bound constraints,” in *Proceedings of IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP)* (Dallas, TX), 766–769.

Chen, H., Kold-Petersen, H., Laher, I., Simonsen, U., and Aalkjaer, C. (2015). Impaired endothelial calcium signaling is responsible for the defective dilation of mesenteric resistance arteries from db/db mice to acetylcholine. *Eur. J. Pharmacol.* 767, 17–23. doi: 10.1016/j.ejphar.2015.09.043

Churchill, G. C., Atkinson, M. M., and Louis, C. F. (1996). Mechanical stimulation initiates cell-to-cell calcium signaling in ovine lens epithelial cells. *J. Cell Sci.* 109, 355–365.

De Pittà, M., Volman, V., Levine, H., and Ben-Jacob, E. (2009). Multimodal encoding in a simplified model of intracellular calcium signaling. *Cogn. Process.* 10(Suppl. 1), S55–S70. doi: 10.1007/s10339-008-0242-y

Dickinson, G. D., and Parker, I. (2013). Temperature dependence of IP3-mediated local and global Ca^{2+} signals. *Biophys. J.* 104, 386–395. doi: 10.1016/j.bpj.2012.12.024

Dupont, G., and Combettes, L. (2009). What can we learn from the irregularity of Ca^{2+} oscillations? *Chaos* 19, 037112. doi: 10.1063/1.3160569

Dupont, G., Combettes, L., Bird, G. S., and Putney, J. W. (2011). Calcium oscillations. *Cold Spring Harb. Perspect. Biol.* 3:a004226. doi: 10.1101/cshperspect.a004226

Ellefsen, K. L., Settle, B., Parker, I., and Smith, I. F. (2014). An algorithm for automated detection, localization and measurement of local calcium signals from camera-based imaging. *Cell Calcium* 56, 147–156. doi: 10.1016/j.ceca.2014.06.003

Frame, M. K., and de Feijter, A. W. (1997). Propagation of mechanically induced intercellular calcium waves via gap junctions and ATP receptors in rat liver epithelial cells. *Exp. Cell Res.* 230, 197–207. doi: 10.1006/excr.1996.3409

Francis, M., Waldrup, J. R., Qian, X., Solodushko, V., Meriwether, J., and Taylor, M. S. (2016). Functional tuning of intrinsic endothelial Ca2+ dynamics in swine coronary arteries. *Circ. Res.* 118, 1078–1090. doi: 10.1161/CIRCRESAHA.115.308141

Fritzsche, M., Fernandes, R. A., Colin-York, H., Santos, A. M., Lee, S. F., Lagerholm, B. C., et al. (2015). CalQuo: automated, simultaneous single-cell and population-level quantification of global intracellular Ca2+ responses. *Sci. Rep.* 5:16487. doi: 10.1038/srep16487

Gallagher, J. A., and Buckley, K. A. (2002). Expression and function of P2 receptors in bone. *J. Musculoskelet. Neuronal Interact.* 2, 432–439.

Gasser, J. A., Hulter, H. N., Imboden, P., and Krapf, R. (2014). Effect of chronic metabolic acidosis on bone density and bone architecture *in vivo* in rats. *Am. J. Physiol. Renal Physiol.* 306, F517–F524. doi: 10.1152/ajprenal.00494.2013

Gerevich, Z., Zadori, Z. S., Köles, L., Kopp, L., Milius, D., Wirkner, K., et al. (2007). Dual effect of acid pH on purinergic P2X3 receptors depends on the histidine 206 residue. *J. Biol. Chem.* 282, 33949–33957. doi: 10.1074/jbc.M705840200

Hannanta-anan, P., and Chow, B. Y. (2016). Optogenetic control of calcium oscillation waveform defines NFAT as an integrator of calcium load. *Cell Syst.* 2, 283–288. doi: 10.1016/j.cels.2016.03.010

Hansen, M., Boitano, S., Dirksen, E. R., and Sanderson, M. J. (1993). Intercellular calcium signaling induced by extracellular adenosine 5′-triphosphate and mechanical stimulation in airway epithelial cells. *J. Cell Sci.* 106(Pt 4), 995–1004.

Isakson, B. E., Evans, W. H., and Boitano, S. (2001). Intercellular Ca2+ signaling in alveolar epithelial cells through gap junctions and by extracellular ATP. *Am. J. Physiol. Lung Cell. Mol. Physiol.* 280, L221–L228.

James, L. R., Andrews, S., Walker, S., de Sousa, P. R., Ray, A., Russell, N. A., et al. (2011). High-throughput analysis of calcium signalling kinetics in astrocytes stimulated with different neurotransmitters. *PLoS ONE* 6:e26889. doi: 10.1371/journal.pone.0026889

Jorgensen, N. R., Geist, S. T., Civitelli, R., and Steinberg, T. H. (1997). ATP- and gap junction-dependent intercellular calcium signaling in osteoblastic cells. *J. Cell Biol.* 139, 497–506. doi: 10.1083/jcb.139.2.497

Jørgensen, N. R., Teilmann, S. C., Henriksen, Z., Civitelli, R., Sørensen, O. H., and Steinberg, T. H. (2003). Activation of L-type calcium channels is required for gap junction-mediated intercellular calcium signaling in osteoblastic cells. *J. Biol. Chem.* 278, 4082–4086. doi: 10.1074/jbc.M205880200

Juhola, M., Penttinen, K., Joutsijoki, H., Varpa, K., Saarikoski, J., Rasku, J., et al. (2015). Signal analysis and classification methods for the calcium transient data of stem cell-derived cardiomyocytes. *Comput. Biol. Med.* 61, 1–7. doi: 10.1016/j.compbiomed.2015.03.016

Kaczmarek-Hájek, K., Lörinczi, E., Hausmann, R., and Nicke, A. (2012). Molecular and functional properties of P2X receptors–recent progress and persisting challenges. *Purinergic Signal.* 8, 375–417. doi: 10.1007/s11302-012-9314-7

Kato, Y., Ozawa, S., Miyamoto, C., Maehata, Y., Suzuki, A., Maeda, T., et al. (2013). Acidic extracellular microenvironment and cancer. *Cancer Cell Int.* 13:89. doi: 10.1186/1475-2867-13-89

King, B. F., Wildman, S. S., Ziganshina, L. E., Pintor, J., and Burnstock, G. (1997). Effects of extracellular pH on agonism and antagonism at a recombinant P2X2 receptor. *Br. J. Pharmacol.* 121, 1445–1453. doi: 10.1038/sj.bjp.0701286

Koshimizu, T., Koshimizu, M., and Stojilkovic, S. S. (1999). Contributions of the C-terminal domain to the control of P2X receptor desensitization. *J. Biol. Chem.* 274, 37651–37657. doi: 10.1074/jbc.274.53.37651

Krieger, N. S., Bushinsky, D. A., and Frick, K. K. (2003). Cellular mechanisms of bone resorption induced by metabolic acidosis. *Semin. Dial.* 16, 463–466. doi: 10.1046/j.1525-139X.2003.16100.x

Krieger, N. S., Frick, K. K., and Bushinsky, D. A. (2004). Mechanism of acid-induced bone resorption. *Curr. Opin. Nephrol. Hypertens.* 13, 423–436. doi: 10.1097/01.mnh.0000133975.32559.6b

Langfelder, A., Okonji, E., Deca, D., Wei, W. C., and Glitsch, M. D. (2015). Extracellular acidosis impairs P2Y receptor-mediated Ca(2+) signalling and migration of microglia. *Cell Calcium* 57, 247–256. doi: 10.1016/j.ceca.2015.01.004

Larsen, A. Z., Olsen, L. F., and Kummer, U. (2004). On the encoding and decoding of calcium signals in hepatocytes. *Biophys. Chem.* 107, 83–99. doi: 10.1016/j.bpc.2003.08.010

Lenertz, L. Y., Baughman, C. J., Waldschmidt, N. V., Thaler, R., and van Wijnen, A. J. (2015). Control of bone development by P2X and P2Y receptors expressed in mesenchymal and hematopoietic cells. *Gene* 570, 1–7. doi: 10.1016/j.gene.2015.06.031

Li, F., Shen, C., Fan, J., and Shen, C. (2007). Image restoration combining a total variational filter and a fourth-order filter. *J. Vis. Commun. Image Represent.* 18, 322–330. doi: 10.1016/j.jvcir.2007.04.005

Lock, J. T., Ellefsen, K. L., Settle, B., Parker, I., and Smith, I. F. (2015). Imaging local Ca^{2+} signals in cultured mammalian cells. *J. Vis. Exp.* e52516. doi: 10.3791/52516

Marquardt, D. W. (1963). An algorithm for least-squares estimation of nonlinear parameters. *J. Soc. Indust. Appl. Math.* 11, 431–441. doi: 10.1137/0111030

Martin, G. R., and Jain, R. K. (1994). Noninvasive measurement of interstitial pH profiles in normal and neoplastic tissue using fluorescence ratio imaging microscopy. *Cancer Res.* 54, 5670–5674.

Massry, S. G., Klin, M., Ni, Z., Tian, J., Kedes, L., and Smogorzewski, M. (1995). Impaired agonist-induced calcium signaling in hepatocytes from chronic renal failure rats. *Kidney Int.* 48, 1324–1331. doi: 10.1038/ki.1995.417

McLachlan, G., and Peel, D. (2005). “ML fitting of mixture models,” in *Finite Mixture Models* (Hoboken, NJ: John Wiley and Sons, Inc.), 40–80. doi: 10.1002/0471721182.ch2

Micali, G., Aquino, G., Richards, D. M., and Endres, R. G. (2015). Accurate encoding and decoding by single cells: amplitude versus frequency modulation. *PLoS Comput. Biol.* 11:e1004222. doi: 10.1371/journal.pcbi.1004222

Miller, P. D. (2012). Unrecognized and unappreciated secondary causes of osteoporosis. *Endocrinol. Metab. Clin. North Am.* 41, 613–628. doi: 10.1016/j.ecl.2012.05.005

Nobile, M., Monaldi, I., Alloisio, S., Cugnoli, C., and Ferroni, S. (2003). ATP-induced, sustained calcium signalling in cultured rat cortical astrocytes: evidence for a non-capacitative, P2X7-like-mediated calcium entry. *FEBS Lett.* 538, 71–76. doi: 10.1016/S0014-5793(03)00129-7

O'Dell, B. L., and Browning, J. D. (2013). Impaired calcium entry into cells is associated with pathological signs of zinc deficiency. *Adv. Nutr.* 4, 287–293. doi: 10.3945/an.112.003624

Oh, S., Woo, H., Yun, S., and Kang, M. (2013). Non-convex hybrid total variation for image denoising. *J. Vis. Commun. Image Represent.* 24, 332–344. doi: 10.1016/j.jvcir.2013.01.010

Palmer, R. K., Yule, D. I., McEwen, E. L., Williams, J. A., and Fisher, S. K. (1994). Agonist-specific calcium signaling and phosphoinositide hydrolysis in human SK-N-MCIXC neuroepithelioma cells. *J. Neurochem.* 63, 2099–2107. doi: 10.1046/j.1471-4159.1994.63062099.x

Patel, T. P., Man, K., Firestein, B. L., and Meaney, D. F. (2015). Automated quantification of neuronal networks and single-cell calcium dynamics using calcium imaging. *J. Neurosci. Methods* 243, 26–38. doi: 10.1016/j.jneumeth.2015.01.020

Picht, E., Zima, A. V., Blatter, L. A., and Bers, D. M. (2007). SparkMaster: automated calcium spark analysis with ImageJ. *Am. J. Physiol. Cell Physiol.* 293, C1073–C1081. doi: 10.1152/ajpcell.00586.2006

Rast, G., Weber, J., Disch, C., Schuck, E., Ittrich, C., and Guth, B. D. (2015). An integrated platform for simultaneous multi-well field potential recording and Fura-2-based calcium transient ratiometry in human induced pluripotent stem cell (hiPSC)-derived cardiomyocytes. *J. Pharmacol. Toxicol. Methods* 75, 91–100. doi: 10.1016/j.vascn.2015.04.005

Rice, S. O. (1973). Efficient evaluation of integrals of analytic functions by the trapezoidal rule. *Bell Syst. Tech. J.* 52, 707–722. doi: 10.1002/j.1538-7305.1973.tb01986.x

Romanello, M., and D'Andrea, P. (2001). Dual mechanism of intercellular communication in HOBIT osteoblastic cells: a role for gap-junctional hemichannels. *J. Bone Miner. Res.* 16, 1465–1476. doi: 10.1359/jbmr.2001.16.8.1465

Ruffinatti, F. A., Lovisolo, D., Distasi, C., Ariano, P., Erriquez, J., and Ferraro, M. (2011). Calcium signals: analysis in time and frequency domains. *J. Neurosci. Methods* 199, 310–320. doi: 10.1016/j.jneumeth.2011.05.009

Shabir, S., and Southgate, J. (2008). Calcium signalling in wound-responsive normal human urothelial cell monolayers. *Cell Calcium* 44, 453–464. doi: 10.1016/j.ceca.2008.02.008

Skupin, A., Kettenmann, H., Winkler, U., Wartenberg, M., Sauer, H., Tovey, S. C., et al. (2008). How does intracellular Ca(2+) oscillate: by chance or by the clock? *Biophys. J.* 94, 2404–2411. doi: 10.1529/biophysj.107.119495

Smedler, E., and Uhlen, P. (2014). Frequency decoding of calcium oscillations. *Biochim. Biophys. Acta* 1840, 964–969. doi: 10.1016/j.bbagen.2013.11.015

Smith, I. F., Wiltgen, S. M., and Parker, I. (2009). Localization of puff sites adjacent to the plasma membrane: functional and spatial characterization of Ca2+ signaling in SH-SY5Y cells utilizing membrane-permeant caged IP3. *Cell Calcium* 45, 65–76. doi: 10.1016/j.ceca.2008.06.001

Steele, E. M., and Steele, D. S. (2014). Automated detection and analysis of Ca(2+) sparks in x-y image stacks using a thresholding algorithm implemented within the open-source image analysis platform ImageJ. *Biophys. J.* 106, 566–576. doi: 10.1016/j.bpj.2013.12.040

Stoehr, A., Neuber, C., Baldauf, C., Vollert, I., Friedrich, F. W., Flenner, F., et al. (2014). Automated analysis of contractile force and Ca2+ transients in engineered heart tissue. *Am. J. Physiol. Heart Circ. Physiol.* 306, H1353–H1363. doi: 10.1152/ajpheart.00705.2013

Sun, F., Berry, D. J., Leong, D. A., and Veldhuis, J. D. (1997). Recruitment of individually (all-or-none) responding cells, rather than amplitude enhancement, is the single-cell mechanism subserving the dose-responsive activation of intracellular calcium second messenger signaling by the human luteinizing-hormone receptor. *Endocrine* 7, 219–226. doi: 10.1007/BF02778144

Moré, J. J. (1978). *The Levenberg-Marquardt Algorithm: Implementation and Theory.* Heidelberg: Springer.

Thurley, K., Tovey, S. C., Moenke, G., Prince, V. L., Meena, A., Thomas, A. P., et al. (2014). Reliable encoding of stimulus intensities within random sequences of intracellular Ca2+ spikes. *Sci. Signal.* 7, ra59. doi: 10.1126/scisignal.2005237

Tibshirani, R., Walther, G., and Hastie, T. (2001). Estimating the number of clusters in a data set via the gap statistic. *J. R. Stat. Soc. B* 63, 411–423. doi: 10.1111/1467-9868.00293

Traub, J. F. (1964). On lagrange-hermite interpolation. *J. Soc. Indust. Appl. Math.* 12, 886–891. doi: 10.1137/0112076

Vogel, C. (2002). *Computational Methods for Inverse Problem*. Philadelphia, PA: Society for Industrial and Applied Mathematics.

Volonté, C., Amadio, S., D'Ambrosi, N., Colpi, M., and Burnstock, G. (2006). P2 receptor web: complexity and fine-tuning. *Pharmacol. Ther.* 112, 264–280. doi: 10.1016/j.pharmthera.2005.04.012

von Kügelgen, I., and Hoffmann, K. (2016). Pharmacology and structure of P2Y receptors. *Neuropharmacology* 104, 50–61. doi: 10.1016/j.neuropharm.2015.10.030

Wildman, S. S. P. (2009). Changes in tubular fluid pH and nucleotide concentration alter P2 receptor-mediated regulation of renal ENaC. *FASEB J.* 23(1 Suppl.), 602–605.

Wong, L. C., Lu, B., Tan, K. W., and Fivaz, M. (2010). Fully-automated image processing software to analyze calcium traces in populations of single cells. *Cell Calcium* 48, 270–274. doi: 10.1016/j.ceca.2010.09.008

Xing, S., Grol, M. W., Grutter, P. H., Dixon, S. J., and Komarova, S. V. (2016). Modeling interactions among individual P2 receptors to explain complex response patterns over a wide range of ATP concentrations. *Front. Physiol.* 7:294. doi: 10.3389/fphys.2016.00294

Yan, Z., Khadra, A., Li, S., Tomic, M., Sherman, A., and Stojilkovic, S. S. (2010). Experimental characterization and mathematical modeling of P2X7 receptor channel gating. *J. Neurosci.* 30, 14213–14224. doi: 10.1523/JNEUROSCI.2390-10.2010

Keywords: algorithm, calcium imaging, kinetics, osteoclast pathophysiology, parameter characterization, purinergic/P2 receptors, real-time imaging

Citation: Mackay L, Mikolajewicz N, Komarova SV and Khadra A (2016) Systematic Characterization of Dynamic Parameters of Intracellular Calcium Signals. *Front. Physiol*. 7:525. doi: 10.3389/fphys.2016.00525

Received: 31 August 2016; Accepted: 24 October 2016;

Published: 10 November 2016.

Edited by:

Krasimira Tsaneva-Atanasova, University of Exeter, UKReviewed by:

Silvina Ponce Dawson, University of Buenos Aires, ArgentinaJames Sneyd, University of Auckland, New Zealand

Martin Falcke, Max Delbrück Center for Molecular Medicine (HZ), Germany

Copyright © 2016 Mackay, Mikolajewicz, Komarova and Khadra. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Anmar Khadra, anmar.khadra@mcgill.ca

^{†}These authors have contributed equally to this work.