Research Article

The Generalized Model N - Pacemaker Curve Phase Response of the Atria, Ventricular Fibliration and AB - Blockade

Sergey Belyakin¹* and Sergey Shuteev²

¹ Department of General Physics, Lomonosov Moscow State University, Russia

² Laboratory of dynamic systems, Lomonosov Moscow State University, Russia

Submission:February 10, 2020; Published:March 13, 2020

*Corresponding author:Sergey Belyakin, Department of General Physics, Lomonosov Moscow State University, Moscow, Russia

How to cite this article:Sergey Belyakin, Sergey Shuteev. The Generalized Model N - Pacemaker Curve Phase Response of the Atria, Ventricular Fibliration and AB - Blockade. Glob J Nano. 2020; 5(2): 555658. DOI: 10.19080/GJN.2020.05.555658

Abstract

In this publication, we generalize the proposed model of two interacting oscillators in the case of a strong difference in their periods (when the pacemaker pulses do not alternate) and propose a General model describing a network of oscillators coupled globally. Our goal is to make the model as simple as possible and enter the minimum number of parameters. Therefore, we will fully characterize the pacemaker of their internal lengths of the cycle and re-present them as pulse oscillators. Interaction of pacemakers is described by PRC.

Keywords: The Function Phase Response Curve; Pacemaker; AB-blockade; The Neural Network Convolution; Magnitude; Multiplicity Phases; Pacemakers

Introduction

The circle map curve the phase response and the Arnold tongues of

Consider some physical quantity ξ which reflects the internal state of the biological oscillator. Let the eigenfrequency of the oscillator be equal T0. Let’s call a marker any event that can be clearly seen in the experiment, which is reached by the value ξ only once per period. Such a marker may be, for example, the beginning of the action potential in the cardiac preparation. Let’s define the oscillator phase as follows. The phase of an arbitrarily selected marking event (for example, the maximum value of ξ) is assumed to be zero. At any next time t, 0 < t < T0, the phase is defined as φ = t ∕ T0 (mod1). Since the rhythm is restored after the perturbation of the system, the introduced phase completely determines the state of the system. Suppose that an external periodic perturbation acts on a nonlinear oscillator. Then each external influence shifts the state of the system to a new state (1):

The function f (φn) is called the phase response curve (PRC) [1] and determines the phase change after the stimulus. It is convenient to represent the points f (φn) of the system state lying on the circle of the unit radius. Then, by iterating the mapping (1), one point of the circle is converted to another point of the same circle. If the circle map is continuous, then it can be characterized by a number called the topological degree sand equal to the number of passes through φn+1 the unit circle during f (φn) the time it passes once. In periodic perturbations of self-oscillations with a stable limit cycle, the dynamics is often described by maps of a circle with a topological degree 0 (when the over-threshold response gives rise to a new cycle) or 1 (which expresses a sub-threshold response to stimulation). The different types of circle maps are shown in Figure 1a-d. Along with the topological degree, an important characteristic of the circle display is the number of rotations. We define it as the time average ratio of the external perturbation period to the period of the perturbed oscillator. If the rotation number is rational, ρ = M/N (here M is the number of cycles of the stimulator, and N is the number of cycles of the nonlinear oscillator), then the dynamics of the system will be periodic with the capture of the multiplicity phase N/M. If the rotation number is irrational, the system demonstrates quasi-periodic or chaotic behavior. In many cases, the disturbance by a single pulse of a spontaneously oscillating system leads to a phase shift of the current rhythm (see, for example, [2] and references there). The magnitude of the shift depends on both the magnitude of the stimulus and its phase in the cycle. The graph of the dependence of the new phase on the previous phase (i.e., PRC) is either a continuous circle map with a topological degree of 1 or 0, or a discontinuous function. Phase shift experiments were performed for a large number of different systems. We are interested, first of all, in the phase response curve, experimentally obtained in the study of cardiac drug.

In [3] the duration of the cycle of spontaneous oscillations of Purkinje fibers after stimulation by short pulses of electric current was measured. The obtained phase response curve of the Biphase form is shown in Figure 2. Based on the study of this experimental material, the following generalizations can be made [3]. After the disturbance, the rhythm is usually restored (after the transition process) with the same frequency and amplitude as before the disturbance, and its phase is shifted. Depending on the phase, a single stimulus can result in either an elongation (early stimulus) or a shortening (late stimulus) of the duration of the perturbed cycle. At some amplitudes of the stimulus, obvious discontinuities are observed. To further study the dynamics of any constructed model, it is necessary, having experimentally obtained PRC, to find a good analytical approximation of this curve. This will allow to investigate the main features of the behavior of the system. The main characteristic of the desired function is the need to directly depend on only two physical parameters: the amplitude of the stimulus and the phase of the applied perturbation. All other (socalled “internal”) parameters describing the course of the curve should (ideally) be reduced to these two. One of the simplest (and coarsest) approximations of a given PRC is the sinusoidal function, which ultimately results in a map of the form (2):

Where a and b are constants. However, despite its simplicity, this approximation correctly reflects the qualitative structure of the phase portrait of the system under study. The analysis of bifurcations of reversible circle maps was undertaken in the last century by A. Poincare and still attracts much -by V. I. Arnold [4] (see also [5] and the references given there). For (Figure 3) the bifurcation diagram of the circle diffeomorphism on the parameter plane (b, a) is shown. This diagram is divided into areas called language (or horns) of Arnold, which correspond to the sustainable capture phase ratio N/M (i.e., N cycles of the stimulator has M cycles of a nonlinear oscillator). Arnold languages exist for all rational relations N/M, where N and M are mutually Prime numbers. This means that there are an infinite number of Arnold languages that correspond to all possible ratios of frequencies of the stimulator and the perturbed oscillator. Between any two languages corresponding to N/M and N*/M* phase captures, there is another capture region corresponding to the capture of multiplicity phases (N+N*)/(M+M*). The structure shown in Figure 3, is the usual behavior for low stimulus amplitudes in simple theoretical models discussed below. However, as the amplitude of the periodic effect increases, this structure collapses. Relaxation model of Poincare oscillator. A widely used idealization of some periodically stimulated oscillators is the Gelfand and cetlin model [6], or the relaxation model [2,7-9].

In this model, the value referred to as activity increases to the upper threshold, leading to some event. Then the activity returns to the lower threshold. If the rates of rise and fall of activity to the thresholds are fixed, and the thresholds are also fixed, then a periodic sequence of events is generated, the frequency of which is easy to determine. Periodic perturbation in relaxation models can be included in the form of threshold modulation, usually sinusoidal. In some works, instead of sinusoidal modulation of the threshold, other functions were considered, for example, Delta function peaks, rectangular and triangular pulses, etc. (see references in) [10]. Arnold in [4] briefly discussed the possibility of using obtained model in the relaxation mapping to study the rhythms Wenkebach. Subsequent researchers found that piecewise linear monotone, discontinuous maps (Figure 1c), similar to those found in the relaxation model, appear in theoretical models of atrial-ventricular communication in AB blockade [11,12]. Such mappings can be experimentally measured and used to predict complex rhythms observed in humans [13].

Despite the wide application, the relaxation model too simplistically describes the interaction of the oscillator and the external perturbation; much more vital is the use of models that consider the individual response of the system to the external perturbation. Cardiologists usually assume that the nonlinear ODE that contains oscillation with a stable limit cycle, represent a suitable model for the generation of periodic activity of the heart [14]. In the case where the limit cycle is quickly achieved after a single stimulus, and the action of a single stimulus is known, it is possible to calculate the effect of periodic stimulation. The prototype of the model with a periodically perturbed limit cycle is the vander Pol equation with a sinusoidal perturbation. Consider the effect of a periodic sequence of short pulses on the oscillations described by a system with a limit cycle (see, for example) [10]. The simplest model is Poincare oscillator. In this model, a stable limit cycle is a circular trajectory. The perturbation is a horizontal displacement of magnitude b, and after stimulation, the system rapidly approaches the limit cycle along its radius. If it is the phase φn immediately preceding the nth stimulus, then the phase preceding (n+1) the nth stimulus is simply φn+1 = τ + g (φn, b) (mod1). Where τ is the time interval between periodic stimuli normalized for the eigenfrequency of the autogenerator, and PRC g(φ) is easily calculated [14-16].

This theoretical model for periodically perturbed limit cycles was independently proposed by several researchers [15-18]. Since for a simple model with a limit cycle RPC is calculated quite easily, it is possible to use analytical and numerical methods to determine the detailed structure of the phase capture zones as a function of the amplitude b and frequency a of the stimulus. In this example, for low stimulus amplitudes (b ≤ 1), the capture zone topology has a classical Arnold structure (Figure 3), and the circle display is a reversible degree 1 display. However, for b > 1, the dynamics is described by displaying a circle of zero topological degree. The extensions of Arnold’s languages have a more complex form. There are bifurcations in the system, leading to doubling of the period and chaos. This model, although different from the exact electrophysiological models, nevertheless, surprisingly well reproduces many features observed in experiments on the study of periodically perturbed aggregates of some heart cells [19]. For example, glass et al. observed period doubling bifurcations and chaotic dynamics in the stimulation of heart fiber aggregates at frequencies slightly lower than the internal frequency at moderate stimulus amplitude. The same behavior was observed in the Poincar \ ‘ e oscillator [15,16,18]. The study of periodic exposure to the limit cycle has a direct application to the study of ventricular, or ventricular parastole. If the ectopic beat occurs outside the period of ventricular refractoriness, it is observed on the electrocardiogram, and the next normal (sinus) beat is blocked.

A similar model can be used to predict the sequences of sinus and ectopic beats in patients with ventricular parasistole (see, for example) [20-24]. In work [24] it is established that the model of pure parasistole (where sinus and ectopic pacemakers coexist without mutual influence on each other) corresponds well to reality. However, the model of the perturbed cycle does not answer the question of what happens at large relaxation times to the limit cycle (this model is built on the assumption of a rapid return of the system to it, which allows only one phase variable to be taken into account in the calculations). In addition, it is not known what happens when the values of the parameters at which the existence of the limit cycle is controversial. Thus, it seems necessary to consider a model based on more General principles that consider the change in the length of the perturbed cycle depending on the type of PRC without any additional conditions imposed on the behavior of the a priori system.

General case of interaction between two pacemakers

Consider two interacting pulse oscillators (or pacemakers) A and B with internal periods of Autonomous beats Ta and Tb, respectively. The interaction between the oscillators is again described using the so-called phase response curves (PRC). Again, we accept the same limitations on the nature of interaction.

Briefly indicate them

a. The phase of the perturbed oscillator changes instantly after exposure.

b. Phase shift depends only on two main parameters: the phase difference of the oscillators and the force of influence. In turn, the force of influence depends on the amplitude and coupling coefficient of the oscillators. In a real system, the coupling coefficient is an average factor that shows how the momentum decreases as it passes from one oscillator to another.

Therefore, the phase shift f, which determines the new phase of the perturbed oscillator with a period T, can be represented as follows: Φ = Δ/Τ ≡ Φ(φ, ε), where Δ is the time shift of the perturbed oscillator, φ is the phase difference of the oscillators, and ε is the force of influence. Pacemakers can be represented as a set of individual peaks in the timeline. Suppose that the moments of the last pulses of oscillators A and B – a and b, respectively (Figure 4). Note that a and b are momentum moments after all previous phase shifts of the oscillators. In other words, you can observe the pulses of the oscillators at these moments. Thus, it is necessary to analyze two cases.

i. b < a. This is the case (1) in Figure 4, i.e. oscillator B excited the medium before oscillator A. let us Follow the dynamics of the system in real time. The nearest event affecting the further behavior of the whole system is the appearance of the impulse A at the time a. let us Stop at this moment and make a forecast. To do this, we define the concept of the moments of the expected pulses of oscillators, ae and be (expected). Suppose we go back in time with respect to moment a. Since oscillator A has not yet acted, we should expect the appearance of the following pulses A and B at moments ae = a and be = b + Tb, respectively, where Tb – period B. Let’s call this situation as “A excites and B is in the expected state” and symbolically denote as (a, be). Now consider moment a. As A pulses, the following expected values aenext and benext can be represented as:

where Δb (φb, εb) - time shift of the oscillator B, due to the action of A. It depends on the phase φb of pacemaker A with respect to B and the force of influence εb. The phase φb can be calculated as follows: φb = (ae - be) ∕ Tb (mod1). Where, the phase φb is a positive value, and it belongs to the segment [0; 1] (negative values in the previous two expressions are excluded by the (mod1) operation). To determine which oscillator will excite the medium next, you need to compare ae next and be next. If, ae next < be next then A pulses and B remains in the expected state until be next, i.e. the system goes into state (a, be)next. Otherwise, if ae next > be next then B excites the medium and A goes into the expected state, and the state of the whole system becomes (ae, b) next.

ii. b>a. This is the case (2) in Figure 4. This inverse situation is similar to the previous case with the difference in reflections due to the pacemaker pulse A to B. Then the expected values of ae and be can be written respectively: ae next = a + Ta and be = b. You can call this case as “ B excites, and A is in the expected state” and denote as (ae, b). The following expected values are given by the expression:

where Δa (φa, εa) - time shift of the oscillator A, due to the action of B. It depends on the phase φb = (ae - be) ∕ Tb (mod1) of pacemaker A with respect toss A and the force of influence εa. Further analysis is also similar to case 1. Namely, if ae next > be next , then B pulsates and A goes into the expected state (ae, b)next , i.e. the state of the system becomes. Otherwise, if ae next < be next, then A excites the environment, and B remains in the expected state, and the system goes into the state (a, be) next. Summing up the above considerations, the model can be represented by the following diagram:

Generalization to N pacemakers

Suppose we have N Autonomous pulse oscillators, or pacemakers. Let’s also assume that all pacemakers are different. This means that each has its own internal cycle length Ti, i = 1, N, and the am plitude of the stimulus. To identify the connection between pacemakers, you must define the topology of the system space. In other words, you need to determine the nearest neighbors of each pacemaker in space. Conversely, it is obvious that determining the relationship between elements of such a spatially discrete system establishes its topology. Let’s assume that all pacemakers interact with each other, i.e. the so-called global coupling is implemented. Let the set of expected pacemaker pulses be located on the time axis see Figure 5. This means that in the absence of coupling, the pacemaker produces pulses at these points in time. In terms of expected values, which we will refer to for convenience as ae ≡ a, be ≡ b , the dynamics of the system can be described by the following difference equation: ssss

where:

To get the next expected values, you must perform the same procedure with the newly received expected pulses. Then the dynamics of the system can be easily represented as the following iterative relation:

where:

Moreover, it is convenient to use functions in the construction of equations {fij(φij n ,εij)} for dimensionless phase differences between pairs of pacemakers.

Analysis of the model of two interacting oscillators

The system consists of two oscillators A and B, connected via PRC Δa(φn a,εa), and Δb(φn b,εa). The map (4) that determines the dynamics of such a system is a special case of the General model (5) for N = 2. We rewrite (4) using the expression for the dimensionless PRC:

Suppose that pacemakers are of the same nature. Hence, we can assume fa(φ,ε) ≡ fb(φ,ε) ≡ f(φ,ε), where f(φ,ε) is a one-parameter function. The parameter ε integrally determines the overall effect of one oscillator on another. In the symmetric case εa ≡ εb ≡ ε. In our further consideration, we do not accept this symmetry.

We introduce the dimensionless phase difference of pacemakers A and B:

The choice Ta as a normalizing multiplier is irrelevant. Once selected, Tb we get similar expressions. Subtracting the second equation of the system (6) from the first and dividing tshe result by, Ta we come to the following expression:

Let’s make a brief analysis of the resulting display. Since in the General case x ∈ (‒ ∞, ∞), equation (7) is a one-dimensional nonlinear mapping of the real axis to itself. Note that this mapping cannot be reduced to a circle mapping by limiting x on the segment [0; 1], as is usually done for two pacemakers interacting via PRC. It is essentially asymmetric with respect to replacing x with - x see Figure 6. If f(x,ε) is a non-monotonic function, then the map (7) is nonlinear and can exhibit different behavior: from complex periodic to chaotic dynamics. Due to the fact that x ∈ (‒ ∞, ∞), in the strict sense it is not the phase difference of pacemakers. In this context, x can be called a generalized phase difference. By analyzing equation (7), you can determine which oscillator (A or B) is pulsating at a given discrete time n. This depends on the sign x: A excites the medium if xn < 0, and B pulsates when xn > 0. This makes it possible to determine the degree of phase capture of pacemakers. However, taking only the values into consideration xn, we will not be able to restore the original time series of events of pacemakers A and B. Using them xn, we can only talk about their phase difference. We will demonstrate how equations (6), (7) can be applied to study the behavior of two interacting pacemakers. Analysis of the real system shows that the function f (x, ε) can take various forms, but, as a rule, it satisfies some General properties that have already been mentioned earlier. For example, f (0, ε) = f (1, ε) =0. It usually has a maximum and a minimum. Sometimes, instead of extremes, it has breaks. Consider f (x, ε) in the elementary frequently used form. Then let f (x, ε) = εsin2πx, and this will lead to a system of dynamic equations:

For Figure 6. examples of direct modeling of the system (4), (5) are presented. The left part of the figure shows some possible phase captures obtained from equation (4). The right and lower parts show the corresponding mappings (5), the periodic orbit defining the values of Lyapunov exponents.

Approximation of the active medium as a lattice of pulse oscillators

In this section, we will demonstrate a way to approximate discrete distributed environments based on the General model of coupled oscillators (5). Looking at the heart pacemaker at a microscopic level, it can be thought of as a large group of cells that generate heart rate and synchronize their action potentials to initiate heart contractions. Thus, instead of considering a single pacemaker, we can construct a lattice of coupled pulse oscillators. In this paper, we have limited ourselves to one-dimensional (chain) and two-dimensional (lattice) cases. Assume that the Autonomous pacemakers are located at the nodes of a two-dimensional square lattice of size (N х M). we denote the lattice element with coordinates (i, j) as Aij, where i=1, N and j=1, M. we restrict ourselves to considering a homogeneous medium and accept some restrictions on anisotropy. This means that the lattice pacemakers are identical, i.e. they have the same cycle length Tij ≡ T , i=1,..., N; j=1,..., M, (however, in reality, the cells on the periphery of the sine pacemaker have the shortest cycle length, although its center acts as the leading pacemaker). This restriction reduces the number of system parameters and therefore makes it easier to study the model. Now we will define the relationship between the elements. In works on lattices of concatenated maps, two main types of coupling are usually considered: nearest neighbor coupling and global coupling. Since in the previous sections we assumed that pacemakers all interact with each other, this time, as an example, we will consider lattices with a connection of the nearest neighbor type: first a two – dimensional lattice, and then a chain of coupled pulse oscillators. In a square grid, each pacemaker Aij interacts with four neighbors in accordance with the schematic (Figure 7a). Considering the restriction of a homogeneous medium, we can assume that all the lattice couplings are identical, i.e. every two adjacent elements interact with each other according to a General law defined by identical PRC f (x,ε). Moreover, suppose that the relationship between a pair of elements is isotropic in the sense that, and is equal to one of the values or, depending on the corresponding arrangement of elements. This means that there is an anisotropy of the force of influence in the vertical and horizontal directions. In other words, if two pacemakers are neighbors in the vertical direction, they interact through f (x, ε1), and if they are horizontal neighbors, they are linked through f (x, ε2).

Note that all restrictions are made to simplify the analytical form of the resulting model. In fact, it is possible to write expressions for a two-dimensional lattice of coupled pulse oscillators without any restrictions. We present equations that determine the iterative dynamics of the expected pacemaker pulses {aij} i=1,N; j=1,M based on the approach presented in paragraph 4. In order to get the (n+1)th value of a single element aij, it is necessary to analyze all the elements of the lattice, since they are connected to each other by local coupling. In other words, the element in question cannot be affected by the others at the (n)th time step, since it is suppressed by the influence of other elements and remains in the expected state until the (n+1)th step. Thus, the dynamics of the model can be described by the following expression:

The constructed model deserves a detailed study based on the approach developed for lattices of concatenated maps [1]. As a second example, consider a chain of identical pulse oscillators connected by the nearest neighbor principle. We will limit ourselves to a homogeneous case with anisotropy of the right and left directions in the force of influence between the nearest neighbors. A schematic picture of the chain is shown in Figure 7b. Similar to the above consideration, you can get:

If, then equations (11) define the so – called open – flow model (159). Since in this section of the thesis we present a General approach to the development of models without a detailed analysis of their behavior, the type of boundary conditions for both lattices was not specified. Therefore, for an analytical or numerical study of such systems, it is necessary to determine both the boundary conditions and the PRC f (x, ε). Usually the boundary conditions are selected as periodic, i.e. ai j+M ≡ aij, ai+Nj for a two-dimensional lattice and for a one-dimensional one ai+N ≡ ai. For an open stream model, the condition of a fixed left border is often accepted. The described models (4), (10), and (11) can be generalized to the natural non-uniform case, assuming the internal lengths of pacemak ers, PRC, and influence forces to be different for different groups of elements. However, considering inhomogeneous anisotropic lattices is an extremely difficult task even for numerical analysis. The first attempts to study inhomogeneous lattices of concatenated maps are described in [25].

References

1. Kaneko K (1989) Spatiotemporal chaos in one-dimensional and two-dimensional coupled map lattices. Physica D: Nonlinear Phenomena 37(1-3): 60-82.
2. Beuter A, Glass L, Mackey MC (2003) Nonlinear Dynamics in Physiology and Medicine. Titcombe MS Eds Springer-Verlag Inc, New York, USA.
3. Guevara MR, Shrier A, Glass L (1986) Phase resetting of spontaneously beating embryonic ventricular heart cell aggregates. Am J Physiol 251(6 Pt 2): H1298-H1305.
4. Arnold VI (1991) Cardiac arrhythmias and circle mappings. Chaos 1(1): 20-24.
5. Kaneko K (1984) Supercritical behavior of disordered orbits of a circle map. Progress of Theoretical Physics 72(6): 1089-1103.
6. Berklinblit MB, Gelfand IM, Garfinkel VS, Fomin SV, et al. (1971) Models of the Structural–Functional Organization of Certain Biological Systems. Cambridge: MIT Press.
7. Glass L, Belair J (1986) Continuation of Arnold tongues in mathematical models of periodically forced biological oscillations. Nonlinear Oscillations in Biology and Chemistry, Lecture Notes in Biomathematics 66 / Ed Othmer HG, Berlin Springer–Verlag pp. 232-243.
8. Glass L, Mackey MC (1979) A simple model for phase locking of biological oscillations. Journal of Mathematical Biology 7: 339-352.
9. Keener JP, Hoppensteadt FC, Rinzel J (1981) Integrate and fire models of nerve membranes response of oscillatory inputs. SIAM J Appl Math 41(3): 503-517.
10. Glass L (1991) Cardiac arrhythmias and circle maps - A classical problem. Chaos 1(1): 13-19.
11. Keener JP (1981) On cardiac arrhythmias: AV conduction block. J Math Biol 12: 215-225.
12. Glass L, Guevara MR, Shrier A (1987) Universal bifurcations and the classification of cardiac arrhythmias. Ann N Y Acad Sci 504: 168-178.
13. Shrier A, Dubarsky H, Rosengarten M, Guevara MR, Nattel S, et al. (1987) Prediction of atrioventricular conduction rhythms in humans using the atrioventricular nodal recovery curve. Circulation 76(6): 1196-1205.
14. Hoppensteadt FC, Keener J (1982) Phase locking of biological clocks. Journal of Mathematical Biology 15: 339-349.
15. Guevara MR, Glass L (1982) Phase locking, period doubling bifurcations and chaos in a mathematical model of periodically driven oscillator: A theory for the entrainment of biological oscillators and the generation of cardiac dysrhythmias. J Math Biol 14(1): 1-23.
16. Keener J, Glass L (1984) Global bifurcations of a periodically forced oscillator. J Math Biol 21(2): 175-190.
17. Zaslavsky GM (1978) The simplest case of a strange attractor. Physics Letters A 69(3):145-147.
18. Ding EJ (1987) Analytic treatment of a driven oscillator with a limit cycle. Phys Rev A Gen Phys 35(6): 2669-2683.
19. Guevara MR, Glass L, Shrier A (1981) Phase locking, period doubling bifurcations, and irregular dynamics in periodically stimulated cardiac cells. Science 214(4527): 1350-1353.
20. Moe GK, Jalife J, Mueller WJ, Moe B (1977) A mathematical model for parasistole and its application to clinical arrhythmias. Circulation 56: 968-979.
21. Jalife J, Antzelevitch C, Moe GK (1982) The case for modulated parasistole. Pacing Clin Electrophysiol 5(6): 911-926.
22. Ikeda N (1982) Model of bidirectional interaction between myocardial pacemakers based on the phase response curve. Biol Cybern 43(3): 157-167.
23. Ikeda N, Yoshizawa S, Sato T (1983) Difference equation model of ventricular parasystole as an interaction between cardiac pacemakers based on the phase response curve. J Theor Biol 103(3): 439-465.
24. Courtemanche M, Glass L, Belair J, Scagliotti D, Gordon D (1989) A circle map in a human heart. Physica D: Nonlinear Phenomena 40(3): 299-310.
25. Loskutov A Yu, Prokhorov AK, Rybalko SD (2002) Analysis of inhomogeneous chains of connected quadratic maps. Theoretical and mathematical physics 132: 105-125.