Systematic Characterization of Dynamic Parameters of Intracellular Calcium Signals

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.
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, 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. consisting of eight subtypes (P2Y [1][2]4,6,[11][12][13][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 signalprocessing 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 ½ 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 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.  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 ≫characterizeDocument("filename.xlsx") The MATLAB code required to execute this command along with examples of data and scripts are provided in Supplmentary Materials (Data Sheet 1).

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-highspeed 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 is the average value of the estimated parameter,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 95% interval of agreement =¯ ± 1.96σ¯ where¯ is the mean percent difference over all recordings and σ¯ is the standard deviation of¯ .
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,û (t), 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 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 nonlinear, 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 = {i = 1, 2, · · · , N |ψ n (t i ) = ∞ } (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 i = 1 − | i |/max i | i | in order to provide an upper bound on the noise of the signal, given by σ n + = j j g j j g j .
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 σ n + . On the other hand, we can also estimate a lower bound on the noise using where ζ n i = F i − (Au n ) 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 0 i = 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 σ n − being smaller than σ 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 = u ′ i /max i u ′ 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 δ n i = 1 − h n i n i . 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 datafidelity 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 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 √ φ, 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 τ n rm , given by is a convergence parameter for the noise rejection method. Positive ζ n i are rejected if they are greater than τ n rm σ n + , while negative ζ n i are rejected if they are less than τ n rm 2 σ 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 = {i = 1, 2, · · · , N | ψ n (t i ) = ∞ } as stated before).

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 i 0 act (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 i max , and (iii) the time at which the deactivation ends and the signal is once again dominated by the drift t i 0 end . 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 0 local = min { i local }. 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 t i 0 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 : q = i min local s , r = i max local s , s ∈ {1, . . . , N local } and select the first local maximum at the location i 0 max = min r : µ r ≥ mean (µ) − 3 · std(µ) and ν r > mean (ν) − std (ν) .
which has an average rate and magnitude within a statistically acceptable range (that excludes small outliers). t i 0 max 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 i 0 act ( Figure 2D). This is because i 0 act depends solely on the derivative around its first local maximum, while i 0 max takes into account an approximation to the average rate of change around all local maxima; the local maximum after i 0 act should only differ from i 0 max 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 i end max 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 i max = min arg max i∈i 0 max ,...,i end max Aû i .
iii. Having identified the time at which activation is likely to begin t i 0 act , 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 i 0 act . 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 i 0 end 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 = 1, . . . , i 0 act i 0 end , . . . , N j (i 0 act and i 0 end 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 θ 0 drift . 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 θ 0 drift 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Ḋ t; θ 0 drift is minimal. To achieve this, we choose w to depend on y = û −Ḋ t; θ 0 drift as follows We then apply the non-linear least squares fitting procedure to minimize the error function and use θ 0 drift as an initial condition for the fitting procedure.

Activation Fitting
For the activation phase of the response, we use the model g act (t; A act , β, n act , m act ) = A act t n act t n act + β n act + m act t 0 x n act x n act + β n act dx.
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 quasilinear 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 t i 0 act , t i max . 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θ act drift 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 D   t; arg min θ ∈ θ drift ,θ act drift S drift (θ )    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 0 max . 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.
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 g de (t; A de , γ , n de ) = A de γ n de t n de + γ n de + m de t 0 γ n de x n de + γ n de dx, 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 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.
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 G resp t; θ resp , θ drift = D t; θ resp + g resp t i ; θ resp if t on ≤ t D t; θ resp otherwise , and minimize the error function (14) to obtain the fitting, where κ is a parameter quantifying the apparent coherence between the drift and reponse models (D and g resp ), given 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 i 0 act 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).

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.
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 halfmax (FWHM) is defined as the time elapsed between the two halfmax 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 t;θ resp , truncated at its half-maximums, whereθ 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 (t peak ,t trough ) of significant peaks and troughs of the deviations, respectively (see Figure 2D). In total, there arẽ N = 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 The weighting functions π peak i , π trough i 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Ġ resp t;θ resp . 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θ resp ), which we take to be the baseline of the MPR-associated deviations. Given the estimate G resp t; θ resp , 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.,t peak ort 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 periodand 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 ofN peak peaks (N trough troughs) occurring at timest peak (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 ℓ osc =t peak N peak −t peak 1 ), 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 nonstationary 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).
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).
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 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. 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. it. In particular, the manual evaluations of AUC values, decay constant and the time of onset, as well as most of the MPRparameter 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 timedomain. 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 largescale 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 multiparametric 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, realtime 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.