The first theoretical study of cardiac arrhythmias was conducted in 1920-ies of the van - der - Pol and van - der - Mark, who proposed an electrical model of the heart [1,2]. They suggested that the activity of the heart can be simulated by three nonlinear oscillators corresponding to ACS, Atria and ventricles. This possibility was confirmed by the fact that the graphical representation of the processes occurring in the system of nonlinear oscillators had a form similar to the action potentials of heart cells. In this system, there was a unidirectional connection between the sinus and atrial oscillators, and the same connection existed between the atrial and ventricular oscillators. Reducing the connection between the latter, they found that it is possible to obtain a number of different rhythms with phase captures, qualitatively corresponding to the AB- blockades of the heart.

After the publication of the work [2], many researchers tried to model the dynamics of heartbeats, believing it to be generated by several coupled oscillators. These developments can be divided into two large groups: analysis based on continuous - time representations (ODE systems) and studies based on discrete - time representations (finite-difference equations). In this section we will consider in detail the main stages of development of these concepts and the corresponding mathematical concepts.

Discrete Models & Of A System of Interacting Pacemaker

Consider some physical quantity ξ, which reflects the internal state of the biological oscillator. Let own oscillator period is equal to 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. Define the phase of the oscillator 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>φ = 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:

                φn+1 = φn + f(φn) (mod1).                                                                                                                           (1)

The function f(φn) is called the phase response curve (PRC) and determines the phase change after the stimulus. The points  f(φn) of the system state are conveniently represented 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 and equal to the number of passes φn+1 on the unit circle for the time in which  f(φn) passes it 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 fig.1.

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 phase of multiplicity N/M. If the rotation number is irrational, the system exhibits quasi-periodic or chaotic behavior.

tn+1 = tn + T1, τn+1 =  τn  + T2 .                                                                          (2)

Figure.8. Phase capture areas of a system of oscillators with two-way coupling (δ = 0,1): (a) ε = 0,1: (b) ε = 0,5.

The seizure phases in the space (ε, γ)

          

                                                                                                                (9)

 where function f is defined by [removed]10):

  f = x + a + εh[a-1 (1 + γh(x) − x)] − γh(x),                                                                                                                 (10)

   where a and q are parameters. The set of parameters here are a, γ, ε and δ.If the map (11) has a cycle p of period t equal to the perturbation period, t = τ,  p = x1, x2,..., xt, the points forming this cycle will obey the following system of equations (11):

                                                                                                                                                                     

                                                             (11)    

To solve the inverse problem, i.e. to find the parameter values at which the map (9) has a given cycle p, it is necessary to express the values from (11) qi. It is clear that not for all possible xi the obtained values of the parameters will satisfy the ratio qi ∈ Q. However, i if this is true for any cycle p = x1, x2,..., xt  you can find the values of the parameters q1,q2,..., qt, for which the perturbed map (9) has such a cycle.

If the multiplier of the cycle  β (p) = Πti=1 f '(x) < 1>xc, the multiplier is always less than one, which ensures stability. In case of unilateral influence of pacemakers on each other, at ε = 0, the display at γ ≥ (√5 ∕ 5) for polynomial PRC taking into account refractoriness (pic.10a). The functions f(x), determined from (10), under the mutual influence of oscillators are shown in (pic.10b) for polynomial PRC.

 As another approximation of the experimental CFO we take the following polynomial function (7):

         h(x) = Cx² (1/2 − x)(1 − x)² .                                                                                                                                (7)

We choose the normalization factor C in such a way that the amplitude of h (x) equals 1, i.e. C = 20√5 (see Fig.2). Then, taking into account the refractoriness, leaving the function continuous, and ε = 0, the map (2) will take the form (8):

                                                                                (8)

where h(x) is defined from (7). In contrast to the sinusoidal approximation, this curve at x = δ touches the abscissa axis. In other words, this means that the polynomial CFO, taking into account the refractoriness on the entire segment [0; 1], is a smooth function. Let us now compare several cases with different values of the refractor time and the amplitude of the influence ε of the first oscillator on the second.

Change of phase diagrams with increasing refractoriness (case of one-way interaction of oscillators)

Consider the system (6,8) without taking into account the mutual influence of oscillators and the period of refractoriness (in other words, ε = δ = 0), i.e. assuming in (10) δ = 0. For fig.5a the regions of phase captures in parametric space (a, γ) obtained as a result of numerical investigation are presented. In this figure, the color gamut is much richer (visualization was carried out using Matlab R2007b, and the equations were solved using Microsoft Visual Studio C++ 2010 Professional (x64)). Stable phase captures up to 10:10 multiplicity captures are shown here. It can be seen that due to the smoothness of the polynomial function, the boundaries of the resonance languages blur significantly less than in the previous phase diagrams.

The Analogy With Abnormal Heart Rhythms

In conclusion, we draw an analogy between the results and pathological conditions of cardiac tissue. With the help of the constructed models it is possible, for example, to describe the interaction of sinus and ectopic pacemakers, PRC and ABU, and the impact of external disturbances on the sinus rhythm. If, for example, consider the first pulse oscillator PRC, and the second – ABU, it can be found that some stable phase captures correspond to the observed pathologies in clinical practice. In this case, among the various grippers are constructed as normal sinus rhythm (capture ratio 1:1) and classical rhythms Wenkebach (captures multiplicity N:(N-1)) and N:1 AB-blockade. If the first pulse count system of the ABU, and the second ACS, then appear inverted rhythms Wenkebach (similar to direct, but that changes the role of the ventricles and the Atria) observed in some patients.

It should be noted that the considered response functions  f₁(x) = γh(x),  f₂(x) = εh(x) for different approximations of the form, the functions h(x) are model. They were taken to analyze the characteristic features of the dynamics of two nonlinear interacting oscillation sources.  In practice, these functions should be chosen taking into account additional physical assumptions about the nature of the interaction and take into account the experimental data on the response of a single oscillatory system to single pulses of external perturbation. For example, in [11] the effect of short pulses on aggregates of spontaneously oscillating cells from the embryo heart was considered. The experimentally obtained phase response curves were approximated by exponential functions, and the "internal" parameters were chosen for the best correspondence of the curve graphs to the experimental points. Their dependence on physical parameters was also chosen in this way. As a result, the phase diagram obtained numerically corresponded well enough to the real dynamics of the system.

Presence of wide areas of phase captures (pic.5 − 9), in such systems, various types of synchronization of two oscillators are possible, which qualitatively correspond to some types of cardiac arrhythmias. The phase diagram allows to reveal under what conditions of interaction (i.e. at what values of parameters a, γ, ε and δ) this or that kind of synchronization is possible. Moreover, all the phase patterns presented in this paper indicate that with increasing nonlinearity (i.e., with the growth of the γ parameter), the regions with different captures begin to overlap. Knowledge of such areas and the dynamics of the system in these areas allows by external perturbation (for example, a series of single pulses) to withdraw the system from the unwanted mode of synchronization to a more favorable mode, which is vital.

Analysis of phase diagrams makes it possible to find ways to control such systems.We consider the effect of additional periodic pulse action on the behavior of interacting oscillatory subsystems. The study of possible modes of behavior of such a system by varying the frequency and amplitude of the external perturbation will lead its dynamics to a predetermined, for example, to complete suppression of the ectopic sine pacemaker. This problem, considered in the next section, is very relevant for the General theory of control of nonlinear dynamic systems and excitable media, in particular, cardiac tissue, which is satisfactorily described in the models [12,13].

Modern methods of removing the heart from the state of fibrillation are very rigid (supply of a short electrical pulse of a huge voltage and a large current). The development of nonlinear dynamics and synergetics made it possible to understand that such a force effect is not necessary. Often enough weak electrical effects directly on the heart muscle. Precisely, if there are spiral waves with opposite directions of rotation in the medium, then, choosing the phase and frequency of external action, it is possible to achieve the movement of the centers of the two waves towards each other and their annihilation. Now the word for careful experimental research. The theory of dynamic systems describes many processes inherent in active media, including some types of arrhythmias [1,3]. Since arrhythmias are caused by certain disorders in the heart muscle and, therefore, are pathological conditions, the modeling of such systems is of great practical interest and can bring closer to solving the problem of the possibility of controlling their behavior through external influences. This, in turn, allows us to come close to the problem of soft withdrawal of active systems from the state of developed space–time chaos that characterizes some types of pathologies [14-17]. A noteworthy work is presented in the publication [18]. In this paper, we have implemented a model of the heart, which describes the real behavior of rhythmic and arrhythmic processes in the heart.