Numerical Analysis of the Coupled Modified Van Der Pol Equations in a Model of the Heart Action

In this paper, a modified van der Pol equations are considered as a description of the heart action. Wide ranges of the model parameters yield interesting qualitative results, e.g. Hopf bifurcation, Bogdanov-Takens bifurcation, transcritical and pitchfork bifurcations but also some stable solutions can be found. The physiological model works in the narrowest range of parameters which allows to obtain a stable behaviour what is important in biological problem. When some kinds of pathologies appear in the heart, it is possible to obtain chaotic behaviour. My aim is to compare the influence of these two types of coupling (unidirectional and bidirectional) on the behaviour of the van der Pol system. The coupling takes place in a system with healthy conductivity, between two nodes: SA and AV, but in some circumstances, a pathological coupling may occur in the heart. The van der Pol oscillator is a type of relaxation oscillator which can be synchronized. In this paper, synchronization properties of such a system are studied as well. For the purpose of a numerical analysis of the system in question, a numerical model was created. Keywords-van der Pol equation; heart action; coupling; synchronization;


I. Introduction
This paper is related to the research on the electrical conduction system of the human heart.In the heart, in addition to ordinarily working fibres, there are pacemaker centres made of special cells that resemble embryonic cells.These are the cells of the electrical conduction system forming the following concentrations: the sino-atrial node (SA) and the atrioventricular node (AV) and His-Purkinje system, [1].The key elements of the conduction system which we consider are the SA node and the AV node.Each of the two nodes is modelled by the modified van der Pol oscillator.This model allows for rendering phenomena observable in clinical experiments, such as Holter recordings.The aim of this work is to create a model which is able to render the behaviour typical of the sinoatrial block.The initial pulse in the heart is usually formed in the SA node, and carried through the atria to the AV node.In the SA block, the electrical impulse is delayed or blocked on the way to the atria, thus delaying atria depolarization.This is different from AV block which occurs in the AV node and delays ventricular depolarization.The SA blocks are categorized into three classes, based on the length of the delay.The first degree SA block is characterized by a prolonged conduction time from the SA node to the surrounding atrial tissue.The second degree block has two types: Wenckebach block and type II.The Wenckebach block shows an irregular rhythm.The pause of the second degree type I is shorter than twice the minimum length of the period.Type II has a regular rhythm, which may be normal or slow.It is followed by a pause, which is a multiple of the period.Conduction across the SA node is normal until the pause, and then it is blocked.The third degree is characterized by lack of atrial activity.Heart rhythm is determined by escape rhytm.In biological systems, the phenomenon of synchronization is extremly important.It is responsible for many periodic processes in the body, e.g. the coupling of the pituitary gland is responsible for production of hormones from the thyroid gland that produces them, without synchronization of the two components, the operation of our endocrine system would not be correct .Also the main oscillator in the human body, i.e. the heart is subject of synchronization.In this paper, I will discuss the impact of different types of couplings connecting the AV node and the SA node, and how they affect the synchronization of the two oscillators.The analysis of synchronization of various modifications of the van der Pol model is the aim of many papers.Synchronization areas near the main parametric resonance and transition conditions from regular to chaotic motion are presented in paper [1].The phenomenon of complete synchronization in a network of four coupled oscillators is described in [2].In paper [3], the authors investigated mechanisms of various bifurcation phenomena observed in the Bonhoffer van der Pol neurons coupled through the characteristics of synaptic transmissions with a time delay.Also synchronization phenomena in van der Pol oscillators coupled by a time-varying resistor is researched in paper [4].However, these articles do not offer any examples of application of this model for recreating pathological behaviour of the electrical-conduction system of the human heart, and therefore the considered ranges of parameters are wider than those applicable for medical applications.In papers [5,6], the authors showed that coupled two van der Pol oscillators modelled behaviour of the heart conduction system, and there are described a heart block as pathologies of coupled van der Pol oscillators.The van der Pol oscillator provides rich dynamical behaviour, which we would like to exploit in the modelling of the heart action [7] also synchronization phenomena.

A. Mathematical model
Because each node is a self-exciting pacemaker, it can be described by a relaxation oscillator, i.e. the van der Pol oscillator.The model by van der Pol and van der Mark was created as a model in the electronic circuit theory in 1927: where f (x) = 1 2 a(x 2 − 1) is a damping coefficient being a function of the x variable, which is negative for |x| < 1 and positive for |x| > 1.The dynamics of Eq. ( 1) is well-known in literature.
The van der Pol model needs some changes in order to reproduce the actual features of the action potential.Postnov [8] introduced modifications that maintain the required structure of the phase space.To be more precise, he substituted the linear term by a nonlinear cubic force called the Duffing term where a, µ, d are positive control parameters.This model can be treated as a SA or AV node model.The mutual interaction of the limit cycle present around an unstable focus with a saddle and a stable node is the main property of a modified relaxation oscillator.As a result, the refraction period and the nonlinear phase sensitivity of the action potential of node cells are reproduced correctly.A solution of this equation in terms of time presents the action potential, whereas a solution in terms of velocity enables us to obtain crucial phase portrait.As we can see, the main qualitative difference between Eqs. ( 1) and ( 2) is the appearance of two additional steady states, i.e. x 2 = −d and x 3 = −2d.As in the previous case, x 1 = 0 is an unstable node or a focus surrounded by a stable unique limit cycle, x 2 = −d is a saddle and x 3 = −2d forms a focus or a node and can be either stable or unstable, depending on the sign of 4d 2 − µ.In the case considered by Postnov [8], the first steady state is an unstable focus, while the third is a stable node which attracts all solutions starting on the right hand-side of the stable manifold of the saddle x 2 .However, in the considered model (2) it is difficult to regulate the location of steady states in the phase space, therefore, a new parameter e is introduced in order to reproduce the heart behaviour.Notice that this modification has no influence on the phase portrait, whereas we have the opportunity to modify the location of steady states.In order to simplify frequency regulation and obtain the proper timescale, the ed factor in the denominator is substituted with independent coefficient f , corresponding to harmonic oscillator's frequency, [9].Below we present the model in its two variable first order form which reads [9] The final system consists of two coupled modified van der Pol oscillators.This model can be treated as the SA and AV node.The final system that we analyze is given in the following form: 3, e 1 = 7, e 2 = 4.5 control parameters.Parameters values for the modified van der Pol model were chosen so that the oscillations frequency correspond to real frequencies of the SA and AV nodes.The selection of appropriate parameters was done after the verification of the model by Grudziński in [9].The aim was to recreate the physiological properties of the biological model using mathematical equations.This information is essential for examining stability of our setup because without such limitations the system could have completely different properties and would not recreate physiological properties.Modification of the e parameter of the node location influences the distance between consecutive potential needles without changing their shape.This means that the mutual position of the saddle and the node influences the time of spontaneous depolarization, which is one of physiological mechanisms of the regulation of the action potential generation frequency.An increase of the value of the parameter e can be interpreted as an increase of the activity of the nervous system.However, the parameter f is the equivalent of the frequency of the harmonic oscillator.The parameter d adjusts position of a fixed point.The system with delayed feedback (sum of delays with feedback) describes various pathologies observed in the heart action, for example, the SA block, which does not conduct the potential in physiological way.There are situations when the output potential from the SA node influences the input and modifies the action of the system, for example, through injury caused by infarction or instrinsic disease in the SA node.

B. Types of coupling
The origins of the phenomenon of self-exciting reconciliation of vibrations of coupled oscillators back to the seventeenth century.Then the Christian Huygens observed that in the clock with two pendulums after sufficient time has always the situation in which both pendulums oscillated with opposite phases occured, regardless of the initial phase difference.The behaviour of cardiac pacemaker cells resembles that relaxation oscillators.A characteristic property of relaxation oscillators is that they may be synchronized by an external signal, if the latter has a periodicity not differing too much from the spontaneous frequency of the oscillator [7].Synchronization which is defined as an adjustment of rhythms due to weak interaction, is one of the most interesting features displayed by coupled oscillators.We investigate a phenomenological model for the heartbeat consisting of two coupled van der Pol oscillators.The coupling between these nodes can be both unidirectional and bidirectional.In addition, feedback may also occur.We know that the sinoatrial node is also reffered to as the pacemaker of the heart.When the impulses generated by the SA node reach the AV node, they are delayed.So we have here unidirectional coupling (s 1 or s 2 is different from zero in Eq.4).However, bidirectional coupling (s 1 and s 2 are different from zero in Eq.4) is also possible in the heart, for example, during the WPW syndrome.Feedback (a part of Eq.4 with k coupling coefficients and with w delays) also appears only in case pathologies, for example, SA and AV blocks.The Pecorra Caroll (PC) theory is considered in this paper.This type of coupling is used, when a state variable from a chaotic system is input into a replica subsystem of the original one,and as a result, both systems can be synchronized identically.
where ẋ1 = (u 1 , v1 ), The drive system is presented as: and response is given as follows: The resulting equation for the PC theory gives the following form:

II. Numerical analysis
For the purpose of numerical analysis of the discussed system, a numerical model was created using Dynamics Solver and a program in C++ was developed.A Dormand Prince 8 integration algorithm was used.This method constitutes a modification of the explicit Runge-Kutta formula with a variable integration step.In this section, we use some numerical simulations in order to illustrate the pathological behaviour described in the Introduction.System with delayed feedback describes various pathologies observed in the heart action, e.g.SA block.When added,the s 2 coupling makes the SA node influence the AV node rhythm.This type of behaviour is of physiological nature.By adding s 1 , we arrive at pathological behaviour, and consequently, a reentry wave in our system.It is typical of the WPW syndrome.Having included feedback and delay, we obtain a time series which resembles the model presented in Figure 1.The red graph presents the initial model without feedback.The blue one presents a modified model with feedback.The parameter values are as follows: k = 1, k 1 = 2.85 and w 1 = 0.75, w 2 = 0.25, and the remaining parameters are the same as in the reference system.In the time series of the modified model with feedback there is a 'delayed impulse'.The period of this oscillator is almost twice as long as in the reference model (the potential period for a single node model of an electrical conduction system with no coupling and feedback is app.1.4 ), similar like in the second type of the SA block.This is one of the mechanisms causing brachycardia.
If we consider the physiological coupling between nodes, then the s 2 coupling is introduced to our system.It means that the SA node directs AV node.With small values of s 2 , the result is similar to the reference one, but with bigger values, (e.g. 10 or 100) some of the amplitudes are synchronized in-phase with x 1 and aperiodic behaviour appears, which is presented in Figure 2.
The addition of coupling to the y 1 term allows us to model the reentry wave, which causes the exceptional situation when AV node is the master for SA node.Such situation takes place in case of the WPW syndrome.Slowing the oscillations down caused by feedback and addition of the s 1 coupling, we obtain the aperiodic behaviour.The big arrhythmia occurs.From the medical point of view, it resembles atrial fibrillation.The oscillator begins to work aperiodically, trying to adjust its frequency to the frequency of the AV node.It has a tendency to shorten the oscillation period despite the lack of periodic behaviour, as presented in Figure 3.
The phase portrait is similar to the previous case, but with s 1 = 3.5, we observe oscillation death, Figure 4.
The amplitude death can be understood as a full heart block.No impulse is conducted to the AV node.As a result, the AV node may take over the function of the pacemaker.
Figures 5 and 6, which also present various plots of synchronization, indicate that synchronization of the system is greater in case of systems with the s 2 coupling than those with the s 1 coupling.We applied the PC theory in our system, in spite of the fact that this theory is typical of chaotic behaviour.By applying the theory, we can observe greater synchronization in phase, especially for s 2 coupling.The model with the s 1 coupling also tries to synchronize in the phase.
In the system with two couplings, already in the case of s 1 = 1, the oscillator corresponding to the SA node tries to adjust its rhythm to the AV node oscillator by shortening its period.With T = 2.3, we obtain biperiodic behaviour, where T = 1.6 and T = 2, Figure 7.
With s 1 = 20 and s 2 = 1, we get aperiodic behaviour-typical arrhythmia.With these parameters, there is no synchronization, whereas with values s 1 = 2 and s 2 = 5, we observe interesting behaviour.Periods of both oscillators are shortened.Oscillator corresponding to the AV node behaves periodically, with its period at the level of 1.6, whereas the one corresponding to the SA node is biperiodic, with periods 1.5 and 1.7.With such system parameters, we observe the antiphase type of synchronization, Figure 8.

III. Conclusion
Although the uncoupled van der Pol equation has quite trivial dynamics as a stable equilibrium, the system with the coupling can be periodic, but also quasi-periodic, and chaotic.Similarly, these couplings of the mathematical system interfere with the work of the heart conduction system (SA block, AV block, bradycardia, WPW syndrome).Synchronization in the discussed cases is rarely in-phase, but often in anti-phase.The PC theory was also applied to a non chaotic system.This unidirectional coupling affects partial synchronization of our system.Bidirectional coupling should be used to describe physiological behaviour of the conduction system, because only this way we can take into account the effect of the AV node as a delay element (coupling s 1 ).In our system without couplings, if we have asystole than we can try to give additional s 1 coupling.With small values of s 1 , we observe asystole, while with greater values, the SA rhythm appears but it is not synchronized.The model offered in this study, is a correct reconstruction of heart action pathologies, such as a SA block or a type of the arrhythmia.

Fig. 1 .
Fig. 1.Time series: red line-without feedback, blue one-with feedback