Computational Analysis of Ca2+ Oscillatory Bio-Signals: Two-Parameter Bifurcation Diagrams

Two types of bifurcation diagrams of cytosolic calcium nonlinear oscillatory systems are presented in rectangular areas determined by two slowly varying parameters. Verification of the periodic dynamics in the two-parameter areas requires solving the underlying model a few hundred thousand or a few million times, depending on the assumed resolution of the desired diagrams (color bifurcation figures). One type of diagram shows period-n oscillations, that is, periodic oscillations having n maximum values in one period. The second type of diagram shows frequency distributions in the rectangular areas. Each of those types of diagrams gives different information regarding the analyzed autonomous systems and they complement each other. In some parts of the considered rectangular areas, the analyzed systems may exhibit non-periodic steady-state solutions, i.e., constant (equilibrium points), oscillatory chaotic or unstable solutions. The identification process distinguishes the later types from the former one (periodic). Our bifurcation diagrams complement other possible two-parameter diagrams one may create for the same autonomous systems, for example, the diagrams of Lyapunov exponents, Ls diagrams for mixed-mode oscillations or the 0–1 test for chaos and sample entropy diagrams. Computing our two-parameter bifurcation diagrams in practice and determining the areas of periodicity is based on using an appropriate numerical solver of the underlying mathematical model (system of differential equations) with an adaptive (or constant) step-size of integration, using parallel computations. The case presented in this paper is illustrated by the diagrams for an autonomous dynamical model for cytosolic calcium oscillations, an interesting nonlinear model with three dynamical variables, sixteen parameters and various nonlinear terms of polynomial and rational types. The identified frequency of oscillations may increase or decrease a few hundred times within the assumed range of parameters, which is a rather unusual property. Such a dynamical model of cytosolic calcium oscillations, with mitochondria included, is an important model in which control of the basic functions of cells is achieved through the Ca2+ signal regulation.


Preliminaries
Periodic cycles-solutions of nonlinear planar (two variables) systems-are covered by the celebrated Poincaré-Bendixson theorem proven (under certain assumptions) with the use of the Green theorem [1]. To verify an absence of limit cycles, one may use the Dulac's criterion [2]. Such an approach becomes problematic for multi-variable nonlinear systems having a wide spectrum of possible steady-state responses that change (bifurcate) with slowly varying parameters. Verification of periodic (or oscillatory chaotic) responses in such cases is done by numerical methods, often with a parallel computing approach [3,4]. Determining bifurcations of nonlinear dynamical systems often require a sophisticated computational approach to visualize the obtained results [5,6].
In this paper, we deal with the computer-aided verification of oscillatory responses and qualitative analysis of nonlinear dynamical systems [7]. An application case study is presented for a rather intriguing cytosolic (with mitochondria) calcium oscillatory dynamical model having numerous nonlinear terms and interesting properties, one of which is a wide change of frequency of oscillations when parameters vary. We compute twoparameter bifurcation diagrams of that system to identify areas of periodic and non-periodic steady-state responses.

The Reliability Issue of the Numerically Computed Periodic and Chaotic Solutions
One of the major reasons of writing this paper follows from the well-documented issue of the reliability of the chaotic and periodic solutions obtained with finite precision computations. The central issue is the fact that obtaining entirely identical long-term chaotic solutions via numerical solvers is practically impossible even when the same algorithm, time step and initial conditions are applied to a particular nonlinear ordinary differential equation (ODE) system. Reports on such an issue with interesting findings for the well-known chaotic systems are discussed in [8]. The concept of critical time, T c , has been introduced, to mark the time instant beyond which two chaotic solutions depart from each other to become completely different. This has been further expanded in [9,10] and the best description of the problem when dealing with chaotic solutions is the following quotation from [9]: ". . . even when the initial conditions are exactly the same, the algorithm is not changed and the step time is kept unchanged, one can find multiple pseudo-orbits simply by changing the way the model is written" and further in the same paper: "The question that arises (. . . ) is: with the same initial condition, same scheme of discretization and equivalent set of differential equations, which orbit is the true one? This is not an easy question to answer since there is almost no option but to use finite precision machine." In [8], the following observation is made: "Parker and Chua [11] pointed out that a 'practical' way of judging the accuracy of numerical results of a non-linear dynamic system is to use two (or more) 'different' routines to integrate the 'same' system: the initial time interval over which the two results agree is then 'assumed' to be accurate and predictable. More precisely speaking, the computed results beyond the critical decoupling time T c are not reliable." The above issues are, to some extent, related to the results presented in the present paper, in which we follow the guidance expressed in the above quotations, and study two different bifurcation 2D diagrams (of n-period and frequency types) for one particular non-linear model with complex Ca 2+ oscillatory solutions. Our goal was to identify and predict periodic and chaotic solutions using those two types of diagrams. Within the n-period type of 2D diagrams, we differentiate our calculations by using three maximum values of n (64, 128 and 182; see Section 3). To further modify our approach and numerical algorithms, we used two other tools to identify periodic and chaotic oscillations, namely the sample entropy and 0-1 test for chaos (Section 5). Throughout the paper, we compare the obtained results (2D diagrams) to conclude that the same types of solutions are obtained no matter which tool (type of diagram) is used. However, we would also like to point out that the issues raised in [8][9][10][11] are not fully relevant to our study. Even when one obtains two different chaotic solutions and attractors beyond the critical time T c , such solutions will still be identified as chaotic. Thus, our diagrams will stay the same, and the color indicating particular response will not change (see Sections 2 and 3). The 0-1 test for chaos (see Section 5) will still result in a number close to 1 indicating chaotic responses.

Autonomous System of Calcium Oscillations and Its Basic Properties
The slowly varying parameters of nonlinear dynamical systems described by ODEs may cause significant changes to the systems' steady-state responses. Typical visualization of those changes when one parameter varies slowly has the form of one-parameter bifurcation diagrams (figures) with the varying parameter representing the horizontal axis while the vertical axis is a certain quantity characterizing the changing responses. That quantity could be the identified maximum values of the response in a chosen time interval, the output of the 0-1 test for chaos (a number between 0 (for periodic responses) and 1 (for chaotic ones)), the frequency of periodic output, entropy or others [4]. In this paper, we illustrate our bifurcation diagrams (mostly with two slowly varying parameters) with the use of an interesting autonomous system described in Appendix A. The system (A3) has been used in the literature for a qualitative analysis of the oscillatory cytosolic calcium responses [12][13][14][15][16][17][18].
As a result of the many nonlinear terms, a computational approach to analyze is (A3), but some analytical analysis is also possible [12]. Taking into account the structure of the right-hand functionsf i , i = 1, 2, 3, in (A3), it can be shown that for the non-negative initial conditions, the three solutions, x i (t) are all non-negative and bounded. Furthermore, by denoting the equilibrium points of the three variables by (x 0 , y 0 , z 0 ), one can compute x 3 fromf 3 , substitute it intof 2 to express x 2 as a function of x 1 and thenf 1 can be transformed into a polynomial equation in x 0 only (with the parameters p). Fixing 15 out of possible 16 parameters and considering the remaining parameter as slowly varying, we can establish a dependence of x 0 on that one varying parameter. Moreover, the linearization of (A3) along x 0 will yield a Jacobian with the characteristic polynomial of the form [13] which can further be used to analyze the stability properties, and, for example, to find the Hopf bifurcation points. Such an analysis is not very relevant to the topic of this paper, therefore it is omitted here. Figure 1 illustrates one of the many one-parameter bifurcation diagrams of (A3) when k ER,ch is the slowly varying parameter and the vertical axis shows the values of maximum points identified in one period of the steady-state solution of the system. The L s notation (with L and s being co-prime integers) indicates a mixed-mode periodic oscillation. L and s are the numbers of large and small maximum values of oscillations in one period. For example, the 2 3 mixed mode oscillation in an interval around k ER,ch = 570 has 2 large maximum values (close to each other) around the value of 0.49 and three small maximum values around 0.33-0.36. The diagram in the marked rectangle in Figure 1 is shown in more detail in Figure 2. Note an interesting property of (A3)-the L s periodic responses form a Farey sequence of co-prime integers. For example, 2 3 is formed from the two neighboring sequences 1 1 and 1 2 , with 2 3 = 1 1 1 2 , where is the Farey addition. In general, L s = L s 1 1 L s 2 2 , with L = L 1 + L 2 and s = s 1 + s 2 . See [14] for more details and the relation of the L s type of oscillations to the Ford circles, Stern-Brockot tree, sequences of firing numbers and the Riemann zeta function. The L s mixed-mode periodic oscillations have also been considered in the context of two-parameter bifurcation diagrams for the modified Chua's circuits in [6]. The intermingling intervals of periodic and chaotic solutions (as in Figures 1 and 2) can also be used to determine fractal dimension of (A3) [14]. We chose to vary the same parameters of (A3) that were used in [12,[15][16][17][18] for their one-parameter bifurcation analysis. The remaining parameters are constant in our 2D bifurcation analysis, as described in Appendix A.   Figure 1). Note that 3 5 = 2 3 1 2 .

Bifurcation Diagrams of Period-n Responses and Frequency Distribution
In the next section, we use the L s type of periodic oscillations as (L+s)-periodic oscillations, indicating that in one period there are L + s maximum values. Thus, the twoparameter diagrams of n maximum values, n = 1, 2, . . . , 64, presented next, indicate that the solutions are periodic, with n = L + s. Thus, no distinction is made to differentiate between large and small maximum values. For this reason, two different L s 1 1 and L s 2 2 oscillations are both of the same period-n type, if L 1 + s 1 = L 2 + s 2 . From the point of view of period-n oscillations, we have a possibility of n + 1 different types of oscillations, called the period-n ones, as, for example, the following six period-5 oscillations may come from the 5 0 , 4 1 , 3 2 , 2 3 , 1 4 and 0 5 mixed-mode oscillations. All these six mixed-mode oscillations are marked by the same color and the number n = 5 in the two-parameter bifurcation diagrams of the number of maximum values in one period.
The above period-n oscillations as illustrated in Figures 1 and 2 (or in the two-parameter diagrams in the next section) do not provide any information as far as the frequency of oscillations is concerned. Once a period-n oscillation is established, an extra computational work can be performed to measure the frequency and, thus, create one-or two-parameter frequency diagrams. Both, the period-n and frequency diagrams complement each other to provide a better understanding of how nonlinear dynamical systems behave when their parameters vary. The frequency diagram may be essential in other kinds of analysis. For examples, it has been established recently that applying the 0-1 test for chaos may require a prior knowledge of the maximum frequency of the continuous spectrum of chaotic signals to properly choose certain parameters used in that test [19,20].
Depending on the ranges of varying parameters, it often happens that the same type of period-n oscillations differ significantly in frequencies, as we have discovered by creating diagrams in the next section. The difference could be by a factor of a few hundred times, as it happens for the system (A3). Further discussion on the range of frequency changes in the calcium oscillatory systems can be found in [21].
The wide range of frequencies of oscillations of (A3) cannot be seen from the diagrams in Figures 1 and 2 or other similar one-parameter diagrams, but it is known that (A3) has that peculiar feature of a wide range of frequency change and rather small amplitude change, even with a small interval of the varying parameters, including k ER,ch ("In some cells, time intervals between two Ca 2+ spikes change from a few seconds to several minutes or even tens of minutes. However, in contrast to the large changes in frequency, the amplitude of Ca 2+ oscillations remains nearly constant in many cell types. This likely allows the cell to encode information in the frequency of Ca 2+ oscillations." See [16]).  Figure 3c is enlarged in Figure 3e, while the diagram in rectangle C in Figure 3d is enlarged in Figure 3f. To create each of the diagrams in Figure 3, we solved (A3) 360,000 times as 600 × 600 discrete points were used. As in the case of one-parameter diagrams in Figures 1 and 2, one cannot deduct any information about the frequency of period-n oscillations from the diagrams in Figure 3. As explained before, each discrete point with n = 1, . . . , 63 maximum values in one period indicates the sum L + s (= n) rather than the L and s separately. Of course, the maximum values of the analyzed variable in one period cannot be retrieved from the diagrams in Figure 3. On the other hand, the periodic oscillations with n ≥ 64 and chaotic areas represented by the white color in Figure 3 are identifiable in the same way one can identify them in one-parameters diagrams. Those are rather narrow intervals with multiple dots in vertical lines distributed between the lower and upper maximum values seen in Figures 1 and 2, for example, for the k ER,ch values less than 566 and around 578 in Figure 2.   Figure 4 presents frequency distribution of the periodic solutions of (A3) in the same rectangular areas as the period-n bifurcation diagrams in Figure 3. Again, the blue color has the same meaning as that explained before for the blue areas in Figure 3. The extra green color in Figure 4 is not associated with any frequency value as it represents the solutions of either period-n oscillations for n ≥ 64 or chaotic solutions for which, obviously, one cannot find frequency. This means that the green areas in Figure 4 correspond to the white areas in Figure 3.

Two-Parameter Frequency Distribution Diagrams
Furthermore, the diagrams in the larger two-parameter rectangles in Figure 5 (with an accompanying Figure 6) and Figure 7 (with an accompanying Figure 8) are the period-n and frequency distribution diagrams, respectively. The smaller rectangles D in those figures are the rectangles in Figures 3a and 4a, respectively. The purpose of creating Figures 5 and 7 is two-fold. First, notice that the central part of the diagram in Figure 5 is dark brown, indicating mostly period-1 and period-2 oscillations.  However, from Figure 7, it follows that some of those oscillations have frequencies that differ significantly. Notice the very dark brown area next to the blue area on the left side in Figure 7. The frequencies there are as low as 0.003 Hz. Similar period-1 frequencies in the very narrow yellow strip next to the blue area in the lower right corner are around 0.6 Hz. Thus, the frequencies of the same type of period-1 oscillations differ by the factor of around 200. Selecting a point on the line a for k ER,ch = 4400 in Figure 3a indicates period-1 oscillations to which we assign the corresponding point on line a in Figure 4a on the narrow yellow band with frequency around 0.6 Hz. Similar significant change in the frequency values can be observed between the yellow and brown areas in the lower right corner areas in Figure 4b-d. This confirms the second reason for the usefulness of the two-parameter diagrams: they give a broader view of the possible changes as far as the types of period-n solutions and their corresponding frequencies are concerned. Figures 6 and 8 show one-parameter diagrams of the n values and frequency, respectively, as functions of k ER,ch for three values of K ch . Those are simply horizontal cross-cuts of Figures 5 and 7. The narrow intervals with spikes at n = 64 in Figure 6 indicate intervals of k ER,ch yielding white bands in Figure 5 and green bands in Figure 7. Similarly, the frequency spikes in Figure 8, reaching the values close to 0.6 Hz, correspond to the points in a narrow yellowish-white band close to the blue area in the lower right corner in Figure 7.
The two types of the two-parameter diagrams for (A3) are not the only ones that can be created for the same autonomous system. It is certainly possible to compute the twoparameter diagrams for Lyapunov exponents [22], stability domains [23], Poincaré return maps [24], entropy [25] and the 0-1 test for chaos [26]. The later two-parameter diagrams for the electric circuits, Lorenz and Rössler chaotic systems, are presented in [20,27,28].
In addition, notice that while the maximum value of n used in the diagrams in Figure 3 is 64, the character of the diagrams is preserved when that value is increased. For example, Figure 9a,b show how the diagram from Figure 3b changes when the maximum value of n is increased to 128 (Figure 9a) and 182 (Figure 9b). The reduction of the amount of pure white color with the increased maximum of n is clearly noticeable. This indicates that for many of those pairs of parameters of (K ch , K ER,pump ), for which n = 128 (the white color in Figure 9a), the color changes to yellow in Figure 9b. There are rather few pairs of parameters (K ch , K ER,pump ) with n values close to 182, as shown in Figure 9b. The sharp jump of n from single digits (i.e., n < 10 represented by the brown color) to high values (i.e., n around 182 represented by the white color) is an indicator of the changes in oscillatory behavior similar to a periodic to chaotic bifurcation or periodic to bursting transition. This happens, for example, around the area of (K ch , k ER,pump ) = (5.5, 17) in Figures 3b and 9a,b.

Details of Computations of the Two-Parameters Bifurcation Diagrams
Generally, if an oscillation system has periodic solution then for any two variables, say x i (t) and x j (t), we have a homoclinic orbit F(x i , x j ) = 0 in the form of closed loops. The number of loops n is equal to the number of local maxima (or minima) of both x i (t) and x j (t) occurring within one period, as shown in Figure 10 for n = 1 and in Figure 11 for n = 5.
To identify the number n, it is first necessary to determine the geometrical center of the innermost loop, denoted by O(x 2m , x 1m ) in Figures 10 and 11 for the variables x 1 and x 2 . It is well-known that for any two-dimensional area A in variables x and y, the center is determined by the double integrals When we consider the most inner loop (the loop of the smallest area), then an approximate method to find the center of that loop uses the four points P 1 , . . . , P 4 defining the quadrangle P 1 P 2 P 3 P 4 , as is shown in Figures 10 and 11. The center of that quadrangle has been determined by the values max(x 2min ), min(x 2max ), max(x 1min ) and min(x 1max ) (see Figure 11), where, in general, max(x jmin ) denotes the largest local minimum of x j (t), while min(x jmax ) denotes the smallest local maximum of x j (t). As the computations showed, the above method of determining the center O proved to be accurate. To avoid any transients, such an analysis was carried out in this paper in the interval of identification [500, 1000] seconds, while the solution of (A3) was computed in the interval 0 ≤ t ≤ 1000. The number of intersections of the graph F(x i , x j ) = 0 with the positive half axis x 2new (the points C i in Figures 10 and 11) is exactly equal to the number of local maxima of x 2 (t) in one period. The solution values of x(t) are computed at discrete points. The intersections values (coordinates of the points C i ) were obtained as a result of linear interpolation, based on the two nearest discrete values (below and above the x 2new axis). The value of the time variable t, for which the intersections with the positive half axis x 2new happen, is determined analogously by linear interpolation. The period T 0 was determined as T 0 = t i+n − t i , where n is the number of different intersections of the positive half axis of x 2new . We assumed in this paper that the intersections of the positive half axis x 2new are considered identical if the values obtained by linear interpolation do not differ by more than the assumed certain error value . We assumed = 0.0001 for all computations of both one-and two-parameter diagrams in this paper. When a chaotic response is obtained, then the calculations of both n and T 0 are not possible, as no homoclinic orbits occur, as illustrated in Figure 12.
The number of intersections of the graph F(x i , x j ) = 0 with the positive half axis x 2new (the points C i in Figures 10 and 11) is exactly equal to the number of local maxima of x 2 (t) in one period. The solution values of x(t) are computed at discrete points. The intersection values (coordinates of the points C i ) were obtained as a result of linear interpolation, based on the two nearest discrete values (below and above the x 2new axis). The value of the time variable t for which the intersections with the positive half axis x 2new happen, is determined analogously by linear interpolation. The period T 0 was determined as T 0 = t i+n − t i , where n is the number of different intersections of the positive half axis of x 2new . We assumed in this paper that the intersections of the positive half axis x 2new are considered identical if the values obtained by linear interpolation do not differ by more than the assumed certain error value . We assumed = 0.0001 in all computations of both one-and two-parameter diagrams in this paper. When a chaotic response is obtained, then the calculations of both n and T 0 are not possible, as no homoclinic orbits occur, as illustrated in Figure 12.   Finally, in the identification interval 500 ≤ t ≤ 1000 used in this paper, we have 57, 205 and 254 points C i in Figures 10-12, respectively. The number n was determined to be respectively 1, 5 and undefined for the three cases. Due to the long computation time required for the creation of two-parameter diagrams, all computations presented in this paper were done by using a 12-core computer and parallel algorithms described in the previous work [28].

Comparison with the Sample Entropy and 0-1 Test for Chaos Diagrams
The diagrams in Figures 3a and 4a are compared with the sample entropy (SE) [29,30] and 0-1 test (T01) [19,20] for chaos diagrams in Figure 13 in the same rectangular areas of the two varying parameters. The diagrams in Figure 13 are of size 1000 × 1000. Thus, system (A3) was solved 10 6 times in the assumed rectangular areas. The sample entropy diagrams in Figure 13a Figure 3a,b, respectively; (c,d) two-parameter 1000 × 1000 diagrams of the 0-1 test for chaos (K values) corresponding to the diagrams in Figure 3a,b, respectively. The Runge-Kutta method of order 4 with constant step-size of integration 0.001 was used to solve (A3). The sample entropy diagrams were obtained with m = 3, r = 0.025 and N = 10, 000 (see [29]). The SE diagrams in (a,b) and the T01 diagrams in (c,d) were identified based on the solutions for 1000 ≤ t ≤ 2000 and 300 ≤ t ≤ 500, respectively. The parameter T = 40 in the 0-1 test method (see [19,28]

Conclusions
Two types of two-parameter bifurcation diagrams (with the period-n responses and frequency distributions) were presented in this paper. Those diagrams provide useful information about the nature of the steady-state solutions of the nonlinear autonomous system of calcium oscillations with the mitochondria included. The same type of diagrams can be, in principle, created for any other nonlinear autonomous system with varying parameters. In general, because the two-parameter diagrams give different information about the changes in the steady-state responses than the one-parameter diagrams do, both kinds of diagrams complement each other. Our computations were done using parallel computing algorithms, since to create one diagram with two varying parameters requires solving the nonlinear autonomous system a few hundred thousands or even millions of times, depending on the two-parameter area and the resolution one desires. The identified possible chaotic oscillations occur in very narrow intervals of the parameter k ER,ch in Figures 1 and 2 and scattered areas of white and green colors in Figures 3 and 4, respectively. Those scatter areas of chaotic solutions are identical for both types of diagrams. Further comparisons with the sample entropy and 0-1 test for chaos diagrams also show a very good agreement of the obtained solutions.
The diagrams presented in this paper complement other types of one-and twoparameter diagrams one can create for nonlinear autonomous systems, such as, for example, those showing the Lyapunov exponents, Poincaré maps [31], the 0-1 test results [28] and entropy distribution [29,30].
where the stands for the time derivative, and J ER,ch = k ER,ch