Next Article in Journal
The Influence of Pressure on the Formation of FM/AF Configurations in LSMO Films: A Monte Carlo Approach
Previous Article in Journal
All-Nitrogen Cages and Molecular Crystals: Topological Rules, Stability, and Pyrolysis Paths
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Dipartimento di Ingegneria Civile Energia Ambiente e Materiali, “Mediterranea” University, Via Graziella Feo di Vito, 89060 Reggio Calabria, Italy
2
Dipartimento di Ingegneria dell’Informazione Infrastrutture Energia Sostenibile, “Mediterranea” University, Via Graziella Feo di Vito, 89122 Reggio Calabria, Italy
*
Authors to whom correspondence should be addressed.
Current address: Cittadella Universitaria, Via Graziella, Feo di Vito, I-89122 Reggio Calabria, Italy.
These authors contributed equally to this work.
Submission received: 13 October 2020 / Revised: 31 October 2020 / Accepted: 3 November 2020 / Published: 6 November 2020
(This article belongs to the Section Computational Biology)

Abstract

:
In this paper, we introduce a new dynamic model of simulation of electrocardiograms ( E C G s) affected by pathologies starting from the well-known McSharry dynamic model for the E C G s without cardiac disorders. In particular, the McSharry model has been generalized (by a linear transformation and a rotation) for simulating E C G s 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 E C G s with or without heart diseases very similar to those achieved directly on the patients. Moreover, verified that the E C G s 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 E C G s 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 E C G s affected by pathologies can be considered as a solid starting point for the development of synthetic pathological E C G s signals.

1. Introduction

Ever since Einthoven [1,2] discovered the possibility of obtaining an electrocardiogram ( E C G ) 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 E C G 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 E C G 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 E C G 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 E C G 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 E C G 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 E C G 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 E C G s, many papers have been published, including [18], in which, to generate E C G s, 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 E C G 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 E C G also turns out to be a valuable aid in the evaluation of the fetal E C G [22,23]. Researchers are also interested in the synthesis of E C G 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 E C G 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 E C G s with heart diseases. From which, it is imperative to make new models for E C G s 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 E C G s 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 E C G 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 E C G (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 E C G 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 E C G 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 E C G reconstructed in the absence and the presence of three specific pathologies are entirely comparable to the E C G obtained directly on patients. These comparisons, taking into account that E C G s 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 E C G (with or without heart diseases) obtained by the proposed model with an E C G (with or without heart diseases) sampled directly on a patient. Based on these results, our model should allow building a complete database of E C G 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 E C G system while Section 3 describes both normal conditions and pathological events in an E C G . After having described the HDV & the 12-Leads E C G 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 E C G 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 E C G s have been achieved (with or without heart diseases), they have been compared, in the fuzzy sense, in Section 10 with E C G s sampled directed on the patients. Numerical results, conclusions and some preliminary perspectives are provided in Section 11 and Section 12, respectively. Finally, mathematical details and proofs of the main results are reported at the end of the paper in appendices.

2. Normal 12-Lead ECG

An E C G 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 E C G ( 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 E C G 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 E C G 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 s c a n n i n g 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].

3. Structure of an ECG : Normal Conditions and Pathological Events

E C G 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 E C G 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). E C G 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 E C G 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 E C G signals, both normal and pathological ones [30,31,32,34].

4. HDV and 12-Leads ECG System: A Brief Overview

Let us consider a region Ω R 3 endowed with a system of Cartesian axes O x y z in which the origin O represents the central area of the heart muscle (see Figure 1a). The unit orthogonal vectors are indicated by i ^ x , i ^ y and i ^ 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 O x y z system is:
d ( t ) = d x ( t ) i ^ x + d y ( t ) i ^ y + d z ( t ) i ^ z .
As t + , d ( t ) varies cyclically in both amplitude and direction, keeping the origin fixed, depicting a 3D spatial loop in O x y z (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 i ^ x + v y i ^ y + v z i ^ z is the direction of the axis of the recording electrode, E C G ( t ) can be written as follows:
E C G ( t ) = d ( t ) · v T = v x d x ( t ) + v y d y ( t ) + v z d z ( t ) .
In accordance with (1) and (2), it is possible to perform a description of the HDV 3D using three E C G 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 E C G system is useful (albeit redundant) for diagnosis.
The 12-lead E C G system records the following 12 linearly independent signals: the three Einthoven points I, I I , I I I ; Wilson three leads a V R , a V L , a V H and the six precordial leads V 1 V 6 [1,3]. However, to get a more detailed description of E C G 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 l e a d can be achieved as follows:
ECG ( t ) 12 l e a d = H · V · d ( t ) + N ( t )
in which N ( t ) , t , represents a vector in which the components are the noise in each of the 12 E C G channels. Equation (3) represents a mathematical model of the 12-lead E C G 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 E C G ( 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 E C G ( t ) in which cardiac pathologies could also be taken into account.

5. The McSharry Dynamical Model

E C G of Patients with Heart Disease: The McSharry Single Channel Modeling Approach

In [24], the P Q R S T complex was obtained by a set of Gaussian functions with different parameters. If s ( t ) is the E C G signal, then:
s ˙ ( t ) = i { P , Q , R , S , T } a i | θ θ i | 2 π e ( | θ θ i | 2 π ) 2 2 b i 2 ( s ( t ) s 0 ( t ) )
where s ˙ ( t ) = d s ( t ) d t , | | 2 π is the module 2 π operator and ( α i , b i , θ i ) is a set of empirical parameters (see Table 2). Finally, θ represents an angular dynamic parameter that influences the deflections.
In [24], (4) was coupled with the following dynamical system:
χ ˙ ( t ) = ( 1 χ 2 ( t ) + γ 2 ( t ) ) χ ( t ) ω γ ( t ) γ ˙ ( t ) = ( 1 χ 2 ( t ) + γ 2 ( t ) ) γ ( t ) + ω χ ( t )
in which ω is a parameter concerning the heart rate, obtaining the following dynamical model:
χ ˙ ( t ) = ( 1 χ 2 ( t ) + γ 2 ( t ) ) χ ( t ) ω γ ( t ) γ ˙ ( t ) = ( 1 χ 2 + γ 2 ( t ) ) γ ( t ) + ω χ ( t ) s ˙ ( t ) = i { P , Q , R , S , T } a i | θ θ i | 2 π e ( | θ θ i | 2 π ) 2 2 b i 2 ( s ( t ) s 0 ( t ) )
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]:
θ = a t a n 2 ( γ ( t ) , χ ( t ) ) = ( a ) a t a n γ ( t ) χ ( t ) , χ ( t ) > 0 ; ( b ) a t a n γ ( t ) χ ( t ) + π , χ ( t ) < 0 γ ( t ) 0 ( c ) a t a n γ ( t ) χ ( t ) π , χ ( t ) < 0 γ ( t ) < 0 ( d ) + π 2 , χ ( t ) = 0 γ ( t ) > 0 ( e ) π 2 , χ ( t ) = 0 γ ( t ) < 0 ( f ) not defined , χ ( t ) = 0 γ ( t ) = 0
so that π a t a n 2 ( γ ( t ) , χ ( t ) ) π .
Remark 3.
| a t a n 2 ( γ ( t ) , χ ( t ) ) θ i | 2 π means that from the division between a t a n 2 ( γ , χ ) θ i and 2 π only the rest must be considered. a t a n 2 ( γ ( 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 a t a n 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:
χ ˙ ( t ) = ( 1 χ 2 ( t ) + γ 2 ( t ) ) χ ( t ) ω γ ( t ) γ ˙ ( t ) = ( 1 χ 2 ( t ) + γ 2 ( t ) ) γ ( t ) + ω χ ( t ) s ˙ ( t ) = i { P , Q , R , S , T } a i | a t a n 2 ( γ ( t ) , χ ( t ) ) θ i | 2 π e ( | a t a n 2 ( γ ( t ) , χ ( t ) ) θ i | 2 π ) 2 2 b i 2 ( s ( t ) s 0 ( t ) )
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 E C G ( 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).

6. On the Stability of the System

6.1. On the Unique Equilibrium Position

The following result holds.
Proposition 1.
(5) admits the following unique equilibrium position:
( χ 0 , γ 0 ) = ( 0 , 0 ) .
Proof. 
The proof is based on the fact that if ( χ 0 , γ 0 ) is an equilibrium position for (5), thus, t , ( χ ( t ) , γ ( t ) ) = ( χ 0 , γ 0 ) ( χ ˙ ( t ) , γ ˙ ( t ) ) = ( 0 , 0 ) . For details, see Appendix A. □

6.2. 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.

7. 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 E C G .
We introduce the following definition.
Definition 1.
f : R n R n defined on X R n , with X = { x x * x x * * } , is said to be uniformly 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 , x + Δ x X } from which { f ( x + Δ x ) f ( x * ) L Δ 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 x X , x + Δ x X , t T t + Δ t T , from which { f ( x + Δ x ) f ( x * ) L Δ x } , then, f ( x , t ) is said to be uniformly Lipschitzian in x X as t varies in T.
The following Lemma holds.
Lemma 1.
Let us consider the following Cauchy problem:
x ˙ ( t ) = f ( x , t ) x ( t 0 ) = x 0
in which f ( x , t ) is continuous as ( x , t ) varies in R n × T , 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 ) .
Remark 6.
(8) can be written as (10). In fact, setting
χ ̲ ( t ) = χ ( t ) γ ( t ) s ( t ) ; χ ̲ ˙ = χ ˙ ( t ) γ ˙ ( t ) s ˙ ( t ) ;
h ̲ ( t , χ ̲ ( t ) ) = h 1 ( t , χ ̲ ( t ) ) h 2 ( t , χ ̲ ( t ) ) h 3 ( t , χ ̲ ( t ) ) = ( 1 χ 2 ( t ) + γ 2 ( t ) ) χ ( t ) ω γ ( t ) ( 1 χ 2 ( t ) + γ 2 ( t ) ) γ ( t ) + ω χ ( t ) i { P , Q , R , S , T } a i | θ θ i | 2 π e ( | θ θ 1 | 2 π 2 b i ) 2 + ( s ( t ) s 0 ( t ) ) ,
we can write:
χ ̲ ˙ ( t ) = h ̲ ( t , χ ̲ ( t ) ) χ ̲ ( 0 ) = χ ̲ 0 .
where χ ̲ ( 0 ) = χ ̲ 0 are the initial conditions. The following results hold.
Proposition 3.
h ̲ ( t , χ ̲ ( 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. □
Proposition 4.
h ̲ ( t , χ ̲ ( t ) ) , is uniformly Lipschitzian.
Proof. 
For details, see Appendix D. □
Theorem 1.
System (8) admits unique solution.
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 a t a n 2 ( γ ( t ) , χ ( t ) ) (see (7)) provide different formulations with consequent different solution. Then, in the following Section the McSharry models obtained for different values of a t a n 2 ( γ ( t ) , χ ( t ) ) are presented and solved.

8. The McSharry Dynamical Model and the Study of the Trajectory: Peculiarities and Weaknesses

8.1. Case (a): χ ( t ) > 0 , a t a n 2 ( γ ( t ) , χ ( t ) ) = a t a n γ ( t ) χ ( t )

Here, taking into account (7), it follows that a t a n 2 ( γ ( t ) , χ ( t ) ) = a t a n γ ( t ) χ ( t ) , and considering the values in Table 2, (8) becomes:
χ ˙ ( t ) = ( 1 χ 2 ( t ) + γ 2 ( t ) ) χ ( t ) ω γ ( t ) γ ˙ ( t ) = ( 1 χ 2 ( t ) + γ 2 ( t ) ) γ ( t ) + ω χ ( t ) s ˙ ( t ) = 5 ( s ( t ) 0.15 ( 2 π f 2 t ) ) 1.1 | a t a n γ ( t ) χ ( t ) + π 3 | 2 π e | a t a n γ ( t ) χ ( t ) + π 3 | 2 π 0.25 2 2 5 | a t a n γ ( t ) χ ( t ) + π 12 | 2 π e | a t a n γ ( t ) χ ( t ) + π 12 | 2 π 0.1 2 2 30 | a t a n γ ( t ) χ ( t ) | 2 π e | a t a n γ ( t ) χ ( t ) | 2 π 0.1 2 2 + + 7.5 | a t a n γ ( t ) χ ( t ) π 12 | 2 π e | a t a n γ ( t ) χ ( t ) π 12 | 2 π 0.1 2 2 + 0.75 | a t a n γ ( t ) χ ( t ) π 2 | 2 π e | a t a n γ ( t ) χ ( t ) π | 2 π 0.4 2 2 .
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 E C G ( t ) : the waves and relative deflections are clearly visible and determined by the periodic terms contained in s ˙ ( t ) .

8.2. Case (b): χ ( t ) < 0 γ ( t ) > 0 , a t a n 2 ( γ ( t ) , χ ( t ) ) = a t a n γ ( t ) χ ( t ) + π

In this case, considering (7), (8) becomes:
χ ˙ ( t ) = ( 1 χ 2 ( t ) + γ 2 ( t ) ) χ ( t ) ω γ γ ˙ ( t ) = ( 1 χ 2 ( t ) + γ 2 ( t ) ) γ + ω χ s ˙ ( t ) = 5 ( s ( t ) 0.15 ( 2 π f 2 t ) ) 1.1 | a t a n γ ( t ) χ ( t ) + 4 π 3 | 2 π e | a t a n γ ( t ) χ ( t ) + 4 π 3 | 2 π 0.25 2 2 5 | a t a n γ ( t ) χ ( t ) + 13 π 12 | 2 π e | a t a n γ ( t ) χ ( t ) + 13 π 12 | 2 π 0.1 2 2 30 | a t a n γ ( t ) χ ( t ) | 2 π e | a t a n γ ( t ) χ ( t ) | 2 π 0.1 2 2 + + 7.5 | a t a n γ ( t ) χ ( t ) 11 π 12 | 2 π e | a t a n γ ( t ) χ ( t ) π 12 | 2 π 0.1 2 2 + 0.75 | a t a n γ ( t ) χ ( t ) π 2 | 2 π e | a t a n γ ( t ) χ ( t ) π | 2 π 0.4 2 2 .
Also in this case the periodic terms contained in s ˙ ( t ) produce the deflections in the corresponding E C G . However, the presence of noise components (mainly due to the particular value of operative parameters in s ˙ ( t ) ) would require the use of suitable filters (see Figure 3b).

8.3. Case (c): χ ( t ) < 0 γ ( t ) < 0 , a t a n 2 ( γ ( t ) , χ ( t ) ) = a t a n γ ( t ) χ ( t ) π

In this case, considering (7), (8) assumes the following form:
χ ˙ ( t ) = ( 1 χ 2 ( t ) + γ 2 ( t ) ) χ ( t ) ω γ γ ˙ ( t ) = ( 1 χ 2 ( t ) + γ 2 ( t ) ) γ + ω χ s ˙ ( t ) = 5 ( s ( t ) 0.15 ( 2 π f 2 t ) ) 1.1 | a t a n γ ( t ) χ ( t ) + 4 π 3 | 2 π e | a t a n γ ( t ) χ ( t ) + 4 π 3 | 2 π 0.25 2 2 5 | a t a n γ ( t ) χ ( t ) + 13 π 12 | 2 π e | a t a n γ ( t ) χ ( t ) + 13 π 12 | 2 π 0.1 2 2 30 | a t a n γ ( t ) χ ( t ) | 2 π e | a t a n γ ( t ) χ ( t ) | 2 π 0.1 2 2 + + 7.5 | a t a n γ ( t ) χ ( t ) 11 π 12 | 2 π e | a t a n γ ( t ) χ ( t ) π 12 | 2 π 0.1 2 2 + 0.75 | a t a n γ ( t ) χ ( t ) π 2 | 2 π e | a t a n γ ( t ) χ ( t ) π | 2 π 0.4 2 2 .
Like the two previous ones, it produces a regular E C G 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).

8.4. Case (d): χ ( t ) = 0 γ ( t ) > 0 , a t a n 2 ( γ ( t ) , χ ( t ) ) = π 2

In this case, system (8) becomes:
γ ˙ ( t ) = ( 1 γ ( t ) ) γ ( t ) s ˙ ( t ) = 8.10 e 21 5 ( z ( t ) 0.15 ( 2 π f 2 t ) ) .
However, (17) is simple to solve analytically by obtaining the following solution (for analytical details, see Appendix E):
γ ( t ) = K 2 e t K 2 e t 1 s ( t ) = e 5 t C + 0.75 1 5 sin ( 2 π f 2 t ) 2 π f 2 25 e 5 t cos ( 2 π f 2 t ) 1 + 2 π f 2 5 2 8.10 e 21 t
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.

8.5. Case (e): χ ( t ) = 0 , γ ( t ) < 0 , a t a n 2 ( γ ( t ) , χ ( t ) ) = π 2

In this case, system (8) becomes:
γ ˙ ( t ) = ( 1 γ ( t ) ) γ ( t ) s ˙ ( t ) = 1.1 e 1 0.25 2 2 ( s ( t ) 0.15 sin ( 2 π f 2 t ) ) .
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.

8.6. Case (f): χ ( t ) = γ ( t ) = 0 , a t a n 2 ( γ ( t ) , χ ( t ) ) Not Defined

We note that if χ ( t ) = γ ( t ) = 0 , it follows that a t a n 2 ( γ ( t ) , χ ( t ) ) is not defined (see option ( f ) in (7)). Moreover, χ ( t ) = γ ( t ) = 0 corresponds to the equilibrium position for (5) that results to be unstable. So that, in (7), position χ ( t ) = 0 γ ( t ) = 0 must be discarded.
Remark 8.
It is worth observing how the E C G s obtained by numerical solution of McSharry models (see Figure 3a,b and Figure 4a) are qualitatively quite similar to the E C G s obtained by clinical investigation of patients not affected by pathologies (a typical E C G related to a subject not affected by pathologies is shown in Figure 4b). However, a quantification of this similarity is required. In this regard, please refer to Section 10.
Remark 9.
As a t a n 2 ( γ ( 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 E C G ( t ) . Excluding the case (f) where a t a n 2 ( γ ( t ) , χ ( t ) ) is not defined, cases (d) and (e) provide physically unreliable results because there a t a n 2 ( γ ( t ) , χ ( t ) ) = ± π 2 . Finally, due to its structure, the McSharry model is unable to simulate E C G ( 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.

9. The Modified Heart Dipole Model (MHDM)

McSharry and co-authors were not able to characterize the distortion of the P Q R S T 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 E C G 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:
χ ^ ( t ) γ ^ ( t ) = cos ϕ sin ϕ sin ϕ cos ϕ r o t a t i o n · k 1 0 0 k 2 · χ ( t ) γ ( t ) l i n e a r t r a n s f o r m a t i o n
so that θ = a t a n 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 ϕ .
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 P Q R S T complex setting k 1 , k 2 and ϕ.
Then, the following result holds.
Proposition 5.
System (20) admits, as unique equilibrium point, the point (9).
Proof. 
For details, see Appendix F. □

MHDM Modelization

Staring from model (8) and considering the transformation (20) we can easily write:
χ ^ ˙ ( t ) = ( 1 ( k 1 χ ( t ) cos ϕ + k 2 γ ( t ) sin ϕ ) 2 + ( k 1 χ ( t ) sin ϕ + k 2 γ ( t ) cos ϕ ) 2 ) · · ( k 1 χ ( t ) cos ϕ + k 2 γ ( t ) sin ϕ ) ω ( k 1 χ ( t ) sin ϕ + k 2 γ ( t ) cos ϕ ) γ ^ ˙ ( t ) = ( 1 ( k 1 χ ( t ) cos ϕ + k 2 γ ( t ) sin ϕ ) 2 + ( k 1 χ ( t ) sin ϕ + k 2 γ ( t ) cos ϕ ) 2 ) · · ( k 1 χ ( t ) sin ϕ + k 2 γ ( t ) cos ϕ ) + ω ( k 1 χ ( t ) cos ϕ + k 2 γ ( t ) sin ϕ ) s ˙ ( t ) = i { P , Q , R , S , T } a i | a t a n 2 ( k 1 χ ( t ) sin ϕ + k 2 γ ( t ) cos ϕ , k 1 χ ( t ) cos ϕ + k 2 γ ( t ) sin ϕ ) θ i | 2 π · · e ( | a t a n 2 ( k 1 χ ( t ) sin ϕ + k 2 γ ( t ) cos ϕ , k 1 χ ( t ) cos ϕ + k 2 γ ( t ) sin ϕ ) θ i | 2 π ) 2 2 b i 2 ( s ( t ) s 0 ( t ) ) .
in which the parameters present have the same meaning as the cases studied above. System (21) admits different formulations as | a t a n 2 ( 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):
  • Case (a):
    k 1 χ ( t ) cos ϕ + k 2 γ ( t ) sin ϕ > 0 ,
    a t a n 2 ( k 1 χ ( t ) sin ϕ + k 2 γ ( t ) cos ϕ , k 1 χ ( t ) cos ϕ + k 2 γ ( t ) sin ϕ ) =
    = a t a n k 1 χ ( t ) sin ϕ + k 2 γ ( t ) cos ϕ k 1 χ ( t ) cos ϕ + k 2 γ ( t ) sin ϕ ;
  • Case (b):
    k 1 χ ( t ) cos ϕ + k 2 γ ( t ) sin ϕ < 0 k 1 χ ( t ) sin ϕ + k 2 γ ( t ) cos ϕ > 0 ,
    a t a n 2 ( k 1 χ ( t ) sin ϕ + k 2 γ ( t ) cos ϕ , k 1 χ ( t ) cos ϕ + k 2 γ ( t ) sin ϕ ) =
    = a t a n k 1 χ ( t ) sin ϕ + k 2 γ ( t ) cos ϕ k 1 χ ( t ) cos ϕ + k 2 γ ( t ) sin ϕ + π ;
  • Case (c):
    k 1 χ ( t ) cos ϕ + k 2 γ ( t ) sin ϕ < 0 k 1 χ ( t ) sin ϕ + k 2 γ ( t ) cos ϕ < 0 ,
    a t a n 2 ( k 1 χ ( t ) sin ϕ + k 2 γ ( t ) cos ϕ , k 1 χ ( t ) cos ϕ + k 2 γ ( t ) sin ϕ ) =
    = a t a n k 1 χ ( t ) sin ϕ + k 2 γ ( t ) cos ϕ k 1 χ ( t ) cos ϕ + k 2 γ ( t ) sin ϕ π ;
  • Case (d):
    k 1 χ ( t ) cos ϕ + k 2 γ ( t ) sin ϕ = 0 k 1 χ ( t ) sin ϕ + k 2 γ ( t ) cos ϕ > 0 ,
    a t a n 2 ( k 1 χ ( t ) sin ϕ + k 2 γ ( t ) cos ϕ , k 1 χ ( t ) cos ϕ + k 2 γ ( t ) sin ϕ ) = π 2 ;
  • Case (e):
    k 1 χ ( t ) cos ϕ + k 2 γ ( t ) sin ϕ = 0 k 1 χ ( t ) sin ϕ + k 2 γ ( t ) cos ϕ < 0 ,
    a t a n 2 ( k 1 χ ( t ) sin ϕ + k 2 γ ( t ) cos ϕ , k 1 χ ( t ) cos ϕ + k 2 γ ( t ) sin ϕ ) = π 2 ;
  • Case (f):
    k 1 χ ( t ) cos ϕ + k 2 γ ( t ) sin ϕ = k 1 χ ( t ) sin ϕ + k 2 γ ( t ) cos ϕ = 0 ,
    a t a n 2 ( k 1 χ ( t ) sin ϕ + k 2 γ ( t ) cos ϕ , k 1 χ ( t ) cos ϕ + k 2 γ ( t ) sin ϕ ) not defined
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), Figure 5, Figure 6a, Figure 7a and Figure 8a display the E C G obtained by varying k 1 and k 2 as per Table 3. E C G ( t ) strongly adhering to real cases are highlighted, for which the methodology proposed to simulate E C G pathologists provides qualitatively reliable results. Similarly, concerning Case (b) and (c) of the proposed methodology we also obtain E C G conforming to real cases (see Figure 9a–d and Figure 10a–d).
Remark 11.
For Case (a) it is useful to observe that the numerically obtained E C G s are completely qualitatively similar to the E C G s sampled directly on patients as evidenced by the comparison between the Figure 4b and Figure 5; between Figure 6a,b; between Figure 7a,b and, finally, between Figure 8a,b. The same comparisons can be made for Case (b): it is sufficient to compare Figure 9a–d, with Figure 4b, Figure 6b, Figure 7b and Figure 8b respectively. Analogously, concerning Case (c), it is sufficient to compare Figure 10a–d with Figure 4b, Figure 6b, Figure 7b and Figure 8b respectively. Also in these cases, for the quantitative comparison, see Section 10.
Remark 12.
It is worth noting that the E C G s 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 E C G s obtained needed filtering, in our work the E C G s obtained did not require this operation.

10. Fuzzy Similarities for Comparing ECG s

A good cardiologist is certainly able, looking at an E C G , 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 E C G with known E C G s (with or without pathologies) in order to evaluate how close (or, dually, how far away it is) an E C G from a class of known E C G s. However, if we did not have the knowledge of an expert, surely we would need a tool that can evaluate how similar an E C G is to another known E C G . Thus, taking into account the fact that the E C G s may be affected by uncertainties and/or inaccuracies, in this paper, to quantitatively evaluate how much one E C G 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].

10.1. Fuzzy Similarity Concept: Aspects for Comparing E C G s

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 E C G Considering A as a freshly sampled E C G and B as an E C G 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 E C G A is next to B. Starting from the assumption that the same pathology produces very similar E C G s, we group the E C G s into classes of pathologies converging in a single fuzzy set that represents the class of E C G s characterized by that pathology. Furthermore, if B is the fuzzy set representing the class of E C G s with a particular pathology, d ( A , B ) provides the measure, in a fuzzy sense, of how much the E C G A is similar to the class B. Furthermore, the presence of a pathology in an E C G is assessed by the reduction of its FS with the class of E C G s 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]:
S : F ( X ) × F ( X ) [ 0 , 1 ]
such that, for A, B, C F ( X ) it follows that
S ( A , B ) = S ( B , A ) to e n s u r e s y m m e t r y ;
i f A B C S ( A , B ) S ( A , C ) S ( B , C ) S ( A , C ) t o e n s u r e m o n o t o n i c i t y
S ( A , A ) = 1 to ensure reflexivity .
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 N P is the healthy E C G s class and K P is the generic E C G s disease class, then the FS of a newly sampled E C G (with unknown pathology) versus N P or a K P , F S ( N P / K P ) , can be calculated by means of a function satisfying the above properties. In this work, we verified the fuzziness of the classes of E C G s using specific entropy indices, such as the linear index ( L I ) [39,40,41]
L I = 2 n i = 1 n min ( μ N P / K P ( · ) ( 1 μ N P / K P ( · ) ) )
and the fuzzy entropy index ( F E )
F E = 1 n i = 1 n min ( μ N P / K P log ( μ N P / K P ( · ) ) ( 1 μ N P / K P ( · ) ) log ( 1 μ N P / K P ( · ) )
and considering that high values of L I and F E indicate high fuzziness in the E C G s belonging to the class, in this paper, four different formulations of F S have been used [39,40,41]:
F S 1 ( N P / K P ) = 1 n i = 1 n min ( μ N P / K P ( · ) μ Y ( · ) ) max ( μ N P / K P ( · ) μ Y ( · ) ) ;
F S 2 ( N P / K P ) = 1 μ N P / K P ( · ) μ Y ( · ) s n ;
F S 3 ( N P / K P ) = 1 i = 1 n μ N P / K P ( · ) μ Y ( · ) s μ N P / K P ( · ) + μ Y ( · ) ;
F S 4 ( N P / K P ) = 1 1 + μ N P / K P ( · ) μ Y ( · ) s ;
in which s { 1 , 2 , } ; n is the number of samples (which is the same for all E C G s); μ N P / K P ( · ) and μ Y ( · ) represent the membership values related to N P , K P and Y, respectively. It is worth noting that, if Y is an E C G without pathologies, then μ Y ( x i ) = μ N P ( x i ) = 1 . Otherwise, if it is related to a certain pathology K P , then μ Y ( x i ) = μ K P ( x i ) < 1 .

10.2. The Exploited Database

In this paper, MIT-BIH Arrhythmia Database by PhysioNet (public domain database of E C G s, 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:
  • NSR sub-directory, which contains E C G s characterized by normal sinus rhythm;
  • AFL sub-directory containing E C G s affected by atrial flutter;
  • AFIB sub-directory in which E C G s with atrial fibrillation are contained;
  • PVC sub-directory containing E C G s with premature ventricle contraction.
Particularly, each recording contains E C G s 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 E C G s are rather noisy, so that they have been suitably filtered. Four classes of E C G s have been set up: the first class is made up of 140 E C G s of subjects without heart disease; the second class, also made up of 140 E C G s, concerns patients with atrial fibrillation; the third class, as for the previous ones formed by 140 E C G s, contains signals from patients suffering from atrial flutter; finally, the last class contains 140 E C G s of patients with premature ventricle contraction.

10.3. Some Statistical Considerations

It is worth noting that the E C G s concerning a particular pathology are qualitatively similar to each other. This has an impact on the fact that the histograms relating to the similar E C G s signals show strong similarities. In parallel, also E C G s signals without pathologies are qualitatively similar to each other. Also in this case, the relative histograms are similar. Usually, E C G s histograms highlight Gaussian trends as depicted in Figure 11a–d where histograms for typical E C G s characterized by normal sinus rhythm, atrial fibrillation, atrial flutter and premature ventricle contraction, respectively, are shown. Therefore, it makes sense to collect E C G s characterized by the same pathology in one single class while all the E C G s without pathologies can be considered as belonging to a single class. In this way, we obtain four classes of E C G s 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 E C G s that show atrial flutter while the last class, the fourth, collects the 140 E C G s, which show premature ventricle contraction.
Thus, it is sufficient to consider mean m i j and standard deviation σ i j to characterize the ith E C G belonging to the jth. Then, different E C G s 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 K P j , j { 1 , 2 , 3 } , is the jth class of E C G s, it follows that
μ K P j = 1 140 m i j 140
that is the mean value of the means and
σ K P j = max { σ i j }
that is the maximum value of all the standard deviation. Similarly, the same procedure has been applied to the class of E C G s without pathologies achieving m N P and σ N P . To build the membership function related to each class K P j and N P we make the relative Gaussian membership functions:
μ K P j = e ( x m K P j ) 2 2 σ K P j 2
and
μ N P = e ( x m N P j ) 2 2 σ N P j 2 .
In this framework, it makes sense to consider the mean of the means of the E C G s 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 E C G s belonging to that class.

10.4. 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 E C G s without pathologies obtained L I = 0.97 and F E = 0.98 , respectively, while the class of the E C G s with atrial fibrillation obtained L I = 0.94 and F E = 0.95 , respectively. Moreover, concerning the class of the E C G s with atrial flutter achieved L I = 0.98 and F E = 0.92 , respectively. Finally, regarding the class of the E C G s with premature ventricle contraction, L I = 0.97 and F E = 0.97 were obtained, respectively. This justifies the use of the fuzzy similarity based approach for the analysis of E C G s.
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 E C G s without heart disease, and the other three related to E C G s affected by atrial fibrillation, atrial flutter and premature ventricle contraction, respectively). Therefore, for each E C G numerically obtained by the McSharry dynamical model (Figure 3a,b and Figure 4a) the four fuzzy similarities, F S 1 , F S 2 , F S 3 and F S 4 , with respect to the four classes of E C G s have been computed and the results have been reported in Table 4 respectively, highlighting as Figure 3a,b and Figure 4a belongs to the class of normal E C G s (i.e., without heart diseases) due to the high fuzzy similarities obtained. This is supported by the fact that the fuzzy similarities values of Figure 3a,b and Figure 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 E C G s obtained with the proposed approach were computed with the representative E C G s 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 E C G s of patients without cardiac pathologies because the F S 1 , F S 2 , F S 3 and F S 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 E C G depicted in Figure 6a belongs to the class of E C G s affected by atrial fibrillation because F S 1 , F S 2 , F S 3 and F S 4 values achieved are very close to unity (in bold). Again, Table 4 proves that the E C G shown in Figure 7a belongs to the class of E C G s affected by atrial flutter ( F S 1 , F S 2 , F S 3 and F S 4 are very close to unity) (in bold). Finally, Table 4 puts in evidence that the E C G depicted in Figure 8a belongs to the class of E C G s affected by premature ventricle contraction ( F S 1 , F S 2 , F S 3 and F S 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 Table 5 and Table 6, respectively.

11. Conclusions

11.1. 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 E C G ( t ) , so that this undesirable instability is, in fact, eliminated.

11.2. 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 a t a n 2 (see (7)), thus that the six specific McSharry systems have been derived, according to the value assumed by a t a n 2 (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 a t a n 2 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® 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 Figure 3a,b and Figure 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.

11.3. 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 E C G 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 E C G s taken directly by patients. Once k 1 and k 2 are set to simulate the presence of atrial fibrillation (see Table 3) the E C G obtained in the three cases considered are shown in the Figure 6a, Figure 9b and Figure 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 E C G 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 Figure 7a, Figure 9c and Figure 10c. Also in these cases, in adherence to the typical clinical E C G detected directly on patients, reliable results are obtained as the obtained E C G 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, Figure 8a, Figure 9d and Figure 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 E C G s 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 E C G 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, F S 1 , F S 2 F S 3 and F S 4 as formulated in (28)–(31) respectively. Each E C G obtained was compared to four classes of E C G s obtained directly on patients (one class for E C G 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 E C G , 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 E C G is marked by comparing the same E C G with a benchmark (that is, an E C G 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 E C G so that the proposed method can be considered to all intents and purposes a classification system of cardiac pathologies present in E C G s If the system detects the presence of a cardiac anomaly in an E C G 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 E C G s 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 E C G 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 E C G s. 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 E C G s because the procedure works well on raw data.

12. 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 E C G with pathologies of different severity levels. Finally, since the E C G 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.

Author Contributions

Conceptualization, M.V. and F.L.F.; methodology, M.V.; software, G.A.; validation, M.V., F.L.F. and G.A.; formal analysis, M.V.; investigation, F.L.F.; resources, G.A.; data curation, G.A.; writing—original draft preparation, M.V.; writing—review and editing, F.L.F.; visualization, M.V.; supervision, F.L.F.; project administration, G.A. All authors have read and agreed to the published version of the manuscript.

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:
E C G Electrocardiogram
AV, AV nodeAtrial-ventricular, atrial-ventricular node
VCGVector Cardiogram
HDV, HDM, MHDMHeart Dipole Vector, Heart Dipole Model, Modified Heart Dipole Model
T, tAbsolute temperature and time variable
d ( t ) Time-varying rotating vector
d x ( t ) , d y ( t ) , d z ( t ) Cartesian component of d ( t )
I, II, IIIStandard leads
aVR, aVL, aVFGoldberger leads
V 1 , V 2 , V 3 , V 4 , V 5 , V 6 Wilson leads
P, Q, R, S, T, UTypical waves in E C G
v direction of the axe of the recording electrode
V , H , N ( t ) Matrices useful to reconstruct E C G
| | 2 π Module 2 π operator
α i , b i , θ i Set of empirical parameters
θ Angular dynamical parameter influencing the deflections in an E C G
γ ( t ) , χ ( t ) Dynamical variables
f 0 , ω Respiratory frequency, angular velocity of the trajectory
FSFuzzy Similarity

Appendix A

Proof of Proposition 1.
From the theory of dynamic systems [38], setting
χ ˙ ( t ) = ( 1 χ 2 ( t ) + γ 2 ( t ) ) χ ( t ) ω γ ( t ) = 0 γ ˙ ( t ) = ( 1 χ 2 + γ 2 ( t ) ) γ ( t ) + ω χ ( t ) = 0
we easily obtain γ 2 ( t ) = χ 2 ( t ) , from which γ ( t ) = χ ( t ) = 0 , so that (9) yields. □

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
u ( t ) = χ ( t ) γ ( t ) , u ˙ ( t ) = d χ ( t ) d t d γ ( t ) d t
and
f ( t , u ( t ) ) = f ¯ ( χ ( t ) , γ ( t ) ) g ¯ ( χ ( t ) , γ ( t ) ) = ( 1 χ 2 ( t ) + γ 2 ( t ) ) χ ω γ ( t ) ( 1 χ 2 ( t ) + γ 2 ( t ) ) γ ( t ) + ω χ ( t ) ,
so that (5) becomes:
u ˙ ( t ) = f ( u ( t ) ) .
To linearize (A4), we set the following change of variables [38]:
χ ( t ) = χ 0 + ϵ ξ ( t ) γ ( t ) = γ 0 + ϵ η ( t )
with ϵ small enough. From (A4), taking into account (A5) and considering that both χ 0 and γ 0 do not depend on t, we can easily write d χ ( t ) d t = ϵ d ξ ( t ) d t = f ¯ ( χ ( t ) , γ ( t ) ) and d γ ( t ) d t = ϵ d η ( t ) d t = g ¯ ( χ ( t ) , γ ( t ) ) . Therefore, developing in Taylor series both f ¯ ( χ ( t ) , γ ( t ) ) and g ¯ ( χ ( t ) , γ ( t ) ) we can easily write:
ϵ d ξ ( t ) d t = f ¯ ( χ 0 + ϵ ξ ( t ) , γ 0 + ϵ η ( t ) ) f ¯ ( χ 0 , γ 0 ) + ϵ f ¯ ( χ 0 , γ 0 ) d χ χ ( t ) + ϵ f ¯ ( χ 0 , γ 0 ) d γ η ( t ) + o ( τ ) ϵ d η ( t ) d t = g ¯ ( χ 0 + ϵ ξ ( t ) , γ 0 + ϵ η ( t ) ) g ¯ ( χ 0 , γ 0 ) + ϵ g ¯ ( χ 0 , γ 0 ) d χ χ ( t ) + ϵ f ¯ ( χ 0 , γ 0 ) d γ η ( t ) + o ( τ )
where τ = χ 2 + η 2 . Moreover, f ¯ ( χ 0 , γ 0 ) = g ¯ ( χ 0 , γ 0 ) = 0 , we get:
d ξ ( t ) d t = f ¯ ( χ 0 , γ 0 ) χ ξ ( t ) + f ¯ ( χ 0 , γ 0 ) γ ξ ( t ) d η ( t ) d t = g ¯ ( χ 0 , γ 0 ) χ ξ ( t ) + g ¯ ( χ 0 , γ 0 ) γ ξ ( t )
which represents a system of first-order ordinary differential equations.
Remark A1.
(A7) admits a matrix representation. Particularly, by placing:
z = ξ ( t ) η ( t ) , z ˙ = d ξ ( t ) d t d η ( t ) d t , A = f ¯ ( χ 0 , γ 0 ) χ f ¯ χ 0 , γ 0 ) γ g ¯ ( χ 0 , γ 0 ) χ g ¯ ( χ 0 , γ 0 ) γ
(A7) becomes:
z ˙ = A z .
Finally, (A9) is the linearized form of the nonlinear system around the equilibrium position.
In our case, from (A3), we can easily obtain f ¯ ( χ 0 , γ 0 ) χ = 1 ; f ¯ ( χ 0 , γ 0 ) γ = ω ; g ¯ ( χ 0 , γ 0 ) χ = ω ; f ¯ ( χ 0 , γ 0 ) γ = 1 so that
A = 1 ω ω 1
and
ξ ˙ ( t ) η ˙ ( t ) = 1 ω ω 1 = ξ ( t ) η ( t ) .
To obtain the eigenvalues for (A10), we solve d e t ( A λ I ) = ( 1 λ ) ( 1 λ ) + ω 2 = 0 from which the eigenvalue, λ 1 and λ 2 , become λ 1 = 1 + j ω and λ 2 = 1 j ω where their real part, for both, is positive. Then, the unique equilibrium position (9) is unstable [38]. This is easily verifiable considering that ξ ( t ) and η ( t ) are writable in the form [38]:
ξ ( t ) = a 1 e λ 1 t + a 2 e λ 2 t = a 1 e t e j ω t + a 2 e t e j ω t η ( t ) = b 1 e λ 1 t + b 2 e λ 2 t = b 1 e t e j ω t + b 2 e t e j ω t
where a 1 , a 2 , b 1 and b 2 constant. Observing that e j ω t and e j ω t are both oscillating terms, as t + , | ξ ( t ) | and | η ( t ) | become unlimited and the equilibrium is unstable. Moreover, deriving the (A12) we easily obtain:
ξ ˙ ( t ) = a 1 e t e j ω t + a 1 ω e t e j ω t + a 2 e t e j ω t a 2 j ω e t e j ω t η ˙ ( t ) = b 1 e t e j ω t + b 1 ω e t e j ω t + b 2 e t e j ω t b 2 j ω e t e j ω t .
In addition, substituting both (A12) and (A13) into (A11), we obtain t = ln a 2 ln a 1 2 j ω and t = ln b 2 ln b 1 2 j ω from which a 2 a 1 = b 2 b 1 so that (A12) becomes:
ξ ( t ) = a 1 e t e j ω t + b 2 b 1 e t e j ω t η ( t ) = b 1 e t e j ω t + b 2 b 1 e t e j ω t
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).
Lemma A1.
Let us consider the autonomous system (A4). For it, let [ 0 ] be an equilibrium point of the system (i.e., f ( [ 0 ] ) = [ 0 ] ) and let f continue in around [ 0 ] and differentiable into [ 0 ] , with Jacobian matrix [ J ( f ( [ 0 ] ) ) ] = A . If [ 0 ] is stable for the linearized system (A9), then [ 0 ] is stable for (A4). Also, if there exists an eigenvalue of A with a positive real part, then [ 0 ] ) is an unstable equilibrium point for (A4).
System (5) can be written in the form (A4). Moreover, (A3) is continued around [ 0 ] and differential into [ 0 ] . Then, exploiting Lemma A1, [ 0 ] is an unstable equilibrium point for (A4). □

Appendix C. Proof of Proposition 3

( 1 χ ( t ) + γ 2 ( t ) ) χ ( t ) ω γ ( t ) and ( 1 χ ( t ) + γ 2 ( t ) ) γ ( t ) + ω χ ( t ) are clearly continuous functions. Moreover, taking into account that θ = a t a n 2 γ ( t ) χ ( t ) (see (7)) we obtain:
i { P , Q , R , S , T } a i | θ θ i | 2 π e | θ θ i | 2 π 2 b i 2 ( s ( t ) s 0 ( t ) ) = a P | θ θ P | 2 π e | θ θ P | 2 π 2 b P 2 ( s ( t ) s 0 ( t ) ) + a Q | θ θ Q | 2 π e | θ θ Q | 2 π 2 b Q 2 ( s ( t ) s 0 ( t ) ) + a R | θ θ R | 2 π e | θ θ R | 2 π 2 b R 2 ( s ( t ) s 0 ( t ) ) + a S | θ θ S | 2 π e | θ θ S | 2 π 2 b S 2 ( s ( t ) s 0 ( t ) ) + a T | θ θ T | 2 π e | θ θ T | 2 π 2 b T 2 ( s ( t ) s 0 ( t ) ) .
Exploiting (7), it follows that a t a n 2 γ ( t ) χ ( t ) is a continuous function, so that h ̲ ( t , χ ̲ ( t ) ) is continuous as ( t , χ ̲ ( t ) ) R 3 × T varies.

Appendix D. Proof of Proposition 4

We need to prove that:
| h 1 ( t 1 , χ ̲ ( t 1 ) ) h 1 ( t 2 , χ ̲ ( t 2 ) ) | L 1 | t 1 t 2 |
| h 2 ( t 1 , χ ̲ ( t 1 ) ) h 2 ( t 2 , χ ̲ ( t 2 ) ) | L 2 | t 1 t 2 |
| h 3 ( t 1 , χ ̲ ( t 1 ) ) h 3 ( t 2 , χ ̲ ( t 2 ) ) | L 3 | t 1 t 2 |
with L 1 , L 2 and L 3 are constant positive and, in addition, h 1 ( t , χ ̲ ( t ) ) , h 2 ( t , χ ̲ ( t ) ) and h 3 ( t , χ ̲ ( t ) ) defined as in (12). From (A16), taking into account the (12), we can write:
| h 1 ( t 1 , χ ̲ ( t 1 ) ) h 1 ( t 2 , χ ̲ ( t 2 ) ) | = = | ( 1 χ 2 ( t 1 ) + γ 2 ( t 1 ) ) χ ( t 1 ) ω γ ( t 1 ) ( 1 χ 2 ( t 2 ) + γ 2 ( t 2 ) ) χ ( t 2 ) ω γ ( t 2 ) | .
Experimentally, χ ( t ) and γ ( t ) can be expressed as [30]:
χ ( t ) = H sin ( 2 π f t )
γ ( t ) = K sin ( 2 π f t )
with H and K constant positive and f a suitable frequency. Then, from (A19) and taking into account both (A20) and (A21), we can also write:
| h 1 ( t 1 , χ ̲ ( t 1 ) ) h 1 ( t 2 , χ ̲ ( t 2 ) ) | = = | H ( sin ( 2 π f t 1 ) sin ( 2 π f t 2 ) ) ω K ( sin ( 2 π f t 1 ) sin ( 2 π f t 2 ) ) H 2 sin 2 ( 2 π f t 1 ) + K 2 sin 2 ( 2 π f t 1 ) H sin ( 2 π f t 2 ) + + H 2 sin 2 ( 2 π f t 2 ) + K 2 sin 2 ( 2 π f t 2 ) H sin ( 2 π f t 2 ) | 2 π f ( H + ω K ) + 2 π H f ) | t 1 t 2 |
from which
| h 1 ( t 1 , χ ̲ ( t 1 ) ) h 1 ( t 2 , χ ̲ ( t 2 ) ) | 2 π f ( H + ω K ) + 2 π H f ) | t 1 t 2 | .
Similarly, (A17) is proved. Regarding the (A18), it is easy to verify that:
| h 3 ( t 1 , χ ̲ ( t 1 ) ) h 3 ( t 2 , χ ̲ ( t 2 ) ) | = 0.75 | sin ( 2 π f t 1 ) sin ( 2 π f t 2 ) | 1.5 π f | t 1 t 2 |
from which | h 3 ( t 1 , χ ̲ ( t 1 ) ) h 3 ( t 2 , χ ̲ ( t 2 ) ) | 1.5 π f | t 1 t 2 | .

Appendix E. System (17) and Its Solution

To obtain (17), we observe that χ ( t ) = 0 and γ ( t ) > 0 . Moreover, a t a n 2 ( γ ( t ) , χ ( t ) ) = π 2 . Then, the first two equations of the system (8) become ω γ ( t ) = 0 and γ ˙ ( t ) = ( 1 γ ( t ) ) γ ( t ) . However, being γ ( t ) > 0 , it follows that ω = 0 . Then (8), taking into account the values contained in Table 2 and observing that | π 2 + π 3 | 2 π = 5 12 | π 2 + π 12 | 2 π = 7 24 ; | π 2 + 0 | 2 π = 1 4 ; | π 2 π 12 | 2 π = 5 24 , becomes:
γ ˙ ( t ) = ( 1 γ ( t ) ) γ ( t ) s ˙ ( t ) = 8.10 e 21 5 ( s ( t ) 0.15 ( 2 π f 2 t ) ) .
From the first equation of the system (A25), we can write d γ ( t ) d t = ( 1 γ ( t ) ) γ ( t ) from which γ ( t ) = K 2 e t K 2 e t 1 . From the second equation in (A25), we can write s ˙ ( t ) = 8.10 e 21 5 s ( t ) + 0.75 sin ( 2 π f 2 t ) from which, by integration, we achieve:
s ( t ) = e 5 t C + 0.75 1 5 sin ( 2 π f 2 t ) 2 π f 2 25 e 5 t cos ( 2 π f 2 t ) 1 + 2 π f 2 5 2 8.10 e 21 t
so that system (A25) becomes:
γ ( t ) = K 2 e t K 2 e t 1 s ( t ) = e 5 t C + 0.75 1 5 sin ( 2 π f 2 t ) 2 π f 2 25 e 5 t cos ( 2 π f 2 t ) 1 + 2 π f 2 5 2 8.10 e 21 t
with the K 2 and C constant to be determined exploiting suitable initial conditions.

Appendix F. Proof of Proposition (5)

From (20) we can write:
χ ^ ( t ) γ ^ ( t ) = c o s ϕ sin ϕ sin ϕ cos ϕ k 1 0 0 k 2 χ ( t ) γ ( t ) = k 1 χ ( t ) cos ϕ + k 2 γ ( t ) sin ϕ k 1 χ ( t ) sin ϕ + k 2 γ ( t ) cos ϕ .
from which
χ ^ 2 ( t ) + γ ^ 2 ( t ) = k 1 2 χ 2 ( t ) + k 2 2 γ 2 ( t ) .
Taking into account (20) and exploiting (A28) and (A29), we can write
χ ^ ˙ ( t ) = ( 1 k 1 2 χ 2 ( t ) + k 2 2 γ 2 ( t ) ) ( k 1 χ ( t ) cos ϕ + k 2 γ ( t ) sin ϕ ) ω ( k 1 χ ( t ) cos ϕ + k 2 γ ( t ) cos ϕ ) = 0 γ ^ ˙ ( t ) = ( 1 k 1 2 χ 2 ( t ) + k 2 2 γ 2 ( t ) ) ( k 1 χ ( t ) sin ϕ + k 2 γ ( t ) cos ϕ ) + + ω ( k 1 χ ( t ) cos ϕ + k 2 γ ( t ) cos ϕ ) = 0
from which
( k 1 χ ( t ) cos ϕ + k 2 γ ( t ) cos ϕ ) 2 = ( k 1 χ ( t ) cos ϕ + k 2 γ ( t ) sin ϕ ) 2
and finally:
k 1 χ ( t ) cos ϕ + k 2 γ ( t ) cos ϕ = ± ( k 1 χ ( t ) cos ϕ + k 2 γ ( t ) sin ϕ ) .
From condition (A31), and considering the sign “−” on the right side, we can write k 2 γ ( t ) cos ϕ = 0 from which
γ ( t ) = 0 .
If in (A32) we consider the sign “+” we easily obtain 2 k 1 χ ( t ) + k 1 γ ( t ) sin ϕ = 0 and taking into account the condition (A33) we achieve:
χ ( t ) = 0 .
Both conditions (A33) and (A34) correspond to the condition (9).

Appendix G. III/IV-Stage Lobatto IIIa Formula

Let us consider the following system of ordinary differential equations (ODEs) [44]:
d y ( r ) d r = F ( r , y ( r ) ) G [ y ( a ) , y ( b ) ] = 0
in which G [ y ( a ) , y ( b ) ] = 0 represents the usual boundary conditions. (A35) can be converted into the following integral equation
y ( x ) = y ( x n ) + x n x F ( r , y ( r ) ) d r .
and replacing y ( x n ) by y n (approximated value), we achieve:
y ( x ) y n + x n x p ( r ) d r ,
where p ( r ) is an interpolation polynomial of degree lower than s interpoling [ x n , i , F ( x n , i ) , y ( x n , i ) ] , i = 1 , 2 , , s , and x n , i = x n + τ i h , i = 1 , , s , 0 τ 1 < < τ s 1 . Exploiting the Lagrange interpolation polynomial technique [25], p ( r ) can be computed as follows:
p ( r ) = j = 1 s F ( x n , j , y ( x n , j ) ) L j ( r ) ,
with L j ( r ) , the fundamental Lagrange polynomials [25,44]. Then, plugging (A38) into (A37) we easily achieve:
y ( x ) y n + j = 1 s F ( x n , j , y ( x n , j ) ) x n x L j ( r ) d r .
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 + j = 1 s F ( x x n , i , y n n , j ) x n x n , i L j ( r ) d r . If τ s = 1 , then y n + 1 = y n , s , otherwise y n + 1 = y n + j = 1 s F ( x n , j , y n , j ) x n x n + 1 L j ( r ) d r .

Appendix G.1. A Brief Overview on the Runge–Kutta Procedure

Runge–Kutta (RK) procedures involve many evaluations of the function F ( x , y ( x ) ) , [ x n , x n + 1 ] . Generally, an RK procedure can be achieved as [25] y n + 1 = y n + h i = 1 s b i k i in which k i = F x n + c i h , y n + h j = 1 s a i j k j , i = 1 , 2 , , s ; (s denotes the number of stage of the procedure). { a i j } , { c i } and { b i } can be collected in the Butcher Tableau as follows [25,44]
c A b T
in which A = ( a i j ) R s × s , b = ( b 1 , , b s ) T R s and c = ( c 1 , , c s ) T R s . IF a i j = 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 d x 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 2 s 2 [25,44].
Definition A1.
(step-size definition)
On the mesh-grid:
0 = a = r 0 < r 1 < < r n = b = R
we define the step-size h m = r m + 1 r m .
Definition A2.
(midpoint and approximation at the midpoint)
( r m , r m 1 ) , we denote the related midpoints by r m + 1 / 2 . Moreover, by y m + 1 / 2 , we denote the approximation of y ( r ) at r m + 1 / 2 .

Appendix G.2. A Brief Discussion on the Order of the Polynomial

Remark A2.
The cubic polynomial p ( r ) satisfies the boundary conditions in (A35) and, moreover, ( r m , r m + 1 ) , the grid (A41) is considered. In addition, p ( r ) is located at the edges of each sub-interval and midpoint as well where p ( r ) is continuous.
This approach is a collocation procedure and, as proved in [25], it is equivalent to the three-stage Lobatto IIIa implicit RK method in which the Butcher tableau is as follows [25,44]
0 0 0 0 1 2 5 24 1 3 1 24 1 1 6 2 3 1 6 1 6 2 3 1 6
Then, the three-stage Lobatto IIIa formula is writable as [44] y m + 1 / 2 = y m + h m 5 24 F ( r m , y m ) + 1 3 F ( r m + 1 / 2 , y m + 1 / 2 ) 1 24 F ( r m + 1 , y m + 1 ) and y m + 1 = y m + h m 1 6 F ( r m , y m ) + 2 3 F ( r m + 1 / 2 , y m + 1 / 2 ) + 1 6 F ( r m + 1 , y m + 1 ) .

Appendix G.3. On the Use of Simpson Quadrature Formula

Exploiting the Simpson quadrature formula, (A36) easily becomes [25] y m + 1 = y m + h m 6 F ( r m , y m ) + F ( r m + 1 , y m + 1 ) + 4 F r m + 1 / 2 , y m + 1 + y m 2 + h m 8 [ F ( r m , F ( r m , y m ) F ( r m + 1 , y m + 1 ) ] .

Appendix G.4. p ( r ) and Its Derivatives

p ( r ) and their derivatives satisfy, r ( a , b ) , p ( l ) ( r ) = y ( l ) ( r ) + o ( h 4 l ) , l = 0 , 1 , 2 , 3 [25]. Moreover, (A35) are satisfied by p ( r ) at each intermediate point and at the midpoint of each interval. Furthermore, we note that p ( r ) is chosen by MatLab® determining unknown parameters, if any. Finally, we can easily write p ( r m ) = F [ r m , p ( r m ) ] , p ( r m + 1 / 2 ) = F [ r m + 1 / 2 , p ( r m + 1 / 2 ) ] , p ( r m + 1 ) = F [ r m + 1 , p ( r m + 1 ) which represent nonlinear equations solvable by a MatLab® solver. In addition, MatLab®, r ( a , b ) , evaluates the cubic polynomial by its special routine bvpval [25].

Appendix G.5. A Guess for the Solution and Initial Mesh

As known, a BVP could have more than one solution [25]. Thus, it is necessary to give a guess for both the initial mess and the solution. The MatLab® solver adapts the mesh achieving a solution by a reduced number of mesh points [25,44]. Usually, a good initial hypothesis is very difficult. Then, the MatLab® solver acts by checking a residue [25] res ( r ) = p ( r ) F [ r , p ( r ) ] while the boundary conditions become g [ p ( a ) , p ( b ) ] . If res ( r ) is small, p ( r ) represents a good solution and, if the problem is well-conditioned, p ( r ) is next to y ( r ) . Here, the MatLab® R2017a bvp5c solver has been used because, on one hand, it implements the collocation procedure by a piece-wise cubic p ( r ) , with coefficients determined requiring that p ( r ) be continuous on ( a , b ) . On the other hand, both mesh and estimation error are based on the evaluation of the residual of p ( r ) in which the control is useful to manage poor or inadequate guesses for both mesh and solution [25,44]. Moreover, this routine is characterized by a very low computational burden to obtain the Jacobian. Finally, bvp4c being a vectorized procedure, it reduces the run-time vectoring F ( r , y ( r ) ) [25,44].

Appendix G.6. Four-Stage Lobatto IIIa Formula

It derives as an implicit RK method with Butcher tableau [25,44]:
0 0 0 0 0 5 5 10 11 + 5 120 25 5 120 25 13 5 120 1 + 5 120 5 + 5 10 11 5 120 25 + 13 5 120 25 + 5 120 1 5 120 1 1 12 5 12 5 12 1 12 1 12 5 12 5 12 1 12

References

  1. Hampton, J. The ECG Made Easy; Elsevier: London, UK; New York, NY, USA; Oxford, UK; Philadelphia, PA, USA, 2019. [Google Scholar]
  2. Gargiulo, G.D.; Bifulco, P.; Cesarelli, M.; McEwan, A.L.; Moeinzadeh, H.; O’Loughlin, A.; Shugman, I.M.; Tapson, J.; Thiagalingam, A. On the Einthoven Triangle: A Critical Analysis of the Single Rotating Dipole Hypothesis. Sensors 2018, 18, 2353. [Google Scholar]
  3. Labatin, R.D.; Munoz, E.; Piuri, V.; Sassi, R.; Scotti, F. Deep-ECG: Convolutional Neural Newtorks for ECG biometric Recognition. Pattern Recognit. Lett. 2019, 20, 168–181. [Google Scholar]
  4. Biglu, M.H.; Ghavami, M.; Biglu, S. Cardiovascular Diseases in the Mirror of Science. J. Cardiovasc. Thorac. Res. 2016, 8, 158–163. [Google Scholar] [PubMed] [Green Version]
  5. La Foresta, F.; Morabito, F.C.; Azzerboni, B.; Ipsale, M. PCA and ICA for the extraction of EEG dominant components in cerebral death assessment. Proc. Int. Jt. Conf. Neural Netw. 2005, 4, 32532–32537. [Google Scholar]
  6. La Foresta, F.; Mammone, N.; Morabito, F.C. Independent component and wavelet analysis for ECG extraction: The ST waveform evaluation. In Proceedings of the 2005 IEEE International Conference on Computational Intelligence for Measurement Systems and Applications, Messina, Italy, 20–22 July 2005; pp. 86–90. [Google Scholar]
  7. Alqudah, A.M.; Albadarneh, A.; Abu-Qasmieh, I.; Alquran, H. Developing of Robust and High Accurate ECG Beat Classification by Combining Gaussian Mixture s and Wavelets Features. Australas. Phys. Eng. Sci. Med. 2019, 42, 149–157. [Google Scholar]
  8. Attia, Z.I.; DeSimone, C.V.; Dillon, J.J.; Sapir, Y.; Somers, V.K.; Dugan, J.L.; Bruce, C.J.; Ackerman, M.J.; Asirvatham, S.J.; Striemer, B.L.; et al. Novel Bloodless Potassium Determination Using a Signal-Processed Single-Lead ECG. J. Am. Heart Assoc. 2016. [Google Scholar] [CrossRef] [Green Version]
  9. Němcová, A.; Smíšek, R.; Maršánová, L.; Smital, L.; Vítek, M. A Comparative Analysis of Methods for Evaluation of ECG Signal Quality after Compression. BioMed Res. Int. 2018, 2018, 1–26. [Google Scholar]
  10. Stauffer, F.; Thielen, M.; Sauter, C.; Chardonnens, S.; Bachmann, S.; Tybrandt, K.; Peters, C.; Hierold, C.; Vörös, J. Skin Conformal Polymer Electrodes for Clinical ECG and EEG Recording. Adv. Healthc. Mater. 2018, 12, 422–441. [Google Scholar]
  11. Manikandan, M.S.; Dandapat, S. Wavelet Energy Based Diagnostic Distortion Measure for ECG. Biomed. Signal Process. Control 2007, 2, 80–96. [Google Scholar]
  12. Strodthoff, N.; Strdothoff, C. Detecting and Interpreting Myocardial Infarction Using Fully Convolutional Neural Networks. Physiol. Meas. 2019, 40, 015001. [Google Scholar]
  13. Selder, J.L.; Breukel, L.; Block, S.; van Rossum, A.C.; Tulevski, I.I.; Allaart, C.P. A Mobile One-Lead ECG Device Incorporated in a Symptom-Driven Remote Arrhythmia Monitoring Programm. The First 5982 Hartwacht ECGs. Neth. Hearth J. 2019, 27, 38–45. [Google Scholar]
  14. Tabassum, T.; Ahmed, M. A Simplified Cardiac Conduction Model and Twelve-Lead ECG Generation. In Proceedings of the IEEE International Conference on Computer, Electrical & Comunication Engineering (ICCECE), Kolkata, India, 17–18 January 2020; pp. 1–5. [Google Scholar]
  15. Hernandez-Caceres, J.L.; Gonzales-Fernandez, R.; Ontivero-Ortega, M.; Molte, G. Introducing BisQ, A Bicoherence-Based Nonlinear Index to Explore the Heart Rhythm. Math. Comput. Appl. 2020, 25, 45. [Google Scholar]
  16. Trayanova, N.A. Whole Heart Modeling: Applications to Cardiac Electrophysiology and Electromechanics. Circ. Res. 2012, 108, 113–128. [Google Scholar]
  17. Vinciguerra, V.; Ambra, E.; Maddiona, L.; Romeo, M.; Mazzillo, M.; Rundo, F.; Fallica, G.; Pompeo, F.; Chiarelli, A.M.; Zappasodi, F.; et al. PPG/ECG Multisite Combo System Based on SiPM Technology. In CNS 2018. Lecture Notes in Electrical Engineering; Springer: Cham, Switzerland, 2018; p. 53. [Google Scholar]
  18. Quiroz-Juárez, M.A.; Jiménez-Ramírez, O.; Vázquez-Medina, R.; Brena-Medina, V.F.; Aragón, J.L.; Barrio, R.A. Generation of ECG Signals from a Reaction-Diffusion Model Spatially Discretized. Sci. Rep. 2019, 9. [Google Scholar] [CrossRef] [Green Version]
  19. Sornmo, L.; Laguna, P. Electrocardiogram (ECG) Signal Processing. In Wiley Encyclopedia of Biomedical Engineering; Wiley: New York, NY, USA, 2006. [Google Scholar]
  20. Rai, H.M.; Trivedi, A.; Shukla, S. ECG Signal Processing for Abnormalities Detection Using Multi-Resolution Wavelet Transform and Artificial Neural Newtork Classifier. Measurement 2013, 46, 3238–3246. [Google Scholar]
  21. Gacek, A. An Introduction to ECG Signal Processing and Analysis; Gacek, A., Pedrycz, W., Eds.; Springer: London, UK, 2012; pp. 21–46. [Google Scholar]
  22. Jaros, R.; Martinek, R.; Kahankova, R. Non-Adaptive Methods for Fetal ECG Signal Processing: A Review and Appraisal. Sensors 2018, 18, 3648. [Google Scholar]
  23. Sameni, R.; Clifford, G.D. A Review of Fetal ECG Signal Processing: Issues and Promising Directions. Open Spacing Electrophysiol. Ther. J. 2010, 1, 4–20. [Google Scholar]
  24. McSharry, P.E.; Clifford, G.D.; Tarassenk, L.; Smith, L.A. A Dynamical Model for Generating Synthetic Electrocardiogram Signals. IEEE Trans. Biomed. Eng. 2003, 50, 289–294. [Google Scholar]
  25. Quarteroni, A.; Sacco, R.; Saleri, F. Numerical Mathematics; Springer: Berlin/Heidelberg, Germany, 2007. [Google Scholar]
  26. Sayadi, O.; Mohammad, B.S.; Gari, D.C. Synthetic ECG generation and Bayesian filtering using a Gaussian wave-based dynamical model. Physiol. Meas. 2010, 31, 1309. [Google Scholar]
  27. Jackson, S.; Cogswell, M.E.; Zhao, L.; Terry, A.L.; Wang, C.Y.; Wright, J.; King, S.M.C.; Bowman, B.; Chen, T.-C.; Merritt, R.; et al. Association Between Urinary Sodium and Potassium Excretion and Blood Pressure Among Adults ioon the United States, National Health and Nutrition Examination Survey. Circulation 2018, 137, 237–246. [Google Scholar]
  28. Aaron, K.J.; Sanders, P.W. Role of Dietary Salt and Potassium Intake in Cardiovascular Health and Disease: A Review of the Evidence. Mayo Clin. Proc. 2014, 11, 9. [Google Scholar]
  29. Vicente, J.; Zusterzeel, R.; Johannesen, L.; Ochoa-Jimenez, R.; Mason, J.W.; Sanabria, C.; Kemp, S.; Sager, P.T.; Patel, V.; Matta, M.K.; et al. Assessment of Multi-Ion Channel Block in a Phase I Randomized Study Design: Results of the CiPA Phase I ECG Biomarker Validation Study. Clin. Pharmacol. Ther. 2018, 105, 111–129. [Google Scholar]
  30. Garcia, M.D.; Tomas, B. 12 Lead ECG: The Art of Interpretation; Pearson: London, UK, 2014. [Google Scholar]
  31. Goldberger, A.L.; Goldberger, Z.D.; Shvilkin, A. Goldberger’s Clinical Electrocardiography; Elsevier: Amsterdam, The Netherlands, 2017. [Google Scholar]
  32. Versaci, M. Soft Computing Approach to Predict Intracranial Pressure Values. Am. J. Appl. Sci. 2014, 11, 844–850. [Google Scholar]
  33. Versaci, M.; Morabito, F.C. Fuzzy Time Series Approach for Disruption Prediction in Tokamak Reactors. IEEE Trans. Magn. 2003, 39, 1503–1506. [Google Scholar]
  34. Cacciola, M.; La Foresta, F.; Morabito, F.C.; Versaci, M. Advanced Use of Soft Computing and Eddy Current Test to Evaluate Mechanical Integrity of Metallic Plates. NdT E Int. 2007, 40, 357–362. [Google Scholar]
  35. Angiulli, G.; Jannelli, A.; Morabito, F.C.; Versaci, M. Reconstructing the Membrane Detection of a 1D Electrostatic-Driven MEMS Device by the Shooting Method: Convergence Analysis and Ghost Solutions Identification. Comput. Appl. Math. 2018, 37, 4484–4498. [Google Scholar]
  36. Schmidt, R.F.; Thews, G. Human Physiology; Springer: Berlin/Heidelberg, Germany, 1989. [Google Scholar]
  37. Agarwal, A.M. Trigonometry; Arihant: Delhi, India, 2018. [Google Scholar]
  38. Barreira, L.; Valls, C. Dynamical Systems; Springer: London, UK, 2013. [Google Scholar]
  39. Cacciola, M.; Calcagno, S.; Morabito, F.C.; Versaci, M. Swarm Optimization for Imaging of Corrosion by Impedance Measurements in Eddy Current Tests. IEEE Trans. Magn. 2007, 43, 1853–1856. [Google Scholar]
  40. Cacciola, M.; Pellicanó, D.; Megali, G.; Lay-Ekuakille, A.; Versaci, M.; Morabito, F.C. Aspects About Air Pollution Prediction on Urban Environment. In Proceedings of the 4th IMEKO TC19 Symposium on Environmental Instrumentation and Measurements 2013: Protection Environment, Climate Changes and Pollution Control, Lecce, Italy, 3–4 June 2013; pp. 15–20. [Google Scholar]
  41. Coletti, G.; Bouchon-Meunier, B. Fuzzy Similarity Measures and Measurement Theory. In Proceedings of the 2019 IEEE International Conference on Fuzzy Systems (FUZZ-IEEE), New Orleans, LA, USA, 23–26 June 2019; pp. 1–8. [Google Scholar]
  42. Moody, G.B.; Mark, R.G. The impact of the MIT-BIH Arrhythmia Database. IEEE Eng. Med. Biol. 2001, 20, 45–50. [Google Scholar]
  43. Goldberger, A.L.; Amaral, L.A.N.; Glass, L.; Hausdorff, J.M.; Ivanov, P.C.; Mark, R.G.; Mietus, J.E.; Moody, G.B.; Peng, C.-K.; Stanley, H.E. PhysioBank, PhysioToolkit, and PhysioNet: Components of a new research resource for complex physiologic signals. Circulation 2000, 101, e215–e220. [Google Scholar]
  44. Iserles, A. A First Course in the Numerical Analysis of Differential Equations; Cambridge University Press: Cambridge, UK, 2009. [Google Scholar]
Figure 1. (a) Sagittal plane ( x z plane), frontal plane ( y z plane) and transverse plane ( x y plane) (public domain image). (b) Typical E C G signals: continuous line represents a normal E C G . Dotted and dashed lines represent pathological E C G s.
Figure 1. (a) Sagittal plane ( x z plane), frontal plane ( y z plane) and transverse plane ( x y plane) (public domain image). (b) Typical E C G signals: continuous line represents a normal E C G . Dotted and dashed lines represent pathological E C G s.
Computation 08 00092 g001
Figure 2. (a) A typical vector cardiogram (VCG) loop for a normal cardiac cycle. (b) A typical VCG loop characterized by a nodal rhythm disease.
Figure 2. (a) A typical vector cardiogram (VCG) loop for a normal cardiac cycle. (b) A typical VCG loop characterized by a nodal rhythm disease.
Computation 08 00092 g002
Figure 3. E C G reconstruction by (a) System (14) (case (a)) and (b) System (15) (case (b)).
Figure 3. E C G reconstruction by (a) System (14) (case (a)) and (b) System (15) (case (b)).
Computation 08 00092 g003
Figure 4. (a) E C G reconstruction by system (16) (case (c)). (b) Typical E C G obtained from a patient not suffering from any cardiac disease.
Figure 4. (a) E C G reconstruction by system (16) (case (c)). (b) Typical E C G obtained from a patient not suffering from any cardiac disease.
Computation 08 00092 g004
Figure 5. Case (a) of the proposed methodology: normal sinus rhythm.
Figure 5. Case (a) of the proposed methodology: normal sinus rhythm.
Computation 08 00092 g005
Figure 6. Case (a) of the proposed methodology: (a) atrial fibrillation and (b) typical E C G with atrial fibrillation.
Figure 6. Case (a) of the proposed methodology: (a) atrial fibrillation and (b) typical E C G with atrial fibrillation.
Computation 08 00092 g006
Figure 7. Case (a) of the proposed methodology: (a) atrial flutter and (b) typical E C G with atrial flutter.
Figure 7. Case (a) of the proposed methodology: (a) atrial flutter and (b) typical E C G with atrial flutter.
Computation 08 00092 g007
Figure 8. Case (a) of the proposed methodology: (a) premature ventricle contraction and (b) typical E C G with premature ventricle contraction.
Figure 8. Case (a) of the proposed methodology: (a) premature ventricle contraction and (b) typical E C G with premature ventricle contraction.
Computation 08 00092 g008
Figure 9. Case (b) of the proposed methodology: (a) normal sinus rhythm (b) atrial fibrillation (c) atrial flutter (d) premature ventricle contraction.
Figure 9. Case (b) of the proposed methodology: (a) normal sinus rhythm (b) atrial fibrillation (c) atrial flutter (d) premature ventricle contraction.
Computation 08 00092 g009
Figure 10. Case (c) of the proposed methodology: (a) normal sinus rhythm (b) atrial fibrillation (c) atrial flutter (d) premature ventricle contraction.
Figure 10. Case (c) of the proposed methodology: (a) normal sinus rhythm (b) atrial fibrillation (c) atrial flutter (d) premature ventricle contraction.
Computation 08 00092 g010
Figure 11. Histograms for typical E C G s characterized by (a) normal sinus rhythm, (b) atrial fibrillation, (c) atrial flutter, (d) premature ventricle contraction.
Figure 11. Histograms for typical E C G s characterized by (a) normal sinus rhythm, (b) atrial fibrillation, (c) atrial flutter, (d) premature ventricle contraction.
Computation 08 00092 g011
Figure 12. Evaluation of linear index (LI) and fuzzy entropy (FE) for each category of E C G s.
Figure 12. Evaluation of linear index (LI) and fuzzy entropy (FE) for each category of E C G s.
Computation 08 00092 g012
Table 1. Characteristics of an E C G .
Table 1. Characteristics of an E C G .
SectionDuration (s)Amplitude (mV)Meaning
P Wave0.07–0.120.2–0.4Atrial Depolarization
QRS Complex0.06–0.101–2Depolarization of the intraventricular septum (Q)
and ventricles (RS)
T Wave0.18–0.200.4–0.5Ventricular re-polarization
U Wave0.08-Re-polarization of the Purkinje System
P-R Interval0.1–0.20-Atrioventricular conduction time
S-T Interval0.30-Duration of ventricular re-polarization
Q-T Interval0.40-Duration of electrical ventricular systole
R-R Interval0.8–0.9-Duration of cardiac cycle
Table 2. Experimental parameters for the McSharry dynamical model.
Table 2. Experimental parameters for the McSharry dynamical model.
iPQRST
time (s)−0.250−0.02500.0250.250
θ i (rad) π / 3 π / 12 0 π / 12 π / 2
a i 1.25−5.0030.00−8.001.00
b i 0.250.10.10.10.5
Table 3. Some synthetic diseases by the modified heart dipole model (MHDM).
Table 3. Some synthetic diseases by the modified heart dipole model (MHDM).
Model ParametersPQRST Complex AlterationCorrelated Disease
k 1 = 1 , k 2 = 1 , ϕ = 0 -normal sinus rhythm
k 1 = 0.3 , k 2 = 0.1 , ϕ = π / 8 rapid ventricularatrial fibrillation
rate normal QRS
k 1 = 0.5 , k 2 = 0.1 , ϕ = π / 4 isoelectric interval TPatrial flutter
k 1 = 0.2 , k 2 = 0.3 , ϕ = 2 3 π abnormal QRSpremature
P wave not associatedventricle contraction
Table 4. F S 1 , F S 2 , F S 3 and F S 4 achieved for E C G s related to Case (a).
Table 4. F S 1 , F S 2 , F S 3 and F S 4 achieved for E C G s related to Case (a).
Figure 3aFigure 3bFigure 4aFigure 5Figure 6aFigure 7aFigure 8a
no pathologies
F S 1 0.960.990.990.970.120.090.15
F S 2 0.940.9710.9890.980.130.110.14
F S 3 0.970.9790.9870.990.1360.120.13
F S 4 0.930.730.9880.980.1340.1290.129
atrial fibrillation
F S 1 0.140.140.1350.010.930.160.09
F S 2 0.160.190.1390.0260.9320.1510.09
F S 3 0.1650.1960.1380.150.9210.1460.09
F S 4 0.170.090.1380.1650.9090.1440.08
atrial flutter
F S 1 0.180.160.140.1640.140.980.07
F S 2 0.160.170.210.1630.1420.990.11
F S 3 0.1560.1880.160.1620.1430.9850.12
F S 4 0.1640.150.040.1610.1450.9870.111
premature ventricle contraction
F S 1 0.20.150.1590.040.1390.040.99
F S 2 0.1750.140.1650.040.1320.0390.97
F S 3 0.1610.1420.1450.0480.1380.0410.98
F S 4 0.1710.030.1610.1450.1370.040.98
Table 5. F S 1 , F S 2 , F S 3 and F S 4 achieved for E C G s related to Case (b).
Table 5. F S 1 , F S 2 , F S 3 and F S 4 achieved for E C G s related to Case (b).
Figure 9aFigure 9bFigure 9cFigure 9d
no pathologies
F S 1 0.980.130.120.16
F S 2 0.980.1440.110.14
F S 3 0.960.1340.130.12
F S 4 0.970.1310.140.11
atrial fibrillation
F S 1 0.130.960.170.12
F S 2 0.1330.9550.1320.11
F S 3 0.1350.9510.1440.15
F S 4 0.420.920.1570.11
atrial flutter
F S 1 0.210.130.990.12
F S 2 0.220.110.990.11
F S 3 0.210.1230.990.13
F S 4 0.210.140.970.11
premature ventricle contraction
F S 1 0.210.160.1690.99
F S 2 0.1180.150.1660.99
F S 3 0.1430.1440.1570.97
F S 4 0.1240.090.1810.968
Table 6. F S 1 , F S 2 , F S 3 and F S 4 achieved for E C G s related to Case (c).
Table 6. F S 1 , F S 2 , F S 3 and F S 4 achieved for E C G s related to Case (c).
Figure 10aFigure 10bFigure 10cFigure 10d
no pathologies
F S 1 0.990.130.120.19
F S 2 0.980.140.130.17
F S 3 0.990.1550.1320.15
F S 4 0.980.1590.1310.137
atrial fibrillation
F S 1 0.120.970.170.11
F S 2 0.1220.9520.170.1
F S 3 0.1390.930.1620.12
F S 4 0.1690.940.1390.13
atrial flutter
F S 1 0.160.150.990.11
F S 2 0.1560.150.990.21
F S 3 0.1460.1770.9840.22
F S 4 0.1540.120.970.211
premature ventricle contraction
F S 1 0.210.160.1620.98
F S 2 0.180.1440.1550.975
F S 3 0.1690.1530.1510.99
F S 4 0.1770.110.1190.99
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Versaci, M.; Angiulli, G.; La Foresta, F. A Modified Heart Dipole Model for the Generation of Pathological ECG Signals. Computation 2020, 8, 92. https://0-doi-org.brum.beds.ac.uk/10.3390/computation8040092

AMA Style

Versaci M, Angiulli G, La Foresta F. A Modified Heart Dipole Model for the Generation of Pathological ECG Signals. Computation. 2020; 8(4):92. https://0-doi-org.brum.beds.ac.uk/10.3390/computation8040092

Chicago/Turabian Style

Versaci, Mario, Giovanni Angiulli, and Fabio La Foresta. 2020. "A Modified Heart Dipole Model for the Generation of Pathological ECG Signals" Computation 8, no. 4: 92. https://0-doi-org.brum.beds.ac.uk/10.3390/computation8040092

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop