Mathematical model of the atrioventricular nodal double response tachycardia and double-fire pathology.

A proposed model consisting of two coupled models (Hodgkin-Huxley and Yanagihara-Noma-Irisawa model) is considered as a description of the heart's action potential. System of ordinary differential equations is used to recreate pathological behaviour in the conducting heart's system such as double fire and the most common tachycardia: atrioventricular nodal reentrant tachycardia (AVNRT). Part of the population has an abnormal accessory pathways: fast and slow (Fujiki, 2008). These pathways in the atrioventricular node (AV node) are anatomical and functional contributions of supraventricular tachycardia. However, the appearance of two pathways in the AV node may be a contribution of arrhythmia, which is caused by coexistent conduction by two pathways (called double fire). The difference in the conduction time between these pathways is the most important factor. This is the reason to introduce three types of couplings and delay to our system in order to reproduce various types of the AVNRT. In our research, introducing the feedback loops and couplings entails the creation of waves which can correspond to the re-entry waves occurring in the AVNRT. Our main aim is to study solutions of the given equations and take into consideration the influence of feedback and delays which occur in these pathological modes. We also present stability analysis for both components, that is Hodgkin-Huxley and Yanagihara-Noma-Irisawa models, as well as for the final double-fire model.

1. Introduction. The paper concerns a research of an electrical conduction system of the human heart. Modeling the formation and conduction of electrical impulses in the heart is one of the most developed areas of mathematical biology. For years the most common models of action potentials that occur in the heart have included the Hodgkin-Huxley [14], van der Pol [7] and the Purkinje cells model [17]. These models allow to reconstruct the dynamic of action potential which occurs in the cardiac conduction system. In this paper we propose ordinary differential system which is based on the Hodgkin-Huxley models and allows to reconstruct pathological behaviors in the conducting heart's system, such as rare arrhythmia and various types of AVNRT. A problem with making the appropriate diagnosis, and therefore the problem of effective treatment of the disease was the motivation of our research. This kind of problem is observed in various types of AVNRT and double-fire. The reasons of those problems lie in not fully known and understood mechanisms of these pathologies. Also the symptoms are often mistakenly taken for other heart disease (double-fire phenomena are often recognized as one of the types of AVNRT, atrial fibrillation or bigeminia). It should be emphasized that the clinical pictures of these diseases are also non-specific. In the last few years we have understood that the construction of the AV node has a multilevel architecture in which there may be many pathways (slow and fast) at different locations in the AV node [16,13]. It helped to recognize many types of AVNRT, which until now were understood as one, although the mechanisms of action were different. In the literature, there are no mathematical models that would specifically model various types of AVNRT. There were attempts of modeling AVNRT only as a single pathology having the pathway of slow and fast type [23,18]. Part of the population has abnormal accessory pathways: fast and slow; cf. [9,15]. These pathways in the AV node are anatomical and functional contributions of the most common supraventricular tachycardia, which is a re-entry tachycardia from the AV node. Clearly, the atrioventricular nodal re-entrant tachycardia is caused by re-entries. AVNRT appears when two electric pathways occur in and around the AV node (for example slow and fast pathways). That gives a way to the occurrence of re-entry. We can distinguish five different forms of the AVNRT; typical: slow/fast, athypical: fast/slow and other forms: slow/slow, more than two re-entry waves, one fast pathway with depolarization of slow pathway [4,8,19]. It depends on multilevel architectonics of the AV node. AVNRT circuit involves larger areas including atrioventricular junction, adjacent atrial structures and in particular so called atrial inputs including at least antero-superior and postero-inferior entries, sometimes also left atrial entry. However, the appearance of two pathways in the AV node may be a contribution of rare arrhythmia, which is caused by coexistent conduction of two pathways. It is called double fire [11] and was first described by Wu et al. in 1975 [21]. It leads to doubling the frequency of ventricle's rhythm. An electrophysiological defect of the AV node is the cause of the doubling conduction and the aforementioned arrhythmia. It is connected with an unidirectional (backward) block in the slow pathway, which prevents backward depolarization. The conduction in the slow pathway is very delayed because the His-Purkinje system must step out with a refraction period after excitation by the fast pathway. The difference in the conduction time between these pathways is the most important factor. Double-fire tachycardia is an uncommon and under-recognized entity with approximately 50 cases reported in the literature [12]. In [11], the author suggests that this is not a rare arrhythmia but rarely recognized and that the main problem lies in the lack of bibliography on the subject, as well as the lack of a precise description of this phenomenon even in text-books of cardiology. Recent articles [6,2,13] about this arrhythmia consider individual clinical cases which are recognized as double-fire mainly on the basis of resistance to pharmacological treatment and for ablation (in the case of double-fire a different type of ablation must be carried out than for atrial fibrillation and other arrhythmias such as AVNRT). Up to our knowledge, there is no mathematical model of the phenomenon of double-fire. In order to better understand the mechanisms that causes these arrhythmias in the paper we propose mathematical models of the AVNRT and double-fire. 2. Mathematical model. In this paper we propose a model consisting of two coupled models. The first model is basic Hodgkin-Huxley model (HH) which reconstructs the pacemaker potential, and the second one is Yanagihara-Noma-Irisawa model (YNI) which reconstructs a sino-atrial node pacemaker potential based on the voltage clamp experiments [22]. We understand that, in our system, YNI model describes the slow pathway which has sinusoidal rhythm (characteristic of this rhythm is similar to the rhythm in the AV node). We also introduce three types of couplings (bidirectional coupling, unidirectional coupling and feedback) and delays in order to reproduce various types of the AVNRT. Introducing feedback loops entails the creation of waves which can correspond to re-entry waves. The aim of this research is to recreate physiological properties of the biological system using differential equations.  [10,5]. Current can be carried through the circuit as ions passing through the membrane or by charging the capacitors of the membrane. From the definition of capacity and the second Kirchoff's law we obtain the first equation of this model: where C m is the membrane capacitance, V is the intracellular potential,V is its derivative with respect to time t, I ion is the net ionic current flowing across the membrane and I ext is the externally applied current.
The ionic current is divided into three currents: sodium, potassium and leakage current. The movement of each of these currents is proportional to the conductance what is presented by the given dependency: I k = G k (V − E k ). The conductance is caused by the opening of many microscopic channels in the membrane. Each individual membrane contains many gates and when all gates are in the permissive state then the channel is considered to be open. The probability of a gate being in the permissive state depends on the current value of the membrane voltage. To develop the differential equations that describe the conductance, the probability of a gate being open is defined as p i for any ion, i. The probability of a gate being close is therefore 1 − p i . From these assumptions we obtain an equation of the first order kinetics which scheme is similar to the kinetics of chemical reactions characterized by rate constants between the given states: where α i and β i are the rate constants dependent on the voltage that describes the transient rates of permissive and non-permissive gates. In this model, we have three types of gates denoted by m, n and h. In each channel we have a specific amount of these gates. Summarizing aforementioned assumptions, we can create the classical Hodgkin-Huxley model, which consists of four equations that reads: where V 1 is the membrane potential, n, m, h are gating variables, α x , β x denote the rate constants of opening and closing of the ionic channel x. The forms of the rate constants were proposed on the basis of experimental studies [10] and read: Parameters used in the classical Hodgkin-Huxley model [10] are the following: In the model we propose in the next section, the HH model is used with the parameters proposed in [1], that is: 2.1.1. Stability analysis of the HH model. There are many papers focused on specific properties of the HH model. Below we present some interesting results concerning Eqs. (1). However, it should be noticed that they were obtained for various parameter values, not necessarily proposed in [10].
In the analysis of differential equations, stability and equilibrium points are one of the most important aspects. We mainly focus on finding equilibria as they may be the final states of the system. For Eqs. (1) equilibria are uniquely defined by zeros of the following function: On the other hand, as in the classical model the applied current depends on time, we can also choose different values. For comparison we take I ext = 6, and then we obtain another equilibrium point (V * 1 , n * , m * , h * ) ≈ (11.7481, 0.9239, 0.9845, 0.0020). From the characteristic equation we are able to calculate eigenvalues for both equilibrium points, which are respectively: For parameter values given in (3), there is again only one equilibrium point of Eqs. (1). Uniqueness of this point with coordinates (V * 1 , n * , m * , h * ) ≈ (−54.6192, 0.4023, 0.09765, 0.405) could be seen from the graph of the function F (V 1 ) (cf. Fig. 1 left), because this function is monotonically decreasing. For this case eigenvalues read: According to the stability theory, the equilibrium point is unstable, and this is the main reason of oscillatory dynamics of the HH model; cf. Fig. 1 right. Now, we focus on another results that could be found in the literature. In [20], the authors used nonlinear systems theory to look for equilibrium points of the HH model and analyse their stability. They focused on the influence of temperature on the model parameters. Varying temperature they also varied parameters, such that they obtained up to 3 equilibrium points (A, B, C), where the point A remains stable when the parameters are varied, B changes stability, while C is always unstable, but can be an unstable node or focus.
The paper [3] discusses the role of change of the external current (I ext ) for an unique steady state voltage in the HH model. It occurs that for the rest of parameters fixed as in [10] this model has a unique steady state voltage V ss for each value of the current I ext . Then the model is linearized around the steady state. The dependence of eigenvalues on I ext is shown in Fig. 2, which is reproduced from [3]. In Figure 2. Stability of steady states and periodic orbits in the HH model. A: Real part of the complex pair of eigenvalues Re(λ 1,2 ) and third negative eigenvalue λ 3 . λ 4 is even more negative such that it is out off the scale of this plot. For I 1 < I ext < I 2 , (λ 1,2 ) is positive (dashed), showing instability of the steady state. B: Values of the leading non-trivial Floquet multiplier along the branch of periodic solutions in log-log axes. The part of the curve with σ 1 > 1 (dashed) indicates instability of the orbit [3].
general the system is stable, but there is a narrow range of the current I 1 < I ext < I 2 for which the steady state is unstable. In this problem there are two Hopf bifurcations; the first bifurcation appears when the steady state loses stability for the value of I ext increasing from minimal values to the left end of the range (I 1 ). The second one takes place when I ext increases through top value of the given range I 2 and the steady state regains stability [3]. The stability of periodic solutions is also described in [3], where it is determined according to the Floquet theory. In this case we also linearize the HH model, but about the periodic solution. The linearized equation has the formẏ = J osc y, where J osc has entries that are T -periodic. The solutions have the form exp(λt)q(t), where q is periodic and exp(λT ) = σ is the Floquet multiplier. If any σ has an absolute value greater than 1 it means that the periodic solution is unstable. In Fig. 2 the non-trivial multiplier σ 1 is plotted vs. I ext along the branch of periodic solutions. In the graph we see two periodic solutions, stable and unstable one. From the analysis presented in [3] we can conclude that the system is bistable, that is the stable steady state and the stable limit cycle coexist.
Stable states indicate that the electrophysiological activity of cell will get to corresponding resting state at last, while periodic phenomena look like a response of pathological cell's action potentials caused by cardiac arrhytmias [24].

2.2.
Yanagihara-Noma-Irisawa model. The next model which is used in our research is Yanagihara-Noma-Irisawa model [22]. The pacemaker activity of the S-A node cell was explained by reconstructing the pacemaker potential using the Hodgkin-Huxley type mathematical model which was based on the reported voltage clamp data. In this model five dynamic currents are taken into account: the sodium current, I N a , the slow inward current, I s , the potassium current, I K , and the hyperpolarization-activated current, I h , including a time-independent leak current, I l . The main difference between this model and the classical HH model is connected with more types of currents flowing through ion channels and, what goes with it, also more types of gates. Hence, the character and the idea of equations are analogous to the case of the reference HH model, although it is particularized, and hence the model consists of seven equations describing membrane potential V 2 , and gating variables d, f , m 1 , h 1 , q, p. As before, we denote by α x and β x the rate constants of opening and closing of the ionic channel (where x = d, f , m 1 , h 1 , q, p), and by I m the total current passing through the unit membrane. With these notations the differential model YNI reads: where: where the currents are treated as functions of V 2 with other variables calculated at the steady state, that is x = αx(V2) αx(V2)+βx(V2) for any x representing gating variable. Fig. 3 shows the graph of G. We see that G is monotonic, and therefore there is only one equilibrium. In the case we study the equilibrium point has coordinates Looking for stability of this equilibrium we obtain the following eigenvalues: implying stability of the equilibrium. Although the equilibrium is stable, we should remember that it is local stability, and we can expect also other model dynamics, like oscillatory behaviour shown in Fig. 4. On the other hand, the YNI model does not possess any stable equilibrium, but have a stable periodic solution for parameters given in [5]. One-parameter bifurcation diagram was created, were I N a is the bifurcation parameter. We reproduce this diagram in Fig. 5. The solid and broken curves in Fig. 5 show stable and unstable Figure 5. One-parameter bifurcation diagram for Eqs. (4). I N a is the bifurcation parameter. In the graph the value of V 2 in the steady state was plotted for each value of I N a [5]. equilibrium points, respectively. For each value of I N a between Hopf bifurcation (HB1) and Hopf bifurcation (HB2), a periodic solution (stable or unstable) exists. Moreover, there can be many bifurcations and chaotic solutions are possible [5]. 3. Construction and analysis of proposed models. In this section we turn to constructing new models for specific pathologies we would like to describe. All the models we propose below consist of two basic models, namely the classic HH model and the YNI model with various types of couplings and feedback loops. Parameters for the part described by the HH model are defined in (3) and for the part described by the YNI model are defined in the subsection above. We use initial conditions in the following form:

3.1.
Model of the double-fire pathology. The first model we would like to propose reproduces the dynamics of action potential occurring in the phenomenon of double-fire. In this case, we have a coexistent conduction by two pathways: slow and fast. In this model we assume that the fast pathway is described by the HH model, and the slow pathway is treated as an action potential in the AV node, so we use the YNI model to describe the pathway. The conduction in the slow pathway is delayed, because the His-Purkinje system must step out with a refraction period after excitation by the fast pathway, such that there is a delay in this part. In this case, we use positive bidirectional coupling with delay in our model, because during the coexistence conduction by these pathways there are the interactions of the two waves. The faster pathway influences the slower pathway and vice versa. But we assume that the coupling coefficient is greater for the faster pathway. Therefore, the considered model reads: where HH(V 1 ) is the classical HH system, Y N I(V 2 ) is the YNI system, k 1 , k 2 are coupling coefficients, τ is the delay and V 1 , V 2 are membrane potentials. In the following we shall use such short description also for other models, but to avoid misunderstanding we write the full model of double-fire below. where , and only the delayed variables are written with their arguments In Fig. 6 we present an action potential which is typical for double-fire pathology. The characteristics of the action potential corresponds to atrial flutters and has a rather irregular rhythm. We obtain doubling the rhythm in relation to the original rhythm, also similar relationships were obtained in the electrophysiological recording. After making accurate invasive tests, e.g. electrophysiology study (EP), such a shape of the potential is diagnosed as double-fire tachycardia.
In the model without delay we obtain the regular rhythm which does not reflect the pathology we would like to describe. The delay added to our couplings can change the amplitude of oscillations and it will be non-physiological behaviour. In fact, in this model we need delay in one type of the pathways to obtain required patterns.

Solutions of this system lie at the intersection of the curves
Typically, studying stability of the system with delay we start from checking the dynamics for τ = 0; cf. [7]. Then we have ordinary differential equation system for which we can calculate eigenvalues. For given equilibrium point we get the following eigenvalues:  All the real parts of λ i , i = 1, 2, ..., 11, are negative which means that the equilibrium point is stable in the case without delay. What is interesting, the full system has the same stability property as the YNI model, and not inherit instability of the HH model. Now, we turn to the analysis for τ > 0. For the specific parameters it is possible to check stability using so-called Mikhailov criterion; cf. [7]. Let W : C → C denote a characteristic quasi-polynomial for the system of delay differential equations (DDEs) with discrete delays. In general W could be expressed as where 0 = h 0 < h 1 < . . . < h m are the delays, A k (z) = n k j=0 a jk z j , a jk ∈ R, are polynomials with real coefficients, n 0 ≥ 1 and n k < n 0 , for k = 1, . . . , m. If W (z) has no zeros on the imaginary axis, then W is stable (meaning that all roots lie in the open left half-plane of the complex plane) if and only if ∆ = n 0 π 2 , where ∆ denotes the change of the argument of W (iω) when ω increases from 0 to ∞. The curve drawn by the vector W (iω) = Re(W (iω)), Im(W (iω)) in the complex plane, when ω increases from 0 to ∞ is called a Mikhailov hodograph.
The characteristic quasi-polynomial for Eqs. (5) is complicated and has the form presented in Appendix. With increasing delay T the shape of Mikhailov hodograps change a little. However, the total change of the argument remains the same and equals 11π/2, independently of the value of the delay; we illustrate it for τ = 10, 20 and 50 in Fig. 8. In our case each hodograph starts on the real axis from negative values, that is the argument for ω = 0 is equal to π, as shown in Fig. 8 top left; here ω ∈ [0, 0.05]. In the consecutive graphs we see next turns around the origin, for ω ∈ [0, 0.192] we have first complete turn, for ω ∈ [0, 1] we have the increase of the argument of more than 3π, for ω ∈ [0, 3] we have more than 4π, and for ω ∈ [0, 22] we have the increase of more than 5π. Eventually, as shown in Fig. 8 bottom right the hodograph remains in the first quadrant yielding the increase of 11π/2. The full hodograph is the curve plotted for ω from 0 to ∞, and for large values of ω the changes near origin are not visible. This is the reason we plotted several parts of the hodograph to show these changes. From this analysis of the Mikhailov hodographs we observe that the number of turns gives stability of the steady state for all considered delays. Again this property is inherited from the YNI model.

3.2.
Model of AVNRT (fast pathway). Next we focus on such type of AVNRT in which only one fast pathway with depolarization of slow pathway occurs. The conduction is only out of the fast pathway, because descending depolarization of slow pathway prevents conduction. To describe that situation, we add only one unidirectional coupling to our model. This coupling does not change the rhythm (cf. Fig. 9). In this case the model could be written in the following form: 3.3. Model of typical and atypical AVNRT. Now we consider situation where the difference in the refractive state of both pathways leads to excitation of one of them causing the re-entry wave [23]. In Fig. 10 left result for a typical AVNRT Figure 9. Action potential of fast pathways described by Eqs. (6) for k = 1 (fast/slow) is presented. We describe this type of AVNRT by the following system: Fig. 10 right shows atypical AVNRT (slow/fast type of AVNRT), which we describe in the following way: In both of those types of AVNRT (fast/slow and slow/fast) we obtain regular fast rhythm, which is a typical behaviour for this kind of pathology (during this type of tachycardia the rhythm of the heart is about 35% more frequent than normal rhythm, which is in accordance with our results). Figure 10. Left: Action potential of fast /slow type of AVNRT for parameters: k = 2, τ = 0.75 (model with feedback and delay, Eqs. (7)). Right: Action potential of atypical AVNRT, that is Eqs. (8), for parameters: k 1 = 3, k 2 = 1, τ = 0.75.

3.4.
Model of AVNRT (slow-slow pathways). The last model is for the situation where there are two re-entry waves but both go through slow pathways (slow/slow type of AVNRT). We reflect it by the following system of equations: In this variant we also observe shortening (like in Fig. 10) of the period of oscillations. In this pathology, the presence of more than two conduction paths is possible which is associated with more re-entry waves. Increasing the number of feedbacks modeling re-entry waves of slow pathways causes a progressive shortening of the period of oscillation, while the rhythm remains regular. For numerical analysis of the discussed systems numerical models were created using Matlab Software. We used dde23 integration algorithm.

4.
Conclusions. The main aim of our paper was to propose the system of differential equations describing the dynamics of action potential that accompanies double-fire tachycardia and five forms of AVNRT. In this work, by using the proposed models we were able to reproduce the most important physiological properties of the discussed pathologies. In the literature, Hodgkin-Huxley type models with couplings and delays have been not considered till now. Also mathematical models of AVNRT and double-fire are not known. Our model could serve as a basis of the study of the influence of treatment. We can add an equation which is responsible for modeling pharmacological treatment of given pathologies. In such a way, we can help to determine optimal treatment. On the other hand, we should keep in mind that this is a phenomenological model, so the results are accurate as far as a simple model can describe the potential found in one of the most complex oscillators appearing in real biological phenomena. Therefore, further validation of this model using medical data (from invasive laboratory EP) is necessary. We hope that such validation will be possible in the future.
Appendix. The exact form of the quasi-polynomial for Eqs. (7)