Improvement of the Cardiac Oscillator Based Model for the Simulation of Bundle Branch Blocks

Abstract

In this paper, we propose an improvement of the cardiac conduction system based on three modified Van der Pol oscillators. Each oscillator represents one of the components of the heart conduction system: Sino-Atrial node (SA), Atrio-Ventricular node (AV) and His–Purkinje system (HP). However, while SA and AV nodes can be modelled through a single oscillator, the modelling of HP by using a single oscillator is a rough simplification of the cardiac behaviour. In fact, the HP bundle is composed of Right (RB) and Left Bundle (LB) branches that serve, respectively, the right and left ventricles. In order to describe the behaviour of each bundle branch, we build a phenomenological model based on four oscillators: SA, AV, RB and LB. For the characterization of the atrial and ventricular muscles, we used the modified FitzHugh–Nagumo (FHN) equations. The numerical simulation of the model has been implemented in Simulink. The simulation results show that the new model is able to reproduce the heart dynamics generating, besides the physiological signal, also the pathological rhythm in case of Right Bundle Branch Block (RBBB) and Left Bundle Branch Block (LBBB). In particular, our model is able to describe the communication interruption of the conduction system, when one of the HP bundle branches is damaged.

1. Introduction

The human heart is a complex electro-pump that, through cycles of depolarization and repolarization, makes the propagation of the action potential, and consequently the contraction of the cardiac muscle tissue, possible. The electrical activity of the heart can be indirectly measured by the electrocardiogram (ECG), a non-invasive method based on 12 electrodes positioned on the body surface [1,2].

Considering the crucial role that the ECG signal plays in the clinical practice, different techniques were presented in the literature in order to model the dynamics of the heartbeat. The electrical and muscular activities of the heart were analyzed through both mathematical modeling and time series analysis [3]. By using mathematical modeling, it is possible to understand the electrical activation of the heart obtaining, by simulation, both normal and abnormal rhythms [4,5,6]. The generation of synthetic ECG signals with a wide range of waveform shapes and heart rates allows the modeling of the characteristics of each subpart composing the cardiac conduction system, and the understanding of its behaviour [7,8,9,10,11]. The nonlinearity and nonstationarity of the cardiovascular system make the use of nonlinear techniques for the modelling of heart activity useful [12,13,14]. Van der Pol (VdP) and Van der Mark described and modeled the behavior of heart using nonlinear relaxation oscillators [15].

FitzHugh [16] proposed an extended version of the Van der Pol equations for the generation of the action potential. In his work, he proposed a model for the emulation of the signal observed in an excitable cell of a living organism. The model is composed of two coupled, nonlinear ordinary differential equations. The first equation describes the evolution of the neuronal membrane voltage (fast action), while the second one the recovery action of the sodium and potassium channels (slow action) [17]. Grudzinski and Zebrowski modified the Van der Pol model to allow the reproduction of the time series of the action potential generated by a natural pacemaker [18,19]. In their model, the intrinsic frequency of the two pacemakers can be changed without the need to change the length of the refractory period. Although the model allows the manipulation of both the diastolic and the refractory period, the frequency obtained is too low if compared with the physiological values. Gois and Savi [20] proposed a mathematical model always based on the Van der Pol equations, but composed of three modified oscillators. This model is able to emulate the cardiac conduction system composed of a Sino-Atrial node (SA), Atrio-Ventricular node (AV) and His–Purkinje system (HS). The oscillators are connected to each other through time-delayed couplings. This model allows for generating the electrical response of the main cardiac pacemakers and obtaining the ECG wave as a composition of these signals. Although this model reproduces normal and pathological ECG signals, it does not model the activity of atrial and ventricular muscles. The characterization of muscle electrical responses proposed by Ryzhii [21] uses a quiescent excitable FitzHugh–Nagumo-type (FHN) oscillator. Successively, in [22,23], the authors improved this model by including the depolarization and repolarization waves of the atrial and ventricules generated by modified FHN systems for each ECG wave. The advanced model was able to generate normal ECG signals and several well-known rhythm disorders. In particular, it reproduces sinus tachycardia, sinus bradycardia, complete SA–AV block and complete AV–HP block [24].

In the models proposed in [20,21,22], the cardiac conduction system was treated as a network of three modified Van der Pol oscillators representing SA, AV nodes and the HP system. While SA and AV nodes can be modelled by using oscillators, the modeling of the HP system by using a single oscillator is a simplification of the cardiac activity as shown in [20,21,22] where the electrical behaviour of right and left bundle branches that composes the His–Purkinje system (Figure 1a) is modelled by using a single oscillator. By using this approach, it is not possible to characterize bundle branches of electrical conduction diseases. In the physiological condition, when the ventricular contraction occurs, the Purkinje fibers bring the impulse from the left and right bundle branches to the ventricles’ myocardium. In this way, the contraction of ventricles muscle tissue and the ejection of blood outside the heart are allowed. However, the degradation of one of the bundle branches causes a defect called Bundle Branch Block (BBB) that implies an alteration of the pathways for the ventricular depolarization [25].

This work proposes a mathematical improvement of the Rizhii model [22] to better describe the heart rhythm by considering four coupled modified Van der Pol oscillators. Besides the use of two oscillators to describe SA and AV nodes, we use two additional oscillators to describe the right and left bundle branches of the His–Purkinje system. In this way, it is possible to reproduce different pathological conditions, such as the bundle branch blocks. Moreover, the FitzHugh–Nagumo (FHN) oscillators are used to describe the electrical process in the cardiac muscles. We underline that the proposed improved model corresponds to a phenomenological model [26]. The paper is organized as follows: in Section 2, the heart conduction system and mathematical models are shown. In Section 3, results of simulations are shown, and, finally, in Section 5, conclusions are drawn.

2. Heart Conduction System and Mathematical Models

Since the dynamic evolution of the heart action potential are close to the dynamical response of the Van der Pol oscillator (VdP) [15], this approach has been widely used in the modelling of the heart conduction system [18,20,21,27,28,29,30,31]. In our model, we modify the heterogeneous oscillator model used for the generation of the ECG signals as proposed by Ryzhii [21]. In addition to coupled VdP oscillators for the SA and AV nodes modeling, we use two additional oscillators to model the right and left bundle branches of the HP system. Furthermore, we modify the excitable FitzHugo–Nagumo equations to take into account the two additional oscillators (right branch (RB) and left branch (LB)). Figure 2 shows the Ryzhii model and the proposed model. In order to better understand the improved model, in the first part of the following section, we provide a theoretical overview of the heart conduction system and a description of the Ryzhii model.

2.1. Conduction System of the Heart

The cardiac conduction system shown in Figure 1a is composed of a group of muscle cells characterized by their own electrical activity that generates the contraction of the cardiac muscles; the Sino-Atrial node (SA), the Atrio-Ventricular node (AV), branches of the right and left beam, and Purkinje fibers.

The SA node is placed in the right atrium and it is the natural cardiac pacemaker. These cells are known as pacemakers because they manifest a spontaneous depolarization. They are able to reach the action potential threshold faster than every other cell. Thanks to this feature, they manage the heart rhythm [32]. The depolarization rises from the SA node, spreads in the whole right and left atria and finally reaches the AV node.

Then, it proceeds along the bundle of his spreading to both left and right bundle branches depolarizing the left and right ventricles. The electrocardiogram (ECG) is a measure of how electrical activity changes during each cardiac cycle [33]. Figure 1b shows a typical ECG waveform.

The atrial depolarization phase is represented by the P wave, the ventricular depolarization by the QRS complex, and the ventricular repolarization by the T wave. In physiological conditions, the SA node works at 60–100 bpm, faster than any other cell; the only conducting path from atria to ventricles is provided by the AV node. If the SA node has a disease and it fails to generate impulses, the cells of the AV node can stimulate the heart at a rate of about 40–60 bpm [32]. If this node also presents damaged cells, the atrial-ventricular conduction block occurs. Depending on the entity of the damage, the block can be complete or partial. In this case, not all the electrical signals generated by the SA node are transmitted to the HP system. If no excitation signal is delivered from SA or AV nodes, the HP cells can fire at a rate of about 20–40 bpm. However, if the right or left bundle branch of the HP system is corrupted, the respective ventricle is not contracted.

This condition is known as Bundle Branch Block (BBB) and it implies an alteration of the pathways for ventricular depolarization. When this condition occurs, the electrical impulse can move through muscle fibers both slowing the electrical activity and changing the propagation direction of the pulses [34].

2.2. Mathematical Model

Ryzhii in [22] describes the three natural pacemakers SA, AV and the HP bundle by a system of modified VdP equations:

$$\mathbf{SA}=\left\{\begin{array}{c}{\dot{x}}_1=y_1,\hfill \\ \begin{array}{c}\hfill {\dot{y}}_1=-a_1{y}_1(x_1-u_{11})(x_1-u_{12})-f_1{x}_1(x_1+d_1)(x_1+e_1),\end{array}\hfill \end{array}\right.$$

(1)

$$\mathbf{AV}=\left\{\begin{array}{c}{\dot{x}}_2=y_2,\hfill \\ \begin{array}{c}\hfill {\dot{y}}_2=-a_2{y}_2(x_2-u_{21})(x_2-u_{22})-f_2{x}_2(x_2+d_2)(x_2+e_2)+K_{SA-AV}(y_1^{{\tau }_{SA-AV}}-y_2),\end{array}\hfill \end{array}\right.$$

(2)

$$\mathbf{HP}=\left\{\begin{array}{c}{\dot{x}}_3=y_3,\hfill \\ \begin{array}{c}\hfill {\dot{y}}_3=-a_3{y}_3(x_3-u_{31})(x_3-u_{32})-f_3{x}_3(x_3+d_3)(x_3+e_3)+K_{AV-HP}(y_2^{{\tau }_{AV-HP}}-y_3),\end{array}\hfill \end{array}\right.$$

(3)

where $x_i\left(t\right)$ and $y_i\left(t\right)$ correspond to the action potential and the transmembrane currents of the heart, $a_i(x_i-u_{ij})(x_i-u_{ik})$ are the damping factors, $f_i{x}_i(x_i-d_{ij})(x_i-e_{ik})$ are the harmonic force terms, $a_i>0$, $u_{ij}$ represent the nonlinear damping force parameters, $f_i$ are the parameters related to the intrinsic frequency of the oscillator, the coupling coefficients $K_{SA-AV}$ and $K_{AV-HP}$ represent the unidirectional coupling between the SA, AV and HP pacemakers, $y_i^{{\tau }_n}\equiv y_i(t-{\tau }_n)$ are the velocity coupling components of the time-delay signal, and ${\tau }_n$ are the time delays [18,20]. The synchronism between the three oscillators depends on the coupling coefficients K [35,36]. A lot of coupling methods can be found in the literature [37,38,39]. However, it is possible to avoid the use of delays ${\tau }_n$ in the coupling terms by choosing the appropriate coupling coefficients $K_{SA-AV}$ and $K_{AV-HP}$. In conclusion, the VdP Equations (1)–(3) become

$$\mathbf{SA}=\left\{\begin{array}{c}{\dot{x}}_1=y_1,\hfill \\ \begin{array}{cc}\hfill {\dot{y}}_1=& -a_1{y}_1(x_1-u_{11})(x_1-u_{12})-f_1{x}_1(x_1+d_1)(x_1+e_1),\hfill \end{array}\hfill \end{array}\right.$$

(4)

$$\mathbf{AV}=\left\{\begin{array}{c}{\dot{x}}_2=y_2,\hfill \\ \begin{array}{cc}\hfill {\dot{y}}_2=& -a_2{y}_2(x_2-u_{21})(x_2-u_{22})-f_2{x}_2(x_2+d_2)(x_2+e_2)+K_{SA-AV}(x_1-x_2),\hfill \end{array}\hfill \end{array}\right.$$

(5)

$$\mathbf{HP}=\left\{\begin{array}{c}{\dot{x}}_3=y_3,\hfill \\ \begin{array}{cc}\hfill {\dot{y}}_3=& -a_3{y}_3(x_3-u_{31})(x_3-u_{32})-f_3{x}_3(x_3+d_3)(x_3+e_3)+K_{AV-HP}(x_2-x_3).\hfill \end{array}\hfill \end{array}\right.$$

(6)

Starting from this model, we propose improvements to describe the phenomenological behaviour of the right and left bundle branch. Leaving the SA and AV node equations (Equations (4) and (5)) unaltered, we describe the right bundle branch (RB) and the left bundle branch (LB) as

$$\mathbf{SA}=\left\{\begin{array}{c}{\dot{x}}_1=y_1,\hfill \\ \begin{array}{cc}\hfill {\dot{y}}_1=& -a_1{y}_1(x_1-u_{11})(x_1-u_{12})-f_1{x}_1(x_1+d_1)(x_1+e_1),\hfill \end{array}\hfill \end{array}\right.$$

(7)

$$\mathbf{AV}=\left\{\begin{array}{c}{\dot{x}}_2=y_2,\hfill \\ \begin{array}{cc}\hfill {\dot{y}}_2=& -a_2{y}_2(x_2-u_{21})(x_2-u_{22})-f_2{x}_2(x_2+d_2)(x_2+e_2)+K_{SA-AV}(x_1-x_2),\hfill \end{array}\hfill \end{array}\right.$$

(8)

$$\mathbf{RB}=\left\{\begin{array}{c}{\dot{x}}_{3_{RB}}=y_{3_{RB}},\hfill \\ \begin{array}{c}\hfill {\dot{y}}_{3_{RB}}=-a_3{y}_{3_{RB}}(x_{3_{RB}}-u_{31})(x_{3_{RB}}-u_{32})-f_3{x}_{3_{RB}}(x_{3_{RB}}+d_3)(x_{3_{RB}}+e_3)+K_{AV-RB}(x_2-x_{3_{RB}}),\end{array}\hfill \end{array}\right.$$

(9)

$$\mathbf{LB}=\left\{\begin{array}{c}{\dot{x}}_{3_{LB}}=y_{3_{LB}},\hfill \\ \begin{array}{c}\hfill {\dot{y}}_{3_{LB}}=-a_3{y}_{3_{LB}}(x_{3_{LB}}-u_{31})(x_{3_{LB}}-u_{32})-f_3{x}_{3_{LB}}(x_{3_{LB}}+d_3)(x_{3_{LB}}+e_3)+K_{AV-LB}(x_2-x_{3_{LB}}).\end{array}\hfill \end{array}\right.$$

(10)

In the physiological condition, the right and left bundle branches are synchronized and the sum of their behaviour has to be the same as the HP oscillator:

$$y_{3_{HP}}={\alpha }_1{y}_{3_{RB}}+{\alpha }_2{y}_{3_{LB}},$$

(11)

where ${\alpha }_1$ and ${\alpha }_2$ are coefficients used to adjust the oscillation amplitudes. In order to obtain the physiological signal, we set ${\alpha }_1=0.5$ and ${\alpha }_2=0.5$. In this manner, each oscillator is weighted at 50% in the generation of the signal. This is because, in a normal rhythm, the contribution of the left and the right bundle branch is the same.

The parameters $a_i$, $u_{ij}$, $f_i$, $d_i$, $e_i$, are chosen to obtain intrinsic oscillation rates of 70 bpm, 50 bpm, 35 bpm and 35 bpm for uncoupled SA, AV, RB and LB, respectively (Figure 3) and with behaviour close to action potentials of real pacemakers: [40,41,42].

In particular, we use the following experimental parameters: $a_1=40$, $a_2=50$, $a_3=50$, $a_1=40$, $u_{11}=0.83$, $u_{21}=0.83$, $u_{31}=0.83$, $u_{12}=-0.83$, $u_{22}=-0.83$, $u_{32}=-0.83$, $f_1=25$, $f_2=8.4$, $f_3=1.5$, $d_1=3$, $d_2=3$, $d_3=3$, $e_1=3.5$, $e_2=5$, $e_3=12$. The coupling coefficients are: $K_{SA-AV}=100$, $K_{AV-RB}=285$, $K_{AV-LB}=285$.

We describe the depolarization and repolarization process for atrial (AT) and ventricular (VN) muscles starting from the modified FitzHugh–Nagumo model [16,43] proposed by Ryzhii in [22]:

$$\mathbf{P}\mathbf{wave}=\{\begin{array}{c}\begin{array}{cc}\hfill {\dot{z}}_1=& k_1(-c_1{z}_1(z_1-w_{11})(z_1-w_{12})-b_1{v}_1-d_1{v}_1{z}_1+I_{A{T}_{De}}),\hfill \end{array}\hfill \\ {\dot{v}}_1=k_1{h}_1(z_1-g_1{v}_1),\hfill \end{array}$$

(12)

$$\mathbf{Ta}\mathbf{wave}=\{\begin{array}{c}\begin{array}{cc}\hfill {\dot{z}}_2=& k_2(-c_2{z}_2(z_2-w_{21})(z_2-w_{22})-b_2{v}_2-d_2{v}_2{z}_2+I_{A{T}_{Re}}),\hfill \end{array}\hfill \\ {\dot{v}}_2=k_2{h}_2(z_2-g_2{v}_2),\hfill \end{array}$$

(13)

$$\mathbf{QRS}=\{\begin{array}{c}\begin{array}{cc}\hfill {\dot{z}}_3=& k_3(-c_3{z}_3(z_3-w_{31})(z_3-w_{32})-b_3{v}_3-d_3{v}_3{z}_3+I_{V{N}_{De}}),\hfill \end{array}\hfill \\ {\dot{v}}_3=k_3{h}_3(z_3-g_3{v}_3),\hfill \end{array}$$

(14)

$$\mathbf{T}\mathbf{wave}=\{\begin{array}{c}\begin{array}{cc}\hfill {\dot{z}}_4=& k_4(-c_4{z}_4(z_4-w_{41})(z_4-w_{42})-b_4{v}_4-d_4{v}_4{z}_4+I_{V{N}_{Re}}),\hfill \end{array}\hfill \\ {\dot{v}}_4=k_4{h}_4(z_4-g_4{v}_4),\hfill \end{array}$$

(15)

where $k_i$ are the scaling coefficients. In our work, we used the experimental parameters: $k_1=2\ast {10}^3$, $k_2=4\ast {10}^4$, $k_3={10}^4$, $k_4=2\ast {10}^3$, $c_1=c_2=0.26$, $c_3=0.3$, $c_4=0.1$, $b_1=b_2=b_4=0$, $b_3=0.015$, $d_{m_1}=d_{m_2}=0.4$, $d_{m_3}=0.09$, $d_{m_4}=0.1$, $h_1=h_2=0.004$, $h_3=h_4=0.008$, $g_1=g_2=g_3=g_4=1$, $w_{11}=0.13$, $w_{12}=1.0$, $w_{22}=1.0$, $w_{21}=0.19$, $w_{31}=0.12$, $w_{32}=1.1$, $w_{41}=0.22$, $w_{42}=0.8$.

The activation currents $I_i$ represent the coupling between the SA and the AT muscles and between the RB–LB pacemakers and the VN muscles. With respect to the Ryzhii model [22], we adjust the activation currents of the QRS complex and T wave in order to have the HP oscillator composed of the RB and LB oscillators:

$$I_{{AT}_{De}}=\{\begin{array}{cc}0,\hfill & \mathrm{for}{y}_1\le 0,\hfill \\ K_{A{T}_{De}}{y}_1,\hfill & \mathrm{for}{y}_1>0,\hfill \end{array}$$

(16)

$$I_{{AT}_{Re}}=\{\begin{array}{cc}-K_{A{T}_{Re}}{y}_1,\hfill & \mathrm{for}{y}_1\le 0,\hfill \\ 0,\hfill & \mathrm{for}{y}_1>0,\hfill \end{array}$$

(17)

$$I_{{VN}_{De}}=\{\begin{array}{cc}0,\hfill & \mathrm{for}{y}_{3_{HP}}\le 0,\hfill \\ K_{V{N}_{De}}(y_{3_{RB}}^{tot}+y_{3_{LB}}^{tot}),\hfill & \mathrm{for}{y}_{3_{HP}}>0,\hfill \end{array}$$

(18)

$$I_{{VN}_{Re}}=\{\begin{array}{cc}-K_{V{N}_{Re}}(y_{3_{RB}}^{tot}+y_{3_{LB}}^{tot}),\hfill & \mathrm{for}{y}_{3_{HP}}\le 0,\hfill \\ 0,\hfill & \mathrm{for}{y}_{3_{HP}}>0.\hfill \end{array}$$

(19)

In particular,

$$\begin{array}{c}{y}_{3_{RB}}^{tot}={\alpha }_1{y}_{3_{RB}}+{\alpha }_3{y}_{3_{LB}}^{{\tau }_{LB}},\hfill \\ y_{3_{LB}}^{tot}={\alpha }_2{y}_{3_{LB}}+{\alpha }_3{y}_{3_{RB}}^{{\tau }_{RB}},\hfill \end{array}$$

(20)

where $y_{3_{HP}}=(y_{3_{RB}}+y_{3_{LB}})$, $K_{A{T}_{De}}=9\ast {10}^{-5}$, $K_{A{T}_{Re}}=1.5\ast {10}^{-5}$, $K_{V{N}_{De}}=9\ast {10}^{-5}$, $K_{V{N}_{Re}}=9\ast {10}^{-5}$, ${\alpha }_1=0.5$, ${\alpha }_2=0.5$ and ${\alpha }_3=0.001$.

In Equation (20), the terms $y_{3_{RB}}$ and $y_{3_{LB}}$ represent the impulses traveling through the right and left ventricles and through the bundle branches. The terms $y_{3_{RB}}^{{\tau }_{RB}}$ and $y_{3_{LB}}^{{\tau }_{LB}}$ represent the impulses traveling through the right and left ventricles through the myocardium cells. These secondary pathways are slower than the main path propagated through the bundle branches. In the physiological case, the contribution of these impulses is negligible with respect to the contribution of the main pathways [12]. On the other hand, in the pathological case, where one branch is damaged, the secondary pathways guarantee the contraction of the ventricle not served by the main pathway. To provide a complete model, we need to take into account the contribution of the secondary pathways. It can be described as an impulse that travels through the bundle branches delayed and attenuated by a multiplicative constant ${\alpha }_3$. We note that, in the actual case, these secondary impulses travel through the myocardium cells and the propagation of electrical activity should be described with a specific mathematical model. However, the scope of our work is to provide a phenomenological model. For this reason, we assumed that the impulses that travel through the myocardium cells can be modeled as an impulse attenuated and delayed in time.

In the model of the left bundle branch block (LBBB), the activation currents $I_{V{N}_{De}}$ and $I_{V{N}_{Re}}$ represent respectively the coupling between the RB and LB pacemakers with the VN muscles; consequently, Equations (18) and (19) become

$$I_{{VN}_{De}}^{LBBB}=\{\begin{array}{cc}0,\hfill & \mathrm{for}{y}_{3_{RB}}\le 0,\hfill \\ K_{V{N}_{De}}({\alpha }_1{y}_{3_{RB}}+{\alpha }_3{y}_{3_{RB}}^{{\tau }_{RB}}),\hfill & \mathrm{for}{y}_{3_{RB}}>0,\hfill \end{array}$$

(21)

$$I_{{VN}_{Re}}^{LBBB}=\{\begin{array}{cc}-K_{V{N}_{Re}}{\alpha }_1{y}_{3_{RB}},\hfill & \mathrm{for}{y}_{3_{RB}}\le 0,\hfill \\ 0,\hfill & \mathrm{for}{y}_{3_{RB}}>0.\hfill \end{array}$$

(22)

For LBBB simulation, we set the experimental parameters: $K_{SA-AV}=100$, $K_{AV-BR}=285$, $K_{AV-LB}=0$, $K_{A{T}_{De}}=9\ast {10}^{-5}$, $K_{A{T}_{Re}}=3\ast {10}^{-5}$, $K_{V{N}_{De}}=9\ast {10}^{-5}$, $K_{V{N}_{Re}}=-8\ast {10}^{-5}$, ${\tau }_{RB}=0.78$, ${\alpha }_1=1$, ${\alpha }_3=1$ and ${\alpha }_2=0$.

In the model of right bundle branch block (RBBB) the activation currents, $I_{V{N}_{De}}$ and $I_{V{N}_{Re}}$ are

$$I_{{VN}_{De}}^{RBBB}=\{\begin{array}{cc}{K}_{V{N}_{De}}({\alpha }_2{y}_{3_{LB}}+{\alpha }_3{y}_{3_{LB}}^{{\tau }_{LB}}),\hfill & \mathrm{for}{y}_{3_{LB}}<0,\hfill \\ K_{V{N}_{De}}{\alpha }_2{y}_{3_{LB}},\hfill & \mathrm{for}{y}_{3_{LB}}>0,\hfill \end{array}$$

(23)

$$I_{{VN}_{Re}}^{RBBB}=\{\begin{array}{cc}-K_{V{N}_{Re}}{\alpha }_2{y}_{3_{LB}},\hfill & \mathrm{for}{y}_{3_{LB}}<0,\hfill \\ 0,\hfill & \mathrm{for}{y}_{3_{LB}}>0.\hfill \end{array}$$

(24)

For RBBB simulation, we set the experimental parameters: $K_{SA-AV}=100$, $K_{AV-BR}=0$, $K_{AV-LB}=285$, $K_{A{T}_{De}}=9\ast {10}^{-5}$, $K_{A{T}_{Re}}=2\ast {10}^{-5}$, $K_{V{N}_{De}}=5.5\ast {10}^{-5}$, $K_{V{N}_{Re}}=7\ast {10}^{-5}$, ${\tau }_{RB}=0.05$, ${\alpha }_1=0$, ${\alpha }_2=1$ and ${\alpha }_3=1$.

Finally, we obtain the synthetic ECG signal as a composition of the AT and VN waveforms

$$ECG=z_0+z_1-z_2+k_R{z}_3+z_4,$$

(25)

where $z_0=0.2$ provides the adjustment of the baseline and $k_R$ is a multiplicative coefficient to modulate the amplitude of the R peak.

3. Results

In order to validate our model, we performed different Simulink simulations reproducing either normal and pathological signals. Simulations are performed using a fixed step Runge–Kutta solver and a fixed step size of 0.0001 seconds. Figure 4 shows the normal condition where all pacemakers are dominated by the SA node. As a consequence, the heart rate follows the SA rate (70 bpm). In this case, both the right and left bundle branches are working, and they give the same contribution. The experimental parameters are the same used in the non-pathological as described in Section 2.

Then, we performed a simulation of a pathological condition of the His bundle. We obtain this condition by changing the physiological parameters with the pathological parameters in order to obtain both complete right bundle branch block (RBBB) and complete left bundle branch block (LBBB). The pathological parameters used are shown in Section 2.

Left Bundle Branch BlockFigure 5 shows the left bundle branch block (LBBB). In this case, the left ventricle is not directly activated by the impulses travelling through the left bundle branch. The right ventricle is normally activated by the right bundle branch signal. These impulses are able to travel through the myocardium of the right ventricle to the left ventricle and depolarize it (yellow arrows in Figure 5). This activation extends the QRS duration to >120 ms [32,33]. The direction of the slow depolarization (from the right to the left) produces tall R-waves in the lateral leads (I, aVL, V5-6) of the electrocardiogram and deep S-waves in the right precordial leads (V1-3). Since the ventricles are activated sequentially (first right then left) rather than simultaneously, this produces a broad or notched (‘M’-shaped) R-wave in the lateral leads (I, aVL, V5-6).

Figure 6a shows simulation results for the LBBB case. Experimental parameters: $K_{SA-AV}=100$, $K_{AV-BR}=285$, $K_{AV-LB}=0$, $K_{A{T}_{De}}=9\ast {10}^{-5}$, $K_{A{T}_{Re}}=3\ast {10}^{-5}$, $K_{V{N}_{De}}=9\ast {10}^{-5}$, $K_{V{N}_{Re}}=-8\ast {10}^{-5}$, ${\tau }_{RB}=0.78$, ${\alpha }_1=1$, ${\alpha }_3=1$ and ${\alpha }_2=0$.

In this case, the term $y_{3_{LB}}$ in Equation (21) is not present because the contribution of the left bundle branch is zero. On the other hand, the impulse generated by the right bundle brunch $y_{3_{RB}}$ is able to both normally depolarize the right ventricle and to travel towards the left ventricle through the myocardium cells. As this second conduction way is slower than the main bundle of His–Purkinje fibres, the term $y_{3_{RB}}$ in Equation (21) (that represents the left ventricle depolarization) is delayed through the coefficient ${\tau }_{RB}$. This produces a broad or notched (‘M’-shaped) R-wave in the ECG signal (Figure 6b). It is possible to compare simulation results with a real patient’s ECG signal affected by LBBB shown in Figure 6c.

Right Bundle Branch Block In the Right Bundle Branch Block (RBBB), the activation of the right ventricle is delayed because the depolarization has to spread from the left ventricle (yellow arrows in Figure 7). The left ventricle is activated normally and the early part of the QRS complex is unchanged. The delayed right ventricular activation produces a secondary R-wave (R’) in the right precordial leads (V1-3) and a wide, slurred S-wave in the lateral leads [32,33] (Figure 7).

In this case, the term $y_{3_{RB}}$ in Equation (24) is not present because the contribution of the right bundle branch is zero. In this case, the impulse is generated by the left bundle brunch $y_{3_{LB}}$ that is able to both normally depolarize the left ventricle and to travel towards the right one through the myocardium cells. Figure 8a shows the simulation results for the RBBB case. For this simulation, we set the following experimental parameters: $K_{SA-AV}=100$, $K_{AV-BR}=0$, $K_{AV-LB}=285$, $K_{A{T}_{De}}=9\ast {10}^{-5}$, $K_{A{T}_{Re}}=2\ast 10$${}^{-5}$, $K_{V{N}_{De}}=5.5\ast {10}^{-5}$, $K_{V{N}_{Re}}=7\ast {10}^{-5}$, ${\tau }_{LB}=0.05$, ${\alpha }_1=0$, ${\alpha }_2=1$ and ${\alpha }_3=1$. In addition, in this case, it is possible to compare simulation results with a real patient’s ECG signal affected by RBBB shown in Figure 8c.

4. Discussion

The model discussed in this paper is an extension of the one proposed in [22]. Both the models (the original proposed by Ryzhii and the one proposed in this paper) are able to generate realistic ECG signals in order to simulate normal and pathological rhythms. The analysis of the synthetic ECG signal generated by the two models can be performed like a real signal of the patient, which is conducted by observing the morphology of different waves and their behavior in time, with specific interest on the QRS complex. This means that, for RBBB and LBBB, as for the other pathologies discussed in this paper and in [22], the anomalies are identified analyzing the shape of the signal and identifying the interest points in which the waves change their morphology. For example, if we want to reproduce the LBBB, we need to modify the shape of the QRS complex, as well as the polarity of the T wave [46]. As shown in Figure 6, selecting proper parameters, our model is able to reproduce the main typical characteristics of the signal related to this pathology.

However, it is important to note that the synthetic signal may not be exactly the same as the patient’s actual ECG signal. To generate a synthetic signal closer, but not equal, to the real one, it is possible to modify, inside a specific operating range, the model parameters to obtain different nuances of the same pathological behaviour. This difference is due to the uniqueness of the signal generated by the cardiac activity of each person. It means that it is possible to compare two signals belonging to two different people by observing the waveforms in time and searching the features of interest, but not comparing the shape of the tracks point by point. In addition to personal traits, there are several factors that make any real signal different from any other—for example, the noise due to the acquisition method, the constitution of patient that influenced the transmission of the electrical signal from heart to body surface, the coexistence of different pathologies, and so on. For this reason, nowadays, the interpretation of an ECG signal is mainly performed by physicians by observing the specific characteristics of tracks.

5. Conclusions

In the paper, we improve the mathematical model of the heart by using four modified Van der Pol oscillators to describe the electrical activity and we modify the FitzHugh–Nagumo equations to reproduce the behaviour of the cardiac muscle cells. Each oscillator represents one of the main natural pacemakers: Sino-Atrial node (SA), Atrio-Ventricular node (AV), Right Bundle Branch (RBB) and Left Bundle Branch (LBB).

Numerical simulations show that the proposed model, based on four oscillators, is able to reproduce the heart behavior by generating the physiological signal like the model proposed by Ryzhii in [22] (based on three oscillators). The main improvement of the proposed model is the description of the His–Purkinje system with two additional oscillators, in order to better describe the behavior of each bundle branch. The model is able to reproduce the communication interruption in the heart electrical conduction system when one of the His bundle branches does not work. In this case, the model reproduces both normal ventricle contraction due to healthy bundle branch and the propagation of action potential through myocardial cells to the ventricle where the bundle branch is pathologically altered.

The waveforms obtained by simulation are comparable with the real ECG signals of patients affected by the LBBB and RBBB pathologies. In other words, the clinical phenotypes observed in the patient’s signals are similar to waveform features obtained with our model. In conclusion, the proposed model makes the study of interactions between the main components of the heart possible. It allows simulation and evaluation of heart activity and dynamics under different types of pacemaker coupling. These aspects are widely useful in clinical practice.