Modified Kuramoto Phase Model for Simulating Cardiac Pacemaker Cell Synchronization

Motohisa Osaka

Department of Basic Science, Nippon Veterinary and Life Science University, Tokyo, Japan.

DOI: 10.4236/am.2017.89092                

Abstract

It has been suggested that sick sinus syndrome, which is due to the dysfunction of the sinus node, may result from the sparser gap junctions and/or lower intrinsic frequencies of pacemaker cells that occur with aging. Hence, in this paper, the synchronization mechanism of pacemaker cells that lie in the sinus node of the heart is examined using a modified Kuramoto phase model. Although each element always interacts with all the others in the Kuramoto phase model, in the proposed model, each element interacts only with the neighbors over a certain time (called the interaction time) during Phase 4 of the action potential. The pacemaker cell elements are arranged on a square lattice, and each element connects with the elements surrounding it. The results indicate that the diversity of intrinsic frequencies of pacemaker cells may be necessary for synchronization. Moreover, increasing the proportion of invalid connections causes the elements to take more time to synchronize until eventually they do not synchronize at all, and decreasing the intrinsic frequencies of the elements prevents them from synchronizing. Probably these elucidate the cause of sick sinus syndrome.

Keywords

Synchronization, Sick Sinus Syndrome, Oscillator, Overdrive Suppression Test, Gap Junction

1. Introduction

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 [1] . 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 [2] . 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 [3] . 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 [4] . 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 [5] . 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 [6] . 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 [7] . 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 [5] . 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 [8] [9] . Although some mechanisms of this dysfunction have been described [10] , 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.

2. Model

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 ( $i=1,2,3,\cdots ,n$ ) is represented as follows:

$\begin{array}{c}\frac{\text{d}{\theta }_i}{\text{d}t}=f_i-\frac{1}{N}\underset{m=1}{\overset{N}{{\displaystyle \sum }}}{K}_{i{i}_m}\sin \left({\theta }_i-{\theta }_{i_m}\right)\\ \times \frac{1-\text{sign}\left(\left|\text{rem}\left({\theta }_i-{\theta }_{i_m},2\text{π}\right)\right|-\text{π}\times G_i\right)}{2}\\ \times \frac{1-\text{sign}\left(\text{rem}\left({\theta }_i,2\text{π}\right)+\text{π}\left(1-\text{sign}\left(\text{rem}\left({\theta }_i,2\text{π}\right)\right)\right)-\text{π}\times H_i\right)}{2};\end{array}$ (1)

In this equation, n is the total number of the elements of the lattice, q_i is the phase of the i-th element, f_i 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 $K_{i{i}_m}$ is the degree of interaction between the i-th and i_m-th elements. For example, i_m for the first element of Figure 1(c) are 2, 3, 4, and 7, and i_m for the fifth element is 2, 4, 6, and 8. The rem function is the remainder operation: $r=\text{rem}\left(a,b\right)$ returns the remainder after the division of a by b and the result has the same sign as dividend a. For example, $\text{rem}\left(-5\text{π},\text{2π}\right)=-\text{π}$ . Function y = sign(x) returns 1 if x > 0, 0 if x = 0, and −1 if x < 0. Moreover, G_i and F_i 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 q_i and ${\theta }_{i_m}$ , for any integer multiple of 2p, is less than p × G_i, 0 if it is larger than p × G_i, and 1/2 if it is equal to p × G_i. However, it is assumed never to equal to p × G_i 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 × G_i is the i-th element influenced by that neighbor (Figure 2). When

Figure 2. Schema of the action potentials of i-thelement and its i_m-thneighbor. Red line, refractory period; Blue line, Phase 4.q_i, phase of the i-th element (black circle); ${\theta }_{i_m}$ , phase of the i_m-th neighbor (red circle). A1, rem(q_i, 2p) ―p × G_i; A2, rem(q_i, 2p) + p × G_i.

$\text{rem}\left({\theta }_i,\text{2π}\right)\ge 0$ , the third line returns 1 if $\text{rem}\left({\theta }_i,\text{2π}\right)$ is less than p × H_i, 0 if it is larger than p × H_i, and 1/2 if it is equal to p × H_i. When $\text{rem}\left({\theta }_i,\text{2π}\right)<0$ , the third line returns 1 if $\text{rem}\left({\theta }_i,\text{2π}\right)+2\text{π}$ is less than p × H_i, 0 if it is larger than p × H_i, and 1/2 if it is equal to p × H_i. It is also assumed neither to be exactly equal to p × H_i 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 × H_i 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 [11] . 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 × G_i, any element influencing the i-th element is restricted to neighbors in or near Phase 4 correctly from $\text{rem}\left({\theta }_i,\text{2π}\right)-\text{π}\times G_i$ to $\text{rem}\left({\theta }_i,\text{2π}\right)+\text{π}\times G_i$ . The interaction time is defined as $\text{π}\times G_i\times \text{2}$ . 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, F_i, the intrinsic frequency of the i-th element ( $i=1,2,3,\cdots ,n$ ), is the sum of a common frequency F_c and random frequency $F_{r_i}$ . Common frequency F_c 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, $F_{r_i}$ is denoted as (0, 0.1) (Table 1, Table 2). Hence, $F_i=0.4+F_{r_i}$ .

Assumption 3: The degree of interaction between two neighbors varies marginally around a certain common degree. Hence, $K_{i{i}_m}$ , the degree of interaction between the i-th element and its i_m-th neighbor, is the sum of a common degree K_c and random degree $K_{r_{i_m}}$ . The random degree is given individually and is generated from a uniform distribution of random numbers between 0 and 1. Then, $K_{r_{i_m}}$ is denoted as (0, 1) (all Tables). Hence, $K_{i{i}_m}=K_c+K_{i{r}_{i_m}}$ .

Assumption 4: The duration of Phase 4 is represented as p×H_i, during which the i-th element interacts with the neighbors. The value of H_i is from 0 to 1. For the sake of model simplicity, it is the same for all elements. Hence, all H_i 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 [11] , H is assumed to be 1. The start

Table 1. Effects of F_r, G_c, 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 K_c, and lattice size on synchronization (SYNC).

cy: cycles/2000 time steps; u: unstable.

point of Phase 4 is defined as 2jp (Figure 2). Because $\text{rem}\left({\theta }_{i_m},2\text{π}\right)$ is between A1 and A2 in this figure, the i_m-th neighbor influences the i-th element. When $\text{rem}\left({\theta }_{i_m},2\text{π}\right)$ is not between A1 and A2, the i_m-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 × H_i, any neighbor of the i-th element is restricted to the neighbors with a phase between ${\theta }_i-\text{π}\times G_i$ and ${\theta }_i+\text{π}\times G_i$ for any integer multiple of 2p. The value $\text{π}\times G_i\times \text{2}$ is called the interaction time. The value of G_i is from 0 to 1, and varies marginally around a certain common value. Hence, G_i is the sum of a common value G_c and a random value $G_{r_i}$ . The random value is given individually and is generated from a uniform distribution of random numbers between 0 and 0.2. Hence, $G_{r_i}$ is denoted as (0, 0.2). Further, $G_i=G_c+G_{r_i}$ .

Assumption 6: The mean M_f and standard deviation SD_f 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 $\text{rem}\left({\theta }_1,\text{2π}\right)$ is employed as a reference peak. One cycle is 2,000/M_f time steps. The difference between the time step at the reference peak and a time step at each of the other peaks of $\text{rem}\left({\theta }_i,\text{2π}\right)$ , i ³ 2, is calculated. The standard deviation of these differences is divided by the time steps of one cycle. This quotient is denoted by SD_p (no units). When $0\le {\text{SD}}_{\text{f}}\le 0.\text{1}\times {\text{M}}_{\text{f}}$ and $0\le {\text{SD}}_{\text{p}}\le 1/8$ , 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 M_f × 0.8 and M_f × 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 M_f, SD_f, and SD_p. M_f and SD_f from 0 to 2000 time steps are 14.2 cy and 0.2 cy, respectively. Because the condition $0\le {\text{SD}}_{\text{f}}\le 0.\text{1}\times {\text{M}}_{\text{f}}$ is not satisfied, frequency-entrainment does not occur from 0 to 2000 time steps. Because M_f and SD_f from 4000 to 6000 time steps are 14.0 cy and 0.1 cy, respectively, the condition: $0\le {\text{SD}}_{\text{f}}\le 0.\text{1}\times {\text{M}}_{\text{f}}$ is satisfied and frequency-entrainment occurs. An arbitrary peak of $\text{rem}\left({\theta }_1,\text{2π}\right)$ 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 $\text{rem}\left({\theta }_i,\text{2π}\right)$ , i ³ 2, was calculated.

Figure 3. Example showing calculation of the mean M_f and standard deviation SD_f of the frequencies of nine elements as well as the standard deviation SD_p of the differences between the phases of the peaks. Here, M_f and SD_f are calculated over 2000 time steps (unit of frequency, cy: cycles/2000 time steps). The parameters are defined in Table 1 (F_c = 0.4, F_r = (0, 0.1), G_c = 0.2, G_r = (0, 0.2), K_c = 1, and K_r = (0, 1)). Each ordinate is $\text{rem}\left({\theta }_i,\text{2π}\right)$ ( $i=1,2,\cdots ,9$ ).

One cycle is 2000/14.0 = 142.9 time steps. Because SD_p = 0.08, the condition $0\le {\text{SD}}_{\text{p}}\le 1/8$ is satisfied. Hence, frequency- and phase-synchronization occurs.

The phases q_i ( $i=1,2,3,\cdots ,n$ ) 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 q_i ( $i=1,2,3,\cdots ,n$ ) 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 F_r, G_c, and lattice size

Table 1 summarizes the effects of F_r, G_c, and lattice size on the results. When all $F_{r_i}$ equal zero (that is, the frequencies of the elements are the same), synchronization of the elements depends on the lattice size. Specifically, when G_c = 0.2, the elements barely synchronize when the lattice size is small and the synchronization is unstable. In contrast, when G_c = 0.8, smaller lattices, even lattices of only nine elements, synchronize stably. Table 3 shows that when all $F_{r_i}$ equal zero, 121-element lattices do not synchronize, even when K_c = 3. In contrast, a 144-element lattice (12 × 12 elements) synchronizes when K_c < 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 F_c = 0.4, F_r = (0, 0.1), K_c = 2, K_r = (0, 1), %Kzero = 10, G_c = 0.2, and G_r = (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, F_i is the sum of common frequency F_c and random frequency $F_{r_i}$ . Several cases in which F_c is varied from 0.4 to 0 are examined (Table 2). When F_c is lower, the elements synchronize less unstably and ultimately do not synchronize at all. In other words, as the ratio of F_c to $F_{r_i}$ is lower, the elements hardly synchronize.

4. Discussion

The findings of the simulation and their interpretations are as follows:

1) When all $F_{r_i}$ equal zero and G_c is small in a small lattice, the elements do not synchronize. In contrast, even if all $F_{r_i}$ equal zero, the elements of a small lattice synchronize with large G_c. Assumption 5 states that the elements influencing the i-th element are restricted to neighbors with a phase between ${\theta }_i-\text{π}\times G_i$ and ${\theta }_i+\text{π}\times G_i$ 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: F_c = 0.4, F_r = (0, 0.1), K_c = 2, K_r = (0, 1), %Kzero = 10, G_c = 0.2, and G_r = (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 G_c 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 G_c 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 F_r does not equal zero or when G_c is larger. When F_r equals zero, a 121-element lattice does not synchronize even when K_c = 3. In contrast, a 144-element lattice synchronizes when K_c < 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 G_c = 0.2, K_c = 2, and %Kzero < 23. It has been reported that gap junctions are sparser in the sinoatrial node than in the surrounding atrial muscle [12] . 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 G_c 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 G_c is necessary for phase-entrainment. These considerations indicate that the conditions F_c = 0.4, F_r = (0, 0.1), K_c = 2, K_r = (0, 1), %Kzero = 10, G_c = 0.2, and G_r = (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 [8] ) 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 [8] [9] . 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 [10] . It has been reported that the remodeling of gap junctions is caused by structural heart diseases, such as ischemia or hypertrophy, and aging [13] . 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, G_c = 0.2, and K_c = 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 [14] . This phenomenon may be simulated by lower F_c. 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 [14] . 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 F_c to $F_{r_i}$ being lower, a high ratio of F_c to $F_{r_i}$ 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.

5. Conclusion

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.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Mirollo, R.E. and Strogatz, S.H. (1990) Synchronization of Pulse-Coupled Biological Oscillators. SIAM Journal on Applied Mathematics, 50, 1645-1662. https://doi.org/10.1137/0150098
[2]	Josephson, B.D. (1962) Possible New Effects in Superconductive Tunnelling. Physics Letters, 1, 251-253. https://doi.org/10.1016/0031-9163(62)91369-0
[3]	Glass, L. (2001) Synchronization and Rhythmic Processes in Physiology. Nature, 410, 277-284. https://doi.org/10.1038/35065745
[4]	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.
[5]	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. https://doi.org/10.1113/jphysiol.1984.sp015461
[6]	Michaels, D.C., Matyas, E.P. and Jalife, J. (1987) Mechanisms of Sinoatrial Pacemaker Synchronization: A New Hypothesis. Circulation Research, 61, 704-714. https://doi.org/10.1161/01.RES.61.5.704
[7]	Strogatz, S. (2000) From Kuramoto to Crawford: Exploring the Onset of Synchronization in Populations of Coupled Oscillators. Physica D: Nonlinear Phenomena, 143, 1-20. https://doi.org/10.1016/S0167-2789(00)00094-4
[8]	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. https://doi.org/10.1161/01.CIR.44.1.59
[9]	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
[10]	Gaffney, B.J., Wasserman, A.G., Rotsztain, A. and Rios, J.C. (1980) Sick Sinus Syndrome: Mechanisms and Management. Cardiovascular Clinics, 11, 7-25.
[11]	Ito, H., Ono, K. and Noma, A (1994) Background Conductance Attributable to Spontaneous Opening of Muscarinic K+ Channels in Rabbit Sino-Atrial Node Cells. Journal of Physiology, 476, 55-68.
[12]	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
[13]	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. https://doi.org/10.1177/002215540004801008
[14]	Josephson, M.E. (1993) Clinical Cardiac Electrophysiology: Techniques and Interpretations. 2nd Edition. Lea and Febeger, Philadelphia, PA.