Relationship between Mathematical Parameters of Modified Van der Pol Oscillator Model and ECG Morphological Features

The mathematical model describes the electrical and mechanical activity of the cardiac conduction system thought set of differential equations. By changing the value of parameters included in these equations, it is possible to change the amplitude and the period of ECG waves. Although this model is a powerful tool for modeling the electrical activity of the heart, its use is often limited to those familiar with the differential equations that describe the system. The purpose of this work is to provide a system that allows generating an ECG signal using Ryzhii model without knowing the details of differential equations. First, we provide the relationships between the ECG wave features and the model parameters; then we generalize them through fitting neural networks. Finally, putting in series fitting neural network and heart model, we provide a system that allows generating a synthetic signal by setting as input only the morphological ECG feature. We computed numerical simulation in Simulink environment and implemented the fitting neural networks in Matlab. Results show that non-linear trends characterize the correlation functions between ECG morphological features and model parameters and that the fitting neural networks can generalized this trend by providing the model parameters given in input the respectively ECG feature. Keywords— heart model; Van der Pol; FitzHugh-Nagumo; relaxation oscillator; fitting neural network; classification.


I. INTRODUCTION
Over the years, the dynamic of the heartbeat was analyzed through both mathematical model and time series analysis. In recent years, further help has come from the spread of various artificial intelligence techniques [1]- [15], many of which are based on neural networks that use the latest technologies [16]- [21]. Despite the wide use of the time series for the study of ECG signal, the complexity, nonlinearity, and nonstationarity of the cardiovascular system make common the use of nonlinear signal analysis for the modeling of heart activity [22]. Mathematical modeling is becoming a powerful tool for improving our understanding of electrical activation in the heart, in particular, the generation and propagation of both normal and abnormal rhythms [23], [24]. Van der Pol and Van der Mark for the first time described and modeled the dynamic behavior of the heart using nonlinear relaxation oscillators [25]. Their work has become the start point of lots of studies on heart behavior modeling. Afterward, a lot of researchers worked on an extended and modified version of Van der Pol model [26]- [28]. In the last years, the characterization of electrical muscle responses was proposed by Ryzhii [29] using a quiescent excitable FitzHugh-Nagumo-type (FHN) oscillators. Successively, Ryzhii improved their model including a description for depolarization and repolarization waves of the atrias and ventricles utilizing different modified FHN systems for each ECG waves [30]. The advanced model was able to generate normal ECG signal and several well-known rhythm disorders. In particular, it reproduces sinus tachycardia, sinus bradycardia, complete SA-AV block, and complete AV -HP block [31]. This model describes the electrical and mechanical activity of the cardiac system thought the set of differential equations. The parameters included in these equations make it possible to describe different cardiac behavior. By modifying the model parameters, it is possible to change the amplitude and the time duration of the ECG waves, reproducing both physiological and pathological behavior. Although this mathematical model allows describing the morphology of the electrical activity, the relationships between model parameters and morphological characteristics of ECG waves are not provided. In particular, it is not possible to set the model parameters in intuitive mode knowing the characteristics of the waves composing ECG signal.
The scope of this work is to provide a system that allows generating synthetic signals by setting in the input only the morphological features of the ECG signal, without knowing the mathematical parameters. In this way, the knowledge and meaning of differential equations that describe the system are not necessary. First, we researched the relationships that link the morphological characteristics of the signal to the variance of the model parameters. In particular, we implemented the mathematical model proposed by Ryzhii [30] in the Simulink environment. Then, we generalized the characteristic relationships using fitting neural networks. We implemented the neural networks in the Matlab environment and, finally, we used the results of the fitting neural networks as the input of the heart model. In this way, we provide a system that, without knowing the meaning of each parameter, allows to generate synthetic ECG waveforms.
The paper is organized as follows: in Sec.2 the principal morphological characteristics of ECG waves are described. In Sec.3 mathematical model to reproduce synthetic ECG signal is explained. In Sect.4 material and method are provided. In Sec.5 results and discussion are shown, and finally, in Sect.6 conclusions are discussed.

A. Electrocardiographic signal
Depolarization and repolarization of cardiac cells can be recorded from the body surface by the ECG. Figure 1 shows a normal ECG signal. The P-wave represents atrial depolarization, the QRS complex represents ventricular depolarization, and the isoelectric period between these two waves represents the AV node and His bundle activation. Finally, the T-wave represents ventricular repolarization [32]. Another important data component is the temporal relationship of these waves. The PR-interval is measured from the beginning of the P -wave to the beginning of the QRS complex and indicates the time interval between the onset of atrial activation and the onset of ventricular activation. The normal PR-interval in adults is less than 0.20 s. The QT -interval is measured from the beginning of the QRS to the terminal portion of the T -wave. It represents a global measure of the ventricular plateau phase duration. Tab.1 shows the characteristic values of ECG waves and their temporal relationship. where x i (t) and y i (t) correspond to the cardiac action potential and transmembrane currents, a i (x i − u ij )(x i − u ik ) are the damping terms, f i x i (x i − d ij )(x i − e ik ) are the harmonic force terms, a i > 0 and u ij are the parameters representing non-linear damping force, fi are the parameters related to the intrinsic frequency of oscillator, the coupling coefficients K SA−AV and K AV−HP represent unidirectional coupling between the pacemakers. The synchronization level between three oscillators depends on the coupling coefficients K ( [33], [34]). Choosing proper coupling coefficients K SA−AV and K AV−HP is possible not to use time delays in coupling terms in order to consider the different frequency rate of the three oscillators.
The parameters a i , u ij , f i , d i , e i , were choosen in order to obtain intrinsic oscillation rates of 70 bpm, 50 bpm, 35 bpm and 35 bpm for uncoupled SA, AV, RB and LB, respectively. In particular, we used the following experimental parameters: a 1 = 40, a 2 = 50, a 3 = 50, a 1 = 40, u 11   The activation currents I i represent the coupling between the SA and AT muscle and between HP pacemaker and VN muscles. Concerning the Ryzhii model [30], we adjusted the activation currents for the QRS complex and T wave to consider HP oscillator composed by the RB and LB oscillators: where K ATDe = 9 * 10−5, K ATRe = 3 * 10−5, K VNDe = 9 * 10−5 and K VNRe = 9 * 10−5. Finally, we obtained the synthetic ECG signal as a composition of AT and VN waveforms as following: ECG = z 0 + z 1 − z 2 + k R z 3 + z 4 (12) where z 0 = 0.2 provides the adjustment of baseline and k R is a multiplicative coefficient to modulate the amplitude of R peak.

II. MATERIAL AND METHOD
The method is composed of the following steps: • To answer this question, we implemented the mathematical model of Ryzhii [30] in the Simulink environment. In the first step, we looked which mathematical constant or variable of the model allows modifying the amplitude and the temporal intervals of the ECG waves in according to the Table.1. In this way, for each morphological feature FECG of the ECG signal corresponds a specific parameter Pm of the model. In the second step, we defined the operating ranges of the model parameters. To evaluate the operating intervals of every P m , we performed a series of simulations in Simulink. In particular, by knowing the specific interval of each ECG wave (Table.1), we have changed the operative value of the matched Pm until it returns the ECG signal with the desired morphological characteristics. In the third step, we have found the correlation function between FECG and the corresponding P m (Fig.4).
We generated the synthetic ECG signal by executing the cardiac model and modifying the parameters within the respective operating ranges. Then, we applied a modified version of the Pan & Tompkins algorithm [35]- [37] to identify the P-wave, QRS-complex, and T-wave of the synthetic signals. Finally, we chart each model parameter according to its corresponding ECG feature.

B. Training, Test, and Validation of the Fitting Neural Networks
In the second phase, we generalized the relationships F ECG -P m defined in the previous paragraph by a neural network. Neural networks provide a tool for the fast solution of repetitive non-linear curve fitting problems. In particular, function fitting is the process of training a neural network on a set of inputs to produce an associated set of target outputs. Fig.5 shows the conceptual scheme.  First, we constructed the network with the desired hidden layers and the training algorithm; then we trained it using the ECG features as a set of training data and the model parameters P m as the training label. Once the neural network adapts to the data, it forms a generalization of the inputoutput relationship. In other words, given the ECG feature as input, the fitting neural network can provide the respectively model parameter as output. Table 2 shows the set-up of fitting neural networks. For each ECG feature, a custom neural network was created.
In particular, the chosen transfer function depending on the relation curve between the ECG wave feature and model parameter. For the training of P -wave and QRS -complex, we used Levenberg-Marquardt function. It is a training function that updates weight and bias values according to Levenberg-Marquardt optimization. Whereas, for T -wave, PR-interval, and RR-interval we used Bayesian Regularization function. It is a training function that updates the weight and bias values according to Levenberg-Marquardt optimization. It minimizes a combination of squared errors and weights and then determines the correct combination to produce a network that generalizes well. The process is called Bayesian regularization.

C. Parameters Obtained through the Neural Network:
Application to the Model Once defined and modeled, with an appropriate neural network, the relationships between F ECG -P m , we have verified the results using the system shown in Fig.6. For example, we want to simulate an ECG signal characterized by a particular amplitude of the P-wave. This feature represents the input of the fitting neural network that returns as output the model parameter K ATde . The parameter K ATde in input to the heart model generates the synthetic ECG signal. To validate that the output signal is correct, the P -wave amplitude (Ḟ ECG ) was extracted using a modified version of the Pan & Tompkins algorithm. Finally, the amplitude of P -wave in the input (F ECG ) and the amplitude of P-wave of the synthetic ECG signal in output (Ḟ ECG ) were compared. The validation test was repeated for all ECG features in their range of applicability. To control the amplitude values of ECG waves (P, T, Rwaves) we worked on the intensity of activation currents I ATDe and I VNRe . These currents represent the coupling between the SA pacemaker and AT muscles in the depolarization phase, and between HP pacemaker and VN muscles in the repolarization phase. To govern duration time and amplitude range of the P and the T-waves we set the constants K ATDe and K VNRe . In fact, they allow modulating the amplitude of transmembrane ionic currents I ATDe and I VNRe which represent the couplings between the SA ad HP pacemakers and AT and VN muscles.  [20,30] To control the time interval between waves (PR, RR, STinterval) we used the coupling coefficients. They provide proper synchronization behavior of pacemakers in a wide range of heart rhythms. In particular, PR-interval represents the time between the atria depolarization, that generates the P-wave, and ventricles depolarization, that generates the QRS complex. To manage this interval range we worked on the K SA−AV coupling coefficient because it provides synchronization in the SA and AV oscillators interaction. It is important to notice that not for all morphological features of ECG signal exist the corresponding model parameter. For example, with this heart model, it is not possible to change the duration of ST and QT-intervals. For this reason, they do not appear in Tab.2 Fig.7 shows the correlation graphs of the related behavior between morphological ECG features and model parameters. From these graphs, it is possible to notice that only k R and R-wave amplitude show linear functional behavior. In this case, increasing k R in the model equation proportionally increases R-amplitude. Whereas, for the other parameters, there are non-linear functional trends. In these cases, it is not possible, by knowing the ECG wave amplitudes or periods, to set the corresponding model parameters. Fig. 7 -Correlation graphs that describe the relation behavior between morphological ECG features and model parameters By implementing the fitting neural networks, we showed how they can predict the model parameters given in entry the ECG wave features. Fig.8 shows the fitting neural networks results. The blue dots represent the training targets, the data presented to the network during the training phase, and the network is adjusted according to its error. The green dots represent the validation targets; the data set used to measure network generalization and to halt training when generalization stops improving. The red dots represent the test targets; these do not affect training and provide an independent measure of network performance during and after training. Finally, the black line represents the fitting curve searched. The blue, green and red crosses represent respectively training, validation and test outputs that are the result of the network for specific ECG features in the input. The error is the difference between targets and outputs. It is possible to notice that in a) and b) cases the error about zero. In the case of c) it is about zero, besides in the part in which the trend changes. It can depend on the training data size. In d) and e) cases the error is more than zero, but the networks can generalize the relationships FECG-Pm.

IV. CONCLUSION
In this paper, we provide a system that allows generating an ECG signal using Ryzhii model [30] without knowing the details of differential equations. Although the mathematical model is a powerful tool for modeling the electrical and mechanical activity of the heart, its use is limited to those familiar with the differential equations that describe the system. We provide the correlation functions that link the ECG features and model parameters and generalized them with fitting neural networks. Finally, we used the output of the neural network as the input of the mathematical heart model. In this way, we generate ECG signal only setting the amplitude and the time duration of waves composing ECG signal. Results show that the relationships between ECG features and model parameters present a non-linear trend and that the fitting neural networks can generalize them with an acceptable error. Considering the enormous spread of IoT devices [38], this type of system, as well as other complex systems useful for health care [39], will necessarily have to be integrated into low energy consumption portable embedded devices [40]-42] which will allow monitoring the quality of patients' health more efficiently.