A Modified Heart Dipole Model for the Generation of Pathological ECG Signals

In this paper, we introduce a new dynamic model of simulation of electrocardiograms (ECGs) affected by pathologies starting from the well-known McSharry dynamic model for the ECGs without cardiac disorders. In particular, the McSharry model has been generalized (by a linear transformation and a rotation) for simulating ECGs affected by heart diseases verifying, from one hand, the existence and uniqueness of the solution and, on the other hand, if it admits instabilities. The results, obtained numerically by a procedure based on a Four Stage Lobatto IIIa formula, show the good performances of the proposed model in producing ECGs with or without heart diseases very similar to those achieved directly on the patients. Moreover, verified that the ECGs signals are affected by uncertainty and/or imprecision through the computation of the linear index and the fuzzy entropy index (whose values obtained are close to unity), these similarities among ECGs signals (with or without heart diseases) have been quantified by a well-established fuzzy approach based on fuzzy similarity computations highlighting that the proposed model to simulate ECGs affected by pathologies can be considered as a solid starting point for the development of synthetic pathological ECGs signals.


Introduction
Ever since Einthoven [1,2] discovered the possibility of obtaining an electrocardiogram (ECG) from cardiac activity, cardiologists have been able to take advantage of this effective and efficient tool for diagnosing cardiac disorders [1,3,4]. Specifically, an ECG is a particular recording of the electrical activity of the heart muscle about the potentials of the body surface [5][6][7][8]. The use of an ECG is nowadays widespread because of its non-invasiveness and easiness of performing as a test. Furthermore, the presence of slight cardiac anomalies produce distortions of the ECG signal with the consequent possibility of accurate diagnosis [9][10][11]. During the years, the scientific community has been actively engaged in the development, implementation and validation of models and algorithms for the analysis and interpretation of ECG signals [12]. In particular, many researchers have been involved in the study of wave morphology and wave complexes, forming a complete cardiac cycle [13]. In addition, to explain the cardiac conductivity, a lot of models of the ECG wave have been proposed. Among them stands [14] in which the ECG model was constructed starting from Gaussian formulations. Moreover, special indexes, such as the bicoherence-derived index, have been proposed for assessing nonlinear processes in cardiac regulation [15]. At the same time, other researchers have focused their attention on the analysis and variations of the models' outputs on a suitable number of heartbeats [12,13,16]. The main reason that guides researchers to implement and validate physical-mathematical cardiac models concerns the impossibility of distorting cardiac activity directly on the patient to study its effects on the ECG signal [17]. Furthermore, the creation of cardiac models has the undoubted advantage of characterizing the entire cardiac cycle, with the possibility of simulating highly dangerous cardiac disorders [17]. With the aim of generating ECGs, many papers have been published, including [18], in which, to generate ECGs, discretized reaction-diffusion systems to produce a set of three nonlinear oscillators simulating the main pacemakers in the heart were exploited. Moreover, the availability of synthetic ECG signals with or without distortions allows the exploitation of consolidated signal processing algorithms [19][20][21]. Finally, the development of models and algorithms for the study of ECG also turns out to be a valuable aid in the evaluation of the fetal ECG [22,23]. Researchers are also interested in the synthesis of ECG in which distortions are present to simulate as many heart diseases as possible. This allows for laboratory analysis of signals without the need to sample signals from patients who are often difficult to transport due to severe heart disease. In this context, McSharry et al. in 2003, propose a dynamic model for the reconstruction of ECG signals combining Gaussian functions with different parameters [24]. In particular, this model was derived from coupling a dynamic first-order model with a differential equation in which the respiratory artifact was taken into account. This model was solved numerically by Runge-Kutta procedures obtaining a typical ECG trajectory without diseases [24]. It is worth noting that the Runge-Kutta methods, in addition to being a valid numerical tool for solving linear differential problems, are able to numerically solve problems characterized by strong nonlinearities when the nonlinear functions are smooth, such as the McSharry model [25]. Although the McSharry model is interesting from a theoretical point of view, it does not supply ECGs with heart diseases. From which, it is imperative to make new models for ECGs where the heart diseases can be easily implemented exploiting suitable numerical procedures. Although many papers have addressed the problem of generalizing the McSharry model to generate ECGs with pathologies (see, for example, [26]), the main objective of this study concerns the generalization of the model [24] showing how it can be extended to simulate the effects on ECG affected by heart disease exploiting a linear transformation and a rotation of two dynamic parameters of the model in [24]. In particular, atrial and ventricular fibrillation and premature ventricle contraction were chosen in this paper as they represent frequent pathologies in patients who are often highly at risk, having no obvious symptoms. From our analysis, six different formulations of the model [24] have been obtained. These have been solved numerically by a procedure based on the Four Stage Lobatto IIIa formula [25] obtaining six different typical trajectories related to normal ECG (without diseases). Although the implicit fourth-order Runge-Kutta methods are suitable for solving nonlinear differential systems, they can be considered equivalent to the four-stage Lobatto IIIa formula both in terms of convergence and stability. In fact, all the above numerical procedures admit the implicit fourth order version and are stable [25]. However, the four-stage Lobatto IIIa formula, unlike the implicit fourth-order Runge-Kutta method, is a numerical collocation method which, when applied to the problem under study, gave slightly better qualitative results than those obtained in [24] who exploited the implicit fourth-order Runge-Kutta procedure. This is evident from the qualitative analysis of the ECGs obtained in our work compared to the ECGs obtained in [24], which requires a, albeit mild, filtering operation. Moreover, in this paper, we have verified whether the model [24] admits existence and uniqueness of the solution (in order to avoid any numerical ghost solutions). Moreover, unlike the model [24], in our work, the result provided by the autonomous system consisting of the first two equations in [24] represents the input for the third equation in [24] where the ECG signal is performed. So, before studying the complete system [24], we need to study the stability of the system constituted by the first two equations in [24]. This is because we have to verify that there were no instability points (which generate unbounded perturbations) on points admitted by the third equation in [24]. In this paper, we have shown that the only position that generates instability, fortunately, is not in the domain of the third equation in [24]. Thus, a matrix procedure here proposed, consisting of a linear transformation and a rotation depending on two characteristic parameters, distorted the ECG signal for implementing the heart diseases. Also, in this case, six different formulations of this new generalized model, named modified heart dipole model (MHDM), have been obtained, reconstructing their typical trajectories, and verifying that, also in this case, the same unstable equilibrium position is obtained, which is discarded since it does not belong to the set of points allowed by the third equation. The results obtained can be considered encouraging because the ECG reconstructed in the absence and the presence of three specific pathologies are entirely comparable to the ECG obtained directly on patients. These comparisons, taking into account that ECGs can be often affected by uncertainty and/or imprecision (due, for example, to artifacts), have also been evaluated by a procedure based on fuzzy similarity measuring how similar is, in a fuzzy sense, an ECG (with or without heart diseases) obtained by the proposed model with an ECG (with or without heart diseases) sampled directly on a patient. Based on these results, our model should allow building a complete database of ECG without the need to carry out a measurement campaign directly on patients (some of which, due to the severity of the pathology they suffer from, are not transportable).
The paper is structured as follows. Section 2 offers a quick description of the 12-lead ECG system while Section 3 describes both normal conditions and pathological events in an ECG. After having described the HDV & the 12-Leads ECG system in Section 4, Section 5 details the McSharry dynamical model. In Section 6 a discussion on the stability of its unique equilibrium position is given. The existence and the uniqueness of the solution for the McSharry model are proved in Section 7, in which the most important peculiarities and weaknesses of the trajectories are also discussed. In Section 8 the six different formulations of the McSharry model are solved numerically. Section 9 introduces our modified heart dipole model (MHDM). This model extends the model described in [24] for the synthesis of ECG signals affected by cardiac disorders. Moreover, we show how an MHDM allows the 3D representation of an HDV, i.e., a vector cardiogram (VCG) through its reconstruction from the 12-leads system. Once the ECGs have been achieved (with or without heart diseases), they have been compared, in the fuzzy sense, in Section 10 with ECGs sampled directed on the patients. Numerical results, conclusions and some preliminary perspectives are provided in Sections 11 and 12, respectively. Finally, mathematical details and proofs of the main results are reported at the end of the paper in appendices.

Normal 12-Lead ECG
An ECG originates following the atrial and ventricular contractions where the concentrations of the sodium and potassium ions are unbalanced so that the heart muscle, being electrically charged, is equivalent to an electric dipole vector with amplitude and variable direction during the depolarization and re-polarization of cells [8,17,[27][28][29]. The electrical behavior of a single pair of ions produces an electrical dipole. All the electric dipoles determine a single resulting vector called Heart Dipole Vector (HDV), d(t), where t is the time variable. Then, an ECG(t) derives from the scalar product between d(t) and the direction of the axis of the recording electrode, v, generating the so-called Heart Dipole Model (HDM). On these premises, an ECG consists of 12 different leads describing the same vector d(t) at the same instant t but from different points of view. d(t) reconstruction requires the analysis of its projection in several different leads, to identify its spatial arrangement within the heart. The electromotive forces are best recorded on the horizontal plane (left to right, parallel to the horizon) and on the vertical (frontal) plane, perpendicular to the shoulder line, from head to toe. Running an ECG includes the following leads ( Figure 1a) [30,31]: (1) six vertical leads (on the frontal plane): (a) bipolar leads of the limbs: I, II, III (Standard leads); (b) uni-polar leads of the augmented limbs: aVR, aVL, aVF (Goldberger leads); (2) six horizontal branches (on the horizontal plane). Uni-polar or chest wall precordial leads: V 1 to V 6 or V 9 (Wilson leads). Uni-polar leads are recorded by coupling a single scanning electrode with a central terminal, while bipolar leads record the potential difference between two electrodes. Leads I, II, II are bipolar; leads aVR, AVL, aVF, and Wilson leads are uni-polar [30][31][32][33].

Structure of an ECG: Normal Conditions and Pathological Events
ECG is the graphic reproduction of the electrical activity of the heart during its operation, recorded at the level of the body surface in which low intensity electric fields are present and can be recorded, which in the individual at rest are mainly due to periodic depolarization and re-polarization of the heart [30,31]. Thanks to the conversion of electrical energy into mechanical energy, electrical variations produce the movement of a writing system. The electrical energy is adequately amplified so as to be able to transcribe sufficiently large excursions to allow the recording of a readable signal. The deflections are imprinted on paper, which moves at a constant speed in contact with the system, which reports on the paper the waves recorded as a function of time. Simultaneously with the vertical oscillation of the lines produced by the variations in potential, the paper flows to the left. This synchronization makes it possible to bring the vertical movement back to a horizontal plane, recording the oscillations in relation to their duration over time. The principle on which the measurement of the electrical activity of the heart is based is physiological: the onset of impulses in the myocardium leads to the generation of potential differences that vary in space and time and which can be acquired through electrodes placed on the body surface. By convention, the ECG trace is reported on graph paper with the time on the abscissa (one second every 25 mm) and the amplitude on the ordinate (one mV every 10 mm). ECG consists of a series of waves [30,31]. The P wave corresponds to the depolarization of the atria and originates from the sinoatrial node. When the electrical impulse leaves the sinus node it produces the depolarization of the neighboring myofibrils which contract, and then continues to propagate in the radial sinus. The PR interval concerns the wave front which, passing through the atria, passes into the atrio-ventricular node (AV) within which the activated cells are few and the dipole generated is too weak to be registered [30][31][32]34,35]. As soon as the depolarization wave reaches the AV node, the electrical conduction slows down until the ventricular conduction system is reached. The QRS complex concerns a set of waves that follow one another, corresponding to the depolarization of the ventricles. The Q wave corresponds to the depolarization of the inter-ventricular septum, the vector produced goes from left to right. The R wave is a very high peak and corresponds to the depolarization of the apical part of the ventricles. The S wave is also small in size and corresponds to the depolarization of the basal and posterior regions of the left ventricle. The ST segment represents the period in which the ventricular cells are all depolarized and therefore no electrical movements are detectable until the beginning of re-polarization. As a result, the ST segment is usually isoelectric, i.e., placed on the base line of the trace from which it can move, upwards or downwards, by no more than 1 mm (equal to 0.1 mV). The T wave is the first wave of the re-polarization of the ventricles. It is not always identifiable because it can be of very small amplitude. The QT interval represents electrical systole, that is the time in which ventricular depolarization and re-polarization occurs. The U wave is a wave that is not always appreciable because it is often of minimal size. It is due to the re-polarization of the papillary muscles, which can be evidenced in the case of myocardial hypertrophy or altered dimensions of the ventricular cavities [30][31][32]34,35]. The most salient characters of the human ECG are shown in Table 1 ( [36]). Excitation initially spreads along the right atrium via the three inter-nodal tracts and reaches the AV node. The resulting atrial vector is directed from the upper right to the lower left, approximately in the same direction as the standard II lead. The PR interval runs from the beginning of the P wave to the beginning of the QRS complex (which represents the depolarization of the ventricles) [30][31][32]34] and indicates the time required for a pulse to travel from the sinus node to the ventricles and includes conduction through the atria and the AV node. This is the last element as certain drugs or pathologies can delay conduction through the atrioventricular node. A shortening of the P-R interval may indicate that conduction occurs through an accessory beam that "skips" the AV node. This phenomenon can be associated with atrial arrhythmia. The normal activation sequence occurs when the impulse runs from the AV node to the branches and to the Purkinje system [30][31][32]34]. The left branch is depolarized a little before the right branch because the initial depolarization vector runs from the left side to the right side of the septum. Initial septal activation produces the initial Q wave in left-directed leads I and VL, as well as in left thoracic leads V 5 and V 6 . Subsequently, various areas of the left/right ventricle are activated simultaneously. When both ventricles are depolarized, there is no vector, and the R wave returns to the baseline. The point where the vector changes direction along the axis of the derivation of the electrocardiography leads is called deflection. The configuration of QRS in the bipolar and uni-polar leads of the limbs (frontal plane) depends on the axis of the heart. In precordial leads (horizontal plane) the amplitude of the R wave increases from V 1 to V 5 while at the same time, the amplitude of the S wave is reduced. The whole QRS complex may be a little smaller between V 4 and V 6 , especially in patients with opiates or with pulmonary emphysema [30][31][32][33][34]. The T wave represents the depolarization of the ventricles, and it is larger than the QRS complex, with slightly asymmetrical branches and a rounded tip, mainly because re-polarization does not occur with the same course as depolarization. Generally, the T wave is positive in those leads where the QRS complex is predominantly positive. Figure 1b depicts ECG signals, both normal and pathological ones [30][31][32]34].

HDV and 12-Leads ECG System: A Brief Overview
Let us consider a region Ω ⊂ R 3 endowed with a system of Cartesian axes Oxyz in which the origin O represents the central area of the heart muscle (see Figure 1a). The unit orthogonal vectors are indicated byî x ,î y andî z . Usually, the cardiac activity, referred to three orthogonal planes (x − z plane (sagittal plane); y − z plane (frontal plane) and x − y plane (transverse plane)) can be simulated by an ad-hoc time-varying rotating vector [5] d(t) = (d x (t), d y (t), d z (t)) in which the Cartesian representation in the Oxyz system is: As t → +∞, d(t) varies cyclically in both amplitude and direction, keeping the origin fixed, depicting a 3D spatial loop in Oxyz (the well-known VCG). As an example, Figure 2a displays a typical VCG loop during a normal cardiac cycle, while Figure 2b depicts the VCG loop with a nodal rhythm disease. If v = v xîx + v yîy + v zîz is the direction of the axis of the recording electrode, ECG(t) can be written as follows: In accordance with (1) and (2), it is possible to perform a description of the HDV 3D using three ECG signals obtaining by the scalar product of d with each of coordinate axes v 1 , v 2 , v 3 .

Remark 1.
The single dipole model is not a perfect modeling of cardiac activity, so that the complete 12-lead ECG system is useful (albeit redundant) for diagnosis.
The 12-lead ECG system records the following 12 linearly independent signals: the three Einthoven points I, I I, I I I; Wilson three leads aVR, aVL, aV H and the six precordial leads V1 − V6 [1,3]. However, to get a more detailed description of ECG when noise, uncertainty and/or imprecision take place, we need to generalize Equation (2) considering a 3 × 3 matrix V, corresponding to the projection of the HDV onto the direction of the recording electrode, and a 12 × 3 matrix, H, that includes the body volume conductor model [1,13,17]. Then, the vector ECG(t) 12−lead can be achieved as follows: in which N(t), ∀t, represents a vector in which the components are the noise in each of the 12 ECG channels. Equation (3) represents a mathematical model of the 12-lead ECG system even when a pathology is in place. In the next section, we present a particular dynamic model in which the parameters are associated with certain cardiac pathologies.

Remark 2.
We observe that (3) is a typical procedure frequently used in dynamic systems (and in particular in signal processing theory) for the reconstruction of a signal in the time domain [25]. This procedure, as seen above, uses the formulation of matrices, H, V and N as holding information contents typical of the system under study (including the components due to noise). Such matrices are usually very difficult to obtain. Therefore, this procedure is not practically exploitable to generate ECG(t) in which cardiac pathologies could also be considered. Thus, for our purposes, we prefer starting from the well-known McSharry et al. dynamical model [24] to obtain a general ECG(t) in which cardiac pathologies could also be taken into account.

ECG of Patients with Heart Disease: The McSharry Single Channel Modeling Approach
In [24], the PQRST complex was obtained by a set of Gaussian functions with different parameters. If s(t) is the ECG signal, then: Table 2). Finally, θ represents an angular dynamic parameter that influences the deflections.
In [24], (4) was coupled with the following dynamical system: in which ω is a parameter concerning the heart rate, obtaining the following dynamical model: where s 0 (t) = A sin(2π f 0 t), with A = 0.15 mV and f 0 (s −1 ) respiratory frequency, is the artifact due to breathing. ω is the angular velocity of the trajectory [24]. Moreover, θ can be evaluated as follows [37]: so that −π ≤ atan2(γ(t), χ(t)) ≤ π. Remark 3. |atan2(γ(t), χ(t)) − θ i | 2π means that from the division between atan2(γ, χ) − θ i and 2π only the rest must be considered. atan2(γ(t), χ(t)) indicates the angle between the positive semi-axis of X of a Cartesian plane and a point of coordinates (γ(t), χ(t)) lying on it. The angle is positive if counterclockwise (half plane of positive ordinates) and negative if clockwise (half plane of negative ordinates). The function is an extension of the atan function since, unlike this, it distinguishes between diametrically opposite angles, not only considering the relationship between the arguments but also their sign.
Finally, the McSharry system (6) can be written in the following form: Remark 4. In [24], (5) was coupled with (4) without evaluating the presence of possibly equilibrium configurations for (5). In our work, unlike [24], we consider (5) as a system in which the solutions are the inputs for the Equation (4) where s(t) represents the ECG(t). Therefore, before studying the complete system (8), we study the stability of (5). This is because we have to verify that there are no traces of instability points (generating unbounded perturbations) on points admitted by (4). The following section proves how (5) admits an unstable equilibrium position, which does not belong to the set of points eligible for (4).

On the Stability of the Unique Equilibrium Position
Proposition 2. (9) is an unstable equilibrium position for (5).
Proof. (5) is nonlinear, so it was linearized (by Taylor series development). Afterward, the only position of equilibrium of the linearized system corresponding to the origin of the Cartesian axes was obtained. Applying the Lyapunov-Poincaré Lemma [38], it has been shown that the origin is also the only configuration of equilibrium for the nonlinearized system. Applying the first Lyapunov criterion, it has been shown that this equilibrium configuration is unstable. Therefore, these trajectories are not contained in closed curves on the phase plane. Hence the confirmation that the origin is an unstable equilibrium point. For details, see Appendix B.

Remark 5.
We note that the only equilibrium configuration admitted by (5), (i.e., (9)) is located in the origin of the Cartesian axes. Furthermore, (9) does not belong to the domain of (7) because, in (9), (7) is not defined. Therefore, there is no instability points for (4).
In [24], (8) was solved using numerical techniques based on Runge-Kutta procedures. However, the solutions obtained could, in some cases, represent ghost solutions. In other words, the numerical procedure still provides a result that may not satisfy any analytical conditions of existence (and possibly uniqueness) of the solution. Thus, we need to verify that (8) still admits a solution and that it is unique. In this paper, we check if this happens.

The McSharry Dynamical Model: Existence and Uniqueness of the Solution
The uniqueness of the solution, from a clinical point of view, is very important as it ensures the one-to-one correspondence between a given pathology with a particular ECG. We introduce the following definition.
Lipschitzian if for it the ratio between the norm of the increases in f and that of the increases in x is maintained above limited. That is, there exists a number L (Lipschitz constant) such that {x ∈ X, The definition remains valid for a function f (x, t) with x ∈ X and with t set in T (real interval), and if then there exists a constant L for which it results that The following Lemma holds.

Lemma 1.
Let us consider the following Cauchy problem: with T real interval, and uniformly Lipschitzian for x ∈ R n with respect to t ∈ T. Then, however chosen (t 0 , x 0 ) ∈ X × T, (10) admits unique solution, x(t).
, as defined in (12), is continuous as (t, χ(t)) varies in R 3 × T, with T real interval.
Proof. The proof is based on the verification of the continuity of each single addendum. For details, see Appendix C.  Proof. It immediately follows from Propositions 3 and 4 and from Lemma 1.

Remark 7.
The resolution of (8) is usually numeric and techniques based on Runge-Kutta procedures are particularly suitable for this goal [24]. However, different values of atan2(γ(t), χ(t)) (see (7)) provide different formulations with consequent different solution. Then, in the following Section the McSharry models obtained for different values of atan2(γ(t), χ(t)) are presented and solved.

The McSharry Dynamical Model and the Study of the Trajectory: Peculiarities and Weaknesses
Here, taking into account (7), it follows that atan2(γ(t), χ(t)) = atan γ(t) χ(t) , and considering the values in Table 2, (8) becomes: The numerical solution of (14), achieved by four-stage Lobatto IIa formula because it is a particularly suitable technique for solving nonlinear differential problems (for detail, see Appendix G), is depicted in Figure 3a, which shows the regular trend of the ECG(t): the waves and relative deflections are clearly visible and determined by the periodic terms contained inṡ(t).
χ(t) + π In this case, considering (7), (8) becomes: Also in this case the periodic terms contained inṡ(t) produce the deflections in the corresponding ECG. However, the presence of noise components (mainly due to the particular value of operative parameters inṡ(t)) would require the use of suitable filters (see Figure 3b).

Case (c)
: (7), (8) assumes the following form: Like the two previous ones, it produces a regular ECG signal in which waves and deflections are clearly visible (see Figure 4a). However, due to the operating parameters used, the signal amplitude appears reduced (albeit slightly).

Case (d)
: In this case, system (8) becomes: However, (17) is simple to solve analytically by obtaining the following solution (for analytical details, see Appendix E): where K 2 and C are constant. We immediately notice that, as t → +∞, form system (18), γ(t) → 1 and s(t) → 0 so that the expected periodicity vanishes. Thus, this case must be discarded.
In this case, system (8) becomes: Here the solution is similar to the previous case therefore the respective system, as the previous case, results decoupled. Furthermore, the periodicity of the solution, excluding the respiratory artifact, has completely disappeared. Thus, also this case must be discarded.

Remark 9.
As atan2(γ(t), χ(t)) varies according to (7), we have six different McSharry systems. However, only three of them provided reliable results from a biomedical point of view. In particular, only (a), (b) and (c) cases in Section 8 provide reliable ECG(t). Excluding the case (f) where atan2(γ(t), χ(t)) is not defined, cases (d) and (e) provide physically unreliable results because there atan2(γ(t), χ(t)) = ± π 2 . Finally, due to its structure, the McSharry model is unable to simulate ECG(t) affected by heart disease. Therefore it is necessary to modify this model to take into account the possible alterations due to cardiac pathologies. Section 9 presents a new approach to modifying McSharry model in order to take into account heart diseases.

The Modified Heart Dipole Model (MHDM)
McSharry and co-authors were not able to characterize the distortion of the PQRST complex [24]. Then, we propose a generalization of the χ − γ limit orbit by modifying (5) so that this limit orbit can be dynamically adapted. To simulate pathological ECG signals, it is essential to choose a set of parameters that act as a controller of the limit orbit generalization. So, we here propose to generalize the limit cycle by an elliptic orbit introducing a linear transformation and a rotation: linear trans f ormation (20) so that θ = atan 2(γ(t),χ(t)), withγ(t) andχ(t) computed by (20). Moreover, we note that for k 1 = k 2 = 1 and φ = 0, (20) does not modify (5). We also note that, as the parameters k 1 , k 2 and φ vary, the PQRST complex undergoes a distortion. Therefore, for certain values of k 1 , k 2 and φ it is possible to simulate specific heart diseases. Table 3 shows some heart diseases associated with specific values of k 1 , k 2 and φ. Table 3. Some synthetic diseases by the modified heart dipole model (MHDM).

Model Parameters PQRST Complex Alteration Correlated Disease
isoelectric interval TP atrial flutter abnormal QRS premature P wave not associated ventricle contraction Remark 10. Also system (20) admits as its unique equilibrium position the origin O (see Proposition 5) verifying that the choice to use a linear transformation and a rotation does not affect the stability of the starting model. In other words, linear transformation and rotation do not generate instabilities different from the one already detected for (5) and therefore can be eliminated since the associated point is outside the domain of (7). Moreover, (20) has a twofold advantage: the computational complexity of the model is almost unvaried and, in addition, it is possible to characterize the distortion of the PQRST complex setting k 1 , k 2 and φ.
Then, the following result holds. Proof. For details, see Appendix F.

MHDM Modelization
Staring from model (8) and considering the transformation (20) we can easily write: in which the parameters present have the same meaning as the cases studied above. System (21) admits different formulations as |atan2(−k 1 χ(t) sin φ + k 2 γ(t) cos φ, k 1 χ(t) cos φ + k 2 γ(t) sin φ) − θ i | varies according the (7). Thus, as obtained before in Section 8, we here achieve six formulations of system (21) depending on (7): 1. Case (a): 3. Case (c): Case (e): from which, by calculations completely similar to those made in Section 8, cases (d), (e) and (f) can be discarded. For each case we obtain a new formulation of the system (21) completely analogous to the formulations (14)-(16), respectively. For cases (a), (b) and (c) values of parameters k 1 and k 2 were considered such as to simulate normal conditions of heart rhythm and pathologies (atrial fibrillation, atrial flutter and premature contraction of the ventricle because they represent frequent pathologies in patients who are often highly at risk, having no obvious symptoms) as described in Table 3. By means of the four-stage Lobatto IIIa formula each system was numerically solved (for details, see Appendix G).
Concerning the Case (a), Figures 5, 6a, 7a and 8a display the ECG obtained by varying k 1 and k 2 as per Table 3. ECG(t) strongly adhering to real cases are highlighted, for which the methodology proposed to simulate ECG pathologists provides qualitatively reliable results. Similarly, concerning Case (b) and (c) of the proposed methodology we also obtain ECG conforming to real cases (see Figures 9a-d and 10a-d).  Figures 4b, 6b, 7b and 8b respectively. Also in these cases, for the quantitative comparison, see Section 10.

Remark 12.
It is worth noting that the ECGs obtained in this work using the four-stage Lobatto IIIa formula are qualitatively better than those obtained in [24], in which an implicit fourth order Runge-Kutta procedure had been exploited. In fact, while in [24] the ECGs obtained needed filtering, in our work the ECGs obtained did not require this operation.

Fuzzy Similarities for Comparing ECGs
A good cardiologist is certainly able, looking at an ECG, to detect if a pathology is in progress or if the signal is devoid of pathological components. Cardio-logical triage in the emergency rooms of the most common hospitals is based on this principle, at least as a first approximation. In other words, in the emergency hospital clinic it is not customary to use comparison procedures between a freshly sampled ECG with known ECGs (with or without pathologies) in order to evaluate how close (or, dually, how far away it is) an ECG from a class of known ECGs. However, if we did not have the knowledge of an expert, surely we would need a tool that can evaluate how similar an ECG is to another known ECG. Thus, taking into account the fact that the ECGs may be affected by uncertainties and/or inaccuracies, in this paper, to quantitatively evaluate how much one ECG is similar to another one, we rely on the computation of Fuzzy Similarity (FS) between signals as this technique, in a fuzzy sense, quantifies a sort of "proximity" between signals which, as in this case, are affected by uncertainty and/or imprecision [39,40].

Fuzzy Similarity Concept: Aspects for Comparing ECGs
As known, a fuzzy set A can be defined by a membership function µ A (·) : A → [0, 1] or considering A as a point in a certain metric space, the distance between two fuzzy sets A and B, indicated as d (A, B), indicates the distance between two points in a certain n-dimensional space. Considering A as an ECG Considering A as a freshly sampled ECG and B as an ECG with known clinical characteristics (i.e., it is known whether it relates to a healthy patient or a known cardiac disease patient), d(A, B) represents how much the ECG A is next to B. Starting from the assumption that the same pathology produces very similar ECGs, we group the ECGs into classes of pathologies converging in a single fuzzy set that represents the class of ECGs characterized by that pathology. Furthermore, if B is the fuzzy set representing the class of ECGs with a particular pathology, d(A, B) provides the measure, in a fuzzy sense, of how much the ECG A is similar to the class B. Furthermore, the presence of a pathology in an ECG is assessed by the reduction of its FS with the class of ECGs related to healthy patients.

Definition 2.
From the formal point of view, on the set of all fuzzy sets, F(X) with X universe of discourse, we define the following normalized function [39][40][41]: such that, for A, B, C ∈ F(X) it follows that Remark 13. It is worth noting that the more a fuzzy set A is similar to a fuzzy set B, the more S(A, B) approaches 1. If A is totally similar to B, then S(A, B) = 1. In this case, A and B, in the fuzzy sense, are to be considered equal.
If NP is the healthy ECGs class and KP is the generic ECGs disease class, then the FS of a newly sampled ECG (with unknown pathology) versus NP or a KP, FS(NP/KP), can be calculated by means of a function satisfying the above properties. In this work, we verified the fuzziness of the classes of ECGs using specific entropy indices, such as the linear index (LI) [39][40][41] and the fuzzy entropy index (FE) and considering that high values of LI and FE indicate high fuzziness in the ECGs belonging to the class, in this paper, four different formulations of FS have been used [39][40][41]: in which s ∈ {1, 2, ∞}; n is the number of samples (which is the same for all ECGs); µ NP/KP (·) and µ Y (·) represent the membership values related to NP, KP and Y, respectively. It is worth noting that, if Y is an ECG without pathologies, then µ Y (x i ) = µ NP (x i ) = 1. Otherwise, if it is related to a certain pathology KP, then µ Y (x i ) = µ KP (x i ) < 1.

The Exploited Database
In this paper, MIT-BIH Arrhythmia Database by PhysioNet (public domain database of ECGs, https://physionet.org/content/mitdb/1.0.0/), studied by the BIH Arrhythmia Laboratory at Beth Israel Deaconess Medical Center, was exploited [42,43]. From the entire Database, the directory labeled by MLII was used. From it the following sub-directories were exploited: 1.
NSR sub-directory, which contains ECGs characterized by normal sinus rhythm; 2.
AFIB sub-directory in which ECGs with atrial fibrillation are contained; 4.
PVC sub-directory containing ECGs with premature ventricle contraction.
Particularly, each recording contains ECGs recorded for 20 s, digitized at 500 Hz with 12-bit resolution over a nominal ±10 mV range. The records were obtained from volunteers both healthy and suffering from heart diseases. In addition, the number of records for each person varies from 2 (collected during one day) to 20 (collected periodically over 6 months). The raw ECGs are rather noisy, so that they have been suitably filtered. Four classes of ECGs have been set up: the first class is made up of 140 ECGs of subjects without heart disease; the second class, also made up of 140 ECGs, concerns patients with atrial fibrillation; the third class, as for the previous ones formed by 140 ECGs, contains signals from patients suffering from atrial flutter; finally, the last class contains 140 ECGs of patients with premature ventricle contraction.

Some Statistical Considerations
It is worth noting that the ECGs concerning a particular pathology are qualitatively similar to each other. This has an impact on the fact that the histograms relating to the similar ECGs signals show strong similarities. In parallel, also ECGs signals without pathologies are qualitatively similar to each other. Also in this case, the relative histograms are similar. Usually, ECGs histograms highlight Gaussian trends as depicted in Figure 11a-d where histograms for typical ECGs characterized by normal sinus rhythm, atrial fibrillation, atrial flutter and premature ventricle contraction, respectively, are shown. Therefore, it makes sense to collect ECGs characterized by the same pathology in one single class while all the ECGs without pathologies can be considered as belonging to a single class. In this way, we obtain four classes of ECGs signals: the first class is related to the 140 signals without pathologies (normal sinus rhythm); the second class relates to the 140 signals affected by atrial fibrillation; the third class is related to the 140 ECGs that show atrial flutter while the last class, the fourth, collects the 140 ECGs, which show premature ventricle contraction. Thus, it is sufficient to consider mean m ij and standard deviation σ ij to characterize the ith ECG belonging to the jth. Then, different ECGs belonging to the same class highlights light variations of both mean and standard deviation. However, we need to consider just only value for both mean and standard deviation for calculating the Gaussian fuzzy membership function related to the jth class. Therefore, if KP j , j ∈ {1, 2, 3}, is the jth class of ECGs, it follows that that is the mean value of the means and that is the maximum value of all the standard deviation. Similarly, the same procedure has been applied to the class of ECGs without pathologies achieving m NP and σ NP . To build the membership function related to each class KP j and NP we make the relative Gaussian membership functions: and In this framework, it makes sense to consider the mean of the means of the ECGs belonging to each class as the mean value representing that class while, for the standard deviation, we consider the maximum value of all the single standard deviations in order to cover all the ECGs belonging to that class.

Results of Interest
To underline the fuzziness content of each class, both fuzziness indexes ( (26) and (27)) for each class were computed: the high values obtained (that is, very close to the unit value) confirmed the high level of fuzziness contained in each class (see Figure 12). In particular, the class of the ECGs without pathologies obtained LI = 0.97 and FE = 0.98, respectively, while the class of the ECGs with atrial fibrillation obtained LI = 0.94 and FE = 0.95, respectively. Moreover, concerning the class of the ECGs with atrial flutter achieved LI = 0.98 and FE = 0.92, respectively. Finally, regarding the class of the ECGs with premature ventricle contraction, LI = 0.97 and FE = 0.97 were obtained, respectively. This justifies the use of the fuzzy similarity based approach for the analysis of ECGs. Thus, a unique fuzzy membership function for each class has been achieved according to Section 10.3, obtaining four membership functions (the first one related to ECGs without heart disease, and the other three related to ECGs affected by atrial fibrillation, atrial flutter and premature ventricle contraction, respectively). Therefore, for each ECG numerically obtained by the McSharry dynamical model (Figures 3a,b and 4a) the four fuzzy similarities, FS 1 , FS 2 , FS 3 and FS 4 , with respect to the four classes of ECGs have been computed and the results have been reported in Table 4 respectively, highlighting as Figures 3a,b and 4a belongs to the class of normal ECGs (i.e., without heart diseases) due to the high fuzzy similarities obtained. This is supported by the fact that the fuzzy similarities values of Figures 3a,b and 4a obtained with the other classes are far from reaching the unit value as shown in Table 4 (in bold). Similarly, the fuzzy similarities of the ECGs obtained with the proposed approach were computed with the representative ECGs of each class. Also in these cases the results confirmed the qualitative analysis. In particular, concerning the Case (a), Figure 5 belongs to the class of the ECGs of patients without cardiac pathologies because the FS 1 , FS 2 , FS 3 and FS 4 values obtained with respect this class are very close to the unit value (in bold), while the values of fiìuzzy similarities obtained with respect to the other classes with pathologies are extremely low (see Table 4). Moreover, Table 4 highlights that ECG depicted in Figure 6a belongs to the class of ECGs affected by atrial fibrillation because FS 1 , FS 2 , FS 3 and FS 4 values achieved are very close to unity (in bold). Again, Table 4 proves that the ECG shown in Figure 7a belongs to the class of ECGs affected by atrial flutter (FS 1 , FS 2 , FS 3 and FS 4 are very close to unity) (in bold). Finally, Table 4 puts in evidence that the ECG depicted in Figure 8a belongs to the class of ECGs affected by premature ventricle contraction (FS 1 , FS 2 , FS 3 and FS 4 are very close to unity) (in bold). Also for Cases (b) and (c) elaborated with the approach proposed in this paper, the quantitative analysis using fuzzy similarities confirmed the qualitative considerations as shown in Tables 5 and 6, respectively.

Some Considerations on the Stability of the Equilibrium Position
In this paper, we have analyzed in-depth the McSharry model [24]. The dynamic system (5) (a part of the McSharry model (6)) was studied, highlighting the behavior of the single equilibrium position (9), which was found to be unstable (Propositions 1 and 2). This unstable equilibrium position, located at the origin of the coordinate reference system, corresponds to the cusp of the T wave (thought as a strong distortion) of the QRST complex. Furthermore, this unstable equilibrium position corresponds to a point not admitted by the generating equation of the ECG(t), so that this undesirable instability is, in fact, eliminated.

On the Uniqueness of the Solution for the McSharry Model and Numerical Approach
Firstly, we observed that the McSharry model depends on the trigonometric function atan2 (see (7)), thus that the six specific McSharry systems have been derived, according to the value assumed by atan2 (systems (14)- (17), (19) and case (f) in Section 8.6). Furthermore, once that the McSharry model was written in the form (13) as a Cauchy problem and verified that h(t, χ(t)) results in a continuous (Proposition 3) and uniformly Lipschitzian (Proposition 4), we have demonstrated that this modified system admits a unique solution (Theorem 1). We observe that the fact that h(t, χ(t)) is not only continuous but also the Lipschitzian suggests that it has limited growth in the sense that the ratio between the variation of the ordinate and the variation of the abscissa can never exceed a fixed value. Lipschitzianity plays a key role in the uniqueness of solutions in Cauchy problems related to ordinary differential equations. It is, in fact, a condition that guarantees the existence and uniqueness of the solution for a certain initial condition. Also in this case the six models obtained as a function of the different values of atan2 admit unique solutions (Theorem 1). However, three of these systems have been discarded as providing solutions that do not adhere to the clinical reality (some terms of the respective solutions cancel out as t increases, see (17)-case (f) in Section 8.6. This is due to the fact that the solutions obtained are functions having terms in the form of negative exponential with a high time constant, so that the transient quickly extinguishes, leaving a term depending on the respiratory artifact. The three remaining systems (14)- (16) have been solved numerically through a procedure based on the four-stage Lobatto IIIa formula implemented in MatLab R R2017a using the bvp5c subroutine, once the starting problem has been transformed into an equivalent boundary value problem [25,44]. The solutions obtained by solving these three systems, as shown in Figures 3a,b and 4a, show a good regularity of the QRST wave complex comparable with those provided by the McSharry model. In these cases, the regularity of the solutions obtained largely depends on the fact that the respiratory artifact is no longer preponderant (as in the excluded cases) while the effects due to the trigonometric terms in which the arguments depend on the dynamic parameters of the system prevail.

On the Proposed Approach for Modifying the McSharry Model
The core of this paper is represented by the modification of the McSharry system (5) to simulate ECG affected by cardiac pathologies. These consist in algebraically modifying the system (5) by applying a linear transformation (with numeric parameters k 1 and k 2 ) and a rotation function of an angle φ, as detailed in (20). Obviously, by varying these parameters, different cardiac pathologies can be simulated since by altering these variables we can alter the PQRST complex as well (we point out that, if k 1 = k 2 = 1 and φ = 0, the normal sinus rhythm takes place). Three heart diseases have been taken into account: atrial fibrillation, atrial flutter and premature ventricle contraction (as they are the most common heart diseases not associated with specific symptoms and therefore at high risk), in which the parameters k 1 , k 2 and φ are specified in Table 3. After verifying that the new system (20) admits the same unstable equilibrium position as (5) (and therefore, as specified above, eliminated) only three cases have been studied (as for the McSharry model (6)) because the remaining three models provide clinically meaningless results. As far as normal sinus rhythm is concerned, all the three considered cases have provided regular trends that can be completely superimposed on ECGs taken directly by patients. Once k 1 and k 2 are set to simulate the presence of atrial fibrillation (see Table 3) the ECG obtained in the three cases considered are shown in the Figures 6a, 9b and 10b, which show how the proposed model is able to reproduce the deformation of the QRST complex related to this pathology. In these cases, according to the ECG traces taken directly on patients affected by atrial fibrillation, an increase in heart rate and an irregular R-R rhythm (at variable distances) are highlighted. The QRS complex appears narrow while the P waves are almost absent and replaced by waves with a different morphology from each other. The P-QRS ratio is not valuable. By setting k 1 , k 2 and φ to simulate the atrial flutter (see Table 3), for the three cases studied, the traces are displayed in Figures 7a, 9c and 10c. Also in these cases, in adherence to the typical clinical ECG detected directly on patients, reliable results are obtained as the obtained ECG trends show a continuous and regular atrial activation with a saw-tooth pattern, evident in the lower leads. Finally, by setting k 1 , k 2 and φ to simulate premature ventricle contraction, Figures 8a, 9d and 10d show, as clinically expected, large QRS complexes, not preceded by the P wave and therefore clearly distinguishable from sinus ones. In this case, retro activation of the atria may or may not occur. The qualitative analysis was also confirmed from the quantitative point of view. In particular, since the ECGs are often affected by uncertainty and/or imprecision due to various types of artifacts (such as breathing, incorrect application of the electrodes), they take on an evident fuzzy connotation and therefore need to be treated in the fuzzy domain. Once the fuzziness content of each ECG has been highlighted using two specific indices formulated in (26) and (27) respectively (in the sense that the high values obtained showed a high content of fuzziness in each signal), this qualitative analysis was quantitatively confirmed through the computation of particular four indices of fuzzy similarities, FS 1 , FS 2 FS 3 and FS 4 as formulated in (28)-(31) respectively. Each ECG obtained was compared to four classes of ECGs obtained directly on patients (one class for ECG without pathologies, three classes for the three pathologies chosen, respectively). Where it is not possible to directly observe the pathology form a direct qualitative analysis of the ECG, the procedure proposed registers the presence of a pathology (far from the fuzzy similarities unitary value) without classifying the type of the pathology because the fuzzy similarity values are dissimilar to the one typical of each class. Moreover, the classification puts in evidence that the detection of the presence of pathologies in an ECG is marked by comparing the same ECG with a benchmark (that is, an ECG without pathologies). In addition, each detection (for each class of pathology) presents ranges of possible values for fuzzy similarities, so that, by comparison, the kind of pathology can be deduced. The results obtained confirm the good quality of the proposed procedure. Specifically, in terms of detection, the system is able to record the presence of a pathology in an ECG so that the proposed method can be considered to all intents and purposes a classification system of cardiac pathologies present in ECGs If the system detects the presence of a cardiac anomaly in an ECG without, however, classifying it, it is possible to create an additional class of pathologies that do not fall into those considered. From the analysis of the results obtained, we deduce that the approach proposed in this paper for the generalization of the McSharry model presented in [24], appears valid in simulating certain cardiac pathologies, highlighting ECGs that are completely consistent with a real clinic. Moreover, the results obtained are encouraging as they emphasize the capacity to develop a database of pathological ECG signals with a sufficient level of differentiation as the pathology changes. In addition, the fuzzy similarity approach obtains relevant results in terms of detection and classification of pathologies in ECGs. Being the procedure characterized by a low computational complexity, it is easy to design a hardware system able to perform the classification really useful for real-time applications. Furthermore, uncertainty and imprecision are well highlighted by computation of suitable indexes of fuzziness and the proposed approach does not require pre-processing of the ECGs because the procedure works well on raw data.

Future Perspectives
However, there is no doubt that any heart disease presents, from a clinical point of view, different levels of severity. Therefore, in the future, it will be desirable to take this important aspect into account by generalizing the proposed procedure in order to simulate ECG with pathologies of different severity levels. Finally, since the ECG is affected by numerous artifacts (not only due to breathing but also, for example, to the incorrect positioning of the electrodes on the skin), it will be desirable, to integrate the model with terms that take into account such eventualities. Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix B
Proof of Proposition 2. To study the stability of (9), we exploit the first Lyapunov criterion based on the linearization of (5) around its equilibrium position [38]. To do this, we set and so that (5) becomes:u (t) = f(u(t)).
Remark A1. (A7) admits a matrix representation. Particularly, by placing: Finally, (A9) is the linearized form of the nonlinear system around the equilibrium position.
In addition, substituting both (A12) and (A13) into (A11), we obtain t = ln a 2 −ln a 1 2jω and t = ln b 2 −ln b 1 2jω from which a 2 a 1 = b 2 b 1 so that (A12) becomes: and finally ξ η = a 1 b 1 which represents, on the plane ξ − η, straight lines passing through the origin and leaving it.
The following Lemma (by Lyapunov-Poincaré) [38] is useful to prove the instability of the only equilibrium position for (5).
Thus, (A39) is forced for all the x n,j , thrust y n,j at collocation node points are obtained, for i = 1, ..., s, by y n,j = y n + ∑ s j=1 F(x x n,i , y n n,j ) x n,i x n L j (r)dr. If τ s = 1, then y n+1 = y n,s , otherwise y n+1 = y n + ∑ s j=1 F(x n,j , y n,j ) x n+1 x n L j (r)dr.
Generally, an RK procedure can be achieved as [25] y n+1 = y n + h ∑ s i=1 b i k i in which k i = F x n + c i h, y n + h ∑ s j=1 a ij k j , i = 1, 2, ..., s; (s denotes the number of stage of the procedure). {a ij }, {ci} and {b i } can be collected in the Butcher Tableau as follows [25,44] c A b T (A40) in which A = (a ij ) ∈ R s×s , b = (b 1 , ..., b s ) T ∈ R s and c = (c 1 , ..., c s ) T ∈ R s . IF a ij = 0 for j ≥ i, with i = 1, 2, ..., s, then ∀k i can be explicitly evaluated using the i − 1 coefficients k 1 ...k i−1 already computed. Therefore, the RK method is called explicit (implicit, otherwise). To compute k i , one must solve an s-dimensional nonlinear system.

The Three-Stage Lobatto IIIa Formula
This approach requires that c i must be chosen as roots of [44] P * s − P * s−2 = d s−2 dx s−2 (x s−1 (x − 1) s−1 ), in which s represents the stage, achieving both c 1 = 0 and c s = 1 ∀s. Thus, the quadrature formula is exact for any polynomial with degree lower than 2s − 2 [25,44].