Synchronization is a universal phenomenon associated with oscillations. The Moon revolves around the Earth with the same period as that of its rotation. The tide-producing force from the Earth to the Moon causes the synchronization of the revolution and rotation of the Moon. Moreover, the spontaneous synchronization of firefly lights is one of several extremely fascinating exhibitions that occur in nature. Some mathematical models have been proposed to simulate this phenomenon  . Each firefly is assumed to be a light-emitting oscillator that generates a limit cycle. Another example is the synchronization of a series of Josephson junctions, each element of which interacts only with its close neighbors  . Similarly, each pacemaker cell of the sinoatrial node (a group of pacemaker cells in the wall of the right atrium of the heart), which is regarded as an electrical oscillator generating a limit cycle, is loosely coupled with its neighboring pacemaker cells so that those cells synchronize. The electrical impulse current of a single pacemaker cell is not sufficient to travel through the impulse-conducting system, so the ventricles of the heart do not contract. In contrast, synchronized pacemaker cells generate an impulse with sufficiently high current to travel through the system so that the ventricles contract and pump out blood normally. Hence, the failure of this synchronization is presumed to cause sinoatrial arrest  . This interruption of the cardiac cycle generally lasts a few seconds before the atrioventricular junction lying in the middle of the impulse-conducting system begins pacing and restores slower ventricular contractions. Sinoatrial arrest is a part of sick sinus syndrome, the symptoms of which include fainting, vertigo, and weakness. Gap junctions in the sinoatrial node have recently been demonstrated to be a key player in the electrical coupling underlying synchronization  . Because gap junctions have low resistance, local circuit currents propagate over short distances. Jalife proposed a very probable hypothesis called the democratic consensus hypothesis, which is based on an experiment using rabbit sinoatrial pacemaker cells  . This hypothesis states that although the individual pacemaker cells in the sinoatrial node beat at slightly different intrinsic frequencies, they interact mutually by electrical coupling to fire at a “consensus” rate. It has been suggested that when two independent groups of fast and slow pacemaker cells are connected through low-resistance junctions, the period resulting from their mutual entrainment should be a function of their respective intrinsic frequencies, their phase relations, and the degree of electrical coupling. Moreover, Jalife et al. elucidated the mechanisms of sinoatrial pacemaker synchronization using a computer simulation of 81 to 225 coupled cells  . Pacemaker activity has been simulated using differential equations that describe transmembrane ionic currents. These results support the hypothesis that sinoatrial node synchronization occurs through a “democratic” process resulting from the phase-dependent interactions of thousands of pacemakers.
The Kuramoto phase model has been used to simulate the synchronization of firefly signals  . Although this firefly signal model is simple, the respective frequency, phase relation, and degree of coupling are included as variables in the model. The responsiveness of pacemaker cells is rather different from that of fireflies in that the action potential of each pacemaker cell has a refractory period, during which it does not respond to external stimuli at all. However, a pacemaker cell can respond to external stimuli during Phase 4 of the action potential (described in detail below). Jalife reported various patterns of interactions between fast and slow pacemaker cells, that is, simple harmonic (e.g., 1:1, 2:1, and 1:2) and more complex (e.g., 3:2 and 5:4) ratios  . Because such complex phenomena have not been reported in the observation of firefly signals, we speculate that these various ratios may be due to the refractory period. Specific ionic currents, over time, slowly depolarize Phase 4 for pacemaking. Hence, it is assumed that each pacemaker cell is influenced by its neighbors, the phases of which are within a certain range of its phase. This range is called the interaction time. We modified the Kuramoto phase model by incorporating a variable that represents the period of Phase 4. In the modified Kuramoto phase model, for each pacemaker cell, the intrinsic frequency, Phase 4 period, phase relation, and degree of coupling between it and the neighboring pacemaker cells can be modulated independently as variables to observe various patterns. The differences between the Kuramoto phase model and the proposed model are that although each element always interacts with all the others in the former, each element interacts only with its neighbors during the Phase 4 period in the latter. We examine how the synchronization of pacemaker cells depends on those variables. Such an examination is presumed to be rather difficult using differential equations describing transmembrane ionic currents.
The repetitive high-frequency stimulation test of the sinoatrial node (called the overdrive suppression test) is used to examine its function clinically. The sinoatrial node comes to a standstill immediately after the repetitive stimulation, then resumes a regular rhythm after a certain pause (called the sinus node recovery time). The length of this pause is presumed to reflect the degree of dysfunction of the sinoatrial node   . Although some mechanisms of this dysfunction have been described  , the factor that determines the length of the sinus node recovery time is unknown. One of the purposes of the present study is to infer that factor using the modified Kuramoto phase model. Additionally, an unexpected finding is that slightly different intrinsic frequencies of the pacemaker cell elements promote their synchronization, although it was expected that the same frequencies would make those elements synchronize more easily.
The simplest model consists of two connected elements, as illustrated in Figure 1(a). Figure 1(b) shows the minimum square lattice in which there is at least one element that is unconnected to each element. For example, the first element does not connect with the fourth element. Figure 1(c) shows a square lattice consisting of nine elements, which is the minimum lattice in which each element has four connections. For example, the elements are numbered in Figure 1(c): the element on the left corner of the first row is element 1, the element of the right corner of the first row is element 3, and the element of the right corner of the last row is element 9. Every element interacts with four elements: every element inside the lattice connects with its immediate neighbors, and everyone on the lattice sides connects with its immediate neighbors and elements on the op-
Figure 1. Arrangement of coupling resistances between simulated pacemaker cells in two-dimensional square lattices. The resistances are expressed as arrows. (a) 2 elements; (b) 2 × 2 elements; and (c) 3 × 3 elements. Every element interacts with four elements, and (c) shows that every element on the sides connects with its immediate neighbors and elements on the opposite side.
posite side. Because every element is equivalent to the others with respect to position, no boundary condition is necessary. Hence, this two-dimensional lattice is assumed to be the surface of a three-dimensional torus made by connecting the top side with the bottom and the left side with the right. The elements directly connected in this manner are regarded as neighbors. In the modified Kuramoto phase model, the dynamics of the i-th cell ( ) is represented as follows:
In this equation, n is the total number of the elements of the lattice, qi is the phase of the i-th element, fi is the intrinsic frequency of the i-th element, N is the number of couplings with the i-th element (N = 1 in Figure 1(a), 2 in Figure 1(b), and 4 in Figure 1(c)), and is the degree of interaction between the i-th and im-th elements. For example, im for the first element of Figure 1(c) are 2, 3, 4, and 7, and im for the fifth element is 2, 4, 6, and 8. The rem function is the remainder operation: returns the remainder after the division of a by b and the result has the same sign as dividend a. For example, . Function y = sign(x) returns 1 if x > 0, 0 if x = 0, and −1 if x < 0. Moreover, Gi and Fi are coefficients between 0 and 2. The phase of every cycle is 2p.
The first line of the equation is just the Kuramoto phase model for the interactions between each element and its neighbors. The second line returns 1 if the difference between qi and , for any integer multiple of 2p, is less than p × Gi, 0 if it is larger than p × Gi, and 1/2 if it is equal to p × Gi. However, it is assumed never to equal to p × Gi exactly. Hence, the second line means that only if the difference between the phase of the i-th element and that of any neighbor is within p × Gi is the i-th element influenced by that neighbor (Figure 2). When
Figure 2. Schema of the action potentials of i-thelement and its im-thneighbor. Red line, refractory period; Blue line, Phase 4.qi, phase of the i-th element (black circle); , phase of the im-th neighbor (red circle). A1, rem(qi, 2p) ―p × Gi; A2, rem(qi, 2p) + p × Gi.
, the third line returns 1 if is less than p × Hi, 0 if it is larger than p × Hi, and 1/2 if it is equal to p × Hi. When , the third line returns 1 if is less than p × Hi, 0 if it is larger than p × Hi, and 1/2 if it is equal to p × Hi. It is also assumed neither to be exactly equal to p × Hi nor to be < 0. Thus, the third line means that only if the phase of the i-th element, for any integer multiple of 2p, is between 0 and p × Hi is the i-th element influenced by the neighbors. Figure 2 shows trains of action potentials for two pacemaker cells. The action potential mainly consists of Phases 0, 3, and 4, because Phases 1 and 2 are small. Phase 0 is the period of rapid depolarization caused by a fast inflow of calcium ions and Phase 3 consists of repolarization caused by a fast outflow of potassium ions. These periods make up the refractory period. Phase 4 is the period during which the inflow of sodium ions begins, and thereafter, the inflow of calcium ion continues before firing beyond the threshold  . Because the i-th element interacts with the neighbors during Phase 4 and the ionic current through the gap junctions depends on the phase, it is presumed that elements interact with each other when the neighbors are in or near Phase 4. Using p × Gi, any element influencing the i-th element is restricted to neighbors in or near Phase 4 correctly from to . The interaction time is defined as . This is the interval between A1 and A2 in Figure 2.
The assumptions are summarized as follows:
Assumption 1: The elements are independent oscillators. In the Kuramoto phase model, each element interacts with all the others. In contrast, in the proposed model, each element interacts only with the connecting elements (neighbors).
Assumption 2: The frequency of each element varies marginally around a certain common frequency. Hence, Fi, the intrinsic frequency of the i-th element ( ), is the sum of a common frequency Fc and random frequency . Common frequency Fc is fixed as 0.4 arbitrarily. The random frequency is given individually and is generated from a uniform distribution of random numbers between 0 and 0.1. Then, is denoted as (0, 0.1) (Table 1, Table 2). Hence, .
Assumption 3: The degree of interaction between two neighbors varies marginally around a certain common degree. Hence, , the degree of interaction between the i-th element and its im-th neighbor, is the sum of a common degree Kc and random degree . The random degree is given individually and is generated from a uniform distribution of random numbers between 0 and 1. Then, is denoted as (0, 1) (all Tables). Hence, .
Assumption 4: The duration of Phase 4 is represented as p×Hi, during which the i-th element interacts with the neighbors. The value of Hi is from 0 to 1. For the sake of model simplicity, it is the same for all elements. Hence, all Hi are expressed as H. The phase of one cycle is 2p. Because the duration of Phase 4 (=p × H) is approximately one half of one cycle  , H is assumed to be 1. The start
Table 1. Effects of Fr, Gc, and lattice size on synchronization (SYNC).
cy: cycles/2000 time steps; u: unstable.
Table 2. Effects of lower common frequency on synchronization (SYNC).
cy: cycles/2000 time steps; u: unstable.
Table 3. Effects of Kc, and lattice size on synchronization (SYNC).
cy: cycles/2000 time steps; u: unstable.
point of Phase 4 is defined as 2jp (Figure 2). Because is between A1 and A2 in this figure, the im-th neighbor influences the i-th element. When is not between A1 and A2, the im-th neighbor does not influence the i-th element.
Assumption 5: When the phase of the i-th element, for any integer multiple of 2p, is between 0 and p × Hi, any neighbor of the i-th element is restricted to the neighbors with a phase between and for any integer multiple of 2p. The value is called the interaction time. The value of Gi is from 0 to 1, and varies marginally around a certain common value. Hence, Gi is the sum of a common value Gc and a random value . The random value is given individually and is generated from a uniform distribution of random numbers between 0 and 0.2. Hence, is denoted as (0, 0.2). Further, .
Assumption 6: The mean Mf and standard deviation SDf of the frequencies of all elements are calculated over a period 2000 time steps in length (unit of frequency, cy: cycles/2000 time steps) (Figure 3). An arbitrary peak of is employed as a reference peak. One cycle is 2,000/Mf time steps. The difference between the time step at the reference peak and a time step at each of the other peaks of , i ³ 2, is calculated. The standard deviation of these differences is divided by the time steps of one cycle. This quotient is denoted by SDp (no units). When and , it is considered that the elements synchronize with respect to frequency and phase (frequency- and phase-synchronization). This means that about 95% of frequencies are between Mf × 0.8 and Mf × 1.2 and about 95% of the peak phases are between the phase of the reference peak −p/2 and the phase of the reference peak +p/2 because the phases are presumed to be distributed normally. Hence, this model is believed to simulate simultaneous firing of pacemaker cells, which generates an impulse traveling through the stimulating conducting system. Figure 3 shows an example of calculation of Mf, SDf, and SDp. Mf and SDf from 0 to 2000 time steps are 14.2 cy and 0.2 cy, respectively. Because the condition is not satisfied, frequency-entrainment does not occur from 0 to 2000 time steps. Because Mf and SDf from 4000 to 6000 time steps are 14.0 cy and 0.1 cy, respectively, the condition: is satisfied and frequency-entrainment occurs. An arbitrary peak of is employed as a reference peak (red circle). The difference between the time step at the reference peak and a time step at each of the other peaks (blue circles) of , i ³ 2, was calculated.
Figure 3. Example showing calculation of the mean Mf and standard deviation SDf of the frequencies of nine elements as well as the standard deviation SDp of the differences between the phases of the peaks. Here, Mf and SDf are calculated over 2000 time steps (unit of frequency, cy: cycles/2000 time steps). The parameters are defined in Table 1 (Fc = 0.4, Fr = (0, 0.1), Gc = 0.2, Gr = (0, 0.2), Kc = 1, and Kr = (0, 1)). Each ordinate is ( ).
One cycle is 2000/14.0 = 142.9 time steps. Because SDp = 0.08, the condition is satisfied. Hence, frequency- and phase-synchronization occurs.
The phases qi ( ) are calculated using MATLAB® to solve the differential equations. The calculation precision depends on the time step. Although 2 time steps, 0.1 and 0.01, have been used preliminarily, both have given the same results to calculate frequency. Hence, time step 0.1 is selected. It is needed that the time span is long enough to examine whether the elements synchronize or not. The time span is from 0 to 2000 (if necessary, 6,000). Hence, the data length of qi ( ) is 20,000 (or 60,000) time steps. The initial condition is a set of uniformly distributed random numbers in the interval (0, 5).
3. Simulation Results
1) Effects of Fr, Gc, and lattice size
Table 1 summarizes the effects of Fr, Gc, and lattice size on the results. When all equal zero (that is, the frequencies of the elements are the same), synchronization of the elements depends on the lattice size. Specifically, when Gc = 0.2, the elements barely synchronize when the lattice size is small and the synchronization is unstable. In contrast, when Gc = 0.8, smaller lattices, even lattices of only nine elements, synchronize stably. Table 3 shows that when all equal zero, 121-element lattices do not synchronize, even when Kc = 3. In contrast, a 144-element lattice (12 × 12 elements) synchronizes when Kc < 3.
2) Effect of fewer connections between neighbors
Each element connects with four neighbors. For example, a 144-element lattice has 290 different connections. Here, the percentage of invalid connections with respect to total connections is denoted as %Kzero. When Fc = 0.4, Fr = (0, 0.1), Kc = 2, Kr = (0, 1), %Kzero = 10, Gc = 0.2, and Gr = (0, 0.2), a 144-element lattice synchronizes with a %Kzero of less than 23 (Figure 4). As %Kzero increases, it takes the elements longer to synchronize (a longer delay).
3) Effects of a lower common frequency on synchronization
Assumption 2 states that the frequency of each element varies marginally around a certain common frequency. Hence, Fi is the sum of common frequency Fc and random frequency . Several cases in which Fc is varied from 0.4 to 0 are examined (Table 2). When Fc is lower, the elements synchronize less unstably and ultimately do not synchronize at all. In other words, as the ratio of Fc to is lower, the elements hardly synchronize.
The findings of the simulation and their interpretations are as follows:
1) When all equal zero and Gc is small in a small lattice, the elements do not synchronize. In contrast, even if all equal zero, the elements of a small lattice synchronize with large Gc. Assumption 5 states that the elements influencing the i-th element are restricted to neighbors with a phase between and for any integer multiple of 2p. The findings indicate
Figure 4. Effect of fewer connections between neighbors. %Kzero: the percentage of invalid connections with respect to total connections. The ordinate is a delay in synchronization. The conditions of parameters: Fc = 0.4, Fr = (0, 0.1), Kc = 2, Kr = (0, 1), %Kzero = 10, Gc = 0.2, and Gr = (0, 0.2). The delay becomes longer exponentially beyond %Kzero = 10. The elements eventually become unable to synchronize at %Kzero = 23.
that elements with the same frequencies barely synchronize when Gc is small. When the lattice size is bigger, the cells synchronize. In the initial condition, the phases of the elements are uniformly random. Because an element interacts only with the neighbors with a similar phase, the probability that it is influenced by neighbors decreases as Gc decreases. Hence, when the frequencies (that is, phases) of the elements are more varied, this probability increases. This suggests that the diversity of frequencies (or phases) of the elements is necessary for the elements to synchronize. This suggestion is supported by the fact that the elements of smaller lattices synchronize easily when Fr does not equal zero or when Gc is larger. When Fr equals zero, a 121-element lattice does not synchronize even when Kc = 3. In contrast, a 144-element lattice synchronizes when Kc < 3. These findings suggest that the diversity of frequencies is more effective for element synchronization than the strength of interaction between neighboring elements.
2) The synchronization of a 144-element lattice occurs when Gc = 0.2, Kc = 2, and %Kzero < 23. It has been reported that gap junctions are sparser in the sinoatrial node than in the surrounding atrial muscle  . Hence, this report suggests that %Kzero, which corresponds to the percentage of invalid gap junctions, is greater than zero. In the present study, the synchronization for generating an impulse traveling through the stimulating conducting system is assumed to be frequency- and phase-entrainment. When Gc is large, each element can be affected by the phases of its neighbors that are rather different from its phase. In this case, although frequency entrainment occurs easily, phase-entrainment does not also occur. Hence, small Gc is necessary for phase-entrainment. These considerations indicate that the conditions Fc = 0.4, Fr = (0, 0.1), Kc = 2, Kr = (0, 1), %Kzero = 10, Gc = 0.2, and Gr = (0, 0.2) for the elements, are likely conditions for simulating pacemaker cells. If %Kzero is increased, it takes the elements more time to synchronize (Figure 4). Since the delay becomes longer exponentially beyond %Kzero = 10, the elements will not synchronize abruptly at %Kzero = 23. It suggests that the maximum sinus node recovery time (clinically about 1.5 s  ) appears around %Kzero = 23.
The overdrive suppression test is used clinically to examine sinoatrial node functionality. In this test, the sinoatrial node comes to a standstill immediately after a high-frequency repetitive stimulation and then restores a regular rhythm after the sinus node recovery time, which indicates the degree of sinoatrial node dysfunction   . Although some mechanisms of its dysfunction have been described, the common final pathway in most instances seems to be chronic fibrosis of the sinus node  . It has been reported that the remodeling of gap junctions is caused by structural heart diseases, such as ischemia or hypertrophy, and aging  . If these causes make gap junctions sparser, sinus node recovery time will lengthen and eventually the sinoatrial node will generate no impulses that can travel through the stimulating conducting system. This phenomenon (sick sinus syndrome) may be postulated from the results that it takes the elements more time to synchronize when %Kzero increases (Figure 4).
3) The results also suggest that if the common frequency is lower, the elements do not synchronize when %Kzero = 10, Gc = 0.2, and Kc = 2. These conditions are selected because the second finding above indicates that they are suitable for simulating pacemaker cells. It has been reported that the intrinsic frequencies of pacemaker cells decreases with aging  . This phenomenon may be simulated by lower Fc. It is presumed to be one of mechanisms that cause the sinoatrial node to be unable to synchronize steadily and could be an age-related cause of sick sinus syndrome. On the other hand, intrinsic heart rate of normal subjects is rather fast  . The reason for this was for a long time unknown. Heart rate is controlled within a normal range of 60 to 100 beats per minute under the autonomic nervous control. Since the elements hardly synchronize with the ratio of Fc to being lower, a high ratio of Fc to is necessary for synchronization. Probably, this is the reason that intrinsic heart rate of normal subjects is rather fast.
4) Study limitations
The largest lattice size is 144 elements (12 × 12 elements). It takes about six hours to solve 144 differential equations using MATLAB®. It takes too long to calculate a lattice with 225 elements (15 × 15 elements) in practice. Hence, the maximum lattice size in the present study is 144 elements. This model is described using some parameters for frequency, the period of Phase 4, and interaction. Although values for these parameters that simulate pacemaker cells were determined from many trial-and-error experiments, the validity of these values should be investigated in the future.
In this study, the Kuramoto phase model was modified by incorporating the interaction time of Phase 4, during which each element can interact with its neighbors, as a variable. The results are as follows: 1) Certain values for the frequency, interaction time, and degree of interaction are found to simulate pacemaker cells; 2) Diversity in the intrinsic frequencies of the elements promotes their synchronization, although the same frequencies should be easier to synchronize; 3) Increasing the proportion of invalid connections in the model (which corresponds physiologically to sparser gap junctions) causes the elements to take longer to synchronize and eventually become unable to synchronize at all; 4) Decreasing the intrinsic frequencies of the elements prevents them from synchronizing. These results indicate a possible mechanism for the age-related causes of sick sinus syndrome.
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.
 Verheule, S., van Kempen, M.J.A., Postma, S., Rook, M.B. and Jongsma, H.J. (2001) Gap Junctions in the Rabbit Sinoatrial Node. American Journal of Physiology—Heart and Circulatory Physiology, 280, H2103-H2115.
 Jalife, J. (1984) Mutual Entrainment and Electrical Coupling as Mechanisms for Synchronous Firing of Rabbit Sino-Atrial Pace-Maker Cells. Journal of Physiology, 356, 221-243.
 Strogatz, S. (2000) From Kuramoto to Crawford: Exploring the Onset of Synchronization in Populations of Coupled Oscillators. Physica D: Nonlinear Phenomena, 143, 1-20.
 Mandel, W., Hayakawa, H., Danzig, R. and Marcus, H.S. (1971) Evaluation of Sino-Atrial Node Function in Man by Overdrive Suppression. Circulation, 44, 59-56.
 Mandel, W., Hayakawa, H., Allen, H.N., Danzig, R. and Kermaier, A.I. (1972) Assessment of Sinus Node Function in Patients with the Sick Sinus Syndrome. Circulation, 46, 761-769. https://doi.org/10.1161/01.CIR.46.4.761
 Masson-Pevet, M.A., Bleeker, W.K. and Gros, D. (1979) The Plasma Membrane of Leading Pacemaker Cells in the Rabbit Sinus Node: A Quantitative Ultrastructural Analysis. Circulation Research, 45, 621-629. https://doi.org/10.1161/01.RES.45.5.621
 Yeh, H.I., Chang, H.M., Lu, W.W., Lee, Y.N., Ko, Y.S., Severs, N.J. and Tsai, C.H. (2000) Age-Related Alteration of Gap Junction Distribution and Connexin Expression in Rat Aortic Endothelium. Journal of Histochemistry & Cytochememistry, 48, 1377-1389.