Cardiac signal estimation based on the arterial and venous pressure signals of a hemodialysis machine

Most patients with chronic kidney disease also suffer from cardiovascular diseases, causing approximately 45% of all deaths among hemodialysis patients—a much larger proportion than in the general population (Shastri and Sarnak 2010). The main causes are coronary artery disease, congestive heart failure, and arrhythmias (Di Lullo et al 2014). Premature ventricular beats and complex ventricular arrhythmias occur more often during hemodialysis than after (Burton et al 2008). Continuous cardiac monitoring during hemo-dialysis improves patient management and helps to reduce the risk of adverse events as it makes manual or automatic modification of the treatment possible (Shamir et al 1999, Selby and McIntyre 2007, Mancini et al 2008, Javed et al 2010, Solem et al 2010, Yang and Chon 2010, Sörnmo et al 2012). Intradialytic cardiac monitoring is, however, not part of today's clinical routine, one reason being the discomfort caused by wearing electrodes for recording the electrocardiogram (ECG) or a sensor for recording the photoplethysmogram (PPG) during several hours. Another reason is the additional work required by the clinical staff to attach the electrodes or sensors. Hence, it is of vital clinical importance to develop a method for intradialytic cardiac monitoring that does not cause discomfort or add to the workload of the staff.

Pressure waves originating from the heart propagate through the body, enter the extracorporeal blood circuit of the hemodialysis machine, and eventually reach the arterial and venous pressure sensors. Blood from the patient flows through the arterial blood lines of the extracorporeal blood circuit at a flow rate determined by the peristaltic pump. The pump pushes the blood through the venous side, including the dialyzer, back to the patient. Due to the occlusive nature of the peristaltic pump, pressure waves are stopped from propagating between the arterial and venous sides. The pump generates strong periodic pressure waves whose amplitude is considerably higher than that of the pressure waves originating from the heart. Hence, the estimation of a cardiac signal represents a challenging signal separation problem since not only is the signal-to-noise ratio (SNR) poor, but the heart rate and the pump rate may coincide.

The estimated cardiac signal may be analyzed for the purpose of detecting different types of events—cardiac arrest and venous needle dislodgement being two particularly important events to detect as they cause a significant number of deaths among hemodialysis patients (Karnik et al 2001, Herzog 2003, Polaschegg 2010, Segelmark et al 2010, Axley et al 2012). Simultaneous estimation of the cardiac signal using both the venous and arterial pressures would allow event-specific alarms; this is not possible with a method confined to using only the venous pressure sensor. The cardiac signal is likely to disappear in both the venous and arterial pressure signals when cardiac arrest occurs, whereas the cardiac signal is likely to disappear only in the venous signal when venous needle dislodgement occurs. When detecting venous needle dislodgement, it is essential that the blood pump is automatically halted so that the patient's life can be saved. When cardiac arrest is detected, it is essential to elicit an alarm so that the clinical staff can initiate cardiopulmonary resuscitation.

It has earlier been demonstrated that cardiac information can be extracted from the signal produced by the extracorporeal venous pressure sensors of a hemodialysis machine (Holmer et al 2015). In that study, however, the venous pressure was analyzed under idealized conditions since only a 20 min data interval was selected in each patient; this interval was excerpted from a part of the treatment where the blood flow remained constant at a normal flow rate. The results showed that the performance deteriorated considerably at lower SNRs, producing an unreliable heart rate estimate.

In the present paper, a novel method for cardiac signal estimation is proposed which rely on information from both the arterial and the venous extracorporeal pressure sensors, hypothesizing that such information improves the quality of the estimated signal. The earlier assumption of all heartbeats having identical morphology is now relaxed so that ventricular ectopic beats can be detected. Rather than modeling the pump pressure wave for a complete pump revolution, separate modeling of the two rollers of the blood pump is introduced so that more accurate characterization can be achieved. The observed signal is processed on a segment-by-segment basis for robust handling of changes in blood flow rate. The performance of the proposed method is evaluated on a set of complete hemodialysis treatments. The main improvements compared to the previous method are: 1. use of signals from both the arterial and venous pressure sensors, enabling the extraction of low cardiac amplitude signals, 2. handling of all possible blood flow rates, even low ones, 3. automatic handling of flow rate changes, and 4. automatic detection and handling of disturbances. Points 2–4 ensure that cardiac signal extraction can be performed throughout the entire treatment session.

From the observed arterial and venous pressure signals, the method separates the heart-related arterial and venous pressure signals from the pump-related arterial and venous pressure signals using an iterative estimation procedure. This procedure alternates between estimating the pump-related signals and estimating the cardiac signals, and terminates when the difference between cardiac signal estimates of successive iterations has become sufficiently small.

The pump profile is an important concept which is introduced for characterizing the typical pump-induced pressure wave. A pump profile is computed for each of the the arterial and the venous pressure signals. The two pump profiles are used to create arterial and venous pump model signals by concatenating the profiles once subjected to suitable time scaling, as described below. The resulting arterial and venous pump model signals are, in turn, subtracted from the corresponding observed pressure signal so that improved estimates of the cardiac signals are obtained.

Table 1 summarizes the notations used in the following to describe the method for cardiac signal estimation.

Table 1. Notations used in the present method.

Variable name	Description
j	Iteration index
q	Subscript denoting roller 1 or 2 of the pump rotor
x	Subscript denoting either arterial (a) or venous (v) pressure
$\hat{a}_{xq,2k}^{\left(\,j\right)}$	Estimate of time scale factor of the kth pressure wave
$\check {{c}}\,_{x}^{\left(\,j\right)}(t)$	Intermediate cardiac signal estimate
$\hat{c}_{x}^{\left(\,j\right)}(t)$	Cardiac signal estimate
${{\hat{c}}^{\left(\,j\right)}}(t)$	Joint cardiac signal estimate
$\Delta \hat{t}$	Estimated time delay between $\hat{c}_{a}^{\left(\,j\right)}(t)$ and $\hat{c}_{v}^{\left(\,j\right)}(t)$
N	Number of complete pump revolutions
$\bar{p}_{xq}^{\left(\,j\right)}(t)$	Pump profile
$\tilde{p}_{x}^{\left(\,j\right)}(t)$	Pump model signal
$\hat{p}_{x}^{\left(\,j\right)}(t)$	Pump signal estimate
${{\hat{t}}_{c,k}}$	Occurrence time of the kth heartbeat
${{\hat{t}}_{x,k}}$	Occurrence time of the kth pressure wave
${{\bar{T}}_{xq}}$	Mean length of one roller period
yx(t)	Observed pressure signal
${{F}_{\text{LP}}}$	Low pass cut-off frequency when filtering ${{\hat{c}}^{\left(\,j\right)}}(t)$
${{F}_{\text{HP}}}$	High pass cut-off frequency when filtering ${{\hat{c}}^{\left(\,j\right)}}(t)$
L	Maximum number of iterations
λ	Max time lag between the heart pulses in $\hat{c}_{a}^{\left(\,j\right)}(t)$ and $\hat{c}_{v}^{\left(\,j\right)}(t)$
α	Threshold used for initial estimate of roller period

2.1. Iterative cardiac signal estimation

An iterative procedure is proposed in which an intermediate cardiac signal estimate $\check {{c}}\,_{x}^{\left(\,j\right)}(t)$ is determined by subtracting a pump model signal $\tilde{p}_{x}^{\left(\,j-1\right)}(t)$ from the observed pressure signal y_x(t), i.e.

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \check {{c}}\,_{x}^{\left(\,j\right)}(t)={{y}_{x}}(t)-\tilde{p}_{x}^{\left(\,j-1\right)}(t),\quad x\in \left\{a,v\right\}. \end{array}\end{eqnarray} \tag{ 1 }$$

A generic subscript x is added to indicate that subtraction is performed for both arterial (a) and venous (v) blood pressure. The procedure is initiated by setting the iteration index j  =  1. The initial pump model signal $\tilde{p}_{x}^{(0)}(t)$ is determined by the procedure described in appendix. Since segment-based processing is employed, the initialization procedure is repeated for each segment. In contrast to (Holmer et al 2015), the present initialization procedure can handle widely different blood flow rates.

The amplitude and time delay $\Delta t$ of the cardiac waves in $y_{v}^{\left(\,j\right)}(t)$ usually differ from the ones in $y_{a}^{\left(\,j\right)}(t)$. This is explained by differences between the arterial and venous sides of the extracorporeal blood circuit. The main differences are in the length of the tubing, the pressure, and the amount of air (due to the air in the drip chamber that is present only on the venous side). On the other hand, the morphology of the cardiac waves in $y_{v}^{\left(\,j\right)}(t)$ and $y_{a}^{\left(\,j\right)}(t)$ are similar, rendering it possible to combine the intermediate venous and arterial cardiac signals once subjected to time alignment.

The time delay $\Delta t$ is estimated by determining the location of the maximum crosscorrelation between $\check {{c}}\,_{a}^{\left(\,j\right)}(t)$ and $\check {{c}}\,_{v}^{\left(\,j\right)}(t)$ (Kay 1993),

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \Delta \hat{t}=\text{arg}\,\underset{| \Delta t|<\lambda}{{\max}}\,{\int}_{0}^{T}\check {{c}}\,_{a}^{\left(\,j\right)}(t)\check {{c}}\,_{v}^{\left(\,j\right)}\left(t+ \Delta t\right)\text{d}t, \end{array}\end{eqnarray} \tag{ 2 }$$

where T denotes the segment length (here set to 1 min) and λ the maximum lag between the cardiac waves of the arterial and venous signals. The delay is assumed to be constant within each segment.

A joint cardiac signal estimate ${{\hat{c}}^{\left(\,j\right)}}(t)$ is then computed from the intermediate cardiac signal estimates in (1) as a weighted average of $\check {{c}}\,_{a}^{\left(\,j\right)}(t)$ and $\check {{c}}\,_{v}^{\left(\,j\right)}(t)$,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{{\hat{c}}}^{\left(\,j\right)}}(t)={{w}_{a}}\check {{c}}\,_{a}^{\left(\,j\right)}(t)+{{w}_{v}}\check {{c}}\,_{v}^{\left(\,j\right)}\left(t+ \Delta \hat{t}\right), \end{array}\end{eqnarray} \tag{ 3 }$$

where two different approaches to determining the weights w_a and w_v are investigated: fixed, equal weights given by arithmetic averaging (AA),

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{w}_{a}}={{w}_{v}}=\frac{1}{2}, \end{array}\end{eqnarray} \tag{ 4 }$$

and data-dependent weights given by principal component analysis (PCA) (Hyvärinen et al 2001) of the data matrix

$$\begin{eqnarray}&&\check {{\mathbf{C}}}\,=\left[\begin{array}{*{35}{l}} {{\check {{\mathbf{c}}}\,}_{a}} \\ {{\check {{\mathbf{c}}}\,}_{v}} \end{array}\right].\end{eqnarray} \tag{ 5 }$$

The row vectors ${{\check {{\mathbf{c}}}\,}_{a}}$ and ${{\check {{\mathbf{c}}}\,}_{v}}$ result from equidistant sampling of $\check {{c}}\,_{a}^{\left(\,j\right)}(t)$ and $\check {{c}}\,_{v}^{\left(\,j\right)}\left(t+ \Delta \hat{t}\right)$, respectively, using a sufficiently high rate. The two eigenvectors of the correlation matrix

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\mathbf{R}}_{c}}=\check {{\mathbf{C}}}\,{{\check {{\mathbf{C}}}\,}^{T}} \end{array}\end{eqnarray} \tag{ 6 }$$

are then determined and used for defining the orthogonal linear transformation that produces the principal components. In order to ensure that the two pressure components remain in phase, the eigenvector of ${{\mathbf{R}}_{c}}$ producing the largest absolute value of the scalar product with the direction $\frac{1}{\sqrt{2}}\left[\begin{array}{*{35}{l}} 1 & 1 \end{array}\right]$ is identified; the two elements of the eigenvector determine the weights w_a and w_v. Selecting the eigenvector with largest absolute value of the scalar product implies that the eigenvector lies in the 1st or the 3d quadrant. Thus the weights w_a and w_v have the same sign, which in turn means that phase between the two signals remains unaltered.

Starting from the joint cardiac signal estimate ${{\hat{c}}^{\left(\,j\right)}}(t)$, the cardiac signal estimates $\hat{c}_{a}^{\left(\,j\right)}(t)$ and $\hat{c}_{v}^{\left(\,j\right)}(t)$ related to the arterial and venous pressure, respectively, are determined by suitable rescaling and realignment,

$$\begin{eqnarray}\begin{array}{*{35}{l}} \hat{c}_{a}^{\left(\,j\right)}(t) & =\frac{{{w}_{a}}}{w_{a}^{2}+w_{v}^{2}}{{{\hat{c}}}^{\left(\,j\right)}}(t), \end{array}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} \hat{c}_{v}^{\left(\,j\right)}(t) & =\frac{{{w}_{v}}}{w_{a}^{2}+w_{v}^{2}}{{{\hat{c}}}^{\left(\,j\right)}}\left(t- \Delta \hat{t}\right). \end{array}\end{eqnarray} \tag{ 8 }$$

These two cardiac signal estimates are used below to compute the pump signal estimates.

2.2. Iterative pump signal estimation

The rotor of the peristaltic pump consists of two rollers, and, accordingly, one complete rotor revolution gives rise to two pump pressure waves, one for each roller actuation. As a result, two subsequent pump pressure waves will differ in morphology, whereas every second wave is generated by the same roller and, therefore, has similar morphology. Hence, to achieve better characterization, each of the two pressure waves has its own pump profile, denoted $\bar{p}_{xq}^{\left(\,j\right)}(t),q=1,2,x\in \left\{a,v\right\}$, leading to that four different pump profiles have to be computed for one complete rotor revolution, see the illustration in figure 1. It should be noted that the waves of the arterial and the venous pump pressures differ with respect to amplitude and morphology since the sensors are located on opposite sides of the two pump rollers, either pulling the blood through the circuit on the arterial side of the pump or pushing the blood on the venous side of the pump.

Zoom In Zoom Out Reset image size

The pump profile $\bar{p}_{xq}^{\left(\,j\right)}(t)$ is determined as the ensemble median resulting from the time-aligned roller periods of the observed signal y_x(t). Prior to time alignment, however, the occurrence time $t_{x,k}^{\left(\,j\right)}$ of each roller period, and the related time scale factor $a_{xq,k}^{\left(\,j\right)}$, need to be estimated; the roller periods are indexed by k, ranging from 1 to 2N where N denotes the number of complete pump revolutions.

The estimation of $t_{x,k}^{\left(\,j\right)}$ takes its starting point in the pump signal estimate $\hat{p}_{x}^{\left(\,j\right)}(t)$, given by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \hat{p}_{x}^{\left(\,j\right)}(t)={{y}_{x}}(t)-\hat{c}_{x}^{\left(\,j\right)}(t) \end{array}\end{eqnarray} \tag{ 9 }$$

which is used for computing the residual

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\epsilon}_{x}}\left(t;t_{x,k}^{\left(\,j\right)}\right)=\hat{p}_{x}^{\left(\,j\right)}(t)-f\left(\bar{p}_{x1}^{\left(\,j-1\right)}(t),\bar{p}_{x2}^{\left(\,j-1\right)}(t);t_{x,k}^{\left(\,j\right)}\right). \end{array}\end{eqnarray} \tag{ 10 }$$

The function $f\left(\centerdot \right)$ determines the pump profile of a complete revolution, defined by the concatenation of $\bar{p}_{x1}^{\left(\,j-1\right)}(t)$ and $\bar{p}_{x2}^{\left(\,j-1\right)}(t)$ once these two profiles have been time-scaled to fit their respective intervals $\left[\hat{t}_{x,k-1}^{\left(\,j\right)},\hat{t}_{x,k}^{\left(\,j\right)}\right]$ and $\left[\hat{t}_{x,k}^{\left(\,j\right)},\hat{t}_{x,k+1}^{\left(\,j-1\right)}\right]$. Note that $\hat{t}_{x,k-1}^{\left(\,j\right)}$ is the preceding occurrence time from the current iteration and $\hat{t}_{x,k+1}^{\left(\,j-1\right)}$ is the succeeding occurrence time that results from the previous iteration. The occurrence time $t_{x,k}^{\left(\,j\right)}$ of the kth roller period is determined by the time instant that minimizes the least square error

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \hat{t}_{x,k}^{\left(\,j\right)}=\text{arg}~\underset{t_{x,k}^{\left(\,j\right)}}{{\min}}\,{\int}^{}\frac{{{\partial}^{2}}\epsilon _{x}^{2}\left(t;t_{x,k}^{\left(\,j\right)}\right)}{\partial {{t}^{2}}}\text{d}t. \end{array}\end{eqnarray} \tag{ 11 }$$

The reason for minimizing the second derivative of the squared residual, rather than the squared residual itself, is to produce a smoother cardiac signal estimate.

For each roller, the mean period length $\bar{T}_{xq}^{\left(\,j\right)}$ is determined by

$$\begin{eqnarray}\begin{array}{*{35}{l}} \bar{T}_{xq}^{\left(\,j\right)} & =\frac{1}{N}\underset{k=1}{\overset{N}{\sum}}\,\left(\hat{t}_{x,2k+1-q}^{\left(\,j\right)}-\hat{t}_{x,2k-q}^{\left(\,j\right)}\right),\quad q=1,2. \end{array}\end{eqnarray} \tag{ 12 }$$

The pump profile $\bar{p}_{xq}^{\left(\,j\right)}(t)$ is computed as the ensemble median at each time instant of all N revolutions, once the duration of each roller period has been scaled in time, i.e.

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \bar{p}_{xq}^{\left(\,j\right)}(t)=\underset{k}{{\text{median}}}\,\,{{y}_{x}}\left(\hat{a}_{xq,2k}^{\left(\,j\right)}t+\hat{t}_{x,2k}^{\left(\,j\right)}\right),\quad 0\leqslant t<\bar{T}_{xq}^{\left(\,j\right)}, \end{array}\end{eqnarray} \tag{ 13 }$$

where the time scale factor $\hat{a}_{xq,k}^{\left(\,j\right)}$ is estimated by

$$\begin{eqnarray}\begin{array}{*{35}{l}} \hat{a}_{xq,2k}^{\left(\,j\right)} & =\frac{\hat{t}_{x,2k+1-q}^{\left(\,j\right)}-\hat{t}_{x,2k-q}^{\left(\,j\right)}}{\bar{T}_{xq}^{\left(\,j\right)}}. \end{array}\end{eqnarray} \tag{ 14 }$$

The pump model signal $\tilde{p}_{x}^{\left(\,j\right)}(t)$ is then produced by alternating the concatenation of the two pump profiles: for every other roller period, the first pump profile is used,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \tilde{{p}}\,_{x}^{\left(\,j\right)}\left(\hat{a}_{x1,k}^{\left(\,j\right)}t+\hat{t}_{x,2k}^{\left(\,j\right)}\right)=\bar{p}_{x1}^{\left(\,j\right)}(t),\quad 0\leqslant t<\bar{T}{{_{x1}^{\left(\,j\right)}}_{{}}}, \end{array}\end{eqnarray} \tag{ 15 }$$

and, otherwise, the second pump profile,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \tilde{{p}}\,_{x}^{\left(\,j\right)}\left(\hat{a}_{x2,k}^{\left(\,j\right)}t+\hat{t}_{x,2k+1}^{\left(\,j\right)}\right)=\bar{p}_{x2}^{\left(\,j\right)}(t),\quad 0\leqslant t<\bar{T}{{_{x2}^{\left(\,j\right)}}_{{}}}. \end{array}\end{eqnarray} \tag{ 16 }$$

2.3. Convergence

The iterative estimation procedure is considered to have converged when the following criterion is fulfilled:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {\int}^{}|{{{\hat{c}}}^{\left(\,j\right)}}(t)-{{{\hat{c}}}^{\left(\,j-1\right)}}(t)|\text{d}t>{\int}^{}|{{{\hat{c}}}^{\left(\,j-1\right)}}(t)-{{{\hat{c}}}^{\left(\,j-2\right)}}(t)|\text{d}t. \end{array}\end{eqnarray} \tag{ 17 }$$

If the procedure has not converged, the iteration index j is incremented by one and the procedure jumps to section 2.1 and continues processing. If convergence has not been reached within L iterations, the procedure is terminated and the estimates of the Lth iteration are used for subsequent analysis.

The block diagram in figure 2 provides an overview of the estimation of a joint cardiac signal from the signals acquired by the arterial and venous pressure sensors.

Zoom In Zoom Out Reset image size

2.4. Estimation of heartbeat occurrence time

Following convergence, the joint cardiac signal estimate ${{\hat{c}}^{(J)}}(t)$ is subjected to bandpass filtering to facilitate the estimation of the heartbeat occurrence times ${{\hat{t}}_{c,k}}$ and to reduce the number of falsely detected beats; J denotes the number of the iterations required for convergence. The cut-off frequencies of the bandpass filter, ${{F}_{\text{LP}}}$ and ${{F}_{\text{HP}}}$, are chosen such that all heart rates can be accommodated in the passband. The occurrence time of a heartbeat is taken as the time of the mid-amplitude of a heartbeat. The mid-amplitude has been successfully used for estimation of the occurrence time in the PPG signal (Lázaro et al 2014).

Figure 3 illustrates the signals produced at different steps of the present method; the final result being a series of heartbeat occurrence times ${{\hat{t}}_{c,1}},{{\hat{t}}_{c,2}},\ldots,{{\hat{t}}_{c,M}}$, where M denotes the number of estimated heartbeats in the segment.

Zoom In Zoom Out Reset image size

2.5. Estimation of heartbeat amplitude

A clinical study was conducted at Skåne University Hospital, Lund, Sweden. The study, adhering to the Declaration of Helsinki, was approved by the local ethical review board, and all patients signed an informed consent form before participating. A total of 8 patients were included, all suffering from chronic renal failure and all undergoing hemodialyis treatment for at least 3 months prior to the study onset. A treatment lasted from 4 to 5 h and was performed according to the ordinary prescription provided by the responsible nephrologist. The blood flow rate varied between 320 and 500 ml/min. Four patients had a history with one or more heart complications. One patient had a graft as vascular access, whereas all others had fistulas. One patient appears twice in the present dataset, and thus a total of 9 treatments were included. Patients undergoing hemodiafiltration treatment and patients participating in other studies were excluded.

All patients were treated with the Gambro AK 200 hemodialysis machine. An external device with pressure sensors identical to the built-in ones was connected to the extracorporeal blood circuit to acquire the arterial and venous pressure signals. The external device was used instead of the built-in pressure sensors in order to avoid the time-consuming work associated with a software update of a commercial machine and subsequent verification and validation. However, the recorded data can still be regarded as derived form the built-in sensors.

A finger pulse oximeter (LifeSense^™) was used for acquiring the PPG signal. The pulse oximeter signal is proportional to changes in blood volume, and therefore deemed appropriate for use as reference signal. One minute signal segments were excluded whenever the dominant spectral peak did not correspond to a heart rate within the interval 30 to 180 beats per minute (bpm); thus, such segments are not part of the results. In addition, PPG segments with motion artifacts were manually excluded, leading to that 4.9% of all 2280 1 min segments were excluded. Recent techniques may be used to lower this percentage, but were not applied in the present study (Chong et al 2014, Salehizadeh et al 2014).

The performance is evaluated with respect to the accuracy of estimating (i) heart rate and (ii) occurrence time. The following three methods are evaluated:

• (i)  
  the joint cardiac signal estimate $\hat{c}_{\text{AA}}^{(J)}(t)$ resulting from arithmetic averaging,
• (ii)  
  the joint cardiac signal estimate $\hat{c}_{\text{PCA}}^{(J)}(t)$ resulting from PCA-based weighting, and
• (iii)  
  the cardiac signal estimate $\hat{c}_{\text{V}}^{(J)}(t)$ resulting from processing of the venous pressure only, i.e. the method described in Holmer et al (2015) but here extended to include the procedure described in appendix so that complete treatments can be analyzed.

The following parameter values are used: ${{F}_{\text{LP}}}=0.5$ Hz, ${{F}_{\text{HP}}}=5.0$ Hz, T  =  1 min, $\lambda =0.5$ s, L  =  50. Estimation of the cardiac signal is performed in consecutive 1 min segments throughout treatment. All signals are processed at a sampling rate of 100 Hz.

The error between the estimated and the reference heart rate, as well as the error between the estimated and the reference occurrence time, is characterized by the median absolute deviation (MAD). Before the MAD was calculated, the cardiac amplitude range was divided into adjacent intervals of equal length (amplitude bins). The heart rates were, for the complete data set, sorted into the amplitude bins, and the MAD was calculated for each bin; the same was done for the occurrence times.

The performance is illustrated by the examples in figure 4 where heart rate is presented for two complete hemodialysis treatments. In these examples, the heart rate is determined from $\hat{c}_{\text{AA}}^{(J)}(t)$ and compared to the heart rate determined from the reference signal. The pump rate harmonics are also displayed, being well-separated from the heart rate in figure 4(a), whereas they overlap with the heart rate in figure 4(b). The estimated heart rate agrees well with the reference heart rate most of the time during treatment. Gaps in the estimated heart rate are due to the exclusion of segments with large disturbances in the observed pressure signal, see appendix. Gaps in the reference heart rate, on the other hand, are due to the exclusion of segments with disturbances in the PPG signal. The exclusion was performed automatically for the estimated heart rate, whereas manually for the reference heart rate.

Zoom In Zoom Out Reset image size

Figure 5 quantifies the performance by presenting the difference between the estimated and the reference heart rate plotted versus the reference heart rate for all treatments (i.e. a Bland–Altman type of plot). Each point refers to the heart rate of a 1 min segment.

Zoom In Zoom Out Reset image size

For the dataset, the mean error of the heart rate was  −0.20 bpm for AA, −0.28 bpm for PCA and  −6.49 bpm for the venous method. The corresponding standard deviation was 3.4, 4.9 and 15 bpm, respectively, whereas the corresponding MAD was 1.2, 1.8 and 10 bpm, respectively. The standard deviation of the occurrence times was 57 ms for AA, 65 ms for PCA and 109 ms for the venous method, whereas the corresponding MAD was 36, 38 and 60 ms, respectively.

The performance of the three methods for extracting heart rate is presented in figure 6(a), where the error between the estimated and the reference heart rate is displayed versus heartbeat amplitude. The MAD is computed for heart rates associated with different amplitudes. The results show that the MAD related to $\hat{c}_{\text{AA}}^{(J)}(t)$ is lower than that related to $\hat{c}_{\text{PCA}}^{(J)}(t)$ when heartbeat amplitudes are below 1 mmHg. For example, at 0.5 mmHg the MAD is 0.6 bpm when using $\hat{c}_{\text{AA}}^{(J)}(t)$, and 0.7 bpm when using $\hat{c}_{\text{PCA}}^{(J)}(t)$. Moreover, the results show that the method based on $\hat{c}_{\text{V}}^{(J)}(t)$, i.e. the cardiac signal estimated from the venous pressure only, performs much worse at low heartbeat amplitudes, e.g. it has a MAD of 1.3 pbm at 0.5 mmHg, thus demonstrating that the use of both arterial and venous pressure is essential for improving performance.

Zoom In Zoom Out Reset image size

For the three methods, the accuracy of the estimated heartbeat occurrence times is presented versus the heartbeat amplitude in figure 6(b). Similar to the results in figure 6(a), a lower MAD is achieved when both the arterial and the venous pressure are analyzed. The performance is similar when the occurrence times are estimated from either $\hat{c}_{\text{AA}}^{(J)}(t)$ or $\hat{c}_{\text{PCA}}^{(J)}(t)$, although the former yields a slightly lower MAD at very low heartbeat amplitudes.

Figure 7 illustrates variations in heartbeat amplitude which are representative during hemodialysis treatment. The gradual decrease occurring in heartbeat amplitude is likely due to that the patient's blood pressure is decreasing, which is often observed during hemodialysis. The sudden drop in heartbeat amplitude of the venous pressure, occurring after about 100 min, is ascribed to the PCA technique which is rather sensitive to disturbances: When incorrect principal components are used, the method cannot separate the pump and heart signals. For weak heartbeat amplitudes, the use of PCA tends to emphasize the pump remainders rather than the cardiac activity, leading to that no improvement can be observed from one iteration to the next.

Zoom In Zoom Out Reset image size

Figure 8 shows a histogram of the heartbeat amplitudes from all treatments, demonstrating the fact that low heartbeat amplitudes are commonly encountered in hemodialysis patients.

Zoom In Zoom Out Reset image size

The main novelty of the present study is the use of both arterial and venous pressures for modeling the pump signal and estimating the cardiac signal. Apart from enabling cardiac signal estimation at low heartbeat amplitudes, averaging of the arterial and venous pressure signals makes the present method much more general than the one in Holmer et al (2015), since the assumption of identical amplitude and wave morphology of the heartbeats is no longer required. The relaxation of this assumption has particular significance with respect to ectopic beats as their morphology in the pressure signal can be quite different from that of normal beats. Based on annotations of the PPG signal, the data set was found to contain 2000 ectopic beats and 140 000 normal beats.

In contrast to Holmer et al (2015), where a 20 min segment with constant flow rate was selected in each patient and subjected to analysis, the present results derive from the analysis of complete hemodialysis treatments. This means that each treatment contains changes in flow rate as well as pressure disturbances, and, consequently, the results reflect performance under clinical conditions.

The heartbeat amplitude histogram in figure 8 shows that weak amplitudes are common, and thus underlines the importance of developing a method which can estimate the cardiac signal at weak amplitudes. Note that the estimated heartbeat amplitudes are much lower than the normal arterial pulse amplitudes (systole–diastole). The fistula, from where the pulsations enter the extracorporeal blood circuit, is a short circuit between an artery and a vein. Depending on how the fistula has matured, the pulsations in it will be more arterial-like (high amplitude) or venous-like (low amplitude). Additionally, the heartbeat amplitude is further weakened on its way through the extracorporeal blood circuit.

The estimated and reference heart rates of the example in figure 4(a) show excellent agreement most of the time during the 4 h treatment. The agreement is slightly less good for the example in figure 4(b) since the heart rate and the blood pump rate coincide.

The error in heart rate exhibits a number of outliers due to a low heartbeat amplitude, see figure 5; all of these outliers originate from 2 of the 9 treatments. Another cause to outliers is the above-mentioned situation with coinciding heart rate and blood pump rate.

Figure 6(a) shows that the error in the estimated heart rate increases rapidly when the heartbeat amplitude drops below 1 mmHg—a property also observed for the estimated heartbeat occurrence times in figure 6(b). If an automatic alert system for sudden cardiac arrest or venous needle dislodgement is implemented, the error rate at low amplitudes may give rise to false alarms. Although it may not be possible to use the proposed method in all patients, it still offers a significant improvement compared to today's clinical practice, as no cardiac monitoring is employed at all.

When calculating the MAD of the heartbeat occurrence times, see figure 6(b), missed and falsely detected beats must first be excluded. On the other hand, such exclusion is not needed when computing the MAD of the heart rate, thus missed and falsely detected beats are included and, accordingly, will influence the result. This observation may explain why, at low heartbeat amplitudes, the difference in MAD between using AA and PCA is less pronounced in figure 6(b) than in figure 6(a). A MAD of 20 mmHg for the heartbeat occurrence times suggest that the accuracy of the timing is not good enough for spectral HRV analysis, but may be good enough for detecting ectopic beats. The pressure signal is obviously not suited for more advanced arrhythmia analysis since it does not convey morphologic information on atrial and ventricular waveforms.

In the present study, the roller period was estimated from the observed pressure signal. When implementing the present method in a hemodialysis machine, it may be preferable to measure the roller period directly from a signal which reflects the rotational speed of the blood pump, possibly replacing the roller period estimation described in appendix.

A limitation of the present study is that the dataset was used for both development and performance evaluation. Thus, another study should be undertaken which involves a larger dataset, so that performance can be better established. None of the parameters which influence performance, i.e. ${{F}_{\text{LP}}},{{F}_{\text{HP}}},L,\lambda$, and α, were, however, subject to optimization, but their values were chosen based on general considerations. For example, L was set to 50 since a higher value would only lead to increased computation time without leading to estimates with better accuracy, and ${{F}_{\text{LP}}}$ and ${{F}_{\text{HP}}}$ were related to the frequency content of interest in the cardiac signal.

Another limitation of the study is that the ECG was not recorded when the database was acquired, although the ECG is a gold standard for analyzing the timing of heartbeats. On the other hand, the PPG was recorded for the purpose of monitoring heart rate and therefore used in the present study; Gil et al have shown that the PPG is a useful surrogate for the ECG when analyzing heart rate variability (Gil et al 2010).

The use of both arterial and venous pressure from the extracorporeal sensors of the hemodialysis machine is found to improve the estimation of heart rate and heartbeat occurrence times at low heartbeat amplitudes, when compared to using the venous pressure only. As a result, cardiac monitoring can be performed in more patients. The new method makes it possible to estimate a cardiac signal also at low blood flow rates, and to automatically handle changes in blood flow, implying that cardiac information can be reliably extracted throughout complete treatments. The results are preliminary, and need to be verified on a separate dataset.

The authors would like to thank Mårten Segelmark, Lena Mattsson, and Sarok Said at Lund University Hospital for all work with the clinical investigation. This work was partially supported by grants from the Swedish Research Council. M Holmer, B Olde, and K Solem are employees of Baxter International Inc.. A patent application related to the work is pending.

Since the flow rate recorded by the machine has insufficient accuracy for the present application, the rotor period is estimated from the observed signals y_x(t). Thus, the procedure for cardiac signal estimation is fully automated.

For each signal segment, initialization starts by determining a rough estimate of the roller period ${{\hat{T}}_{x}}$ using the entire segment. Then, based on ${{\hat{T}}_{x}}$, estimates of each roller period ${{\hat{T}}_{x,k}}$ is determined in order to reflect the fact that the period varies over time.

The rough estimate of the roller period ${{\hat{T}}_{x}}$ is determined by finding the peak in the autocorrelation function of y_x(t) which corresponds to the shortest roller period in the interval ${{\tau}_{1}}<{{T}_{x}}<{{\tau}_{2}}$, where ${{\tau}_{1}}$ and ${{\tau}_{2}}$ are defined by the flow range of the pump. This interval is chosen such that all adequate blood flow rates are covered. At low flow rates, small peaks may occur in the autocorrelation function at shorter lags than the roller period and, therefore, a threshold α is used to differentiate such peaks from the peak corresponding to the roller period. If no peak exceeds α, a blood flow rate change or a large pressure disturbance has likely occurred, and, therefore, the analyzed signal segment is discarded from subsequent analysis.

The estimation of each separate roller period ${{\hat{T}}_{x,k}}$, starts by estimating the occurrence time of the first roller actuation time,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \hat{t}_{x,0}^{(0)}=\text{arg}~\underset{0<t<\left(1+\delta \right){{T}_{x}}}{{\max}}\,|\frac{\text{d}{{y}_{x}}(t)}{\text{d}t}|, \end{array}\end{eqnarray} \tag{ A.1 }$$

where δ defines the maximal variation allowed in period time between consecutive roller waves. The absolute value in (A.1) ensures that maximization can be performed on both arterial and venous signals, since the maximum occurs at positive derivative for arterial signals, whereas negative for venous signals. Starting from $\hat{t}_{x,0}^{(0)}$, the roller period $\hat{T}_{x,k}^{(0)}$ for all revolutions in the segment are estimated as

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \hat{T}_{x,k}^{(0)}=\text{arg}~\underset{\left(1-\delta \right){{T}_{x}}<T<\left(1+\delta \right){{T}_{x}}}{{\max}}\,|\frac{\text{d}y\left(\hat{t}_{x,k-1}^{(0)}+T\right)}{\text{d}t}|, \end{array}\end{eqnarray} \tag{ A.2 }$$

where

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \hat{t}_{x,k}^{(0)}=\underset{l=1}{\overset{k}{\sum}}\,\hat{T}_{x,l}^{(0)}+\hat{t}_{x,0}^{(0)},\quad k=1,\ldots,2N. \end{array}\end{eqnarray} \tag{ A.3 }$$

The pump model signal $\tilde{p}_{x}^{(0)}(t)$ is then produced by the procedure described in section 2.2, except that the estimated occurrence times described in Appendix should be inserted in equations (13)–(16), using j  =  0.

The following parameter values were used for initialization: ${{\tau}_{1}}=0.5$ s, ${{\tau}_{2}}=1.3$ s, $\alpha =0.04$, and $\delta =0.25$.