Glycolytic Synchronization in Yeast Cells via ATP and Other Metabolites : Mathematical Analyses by Two-Dimensional Reaction-Diffusion Models

Possibilities of synchronized oscillations in glycolysis mediated by various extracellular metabolites are investigated theoretically using two-dimensional reaction-diffusion systems, which originate from the existing seven-variable model. Our simulation results indicate the existence of alternative mediators such as ATP and 1,3-bisphosphoglycerate, in addition to already known acetaldehyde or pyruvate. Further, it is also suggested that the alternative intercellular communicator plays a more important role in the respect that these can synchronize oscillations instantaneously not only with difference phases but also with different periods. Relations between intercellular coupling and synchronization mechanisms are also analyzed and discussed by changing the values of parameters such as the diffusion coefficient and the cell density that can reflect intercellular coupling strength.


Introduction
Synchronization in biological systems is widely observed in the natural world.This phenomenon is one of the collective behaviors thought to have a crucial role in maintaining the individual life or in giving benefits for communities.One of the well-known examples is synchronized flashing of male fireflies [1] [2].
Glycolysis is a biological mechanism to decompose glucose and to store energy in the form of ATP.In this chemical process, synchronized variations can be seen for the concentrations of various metabolites [3]- [8].For example, yeast cells exhibit glycolytic synchronized oscillations under certain conditions.A suspension of yeast cells with high population densities shows synchronization, whereas that with low densities does not.The emergence of the collective behavior such as synchronization may occur above a critical cell density.
The theoretical studies of glycolytic oscillations in yeast cells using differential equations started substantially when Sel'kov presented a simple kinetic model of an enzyme reaction with substrate inhibition and product activation [9].His ordinary differential equation system with two variables showed the occurrence of periodic selfoscillations in glycolysis.Another starting point is the allosteric enzyme model by Goldbeter and Lefever and Goldbeter, which also consists of two variables that referred to ATP and ADP, respectively [10] [11].They adopted a partial differential equation system, which afterward facilitated a large extent of investigations for pattern formation and spatiotemporal structures.As these examples, we exemplify the studies of target pattern formation and spatiotemporal chaos by Zhang et al. [12] and of inward rotating spiral waves by Straube et al. [13].Meanwhile, the Sel'kov model is succeeded by Lavlova et al. in the study of inward and outward wave propagation [14].Besides these two streams, Bier et al. proposed a simple conceptual model consisting of glucose and ATP for explanation of glycolytic synchronized oscillations [15].
The epoch-making study was performed by Wolf and Heinrich in order to elucidate the mechanisms of synchronous behaviors in glycolytic oscillations [16].Based on the detailed description of glycolytic reaction processes, they constructed a seven-variable model, by which the effect of intercellular coupling on oscillatory dynamics was theoretically analyzed.Further, Henson et al. connected this model with the cell-ensemble modelling technique and showed that a large ensemble of about 1000 population cells were required to adequately capture complex dynamic behaviors in glycolytic oscillations [17].
However, it seems that no model has succeeded in demonstrating perfect synchronization of a large number of yeast cells that oscillated with different phases and periods.We guess that one of the reasons is that these models assumed acetaldehyde or pyruvate to be an intercellular communication substance for glycolytic synchronization.
It is known that acetaldehyde mediates the synchronization of glycolytic oscillations [6].Various experiments demonstrated that acetaldehyde had a very strong synchronization effect on a suspension of yeast cells [8].However, these studies do not exclude the possibility that extracellular nucleotides such as ATP and ADP can play an important role in synchronization of glycolytic oscillations in yeast cells.These are signaling molecules contained in all tissues [18].Other studies also showed that Saccharomyces cerevisiae released ATP to extracellular solutions, suggesting that extracellular ATP, ADP, AMP or adenosine played a role in yeast physiology as intercellular communicators [19] [20].
Then, we attempted to exchange acetaldehyde or pyruvate to ATP or other metabolites in the original model by Wolf and Heinrich [16].Our present study using these modified seven-variable models shows that glycolytic synchronization in yeast cells may occur via intercellular mediators such as ATP and other metabolites, and that ATP is the most effective synchronizer even under the condition where individual cells have quite different properties.

Mathematical Models
Our mathematical models are derived from the glycolytic oscillation model originally presented by Wolf and Heinrich [16].Major modifications are the following three: 1) Extension to two-dimensional partial differential equation systems with the diffusion term, 2) Exchange of the intercellular coupling substance and, 3) Non-dimensionalization for variables and parameters.
The first modification enables us to observe directly various kinds of oscillations including the synchronous one, and the possibilities of synchronization significantly increase by the second modification.Moreover, the third modification, which leads to the reduction of parameter numbers, eases mathematical analyses of simulation models.
The schematic diagram of the model is sketched in Figure 1, where the simulation area consists of N × N partitioned square compartments divided at regular intervals.One cell is embedded within each compartment [16] [21], whose volumetric ratio φ to the compartment is one of the key control parameters of our models.In this ar-ticle, we adopt N = 11 in all simulations, resulting in the total cell number is N × N (=121).However, the number of N can be decided arbitrarily.
Six major substances are contained within a cell, which are denoted as S 1 , S 2 , S 3 , S 4 , N 2 and A 3 .These are regarded as independent variables, the meanings of which are explained in Table 1.In addition to these six, two more substances, N 1 and A 2 , are also contained, which are connected by the relations, N 1 + N 2 = N and A 2 + A 3 = A, respectively.As N and A are constant, two variables, N 1 and A 2 , are not independent, whose numbers are automatically identified.In the original model, it is assumed that only S 4 can permeate through the membrane, then, the substance in the external solution is labeled as 4 ex S [16].The model of this type, via the S 4 intercellular mediator, is also the starting point of our studies.The transmembrane coupling metabolite is later altered to the other one such as A 3 , S 3 , S 2 and N 2 ; thereafter, the variations of system behaviors are investigated.
The meanings of parameters used in this article are also explained in Table 1, most of which are recycles of those in the original model [16].The dimensionless parameter values are listed in  In Model I, the transmembrane substance is S 4 , whereas in Models II, III, IV and V, this is altered to A 3 , S 3 , S 2 and N 2 , respectively.Intercellular coupling metabolites diffuse through the boundary between adjacent compartments, which is expected to be feasible mechanisms to induce synchronized glycolytic oscillations.As for identification of S 1 , S 2 , S 3 , S 4 , N 1 , N 2 , A 2 and A 3 , see Table 1.Numerical analyses are performed using these parameter values.A parameter k7 is newly introduced for stabilization of the system.The diffusion coefficient dF and the parameter f0 that defines the amplitude of randomization are also employed.When f0 = 0.1, for example, randomized items such as initial values and parameter values are scattered within the range from 90% to 110% around the above values.As for two control parameters, φ and dF, two different values are provided in Model II, III and IV, while only dF is varied in Model I. Reference values are converted by non-dimensionalization procedures from those in the original model by Wolf and Heinrich [16].
Two parameters, k 1 and K I , are eliminated by these processes [21].Besides, the original parameter k is renamed r, and a new parameter k 7 is employed.The system is stabilized by the addition of k 7 , especially when the transmembrane substance is changed to A 3 , S 3 , S 2 or N 2 .As a result, our seven-variable, two-dimensional reaction-diffusion model with S 4 being the intercellular mediator reads as follows.
Here, the diffusion term is incorporated in the last equation, which describes the intercellular coupling via S 4 .This basic model described by Equations ( 1) is referred to as Model I in this article.
Next, we exchange the coupling substance in the external medium in order to examine other candidates that can produce synchronous oscillations.It should be noted that a new parameter k 7 takes an important role as a stabilizer in these models.For example, when the substance that functions as the communicator is A 3 , the mathematical model is described such as To be accurate, one more variable A 2 or 2 ex A should be added besides 3 ex A , however, it is confirmed that the addition of either variable hardly influences the simulation result.The system described by Equations ( 2) is referred to as Model II in this article.
Similarly, we can construct the system with the intercellular communicator S 3 (Equations ( 3)) and the system with the intercellular communicator S 2 (Equations ( 4)).These are referred to as Model III and Model IV, respectively.

( ) (
) ( ) In above formulations of Equations ( 3) and Equations ( 4), several equations are omitted, which are the same as those in Equations ( 1).Moreover, we can construct Model V, where the intercellular mediator is N 2 , in the same manner, formulations of which are not shown explicitly.
Although two parameters are eliminated by non-dimensionalization procedures, there are still twelve parameters in our models except for two control parameters d F and φ.In principle, we try to use parameter values of the reference state converted from the original model [16], which are listed in Table 2.If it is impossible to make stable limit cycle oscillations using these values, we adjust the values of such parameters as J 0 , r and κ, which are thought to be comparatively easy to manipulate.The exception is Model V, where the value of k 5 is varied.Besides the reference [16], the parameter values in the reference [17] are also referred to.
It is also important to choose properly the initial value of each variable.We determine these values in reference to fixed points, which are calculated in advance.Table 3 is a list of initial values used to create synchronous oscillations.In the case of Model II, for example, the initial value of S 1 is fluctuated within the range from 4.68 (= 5.2 × 0.9) to (5.72 = 5.2 × 1.1), S 2 = 0.43 (= 4.3 × 0.1), and so on.We also examine the case where not only initial S 1 values but also values of nine parameters, J 0 , k 2 , k 3 , k 4 , k 5 , k 6 , k 7 , r and κ, are randomized.In this case, oscillation periods are also disturbed as well as phases.However, it seems impossible to annihilate these differences, which leads to the conclusion that oscillations in Model I are asynchronous for randomization of both initial S 1 values and nine parameter values, i.e., for disturbances of both oscillation phases and periods.

Synchronization and Desynchronization in Model I
Assuming that the diffusion coefficient d F = ∞, our Model I corresponds exactly with the original model by Wolf and Heinrich [16].However, the simulation results are not altered significantly compared with the case of  d F = 1.0, thus, we can expect almost the same results for simulations in the range of d F ≥ 1.0.These situations are also true for Models II, III, IV and V as well.

Synchronization and Desynchronization in Models II, III, IV and V
Despite almost complete synchronization for randomization of initial S 1 values, we did not succeed in detecting any synchronous oscillation for randomization of parameter values.Then, we exchange the intercellular mediator S 4 to other substances in an attempt to synchronize oscillations with different periods.In the end, we examine Model V with the N 2 intercellular mediator.Despite a large extent of surveys, we failed in finding any parameter set that realized synchronous oscillations.that N 2 is concerned with synchronization.
These simulation results in five models are summarized in Table 4.The oscillation period of Model II is extremely short and almost half compared with those of other models.

Two Levels of Synchronization
It seems to be two levels in synchronization of oscillations, namely, synchronization of phases and that of periods.In general, perturbations of initial values cause merely phase shifts.Meanwhile, perturbations of parameter values are at least required to cause the difference of oscillation periods.Thus, it is likely that synchronization of different periods is more difficult and essential than that of phase differences.Taking this inference into consideration, we cannot say that synchronization is completed or perfect until oscillations become in phase even for randomization of parameters.In this sense, synchronization in Model I (Figure 2 Comparing three kinds of synchronization in Models II, III and IV, it seems that overlapping of oscillatory patterns in Model II is more outstanding than the other two.Thus, it could be speculated that synchronization via  the A 3 intercellular metabolite is the most perfect, indicating that the exchange of A 3 is the most feasible mechanism for glycolytic synchronization in yeast cells.

Quantification of Synchronization
In order to certify above inference numerically, the standard deviations of N 2 concentrations are calculated for total 121 cells in Models II, III, IV and V, among which the former three are the synchronous cases and the last is the asynchronous case.The abscissae of It is supposed that these normalized standard deviations reflect the degree of synchronization, and that the smaller value, the more perfect is synchronization.Among three synchronous cases in Figure 6, the values of Figure 6(a) are smaller than those of Figures 6(b) and 6(c), indicating that synchronization in Model II is more perfect than that in Model III or IV.This is the reason why we predict that A 3 is the most probable metabolite involved in synchronous glycolytic oscillations in yeast cells.
Meanwhile, Figure 6(d) is the asynchronous case where normalized standard deviations keep high values.Moreover, any periodic structure is not detected, which is thought to be a typical characteristic of asynchronous oscillations.

Dependencies on Diffusion Coefficient and Cell Density
It is well known that glycolytic oscillations initiate synchronization under the densely populated condition [3]- [8].Considering that high population densities mean intense coupling between cells, it is thought that collective properties of synchronization in glycolytic oscillations is triggered by the increase in such parameters as the diffusion coefficient d F , the volumetric ratio of the cell to the compartment φ, the kinetic constant relating to the permeability κ, and so on.These are connected with such parameters as the compartment size L, the compartment volume V, the cellular surface area A c , the cellular volume V c and the permeability of the cellular membrane P, by the following relations [16] [21].

,
, .In particular, we focus on d F and φ among these parameters in this article.According to our simulation results of Models II, III and IV, the increase in d F gives rise to the transition from the asynchronous oscillation (for example, Figure 4(c) and Figure 4(d)) to the synchronous one (for example, Figure 4(a) and Figure 4(b)).This phenomenon can be called the Kuramoto transition [22].On the other hand, the direct transition from the convergence to the fixed point (for example, ) is induced with the increase in φ in Models II and III.The characteristic of this phenomenon is that the shift to the synchronous oscillation proceeds directly without passing through the intermediate asynchronous oscillation.Thus, it could be speculated that the shifts to synchronization with the increases in the density parameter φ are due to mechanisms of the dynamical quorum sensing [23]- [25].

Conclusions
1) Glycolytic oscillations with different phases can be synchronized by means of intercellular coupling via such substances as S 4 , A 3 , S 3 and S 2 , as demonstrated in Models I, II, III and IV.Meanwhile, synchronization of oscillations with different periods can be mediated by intercellular coupling substances such as A 3 , S 3 and S 2 , as demonstrated in Models II, III and IV.The latter synchronization is characterized by speedy convergence to the synchronous oscillatory state.
2) Among three candidates that can induce synchronization for different periods, A 3 could be the most responsible for the phenomenon, because the normalized standard deviations of N 2 concentrations in Model II are the smallest compared with those in Model III or IV.
3) The transition from the asynchronous to the synchronous oscillation is observed with the increase in the diffusion coefficient d F , which could be referred to the Kuramoto transition mechanism.On the other hand, the direct transition from the convergence to the fixed point to the synchronous oscillation is observed with the increase in the density parameter φ, which could be referred to the dynamical quorum sensing.Here, the diffusion term is neglected.
For another example, the fixed point of Model II is computed, as follows.Similarly, we can identify the coordinates of fixed points for Models III, IV and V as well.We would like to stress that coordinates of fixed points for each parameter value are calculated or identified precisely in all models presented in this article.

Figure 1 .
Figure 1.Schematic diagram of Models I, II, III, IV and V. (a) The simulation area is composed of N × N square compartments, in which a single glycolytic cell is embedded [16][21].In all simulations, the number of N is fixed at 11, thus, the total number of cells is N × N (=121).The temporal behaviors of five gray cells, one central cell plus four nearest neighbors, are selectively monitored in Figure2, Figure4and Figure5(figures on the left side).(b) Each cell contains eight metabolites, S 1 , S 2 , S 3 , S 4 , N 1 , N 2 , A 2 and A 3 , among which non-independent N 1 and A 2 are not shown.In Model I, the transmembrane substance is S 4 , whereas in Models II, III, IV and V, this is altered to A 3 , S 3 , S 2 and N 2 , respectively.Intercellular coupling metabolites diffuse through the boundary between adjacent compartments, which is expected to be feasible mechanisms to induce synchronized glycolytic oscillations.As for identification of S 1 , S 2 , S 3 , S 4 , N 1 , N 2 , A 2 and A 3 , see Table1.

Figure 2
Figure 2 displays the temporal changes in N 2 concentrations from dimensionless time t = 7500 to t = 8000 in Model I, where the transmembrane metabolite is assumed to be S 4 , and only initial values of S 1 are randomized within the fluctuation range from 75% to 125% around the main value.It should be noted that the perturbation of initial values causes the difference of phases only.Two figures (a) and (b) are the simulation results for d F = 1.0 and φ = 0.1, and (c) and (d) are those for d F = 0.01 and φ = 0.1.Meanwhile, (a) and (c) exhibit the temporal changes in central five cells, one central cell plus four nearest neighbors, and (b) and (d) do those in mean N 2 concentrations for total 121 cells.In the case of d F = 1.0 and φ = 0.1, as shown in (a) and (b), the tightly synchronized oscillation is observed, which is strongly certified by the fact that wave patterns of five central cells are nearly overlapped and almost equal with the averaged one.In the case of d F = 0.01 and φ = 0.1, on the other hand, it seems that oscillations continue to be out of phase, as shown in (c), thus, Model I does not generate synchronous oscillations in this case.The time series of two-dimensional distributions of N 2 concentrations are illustrated in Figure 3. Upper three figures (a), (b) and (c) exhibit the synchronous oscillation, where d F = 1.0 and φ = 0.1, while lower three (d), (e) and (f) do the asynchronous oscillation, where d F = 0.01 and φ = 0.1.For each case, three drawing times are chosen so that the central cell at (5, 5) just passes the mean N 2 value of the amplitude.Thus, they are slightly proceeded compared with the exact times, t = 2000, 4000 and 8000.As a result, it is confirmed that six central cells in Figure 3 are drawn by the same color.It should be noted that the synchronization processes in Model I advance comparatively at a slow pace.We also examine the case where not only initial S 1 values but also values of nine parameters, J 0 , k 2 , k 3 , k 4 , k 5 , k 6 , k 7 , r and κ, are randomized.In this case, oscillation periods are also disturbed as well as phases.However, it seems impossible to annihilate these differences, which leads to the conclusion that oscillations in Model I are asynchronous for randomization of both initial S 1 values and nine parameter values, i.e., for disturbances of both oscillation phases and periods.Assuming that the diffusion coefficient d F = ∞, our Model I corresponds exactly with the original model by Wolf and Heinrich[16].However, the simulation results are not altered significantly compared with the case of

Figure 2 .
Figure 2. Temporal changes in N 2 concentrations in Model I, t = 7500 -8000.Two figures on the left side show the temporal changes of central five cells specified in Figure 1 (a), while two on the right side show those averaged for total 121 cells.(a) and (b) show the synchronous oscillation, where d F = 1.0 and φ = 0.1.Meanwhile, (c) and (d) show the asynchronous oscillation, where d F = 0.01 and φ = 0.1.In both cases, only initial values of S 1 are randomized.Other parameter values are fixed in accordance with Table 2.It is recognized that the oscillations of five central cells are overlapped with each other, as shown in (a), and extremely similar with the averaged one, as shown in (b), indicating that almost all the cells oscillate with the same phases and periods.

Figure 3 .
Figure 3.Time series of two-dimensional distributions of N 2 concentrations in Model I. Parameter settings are the same as those in Figure 2. (a), (b) and (c) show the evolution to the synchronous oscillation, where d F = 1.0, φ = 0.1 and (a) t~2000, (b) t~4000, (c) t~8000, respectively.Meanwhile, (d), (e) and (f) show the succession of the asynchronous oscillation, where d F = 0.01, φ = 0.1 and (d) t~2000, (e) t~4000, (f) t~8000, respectively.Elapsed time is slightly adjusted such that the central cell at (5, 5) takes the mean N 2 value of the oscillation amplitude.

Figure 4
exhibits the simulation results of Model II, where the transmembrane substance is A 3 , and both initial S 1 values and nine parameter values of J 0 , k 2 , k 3 , k 4 , k 5 , k 6 , k 7 , r, and κ are disturbed at the same time.Two figures (a) and (b) are the simulation results for d F = 1.0 and φ = 0.5, (c) and (d) are those for d F = 0.01and φ = 0.5 and (e) and (f) are those for d F = 1.0 and φ = 0.1, respectively.Meanwhile, three figures (a), (c) and (e) on the left side represent the temporal changes in central five cells, while (b), (d) and (f) on the right side do those in mean N 2 concentrations for total cells.Different from the case of Model I, it is clear that the differences of both phases and periods disappear, and that synchronization occurs even for oscillations with different periods when d F = 1.0 and φ = 0.5, as shown in (a) and (b).Moreover, synchronization is fully completed until t~500, which is much faster than in Model I (Figure2(a) and Figure 2(b)).However, oscillations continue to be asynchronous when d F = 0.01 and φ = 0.5, as shown in (c) and (d), and the system behavior results in the convergence to the fixed point when d F = 1.0 and φ = 0.1, as shown in (e) and (f).Simulation results of Model III are also presented in Figure 5(a) and Figure 5(b), where the transmembrane substance is S 3 , and both initial S 1 values and nine parameter values are simultaneously randomized.Temporal behaviors are also classified into three modes as well as in the case of Model II, namely, synchronous oscillations for d F = 1.0 and φ = 0.5 (Figure 5(a) and Figure 5(b)),asynchronous oscillations for d F = 0.01 and φ = 0.5 (not shown) and the convergence to the fixed point for d F = 1.0 and φ = 0.1 (not shown).Meanwhile in the case of Model IV, in which the coupling substance is S 2 , we confirmed no more than two modes of temporal behaviors.These are synchronous oscillations for d F = 1.0 and φ = 0.5 (Figure 5(c) and Figure 5(d)), andasynchronous oscillations for d F = 0.01 and φ = 0.5 (not shown).

Figure 5 (
e) and Figure 5(f) exhibit one of the examples that display asynchronous oscillations.Both initial S 1 values and nine parameter values are randomized as well as in Models II, III and IV.Needless to say, this does not necessarily exclude the possibility

Figure 4 .
Figure 4. Temporal changes in N 2 concentrations in Model II.Three figures on the left side show the temporal changes of central five cells, while three on the right side show those averaged for all cells.(a) and (b) show the synchronous oscillation, where d F = 1.0, φ = 0.5 and t = 500 -1000.Meanwhile, (c) and (d) show the asynchronous oscillation, where d F = 0.01, φ = 0.5 and t = 7500 -8000.Further, (e) and (f) show the convergence to the fixed point, where d F = 1.0, φ = 0.1 and t = 7500 -8000.The initial values of S 1 and the values of nine parameters, J 0 , k 2 , k 3 , k 4 , k 5 , k 6 , k 7 , r, and κ, are simultaneously randomized.Other parameter values are fixed in accordance with Table 2.It should be noted that the oscillations in Model II synchronize more rapidly, as shown in (a) and (b), compared with that in Model I (Figure 2).
(a) and Figure 2(b)) via the S 4 communicator is imperfect, while those of Models II (Figure 4(a) and Figure 4(b)), III (Figure 5(a) and Figure 5(b)) and IV (Figure 5(c) and Figure 5(d)) via the A 3 , S 3 and S 2 communicators are perfect.

Figure 5 .
Figure 5. Temporal changes in N 2 concentrations in Models III, IV and V. Three figures on the left side show the temporal changes of central five cells, while three on the right side show those averaged for all cells.(a) and (b) show the synchronous oscillation in Model III, where d F = 1.0, φ = 0.5 and t = 500 -1000.Moreover, (c) and (d) show the synchronous oscillation in Model IV, where d F = 1.0, φ = 0.5 and t = 500 -1000.Meanwhile, (e) and (f) show the asynchronous oscillation in Model V, where d F = 1.0, φ = 0.5 and t = 7500 -8000.The initial values of S 1 and the values of nine parameters, J 0 , k 2 , k 3 , k 4 , k 5 , k 6 , k 7 , r, and κ, are simultaneously randomized.Other parameter values are fixed in accordance with Table2.

Figures 6 (
a)-(d) are exactly the same as those of Figure 4(b), Figure 5(b), Figure 5(d) and Figure 5(f), respectively.Standard deviations of N 2 concentrations are divided by the mean value at each time for normalization, which guarantees the proper comparison between different models.

Figure 6 .
Figure 6.Temporal changes in normalized standard deviations of N 2 concentrations in Models II, III, IV and V.The standard deviations of N 2 concentrations are divided by the mean value at each time.Four figures correspond to Figure 4(b), Figure 5(b), Figure 5(d) and Figure 5(f), respectively.These normalized standard deviations are introduced to estimate the degree of synchronization, and the smaller value means more perfect synchronization.It is clear that synchronization in Model II is the most perfect, as shown in (a), where A 3 is the intercellular communicator.The last figure (d) shows the case of the asynchronous oscillation, where not only the values of standard deviations are high, but also no periodic structure is recognized.

Figure 4 (
e) and Figure 4(f)) to the synchronous oscillation (for example, Figure 4(a) and Figure 4(b)

Table 2 ,
which are obtained by non-dimensionalization processes such as

Table 1 .
Variables and parameters in Models I, II, III, IV and V.

Table 2 .
Parameter settings in Models I, II, III, IV and V.

Table 3 .
Fixed points and initial values of variables corresponding to synchronous oscillations in Models I, II, III and IV.
"FP" denotes the fixed point.For detailed calculations, see Appendix.

Table 2 .Table 4 .
Correlations among transmembrane communicators, randomized items and oscillation modes.oscillations are observed for randomization of initial S1 values in Model I. Meanwhile in Models II, III and IV, synchronous oscillations can also take place for randomization of both initial S1 values and values of nine parameters J0, k2, k3, k4, k5, k6, k7, r and κ.However, any synchronous oscillation is not detected in Model V when N2 is an intercellular mediator.Dimensionless time periods corresponding to synchronous oscillations are listed for reference. Synchronous