Original article Biomath 3 (2014), 1312281, 1–7

B f

Volume ░, Number ░, 20░░ 

BIOMATH

 ISSN 1314-684X

Editor–in–Chief: Roumen Anguelov  

B f

BIOMATH
h t t p : / / w w w . b i o m a t h f o r u m . o r g / b i o m a t h / i n d e x . p h p / b i o m a t h / Biomath Forum

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

Heart Action
Beata Zduniak

Faculty of Applied Informatics and Mathematics
Warsaw University of Life Sciences (SGGW)

Warsaw, Poland
Email: beata.zduniak@wp.pl

Received: 11 October 2013, accepted: 28 December 2013, published: 21 May 2014

Abstract—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 bifur-
cation, Bogdanov-Takens bifurcation, transcritical
and pitchfork bifurcations but also some stable so-
lutions 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 con-
duction system which we consider are the SA node
and the AV node. Each of the two nodes is mod-
elled 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 sinoa-
trial 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

Citation: Beata Zduniak, Numerical Analysis of the Coupled Modified Van Der Pol Equations in a
Model of the Heart Action, Biomath 3 (2014), 1312281,
http://dx.doi.org/10.11145/j.biomath.2013.12.281

Page 1 of 7

http://www.biomathforum.org/biomath/index.php/biomath
http://dx.doi.org/10.11145/j.biomath.2013.12.281


B Zduniak, Numerical Analysis of the Coupled Modified Van Der Pol Equations...

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 synchro-
nization 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 transi-
tion 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 char-
acteristics 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 parame-
ters 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 phe-
nomena.

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:

ẍ + 2 f (x)ẋ + x = 0, µ > 0, (1)

where f (x) = 12 a(x
2 − 1) is a damping coefficient

being a function of the x variable, which is neg-
ative 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

ẍ + a(x2 −µ)ẋ +
x(x + d)(x + 2d)

d2
= 0, (2)

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 cor-
rectly. A solution of this equation in terms of

Biomath 3 (2014), 1312281, http://dx.doi.org/10.11145/j.biomath.2013.12.281 Page 2 of 7

http://dx.doi.org/10.11145/j.biomath.2013.12.281


B Zduniak, Numerical Analysis of the Coupled Modified Van Der Pol Equations...

time presents the action potential, whereas a so-
lution 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.
x2 = −d and x3 = −2d. As in the previous case,
x1 = 0 is an unstable node or a focus surrounded
by a stable unique limit cycle, x2 = −d is a saddle
and x3 = −2d forms a focus or a node and can be
either stable or unstable, depending on the sign
of 4d2 − µ. 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 x2. 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 fre-
quency 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]

ẋ = y,
ẏ = −a(x2 − 1)y − f x(x + d)(x + e).

(3)

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:

ẋ1 = y1 + (k − k1)x1[t − w1] − kx1 + k1 x1[t − w2],
ẏ1 = −a1(x12 − 1)y1 − f1 x1(x1 + d1)(x1 + e1)+
+s1(x2 − x1),
ẋ2 = y2,
ẏ2 = −a2(x22 − 1)y2 − f2 x2(x2 + d2)(x2 + e2)+
+s2(x1 − x2),

(4)
where k, k1 coupled coefficients, s1, s2 coupled
coefficients, w1, w2 delays, a1 = a2 = 5, f1 =
f2 = 3, d1 = d2 = 3, e1 = 7, e2 = 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 param-
eters 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 fre-
quency. 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 os-
cillator. The parameter d adjusts position of a fixed
point. The system with delayed feedback (sum of
delays with feedback) describes various patholo-
gies 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 pace-
maker 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

Biomath 3 (2014), 1312281, http://dx.doi.org/10.11145/j.biomath.2013.12.281 Page 3 of 7

http://dx.doi.org/10.11145/j.biomath.2013.12.281


B Zduniak, Numerical Analysis of the Coupled Modified Van Der Pol Equations...

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 phenomeno-
logical 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 (s1 or s2 is different from
zero in Eq.4). However, bidirectional coupling (s1
and s2 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.

ẋ1 = f (x1),
ẋ2 = f (x2),

(5)

where ẋ1 = (u̇1, v̇1), x1 ∈ Rn, u1 ∈ Rp, v1 ∈ Rq. The
drive system is presented as:

u̇1 = g(u1, v1),
v̇1 = h(u1, v1),

(6)

and response is given as follows:v̇2 = h(u1, v2).
The resulting equation for the PC theory gives the
following form:

ẋ1 = y1 + (k − k1)x1[t − w1] − kx1+
+k1 x1[t − w2],
ẏ1 = −a1(x12 − 1)y1 − f1 x1(x1 + d1)(x1 + e1)+
+s1(x2 − x1),
ẋ2 = y2,
ẏ2 = −a2(x12 − 1)y2 − f2 x1(x1 + d2)(x1 + e2)+
+s2(x1 − x2),

(7)

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

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 s2 coupling
makes the SA node influence the AV node rhythm.
This type of behaviour is of physiological nature.
By adding s1, 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, k1 = 2.85 and w1 = 0.75, w2 = 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 be-
tween nodes, then the s2 coupling is introduced
to our system. It means that the SA node directs

Biomath 3 (2014), 1312281, http://dx.doi.org/10.11145/j.biomath.2013.12.281 Page 4 of 7

http://dx.doi.org/10.11145/j.biomath.2013.12.281


B Zduniak, Numerical Analysis of the Coupled Modified Van Der Pol Equations...

Fig. 2. Time series: physiological coupling between nodes

Fig. 3. Time series: pathological coupling between nodes

AV node. With small values of s2, 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 x1 and aperiodic
behaviour appears, which is presented in Figure 2.

The addition of coupling to the y1 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 s1
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 s1 = 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 s2 coupling than those with the s1 coupling.

Fig. 4. Phase portrait: pathological coupling between nodes

Fig. 5. Synchronization plot for unidirectional case: s1 = 5

- 2 - 1 0 1

- 2

- 1

0

1

x

w

Fig. 6. Synchronization plot for unidirectional case: s2 = 10

Biomath 3 (2014), 1312281, http://dx.doi.org/10.11145/j.biomath.2013.12.281 Page 5 of 7

http://dx.doi.org/10.11145/j.biomath.2013.12.281


B Zduniak, Numerical Analysis of the Coupled Modified Van Der Pol Equations...

Fig. 7. Time series:bidirectional coupling between nodes

Fig. 8. Synchronization:bidirectional coupling between
nodes

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 s2
coupling. The model with the s1 coupling also tries
to synchronize in the phase.

In the system with two couplings, already in the
case of s1 = 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 s1 = 20 and s2 = 1, we get aperi-
odic behaviour- typical arrhythmia. With these
parameters, there is no synchronization, whereas
with values s1 = 2 and s2 = 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 anti-
phase 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 synchroniza-
tion 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 s1). In our system
without couplings, if we have asystole than we
can try to give additional s1 coupling. With small
values of s1, 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.

References

[1] J. Warmiński, Synchronisation effects and chaos in the
van der Pol-Mathieu oscillator, Journal of Theoretical and
Aapplied Mechanics, 4, 39, 2001.

[2] P. Perlikowski, A. Stefanski, T. Kapitaniak, Discontin-
uous synchrony in an array of Van der Pol oscillators,
International Journal of Non-Linear Mechanics, 895–901,
45, 2010.

[3] K. Tsumoto, T. Yoshinaga, H. Kawakami Bifurcations of
synchronized responses in synaptically coupled Bonhof-
fervan der Pol neurons, Physical Review E, 65, 2002.

[4] Y. Uwate, Y. Nishio Synchronization phenomena in van
der Pol oscillators coupled by a time-varying resistor,
International Journal of Bifurcation and Chaos, 17(10),
35653569, 2007.

[5] J. J. Źebrowski, K. Grudziński, T. Buchner, P. Kuklik, J.
Gac, G. Gielerak Nonlinear oscillator model reproducing
various phenomena in the dynamics of the conduction
system of the heart, Chaos, 17, 2007.

[6] A. M. dos Santos, S. R. Lopes, R. L. Viana Rhythm
synchronization and chaotic modulation of coupled Van
der Pol oscillators in a model for the heartbeat, Physica
A, 338, 335 355, 2004.
http://dx.doi.org/10.1016/j.physa.2004.02.058

Biomath 3 (2014), 1312281, http://dx.doi.org/10.11145/j.biomath.2013.12.281 Page 6 of 7

http://dx.doi.org/10.1016/j.physa.2004.02.058
http://dx.doi.org/10.11145/j.biomath.2013.12.281


B Zduniak, Numerical Analysis of the Coupled Modified Van Der Pol Equations...

[7] L. Henk van der Tweel, F. L. Meijler, F. J. L. van
Capelle Synchronization of the heart, Journal of Applied
Physiology, 34(2), 1973.

[8] D. Postnov, S.K. Han, and H. Kook Synchronization
of diffusively coupled oscillators near the homoclinic
bifurcation, Physical Review E, 60(3), 1999.

http://dx.doi.org/10.1103/PhysRevE.60.2799
[9] K. Grudziński, Modelowanie czynności elektrycznej

ukÅadu przewodzenia serca, Rozprawa Doktorska,
WydziaÅ Fizyki Politechniki Warszawskiej, 2007.

Biomath 3 (2014), 1312281, http://dx.doi.org/10.11145/j.biomath.2013.12.281 Page 7 of 7

http://dx.doi.org/10.1103/PhysRevE.60.2799
http://dx.doi.org/10.11145/j.biomath.2013.12.281

	Introduction
	Mathematical model
	Types of coupling

	Numerical analysis
	Conclusion
	References