The dynamical model

The model generates a trajectory in a three-dimensional state space with co-ordinates $(x,y,z)$. Quasi-periodicity of the ECG is reflected by the movement of the trajectory around an attracting limit cycle of unit radius in the $(x,y)$-plane. Each revolution on this circle corresponds to one RR-interval or heart beat. Inter-beat variation in the ECG is reproduced using the motion of the trajectory in the $z$-direction. Distinct points on the ECG, such as the P,Q,R,S and T are described by events corresponding to negative and positive attractors/repellors in the $z$-direction. These events are placed at fixed angles along the unit circle given by $\theta_P$, $\theta_Q$,$\theta_R$,$\theta_S$ and $\theta_T$ (see Fig. 2). When the trajectory approaches one of these events, it is pushed upwards or downwards away from the limit cycle, and then as it moves away it is pulled back towards the limit cycle. The dynamical equations of motion are given by a set of three ordinary differential equations

$${\dot x} = \alpha x - \omega y,$$

$${\dot y} = \alpha y + \omega x,$$

$${\dot z} = - \!\!\!\!\!\! \sum_{i \in \{P,Q,R,S,T\}} \!\!\!\!\!\! a_i \Delta \theta_i \exp(-\Delta \theta_i^2 / 2 b_i^2) - (z - z_0),$$ (1)

where $\alpha = 1 - \sqrt{x^2 + y^2}$, $\Delta \theta_i = (\theta - \theta_i) \ {\rm mod} \ 2 \pi$, $\theta = {\rm atan2}(y,x)$ and $\omega$ is the angular velocity of the trajectory as it moves around the limit cycle. Baseline wander was introduced by coupling the baseline value $z_0$ in (1) to the respiratory frequency $f_2$ using

$$z_0(t) = A \sin(2 \pi f_2 t),$$ (2)

where $A = 0.15$ mV. These equations of motion given by (1) were integrated numerically using a fourth order Runge-Kutta method [15] with a fixed time step $\Delta t = 1/f_s$ where $f_s$ is the sampling frequency. Visual analysis of a section of typical ECG from a normal subject was used to suggest suitable times (and therefore angles $\theta_i$) and values of $a_i$ and $b_i$ for the PQRST points. The times and angles are specified relative to the position of the R-peak as shown in Table I. A trajectory generated by equation (1) in three-dimensions corresponding to $(x,y,z)$ is illustrated in Fig. 2. This demonstrates how the positions of the events $P,Q,R,S,T$ act on the trajectory in the $z$-direction as it precesses around the unit circle in the $(x,y)$-plane. The $z$ variable from the three-dimensional system (1) yields a synthetic ECG with realistic PQRST morphology (Fig. 3). The similarity between the synthetic ECG and the real ECG may be seen by comparing Fig. 3 with Fig. 1. Note that noise has not been added to the model at this point.

Table I: Parameters of the ECG model given by (1)

Index (i)	P	Q	R	S	T
Time (secs)	-0.2	-0.05	0	0.05	0.3
$\theta_i$ (radians)	$-\frac{1}{3}\pi$	$-\frac{1}{12}\pi$	0	$\frac{1}{12}\pi$	$\frac{1}{2}\pi$
$a_i$	1.2	-5.0	30.0	-7.5	0.75
$b_i$	0.25	0.1	0.1	0.1	0.4

Figure 2: A typical trajectory generated by the dynamical model (1) in the three-dimensional space given by $(x,y,z)$. The dashed line reflects the limit cycle of unit radius while the small circles show the positions of the P,Q,R,S,T events.

Figure 3: Morphology of one PQRST-complex of the ECG.

By contrasting the dynamical model (1) with the mechanisms underlying the cardiac cycle, it is obvious that the time required to complete one lap of the limit cycle is equal to the RR-interval of the synthetic ECG signal. Variations in the length of the RR-intervals can be incorporated by varying the angular velocity $\omega$. The effects of both RSA and Mayer waves in the power spectrum $S(f)$ of the RR-intervals are incorporated by generating RR-intervals which have a bimodal power spectrum consisting of the sum of two Gaussian distributions,

$$S(f) = \frac{\sigma_1^2}{\sqrt{2 \pi c_1^2}} \exp \left( ... ...c_2^2}} \exp \left( \frac{(f - f_2)^2}{2 c_2^2} \right),$$ (3)

with means $f_1,f_2$ and standard deviations $c_1,c_2$. Power in the LF and HF bands are given by $\sigma_1^2$ and $\sigma_2^2$ respectively whereas the variance equals the total area $\sigma^2 = \sigma^2_1+\sigma^2_2$, yielding an LF/HF ratio of $\sigma^2_1/\sigma^2_2$. Fig. 4 shows the power spectrum $S(f)$ given by $f_1 = 0.1$, $f_2 = 0.25$, $c_1 = 0.01$, $c_2 = 0.01$ and $\sigma^2_1/\sigma^2_2 = 0.5$. The Gaussian frequency distribution is motivated by the typical power spectrum of a real RR tachogram [7].

Figure 4: Power spectrum $S(f)$ of the RR-interval process with a LF/HF ratio of $\sigma _1^2/\sigma _2^2 = 0.5$.

A RR-interval time series $T(t)$ with power spectrum $S(f)$ is generated by taking the inverse Fourier transform of a sequence of complex numbers with amplitudes $\sqrt{S(f)}$ and phases which are randomly distributed between 0 and $2 \pi$. By multiplying this time series by an appropriate scaling constant and adding an offset value, the resulting time series can be given any required mean and standard deviation. Suppose that $T(t)$ represents the time series generated by the RR-process with power spectrum $S(f)$. The time-dependent angular velocity $\omega(t)$ of motion around the limit cycle is then given by

$$\omega(t) = \frac{2 \pi}{T(t)}.$$ (4)

In this way the series of RR-intervals of the resultant synthetic ECG will also have a power spectrum equal to $S(f)$; this will be demonstrated in the next section.