On the Heart Rhythm Analysis by an Intelligent Control Method Using Artificial Neural Networks: Analytical Study and Electronic Simulation

Mathematical modelling

As is well known, two pumps - the left and right ventricles - are essential parts of the heart. They are contracted by the electrical impulse produced by the heart through its electrical system, known as the electrical conduction system of the heart. Its role is to ensure a regular heartbeat. For good functioning of the electrical conduction system, the amount of blood needed to oxygenate the lungs is pumped by the right ventricle, and the amount of blood needed to oxygenate the whole body is pumped by the left ventricle. This analysis shows that the greatest mechanical work is done by the left ventricle. All this is due to the regularity of the cardiac conduction system It is composed of the sinus node, the atrioventricular node and Purkinje's His bundles. The electrical impulse which causes the contraction of the atria and ventricles is initiated by the sinus node (the heart's natural pacemaker). After receiving information from the brain, the SA node adapts the heart rate in accordance with the body's energy requirements. In its normal state, heart rate varies between 1 and 1.67 Hz. The electrical impulse initiated by the SA node enables the blood to reach the ventricles after having passed through the two auricles by contracting them. The atrioventricular node in turn receives this electricity, filters it, slowing it down, if necessary, in cases of rhythm disorders, and then directs it through Purkinje's His bundles to the ventricles, also enabling the ventricles to contract. This is how electricity travels through the heart. Each heartbeat forming the electrocardiogram (ECG) signal is caused by each electrical current pathway from the sinus node to the ventricles.

Recent work has shown that heartbeat modeling is a crucial process for the automatic identification of characteristic waves [10,11]. The objective is to propose the simplest, most compact mathematical model capable of qualitatively and quantitatively simulating the shape of the five main waves (P,Q,R,S and T waves) constitutive of heartbeats. Indeed, the most natural representation of the waves would be to describe the signal by its amplitude at each instant; this representation would therefore be a vector in a space whose dimension would be equal to a few hundred. From these analyses, cardiac physiology has been modeled as a system of three oscillators subjected to external stimuli and coupled together by delayed connections [36]. These stimuli have the effect of increasing the dimension of the reduced-order model system and also represent spatio-temporal stimuli.

In this subsection, the mathematical model is that proposed by Fonkou and Savi [36]. This is a system of three nonlinear oscillators coupled by delayed Duffing-type connections incorporating cubic non-linear terms. These three oscillators consist of the SA node, the AV node and the Purkinje fibers linked by Duffing-type connections with delay and non-delay. For this work, in addition to the external stimuli to which each oscillator is subjected, the SA oscillator is subjected to the controller u as presented by its conceptual model of general functioning given in figure 2. The mathematical model describing them is given by equation (1).

Figure 2: Conceptual model of general functioning of the heart with controller.

$\begin{array}{l}\begin{array}{rr}\hfill {\dot{x}}_1=x_2& \hfill \\ \hfill {\dot{x}}_2=& \hfill F_{SA}(t)-{\alpha }_{SA}{x}_2\left(x_1-v_{S{A}_1}\right)\left(x_2-v_{S{A}_1}\right)-\omega \frac{x_1\left(x_1-d_{SA}\right)\left(x_1-e_{SA}\right)}{d_{SA}{e}_{SA}}-k_{AV-SA}\left(x_1-x_1^3\right)-k_{HP-SA}\left(x_1-x_1^3\right)\\ \hfill & \hfill +k_{AV-SA}^{\tau }\left(x_3^{{\tau }_{AV-SA}}+{\left(x_3^{{\tau }_{AV-SA}}\right)}^3\right)+k_{HP-SA}^{\tau }\left(x_5^{{\tau }_{HP-SA}}+{\left(x_5^{{\tau }_{HP-SA}}\right)}^3\right)+u\end{array}\text{ (1)}\\ {\dot{x}}_3=x_4\\ \begin{array}{r}\hfill \\ \hfill {\dot{x}}_4=F_{AV}(t)-{\alpha }_{AV}{x}_4\left(x_3-v_{A{V}_1}\right)\left(x_3-v_{A{V}_1}\right)-\frac{x_3\left(x_3+d_{AV}\right)\left(x_3+e_{AV}\right)}{d_{AV}{e}_{AV}}-k_{SA-AV}\left(x_3-x_3^3\right)-k_{HP-AV}\left(x_3-x_3^3\right)\\ \hfill +k_{SA-AV}^{\tau }\left(x_1^{{\tau }_{SA-AV}}+{\left(x_1^{{\tau }_{SA-AV}}\right)}^3\right)+k_{HP-SA}^{\tau }\left(x_5^{{\tau }_{HP-AV}}+{\left(x_5^{{\tau }_{HP-AV}}\right)}^3\right)\end{array}\\ {\dot{x}}_5=x_6\\ \begin{array}{r}\hfill \\ \hfill {\dot{x}}_6=F_{HP}(t)-{\alpha }_{HP}{x}_6\left(x_5-v_{H{P}_1}\right)\left(x_5-v_{H{P}_1}\right)-\frac{x_5\left(x_5-d_{HP}\right)\left(x_5-e_{HP}\right)}{d_{HP}{e}_{HP}}-k_{SA-HP}\left(x_5-x_5^3\right)-k_{AV-HP}\left(x_5-x_5^3\right)\\ \hfill +k_{SA-HP}^{\tau }\left(x_1^{{\tau }_{SA-HP}}+{\left(x_1^{{\tau }_{SA-HP}}\right)}^3\right)+k_{AV-HP}^{\tau }\left(x_3^{{\tau }_{AV-HP}}+{\left(x_3^{{\tau }_{AV-HP}}\right)}^3\right)\end{array}\end{array}$

Where u is the controller. The SA, AV and HP nodes are described respectively by the pairs $(x_1,x_2),(x_3,x_4)and(x_5,x_6)$ . Coefficients between two nodes m and $n(m\ne n)$  without delay and with delay are given by $k_{m-n}$  and $k_{m-n}^{\tau }$ (m and n can express the SA, AV or HP node as well). The linear and cubic terms with delay are: $x_j^{{\tau }_{m-n}}=x_j\left(t-{\tau }_{m-n}\right)$  and $x_j$ , and those without delay are $x_j$  and $x_j{ }^3(j=1,2,\dots ,6).F_n(t)={\rho }_n\text{sin}\left({\omega }_{n}t\right)$  is the external excitation making the system temporarily dependent and increasing its dimension. $\omega$ is the transient frequency disturbing the SA node. It is a function of time and is described by

$\begin{array}{c}\omega =\omega (t)=\frac{{\omega }_a+{\omega }_b}{2}+\frac{{\omega }_a-{\omega }_b}{2}tanh\left[\gamma \left(t-t_0\right)\right]\end{array}\text{ (2)}$

$\frac{{\omega }_a-{\omega }_b}{2}$ is the transitional phase between ${\omega }_a$ and ${\omega }_b,t_0$ the transition time. In this work its influence will be neglected, it will take the values ${\omega }_a={\omega }_b=1$

Using a combination of three nodes (SA, AV and HP), the mathematical expression of the ECG signal is given by

$\begin{array}{c}X=ECG=X_1+X_3+X_5\end{array}\text{ (3)}$

With

$\begin{array}{c}\{\begin{array}{l}{X}_1=\frac{a_0}{3}+a_1{x}_1\hfill \\ X_3=\frac{a_0}{3}+a_2{x}_3\hfill \\ X_5=\frac{a_0}{3}+a_3{x}_5\hfill \end{array}\end{array}\text{ (4)}$

Where $a_0,a_1,a_2$ and $a_3$ are constant parameters. In addition, it is possible to write:

$\begin{array}{c}\dot{x}=\frac{d(ECG)}{dt}=a_1{x}_2+a_2{x}_4+a_3{x}_6\end{array}\text{ (5)}$

Considering the equations of the system (1) adimensioned, the parameters to conform with them are given by: $\overline{t}[s]:\overline{t}=a_{\tau }t$ where $\left[a_{\tau }\right]=s$ . It can also be considered as the ratio between the experimental value and the numerical value $:\frac{\text{mean}\left(R{R}_{exp }\right)}{\text{mean}\left(R{R}_{num }\right)}$ . Here the external excitations are considered as harmonics: $F_{SA}(t)={\rho }_{SA}\text{sin}\left({\omega }_{SA}t\right),F_{AV}(t)={\rho }_{AV}\text{sin}\left({\omega }_{AV}t\right)$ and $F_{HP}(t)={\rho }_{HP}\text{sin}\left({\omega }_{HP}t\right)$ . All simulations consider $a_{\tau }=1.0mV,a_1=2.0mV,a_2=0.01mV$ and $a_3=0.15mV$ , for all the simulations, we will take as initial conditions : $x_{1,3,5}=0.01$ and ${\dot{x}}_{1,3,5}=0.04$ .

Note that the model parameters are estimated from the same real ECG data and can be either deterministic or non-deterministic. In deterministic situations, each rhythm is considered to be associated with a set of parameters, and a change in that rhythm is attributed to a change in those parameters. In other words, a pathological rhythm can be understood as resulting from a change in these parameters. In contrast, nondeterministic situations involve random variations in the parameters that can lead to pathological rhythms. t should also be noted that the constants of the coupled Duffing model are physiologically grounded: the natural frequencies ${\omega }_i$ reproduce the intrinsic automaticity of the SA node ( 1 Hz ), the AV node ( 0.7 Hz ), and the His-Purkinje system (0.4" " Hz) [58]. The delays τ_i correspond to the PR interval of 0.1 s and the QRS interval of 0.08 s measured on a 12 -lead ECG [59]. The coupling coefficients $k_{i-j}$ model the conductances of connexin-mediated junctions, the reduction of which explains AV block [60]. The Duffing nonlinearities α_i and $a_i(a_0,a_1,a_2$ and $a_3)$ are justified by the nonlinear current-voltage relationship of the voltage-dependent ${\text{Na}}^{+},{\text{Ca}}^{+},{\text{K}}^{+}$ channels [28]. The synthetic ECG generated is therefore the weighted sum of the three oscillators, in accordance with the principle of superposition of cardiac dipoles [61].

It should also be noted that changing the value of a constant alters the nature of the ECG signal in accordance with the clinical case. For example, a decrease in one of the $k_{i-j}$ values can lead to a prolongation of the PR interval, which would imply the onset of AV block. An increase or decrease in ${\alpha }_i$ can alter the T wave and, in some cases, lead to hyperkalemia.

All of these analyses were presented in the manuscript, and the parameter changes were provided in Table 2 of the manuscript, yielding the curves shown in Figures 7-10.

This transforms the mathematical model into a physiological model. Furthermore, considering the equation for a cardiac oscillator given by:

${\dot{x}}_{2i}=F_i(t)-{\alpha }_i{x}_{2i}\left(x_{1i}-v_{i_1}\right)\left(x_{2i}-v_{i_2}\right)-\omega \frac{x_{1i}\left(x_{1i}-d_i\right)\left(x_1-e_i\right)}{d_i{e}_i}$

${\alpha }_i$ modifies the shape of the pulse, resulting in a change in the refractory period. e_i Controls the diastolic period; the self-oscillatory nature is maintained whenr $v_{i1}{v}_{i2}<0$ . Varying $v_{i1}$ and $v_{i2}$ allows the resting potential to be lowered or raised, which decreases or increases the frequency of action potential generation. This occurs without affecting the amplitude of the oscillations [62] $F_i(t)$ is the external stimulus that originates from a spatiotemporal stimulus and is therefore considered a reduced-order representation of spatiotemporal aspects

Fixed points analysis

The search and analysis of equilibrium points is carried out by considering the model in its unidirectional connection state. In this case, its conceptual model is given in figure 3, the mathematical model characterizing it is given by equation (6) and the parameters for obtaining normal dynamics in Table 1, Figure 3.

Figure 3: Conceptual model of normal functioning of the heart.

$\begin{array}{c}\{\begin{array}{l}\begin{array}{c}{\dot{x}}_1=x_2\\ {\dot{x}}_2=-{\alpha }_{SA}{x}_2\left(x_1-v_{S{A}_1}\right)\left(x_2-v_{S{A}_1}\right)-\omega \frac{x_1\left(x_1-d_{SA}\right)\left(x_1-e_{SA}\right)}{d_{SA}{e}_{SA}}\end{array}\hfill \\ \begin{array}{r}\hfill {\dot{x}}_3=x_4\\ \hfill \begin{array}{c}{\dot{x}}_4=-{\alpha }_{AV}{x}_4\left(x_3-v_{A{V}_1}\right)\left(x_3-v_{A{V}_1}\right)-\frac{x_3\left(x_3-d_{AV}\right)\left(x_3-e_{AV}\right)}{d_{AV}{e}_{AV}}-k_{SA-AV}\left(x_3-x_3^3\right)\\ +k_{SA-AV}^{\tau }\left(x_1^{{\tau }_{SA-AV}}+{\left(x_1^{{\tau }_{SA-AV}}\right)}^3\right)\end{array}\\ \hfill \begin{array}{c}{\dot{x}}_5=x_6\\ {\dot{x}}_6=-{\alpha }_{HP}{x}_6\left(x_5-v_{H{P}_1}\right)\left(x_5-v_{H{P}_1}\right)-\frac{x_5\left(x_5-d_{HP}\right)\left(x_5-e_{HP}\right)}{d_{HP}{e}_{HP}}-k_{AV-HP}\left(x_5-x_5^3\right)\\ +k_{AV-HP}^{\tau }\left(x_3^{{\tau }_{AV-HP}}+{\left(x_3^{{\tau }_{AV-HP}}\right)}^3\right)\end{array}\end{array}\hfill \end{array}\end{array}\text{ (6)}$

Fixed points

The fixed points or equilibrium points of this system are obtained by solving the system (7).

$\{\begin{array}{l}{\dot{x}}_1=0\hfill \\ {\dot{x}}_2=0\hfill \\ {\dot{x}}_3=0\hfill \\ {\dot{x}}_4=0\hfill \\ {\dot{x}}_5=0\hfill \\ {\dot{x}}_6=0\hfill \end{array}\text{ (7)}$

(7) To facilitate the search for fixed points, a Taylor expansion [63] is proposed. The delay terms make the system more complex, but this allows for more manageable equations. This expression is given in equation (8) and the new system in equation (9).

$\begin{array}{l}x(t-T)=x^T=x(t)-T\dot{x}(t)\text{ (8)}\\ x_2=0\end{array}$

$\begin{array}{c}\{\begin{array}{l}{x}_1\left(x_1+d_{SA}\right)\left(x_1+e_{SA}\right)=0\hfill \\ x_4=0\hfill \\ a_{AV}{x}_3^3+b_{AV}{x}_3^2+c_{AV}{x}_3+k_{AV}=0\hfill \\ x_6=0\hfill \\ a_{HP}{x}_5^3+b_{HP}{x}_5^2+c_{HP}{x}_5+k_{HP}=0\hfill \end{array}\end{array}\text{ (9)}$

With

${\omega }_{SA}\text{ (10)}$

A detailed, step-by-step derivation of equations (9) and (10) is provided in Appendix A Using Cardran's solution method [64], we obtain 27 equilibrium points, 5 of which are actual equilibrium points, namely:

$\begin{array}{c}{P}_1=(0,0,0,0,0,0,);P_4=(0,0,-1.48,0,-1.62,0,);P_7=(0,0,-0.27,0,-0.18,0,)\\ P_{10}=(-2,0,-2.76,0,-1.27,0,);P_{19}=(-3,0,-3.6,0,-3.8,0) And 22\end{array}$

complex equilibrium points, namely:

$P_2=(0,0,0,0,-0.46+1.22i,0,);P_3=(0,0,0,0,-0.46-1.22i,0,);P_5=(0,0,-1.48,0,0.35+\text{1}.\text{64i},0,)$

$P_6=(0,0,-1.48,0,0.35-1.64i,0,);P_8=(0,0,-0.27,0,-0.37+1.2i,0,);P_9=(0,0,-0.27,0,-0.37-\text{1}.\text{2i},0,)$

$P_{11}=(-2,0,-2.76,0,0.17+1.45i,0,);P_{12}=(-2,0,-2.76,0,0.17-1.45i,0,);P_{13}=(-2,0,0.5+1.7i0,0.1+1.7i,0$

$P_{14}=(-2,0,0.5+1.7i,0,-1.5-0.33i,0,);P_{15}=(-2,0,0.5+1.7i,0,-0.5-1.4i,0,$

$P_{16}=(-2,0,0.5-1.7i,0,0.1-1.7i,0,);P_{17}=(-2,0,0.5-1.7i,0,-1.5+0.33i,0,$

$\begin{array}{c}{P}_{18}=(-2,0,0.5-1.7i,0,-0.5+1.4i,0,);P_{20}=(-3,0,-3.6,0,1.46+3.3i,0,);P_{21}\\ =(-3,0,-3.6,0,1.46-3.3i,0,)\end{array}$

$P_{22}=(-3,0,0.94+2.5i,0,0.58+2.5i,0,);P_{23}=(-3,0,0.94+2.5i,0,-2.6-0.37i,0,$

$P_{24}=(-3,0,0.94+2.5i,0,1.14-2.15i,0,);P_{25}=(-3,0,0.94-2.5i,0,0.58-2.5i,0,$

$P_{26}=(-3,0,0.94-2.5i,0,-2.6+0.37i,0,);P_{27}=(-3,0,0.94-2.5i,0,1.14+2.15i,0,)$

The stability analysis is made by computing the Jacobean matrix of the system (6). This matrix is given by equation (11).

$\begin{array}{c}J=\left[\begin{array}{cccccc}0& 1& 0& 0& 0& 0\\ B_1& B_2& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ B_3& B_4& B_5& B_6& 0& 0\\ 0& 0& 0& 0& 0& 1\\ 0& 0& B_7& B_8& B_9& B_{10}\end{array}\right]\end{array}\text{ (11)}$

With

$\left\{\begin{array}{l}{B}_1=-2{\alpha }_{SA}{x}_1{x}_2+{\alpha }_{SA}\left(v_{S{A}_1}+v_{S{A}_2}\right)x_2-\frac{3}{d_{SA}{e}_{SA}}{x}_1^2-\frac{2}{d_{SA}{e}_{SA}}\left(d_{SA}+e_{SA}\right)x_1-1\\ B_2=-{\alpha }_{SA}{x}_1^2+{\alpha }_{SA}\left(v_{S{A}_1}+v_{S{A}_2}\right)x_1-{\alpha }_{SA}{v}_{S{A}_1}{v}_{S{A}_2}\\ B_3=k_{SA-AV}^{\tau }\left(1+3x_1^2-6{\tau }_{SA-AV}{x}_1{x}_2+3{\tau }_{SA-AV}^2{x}_2^2\right)\\ B_4=k_{SA-AV}^{\tau }\left(-{\tau }_{SA-AV}-3{\tau }_{SA-AV}{x}_1^2+6{\tau }_{SA-AV}^2{x}_1{x}_2-3{\tau }_{SA-AV}^3{x}_2^2\right)\\ B_5=-2{\alpha }_{AV}{x}_3{x}_4+{\alpha }_{AV}\left(v_{A{V}_1}+v_{A{V}_2}\right)x_4-\frac{3}{d_{AV}{e}_{AV}}{x}_3^2-\frac{2}{d_{AV}{e}_{AV}}\left(d_{AV}+e_{AV}\right)x_3-1-k_{SA-AV}\left(1+3x_3^2\right)\\ B_6=-{\alpha }_{AV}{x}_3^2+{\alpha }_{AV}\left(v_{A{V}_1}+v_{A{V}_2}\right)x_3-{\alpha }_{AV}{v}_{A{V}_1}{v}_{A{V}_2}\\ B_7=k_{AV-HP}^{\tau }\left(1+3x_3^2-6{\tau }_{AV-HP}{x}_3{x}_4+3{\tau }_{AV-HP}^2{x}_4^2\right)\\ B_8=k_{AV-HP}^{\tau }\left(-{\tau }_{AV-HP}-3{\tau }_{AV-HP}{x}_3^2+6{\tau }_{AV-HP}^2{x}_3{x}_4-3{\tau }_{AV-HP}^3{x}_4^2\right)\\ B_9=-2{\alpha }_{HP}{x}_5{x}_6+{\alpha }_{HP}\left(v_{H{P}_1}+v_{H{P}_2}\right)x_6-\frac{3}{d_{HP}{e}_{HP}}{x}_5^2-\frac{2}{d_{HP}{e}_{HP}}\left(d_{HP}+e_{HP}\right)x_5-1-k_{AV-HP}\left(1+3x_5^2\right)\\ B_{10}=-{\alpha }_{HP}{x}_5^2+{\alpha }_{HP}\left(v_{H{P}_1}+v_{H{P}_2}\right)x_5-{\alpha }_{HP}{v}_{H{P}_1}{v}_{H{P}_2}\end{array}\right\}\text{ (12)}$

The verification of the Jacobian matrix reductions is provided in Appendix B. By resolving $\text{Det}\left(J-\lambda I_6\right)=0$ , we obtain the characteristic equation (13).

$\begin{array}{c}{b}_6{\lambda }^6+b_5{\lambda }^5+b_4{\lambda }^4+b_3{\lambda }^3+b_2{\lambda }^2+b_1\lambda +b_0=0\\ \begin{array}{llllll}1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}\\ I_6=\begin{array}{llllll}0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 1\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill \end{array} is the identity matrix of order 6, and b_1=-B_1{B}_5{B}_{10}-B_1{B}_5-B_5{B}_2\\ b_6=1;b_5=-B_2-B_6-B_{10};b_4=B_2{B}_6-B_5+B_2{B}_{10}+B_6{B}_{10}-B_9;b_0=-B_1{B}_5{B}_9\\ b_2=B_1{B}_5+B_5{B}_9-B_2{B}_6{B}_9-B_1{B}_6{B}_{10}-B_2{B}_5{B}_{10}\\ b_3=B_1{B}_6+B_2{B}_5+B_2{B}_9+B_6{B}_9-B_2{B}_6{B}_{10}+B_5{B}_{10}\end{array}\text{ (13)}$

Using the Ruth-Hurwitz stability criterion, equation (13) gives us 4 stable fixed points and 23 unstable fixed points. The stable fixed points are: $P_{19},P_{20},P_{21}$ et $P_{23}.P_{20}$ and $P_{21}$ are asymptotically stable, $P_{19}$ and $P_{23}$ are marginally stable.

For the fixed point $P_{19}$ , we have the following eigenvalues: ${\lambda }_1=-1.4090;{\lambda }_2=-1.2390;{\lambda }_3=-1.1526;{\lambda }_4=-0.0077;{\lambda }_5=0.0002\text{i};{\lambda }_6=-0.0003\text{i}$

For the fixed point $P_{20}$ , we have the following eigenvalues: ${\lambda }_1=-0.8306-1.6664\text{i};{\lambda }_2=-1.4094;{\lambda }_3=-1.2376;{\lambda }_4=-0.0062-0.0002\text{i};{\lambda }_5=-0.0002-0.0003\text{i};{\lambda }_6=-0.0003+$ 0.0002i

For the case of a fixed point $P_{21}$ , the eigenvalues are: ${\lambda }_1=-0.8306+1.6664\text{i};{\lambda }_2=-1.4094;{\lambda }_3=-1.2376;{\lambda }_4=-0.0062+0.0002\text{i};{\lambda }_5=-0.0002+0.0003\text{i};{\lambda }_6=-0.0003-$ 0.0002i

On the other hand, the fixed point $P_{23}$ , we have the following eigenvalues: ${\lambda }_1=-1.2380;{\lambda }_2=-0.6149-0.6866\text{i};{\lambda }_3=-0.3885-0.1782\text{i};{\lambda }_4=-0.0083+0.0009\text{i};{\lambda }_5=-0.0002\text{i};{\lambda }_6=0.0003\text{i}$

Normal synthetic ECG

Considering the synthetic model of heart function (Figure 3) and the parameters for obtaining normal rhythms given in Table 1, the synthetic ECG is observed in Figure 4.

Figure 4 shows the synthetic ECG signal obtained by digital simulation. The presence of the five main ECG waves is observed, namely the P waves, the QRS complex and the T wave. A comparison with a signal recovered from a real recording is made [65]. The different dynamics of the SA, AV and HP nodes are given in Figure 5.

Figure 4: Normal ECG. (a) Experimental recorded ECG signal [60] (b) Synthetic ECG time series; (c) Synthetic ECG state space.

Figure 5: Normal synthetic ECG represented by each one of the oscillator responses in the form of time series and state space. (a) SA node; (b) AV node; (c) HP complex.

Pathological synthetic ECGs

In this section, the presentation of cardiac pathologies obtained when only the sinus node is disturbed in the absence of any control signal (u = 0) is described following the mathematical model. None of the parameters of the main system in its normal state are modified. Only external stimulus parameters related to the sinus node are modified. In this case, the system is disturbed only at the level of its sinus node and the other disturbances are not considered since we only wanted to control the heart by following the trajectory of the SA node. According to the values taken by the parameters ${\rho }_{SA}$ and ${\omega }_{SA}$ , four types of pathologies dangerous for the normal functioning of the heart are highlighted in order to observe the power of our controller. These are ventricular flutter, ventricular fibrillation, atrial fibrillation and atrial flutter. Their conceptual model is given in Figure (6) and the parameters of obtaining these different pathologies are in table 2,

Figure 6: Conceptual model of heart functioning in its pathological state.

Case of ventricular flutter

Ventricular flutter is a highly intolerable arrhythmia. In this arrhythmia, the ECG shows an equalsized undulation of the QRS complexes, with no isoelectric line between them and no T wave. There is no evidence of any separation between the QRS and the ST-T. During the onset of ventricular flutter, a high heart rate of between 250 and 300/min is noted. As in the case of ventricular fibrillation, an emergency response is required to avoid recurrence, and a therapeutic approach is necessary. Although it can occur spontaneously, supraventricular tachycardia (SVT) is usually the cause of a heart attack.

It is also manifested by a regular appearance of the cardiac rhythm, absence of the P wave, no true QRS complexes, they are widened (> 0.10 seconds) and of sinusoidal shape, a ratio P wave complex QRS not measurable, a PR interval also not measurable [66]. The presence of influx formation due to an ectopic focus on the left and right branches of the His bundle or the Purkinje network (tertiary automatism center) is also observed. In addition, we also note incessant re-entry of the influx at the level of the branches and network of Purkinje or sometimes in the myocardial fibers.

The dynamics of the ECG are given in Figure 7. Figure 7b shows the time series and Figure 7c the state space. The signal of an ECG obtained from a real recording is shown in Figure 7a [67]. It allows a qualitative comparison with the digitally obtained ECG (Figure 7b).

Figure 7: Ventricular flutter ECG signal. (a) Recorded experimental ECG signal [67]; (b) Synthetic ECG time series; (c) Synthetic ECG state space.

Case of atrial fibrillation

This is a cardiac arrhythmia characterized by irregular heartbeats due to abnormal acceleration of the heart caused by rhythm disorders [68]. Old age is also a factor favoring the onset of this arrhythmia. The presence of cardiac pathology (hypertension, heart valve disease, etc.), obesity, or obstructive sleep apnea syndrome is also noted. When atrial fibrillation occurs, the atria are the focus of frequent, ineffective electrical impulses, which can lead to:

• very rapid, jerky contractions (400 to 600 beats / minute) of the atria, resulting in immobility of this part of the heart.
• Poor contraction of the atria due to blood stagnation. This is the risk of a cerebral vascular accident due to the propulsion of clots or thrombus in the artery.
• Rapid, irregular beating of the ventricles at around 150 beats per minute (known as tachyarrhythmia) caused by accelerated contraction of the ventricles, located below the atria. Heart failure can occur.

Arrhythmia can also occur in episodes lasting several days, between which a regular heartbeat is observed. It can also be permanent (in which case the irregularity of cardiac contractions is constant).

The ECG dynamics are given in Figure 8. The real ECG recording is observed in figure 8a [67], the time series and the state space of the ECG obtained digitally in figures 8b and 8c. A qualitative comparison is made between Figures 8a and 8b.

Figure 8: Atrial fibrillation ECG signal. (a) Recorded experimental ECG signal [67]; (b) Synthetic ECG time series. (c) Synthetic ECG state space.

Case of atrial flutter

This is a cardiac arrhythmia marked by the presence of a regular rapid rhythm secondary to a macrocircuit of intra-atrial reentry [69]. When palpitations and sometimes asthenia, exertional intolerance, dyspnea, and presyncope are present, these are referred to as symptoms of atrial fibrillation. An ECG helps to better diagnose the condition. This arrhythmia affects people above the age of 60. It accounts for some 60-65% of adult atrial tachycardias (excluding atrial fibrillation). Approximately 1% of people over 80 account for the number of cases of atrial fibrillation.

In the right atrium, there is a macroreentry circuit. Depolarization of the atria is also observed, with a rate of 250 to 300/min (usually 300/min). This frequency does not affect the AV node, and the ventricles are affected by every other atrial activation (2:1 block), inducing a regular ventricular frequency of 150 beats/minute. Occasionally, in AV block, an irregular ventricular rhythm is observed due to beat-to-beat variation. A fixed 3:1, 4:1 or 5:1 block may also be present, although this is generally less frequent.

By medication, prevention of thrombo-embolic events by anticoagulants, etc., the heart rate can be controlled.

For this case, Figure 9 shows the ECG dynamics: Figure 9a we have the ECG signal obtained from a real recording [67], Figure 9b the time series and Figure 9c the state space obtained digitally. A qualitative comparison is made between the digital ECG and the recorded ECG.

Figure 9: Atrial flutter ECG signal. (a) Recorded experimental ECG signal [67]; (b) Synthetic ECG time series. (c) Synthetic ECG state space.

Case of ventricular fibrillation

Considered the most severe arrhythmia, ventricular fibrillation is an uncoordinated, potentially lethal series of very rapid, ineffective contractions of the ventricles (lower chambers of the heart), caused by numerous chaotic electrical impulses [70]. In such arrhythmias, one observes very weak or absence of contraction on the part of the ventricles, resulting in very rapid, asynchronous pulsations. One of the consequences is cardiac arrest due to the complete cessation of blood flow. The heart frequency is constantly above 300, with discrete QRS complexes on the electrocardiogram.

Age is an important factor in the development of ventricular fibrillation. Other factors can also be a cause of ventricular fibrillation, including diabetes, hypertension, and underlying heart disease. Major consequences include heart failure and CVA.

For this case, Figure 10 shows the ECG dynamics: Figure 10a present the ECG signal obtained from a real recording [67], Figure 10b the time series and Figure 10c the state space obtained digitally. A qualitative comparison is made between the digital ECG and the recorded ECG.

Figure 10: Ventricular fibrillation ECG signal. (a) Recorded experimental ECG signal [67]; (b) Synthetic ECG time series; (c) Synthetic ECG state space.

Controller

Recent work on control strategies shows a wide variety of intelligent control systems, used in the control of chaotic systems [71,72], robotic systems [73], under-actuated systems [69], also for the control of unexpected dynamics in non-linear systems [74]. In cardiac systems, their particularity lies in their ability to pick up cardiac signals, analyze them and then ask the controller to act [41]. In order to prevent the heart from being blocked, the systems used as controllers do not apply permanent control. They act only when they observe an abnormal rhythm, and then cease to act once the heart has returned to its normal operating state.

In order to return the heart to its normal functioning state when confronted with external disturbances that could lead it to a pathological operating state, an intelligent controller is implemented. It consists in tracking the trajectory of the SA node signals (in order to stabilize the ECG), the general operation depending only on the state to be controlled and its derivatives (generally the first derivative), and the calculation of the tracking error. The mathematical expression of this controller is given below by equation (14):

$\begin{array}{c}\ddot{x}SA=g+b_{-}estu+\rho \end{array}\text{ (14)}$

With u the control signal likely to act on the sinus node in the event of an abnormal heart functioning. The vector field part of the system (1) is represented by g. The hidden dynamics and other specifics of the heart rhythm are represented by p. The controlled states are ( $x_{SA},{\dot{x}}_{SA}$ ). Using the feedback linearization method [41], the control signal obtained by a combination of feedback errors, desired states, and estimation of the unknown dynamics is given by:

$\begin{array}{c}u=\left(- g\_est +\ddot{x}SAd-2\lambda es\dot{A}-\lambda 2esA-{\rho }_{\_est}\right)/b_{\_}est\end{array}\text{ (15)}$

Where $g_{\_est },b_{\_est }$ and p_ est are the estimates of g,b and p respectively, the feedback error is $e_{SA}=x_{SA}-x_{SAd }$ , its first derivative is $e_{S\dot{A}},x_{SAd }$ is the desired state of the sinus node, and  being the control parameter. Replacing u by its expression in equation (14) gives

$\begin{array}{c}{e}_{S\ddot{A}}+2\lambda e_{S\dot{A}}+{\lambda }^2{e}_{SA}=\tilde{p}\end{array}\text{ (16)}$

Equation (16) is obtained by considering all modeling uncertainties are correctly represented by ρ, and also assuming that $g=g_{\_est }.\tilde{\rho }=\rho -{\rho }_{\_est }$ is the approximation error.

Using the principle of the sliding mode method, the linear combination of feedback errors is given by

$\begin{array}{c}s\left(x_{SA},{\dot{x}}_{SA}\right)=e_{S\dot{A}} +\lambda e_{SA}\end{array}\text{ (17)}$

It is found that if s = 0, a first-order differential equation whose states $e_{SA}$ and $e_{S\dot{A}}$ will converge to zero with a speed co-controlled by λ is obtained, i.e. the more lambda increases, the more the feedback errors converge to zero. From equation (17), the closed-loop dynamics (Equation (16)) becomes:

$\begin{array}{c}\dot{s}+\lambda s=\tilde{p}\end{array}\text{ (18)}$

It is known that any approximation can be performed using radial basis function (RBF) networks. Using the desired degree of accuracy, an approximation of p is given by $\rho \left(x_{SA},{\dot{x}}_{SA},t\right)={\rho }_{\_est }^{\text{*}}\left(s\left(x_{SA},{\dot{x}}_{SA}\right)\right)+ϵ(t)$ , with ${\rho }_{\_est }^{\text{*}}$ optimal estimation and $|ϵ(t)|\le$

Thus, the calculation of uncertainty and estimated unmodeled dynamics is done using the RBF in Figure 11. Its expression is:

Figure 11: Radial basis function network with one input, n neurons in the hidden layer and one output.

${\rho }_{\_est}=\text{dot}(W,\varpi )$

Where $W=\left[W_1\dots W_n\right]$ is the weight vector and $\varpi =\left[{\varpi }_1\dots {\varpi }_n\right]$ represents the vector of activation functions ${\varpi }_i,n$ being the number of neurons in the hidden layer. The dot function calculates uncertainty and estimated unmodeled dynamics: an internal product of the weight neural network and the activation function

The activation function at each neuron in the neural network is similar to a Gaussian distribution and given by the equation (19)

$\begin{array}{c}\varpi (i)=e^{-0.5{\left(\frac{s-C(i)}{\varphi (i)}\right)}^2}\end{array}\text{ (19)}$

Where $C(i)$ and $\varphi (i)$ respectively represent the centers and widths of the activation function. By applying s to the input, the order of complexity of the neural network is reduced and made sufficiently flexible, which allows it to be easily inserted into intelligent systems for cardiac rhythm control. A study of the bounding and convergence of the s and W signals is carried out using Lyapunov-type stability. For this purpose, a difference between the actual weight vector and the optimal weight vector, respectively, W and W⁺ is define by $\delta =W-W^{\text{*}}$

Knowing that $\tilde{\rho }=\rho -{\rho }_{\_est}$ , with $\rho ={\rho }_{\_est }^{\text{*}}+ϵ$ , and ${\rho }_{\_est }=\text{dot}(W,\varpi )=W^T\varpi$ . Thus equation (18) takes the form :

$\begin{array}{l}\dot{s}+\lambda s=\rho -{\rho }_{\_est}\\ {\rho }_{-}\text{*}est+ϵ-{\rho }_{\_}est\\ =-{\left[W-W^{\text{*}}\right]}^T\varpi +ϵ\text{(20)}\\ =-{\delta }^T\varpi +ϵ\end{array}$

From equation (20), we obtain the differential equation (21) given by:

$\begin{array}{c}\dot{s}=-\lambda s-{\delta }^T\varpi +ϵ\end{array}\text{ (21)}$

Equation (21) confirms the influence of approximation errors on the dynamic behavior of the system when operating in a closed loop. Hence the interest in integrating a neural network capable of making a universal estimate [71] for which the weights can be known on-line.

Suppose a positive function U defined by the equation (22)

$\begin{array}{c}U(t)=\frac{1}{2}{s}^2+\frac{1}{\mu }{\delta }^T\delta \end{array}\text{ (22)}$

Where µ is the learning rate. For µ = 0, neural network learning is disabled and only the FBL scheme is applied. For $U>0,s\ne 0$ and $W>W^{\text{*}}$ . For $U=0,s=0$ and $W=W^{\text{*}}$ Assuming $W^{\text{*}}=C$ ste, $\delta =\dot{W}$ and U is given by:

$\begin{array}{cc}U(t)=s\dot{s}+{\mu }^{-1}{\delta }^{T}W& \\ =-s\left[\lambda s+d_t{\delta }^T\varpi -ϵ\right]+{\mu }^{-1}{\delta }^{T}W& \\ =-s[\lambda s-ϵ]+{\mu }^{-1}{\delta }^T[W-\mu s\varpi ]& Initializing \end{array}$

W for $W=\mu s\varpi$ , we have:

$\begin{array}{c}V(t)=-[\lambda s-ϵ]s\le -[\lambda |s|-ϵ]\left|s\right|\end{array}\text{ (23)}$

From equation (23) it can be seen that W is only limited for $\left|s\right|\le /\lambda$ . Using the projection algorithm studied in [72], the limit problem of is guaranteed and its domain of definition is taken equal to $Df=W$

$\begin{array}{rrr}\hfill & \hfill & \hfill \left\{W\in I{R}^n:\text{dot}(W,W)\le \eta \right\}:\\ \hfill & \hfill \dot{W}=\{\begin{array}{cc}\mu s\varpi & if \Vert W\Vert <\eta or \Vert W\Vert =\eta and \mu \text{sdot}(W,\varpi )\\ \left(I-\frac{W{W}^T}{\text{dot}(W,\varpi )}\right)\mu s\varpi & otherwise \end{array}& \hfill \end{array}\text{ (24)}$

Where $\eta$ is the desired upper limit. It is used to control the limit of increase of the weight neural network. So any value of. So any value of $\Vert W\Vert$ , it will always be less than $\eta$ , hence $\left|s\right|\le /\lambda$ , which implies $-/\lambda \le s\le /\lambda$ . By replacing s by its expression taken from equation (9), we obtain equation

$\begin{array}{c}-/\lambda \le e_{S\dot{A}} +\lambda e_{SA}\le /\lambda \end{array}\text{ (25)}$

The multiplication of equation (25) by $e^{\lambda t}$ gives the equation

$\begin{array}{c}-e^{\lambda t}(/\lambda )\le \frac{d}{dt}\left(e_{SA}{e}^{\lambda t}\right)\le e^{\lambda t}(/\lambda )\end{array}$

By integrating (26) between 0 and t and multiplying the result $e^{-\lambda t}$ , we obtain

$\begin{array}{llllllll}{\lambda }_2\hfill & SA\hfill & {\lambda }_2\hfill & SA\hfill & {\lambda }_2\hfill & SA\hfill & {\lambda }_2\hfill & (27) - — +[e\hfill \end{array}$

For $t\to \infty$ equation (27) results in

$\begin{array}{rrr}\hfill e\le -+[e& \hfill & \hfill (0)--]e-\lambda t\\ \hfill & \hfill -\overline{{\lambda }^2}\le e_{SA}\le \frac{\overline{{\lambda }^2}}{ }& \hfill \end{array}\text{ (28)}$

Inserting (28) into (25) gives :

$\begin{array}{c}-\frac{2}{\lambda }\le {\dot{e}}_{SA}\le \frac{2}{\lambda }\end{array}\text{ (29)}$

Equation (29) leads to the conclusion that the exponential convergence of the tracking error to the closed region $\varphi =\{\left(e_{SA},{\dot{e}}_{SA}\right)\in I{R}^2:\left|e_{SA}\right|\le \frac{\epsilon }{{\lambda }^2}$ et $\left|{\dot{e}}_{SA}\right|\le \frac{2\epsilon }{\lambda }\}$ is guaranteed by the controller.

Controlled dynamics

In this section, the results obtained digitally after implementation of the control are presented. The sampling frequency of the signals obtained is taken to be 100 Hz. To obtain the numerical results, the controller is programmed to follow the disturbed trajectories of the SA node, since in its disturbed state, the heart operates in pathological behavior (Figure 5). To normalize the system, a text file with the desired SA node signals extracted from a normal rhythm (from previous simulations) is read in. Since the sinus node is the seat of the electrical impulse, because it takes information from the brain and adapts the heart rate according to the body's needs, the electrical impulse from the sinus node will in turn travel through the two atria, contracting them and encouraging the passage of blood to the ventricles. The electrical impulse then reaches the AV node, where it is filtered and slowed down in the event of disruption of the SA node, before reaching the ventricles. The electrical impulse reaching the ventricles contracts them, enabling the right ventricle to propel blood to the lungs and the left ventricle to the rest of the body. By disrupting the sinus node, the course and nature of the electrical impulse through the heart is also disrupted, resulting in cardiac arrhythmias in the heart. In the present case, the pathological rhythms revealed are: ventricular flutter, ventricular fibrillation, atrial flutter or atrial fibrillation.

The learning rate and control parameter are $\mu =1500$ and $\lambda =50$ respectively. Taking $g_{\_est }=$ 0, it is assumed that the control system is known without any knowledge about the cardiac model. Thus, to handle all neglected dynamic effects ${\rho }_{\_est }$ is studied. The centers and widths of the activation function defined in equation (11) are given by $C=[-K/2.5,-K/5,-K/10,K/10,K/5,K/2.5]$ et, $\varphi =[K/2,K/4,K/8,K/8,K/4,K/2].\text{K}$ is the proportionality constant for distributing the neural network activation functions and is taken to be K = 40.0. The weight vector is initialized to W = 0.

His analyses and the results obtained from numerical simulations show that, when the heart is in a situation of ventricular flutter, ventricular fibrillation, atrial fibrillation, or atrial flutter, once the command has been applied, the ECG dynamics show a considerable approximation to normal rhythms. Abnormal rhythms are thus avoided, enabling the heart to pump enough blood to satisfy the body's needs, while stopping the sudden onset of a heart attack.

Figure 12 shows the conceptual model of the heart in its controlled state.

Figure 12: Conceptual model of heart functioning in its pathological state.

When the control signal u is applied to the system (1) in its sick state, the controlled ECG dynamics are shown in Figures 13-16 respectively. Figure 13 shows those obtained in the case of ventricular flutter. For ventricular fibrillation, they are shown in figure 14. Figure 15 shows the results for atrial fibrillation, and figure 16 for atrial flutter.

Figure 13: (a) Control signal. (b) Uncontrolled synthetic ECG time series. (c) Uncontrolled synthetic ECG state space. (d) Controller synthetic ECG time series. (e) Controller synthetic ECG state space.

Figure 14: Control signal. (b) Uncontrolled synthetic ECG time series. (c) Uncontrolled synthetic ECG state space. (d) Controller synthetic ECG time series. (e) Controller synthetic ECG state space.

Figure 15: Control signal. (b) Uncontrolled synthetic ECG time series. (c) Uncontrolled synthetic ECG state space. (d) Controller synthetic ECG time series. (e) Controller synthetic ECG state space.

Figure 16: Control signal. (b) Uncontrolled synthetic ECG time series. (c) Uncontrolled synthetic ECG state space. (d) Controller synthetic ECG time series. (e) Controller synthetic ECG state space.

Figure 16. Control signal. (b) Uncontrolled synthetic ECG time series. (c) Uncontrolled synthetic ECG state space. (d) Controller synthetic ECG time series. (e) Controller synthetic ECG state space Figures 13-16 show the action of the controller in cardiac arrhythmia situations. The synthetic ECGs show a qualitative and quantitative approximation between controlled and normal dynamics. Both the phase spaces and the time traces clearly show a very strong improvement in the performance of the highlighted control scheme, as well as its ability to replace a pathological cardiac rhythm with a normal one. The proposed intelligent control scheme shows its ability to make abnormal dynamics disappear.

It can be seen from the curves in figures 13-16 that the action of the controller starts at 15 s. Before 15 s, the controller has no effect on the heart. In the case of ventricular flutter, the control signal, once engaged, evolves at a frequency of 2.4 Hz. In ventricular fibrillation, the control signal has a frequency of 7.2 Hz. In atrial fibrillation, the control signal has a frequency of 8.2 Hz. In atrial flutter, a frequency of 0.6 Hz is observed. Whether in ventricular flutter, ventricular fibrillation, atrial fibrillation or atrial flutter, as observed in normal rhythm, the controlled dynamics clearly show the presence of the QRS complex and P and T waves. We also note the disappearance of the R and S peaks, reflecting the controller's skill in solving this problem

Thus, the feedback linearization control (FBL) method with the integration of an artificial neural network is very effective and robust for the control of pathological heartbeats.

Electronic implementation

In this subsection, we present the results obtained by electronic simulation using OrCAD-PSpice software. This is one of the most widely used simulations in electronic product design, whether for PCB design, IC design or other electrical systems. For the analysis of the most important electronic circuits, OrCAD-PSpice seems to be the most suitable simulation. Wired assemblies in this software allow the use of components with values that reproduce real components as effectively as possible. Once the device has been wired, the results can be analyzed in all their facets, using tools as varied as they are sophisticated. A large number of research studies have shown that OrCAD-PSpice is a world reference in the field of circuit design and simulation [4,27], as it possesses some of the most interesting features, including: the possibility of creating virtual components, using manufacturers' data sheets, downloading component libraries.

These include: reducing time wasted in choosing component values before the practical implementation of the system and exploring the advantages of different practical implementation options. Faced with the complexity of creating electronic assemblies whose results can be used to interpret real phenomena, OrCAD-PSpice software helps designers to solve problems at different levels of abstraction, to better understanding the interactions between the various components making up the system.

To implement the mathematical system on the OrCAD-Pspice software, Kirchhoff's and Ohm's laws are used and the system is modeled in the form of differential equations [27]. In this subsection, the dynamics of the system in its state of normal functioning are presented. The system used is the one described by the functional model in Figure 3.

Using electronic equivalence techniques and operational amplifier integration (AOP) principles, the OrCAD-PSpice assembly of the SA oscillator, AV oscillator, and HP oscillator coupled together is shown in Figures 17-19 respectively.

Figure 17: Electronic equivalence of the circuit representing the SA oscillator mounted on OrCAD-PSpice software.

Figure 18: Electronic equivalence of the circuit representing the AV oscillator mounted on OrCAD-PSpice software.

Figure 19: Electronic equivalence of the circuit representing the HP oscillator mounted on OrCAD-PSpice software.

The above circuit implements the analog equivalence of the SA oscillator. It is composed of nine (9) resistors ( ${\text{R}}_1\text{SA},{\text{R}}_2\text{SA},{\text{R}}_3\text{SA},{\text{R}}_4\text{SA},{\text{R}}_5\text{SA},{\text{R}}_6\text{SA},\dots \dots \dots$ , and ${\text{R}}_9\text{SA}$ ), 2 capacitors ( ${\text{C}}_1\text{SA}$ and ${\text{C}}_2\text{SA}$ ), 4 operational amplifiers (TL082), 4 multipliers (AD633/AD), DC voltage sources to power integrated circuits (AOP and multipliers) for some and for others they are used as parameters of the equivalent circuit. We also note an AC power source ( ${\text{f}}_{\text{SA}}(\text{t})$ ) used to excite the circuit. Using the principle of Kirchhoff's laws to Ohm's law, the relationship between the inputs is given below:

${\text{v}}_2=\frac{1}{10}\left({\text{x}}_1-{\text{V}}_{1\text{SA}}\right)\left({\text{x}}_1-{\text{V}}_{2\text{SA}}\right)\text{ (30)}$

${\text{B}}_2=\frac{1}{100}{\text{x}}_2\left({\text{x}}_1-{\text{V}}_{1\text{SA}}\right)\left({\text{x}}_1-{\text{V}}_{2\text{SA}}\right)\text{ (31)}$

${\text{ v}}_3=\frac{1}{10}\left({\text{x}}_1+{\text{d}}_{\text{SA}}\right)\left({\text{x}}_1+{\text{e}}_{\text{SA}}\right)\text{ (32)}$

${\text{B}}_3=\frac{1}{100}{\text{x}}_1\left({\text{x}}_1+{\text{d}}_{\text{SA}}\right)\left({\text{x}}_1+{\text{e}}_{\text{SA}}\right)\text{ (33)}$

${\text{B}}_1={\text{C}}_{1\text{SA}}{\text{R}}_{1\text{SA}}\frac{{\text{dx}}_2}{\text{dt}} \text{and} x_2={\text{C}}_{2\text{SA}}{\text{R}}_{2\text{SA}}\frac{{\text{dx}}_1}{\text{dt}}\text{ (34)}$

$\frac{\frac{{\text{B}}_1}{{\text{R}}_{7\text{SA}}}+\frac{{\text{B}}_2}{{\text{R}}_{5\text{SA}}}+\frac{{\text{B}}_3}{{\text{R}}_{6\text{SA}}}}{\frac{1}{{\text{R}}_{7\text{SA}}}+\frac{1}{{\text{R}}_{5\text{SA}}}+\frac{1}{{\text{R}}_{6\text{SA}}}}=\frac{{\text{R}}_{9\text{SA}}{\text{f}}_{\text{SA}}(\text{t})}{{\text{R}}_{9\text{SA}}+{\text{R}}_{8\text{SA}}}\text{ (35)}$

From equation (35), we obtain equation (36) below.

${\text{B}}_1=\frac{{\text{R}}_{7\text{SA}}{\text{R}}_{9\text{SA}}}{{\text{R}}_{9\text{SA}}+{\text{R}}_{8\text{SA}}}{\text{f}}_{\text{SA}}\left(\text{t}\right)-\frac{{\text{R}}_{7\text{SA}}}{{\text{R}}_{5\text{SA}}}{\text{B}}_2-\frac{{\text{R}}_{7\text{SA}}}{{\text{R}}_{6\text{SA}}}{\text{B}}_3\text{ (36)}$

The substitution of equations (31) and (33) in (35) gives us equation (37).

${\text{B}}_1=\frac{{\text{R}}_{7\text{SA}}{\text{R}}_{9\text{SA}}}{{\text{R}}_{9\text{SA}}+{\text{R}}_{8\text{SA}}}{\text{f}}_{\text{SA}}\left(\text{t}\right)-\frac{{\text{R}}_{7\text{SA}}}{100{\text{R}}_{5\text{SA}}}{\text{x}}_2\left({\text{x}}_1-{\text{V}}_{1\text{SA}}\right)\left({\text{x}}_1-{\text{V}}_{2\text{SA}}\right)-\frac{{\text{R}}_{7\text{SA}}}{100{\text{R}}_{6\text{SA}}}{\text{x}}_1\left({\text{x}}_1+{\text{d}}_{\text{SA}}\right)\left({\text{x}}_1+{\text{e}}_{\text{SA}}\right)\text{ (37)}$

The relationships between the parameters of the SA node and the data in equation (37) allow us to establish the following correspondences given by equation (38).

$\frac{{\text{R}}_{7\text{SA}}{\text{R}}_{9\text{SA}}}{{\text{R}}_{9\text{SA}}+{\text{R}}_{8\text{SA}}}=1;{\text{ C}}_{1\text{SA}}{\text{R}}_{1\text{SA}}={\text{C}}_{2\text{SA}}{\text{R}}_{2\text{SA}}=1;\frac{{\text{R}}_{7\text{SA}}}{100{\text{R}}_{5\text{SA}}}={\alpha }_{\text{SA}}\text{ and }\frac{{\text{R}}_{7\text{SA}}}{100{\text{R}}_{6\text{SA}}}=\frac{1}{{\text{e}}_{\text{SA}}{\text{d}}_{\text{SA}}}\text{ (38)}$