Modified Kuramoto Phase Model for Simulating Cardiac Pacemaker Cell Synchronization

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.


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 be-M.Osaka DOI: 10.4236/am.2017.890921229 Applied Mathematics tween 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.

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- ( ( ) In this equation, n is the total number of the elements of the lattice, θ 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 θ , for any integer multiple of 2π, is less than π × G i , 0 if it is larger than π × G i , and 1/2 if it is equal to π × G i .However, it is assumed never to equal to π × 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 π × G i is the i-th element influenced by that neighbor (Figure 2).When is less than π × H i , 0 if it is larger than π × H i , and 1/2 if it is equal to π × H i .It is also assumed neither to be exactly equal to π × 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 2π, is between 0 and π × 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 π × G i , any element influencing the i-th element is restricted to neighbors in or near Phase 4 correctly from ( ) The interaction time is defined as π 2 i G × × .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 F is denoted as (0, 0.1) (Table 1, Table 2).Hence, 0.
Assumption 3: The degree of interaction between two neighbors varies marginally around a certain common degree.Hence, m ii K , 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 i m r K .The random degree is given individually and is generated from a uniform distribution of random numbers between 0 and 1.
Then, i m r K is denoted as (0, 1) (all Tables).Hence, Assumption 4: The duration of Phase 4 is represented as π×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 2π.Because the duration of Phase 4 (=π × H) is approximately one half of one cycle [11], H is assumed to be 1.The start  point of Phase 4 is defined as 2jπ (Figure 2).Because 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 2π, is between 0 π × H i , any neighbor of the i-th element is restricted to the neighbors with a phase between π for any integer multiple of 2π.The value π 2 i G × × 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 i r G .The random value is given individually and is generated from a uniform distribution of random numbers between 0 and 0.2.Hence, i r G is denoted as (0, 0.2).Further, 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 ( ) rem , 2π θ 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 ( ) 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 , 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 −π/2 and the phase of the reference peak +π/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 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: is satisfied and frequency-entrainment occurs.An arbitrary peak of ( ) rem , 2π θ 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 ( ) One cycle is 2000/14.0= 142.9 time steps.Because SD p = 0.08, the condition p 0 SD 1 8 ≤ ≤ is satisfied.Hence, frequency-and phase-synchronization occurs.
The phases θ i ( 1, 2, 3, , i n =  ) are calculated using MATLAB ® to solve the diffe- rential 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 θ i ( 1, 2, 3, , i n =  ) is 20,000 (or 60,000) time steps.The initial condi- tion is a set of uniformly distributed random numbers in the interval (0, 5).

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 i r F equal zero (that is, the frequencies of the elements are the same), syn- M. Osaka DOI: 10.4236/am.2017.890921235 Applied Mathematics chronization 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 i r F 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.
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 i r F .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 i r F is lower, the elements hardly synchronize.

Discussion
The findings of the simulation and their interpretations are as follows: 1) When all 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 i r F being lower, a high ratio of F c to i r F 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.

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.

Figure 1 .
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.
(a), 2 in Figure1(b), and 4 in Figure1(c)), and m ii K is the degree of interaction between the i-th and i m -th elements.For example, i m for the first element of Figure1(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: 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, G i and F i are coefficients between 0 and 2. The phase of every cycle is 2π.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 θ i and m i

( 1 F
is the sum of a common frequency F c and random frequency i r .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, i r in this figure, the i m -th neighbor influences the i-th element.When

Figure 3 .
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 ( ) rem , 2π i θ ( 1, 2, , 9 i =  ).

Table 1 .
Effects of F r , G c , and lattice size on synchronization (SYNC).

Table 2 .
Effects of lower common frequency on synchronization (SYNC).

Table 3 .
Effects of K c , and lattice size on synchronization (SYNC).