A novel particle filtering method for estimation of pulse pressure variation during spontaneous breathing

Background We describe the first automatic algorithm designed to estimate the pulse pressure variation (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {PPV}$$\end{document}PPV) from arterial blood pressure (ABP) signals under spontaneous breathing conditions. While currently there are a few publicly available algorithms to automatically estimate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {PPV}$$\end{document}PPV accurately and reliably in mechanically ventilated subjects, at the moment there is no automatic algorithm for estimating \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {PPV}$$\end{document}PPV on spontaneously breathing subjects. The algorithm utilizes our recently developed sequential Monte Carlo method (SMCM), which is called a maximum a-posteriori adaptive marginalized particle filter (MAM-PF). We report the performance assessment results of the proposed algorithm on real ABP signals from spontaneously breathing subjects. Results Our assessment results indicate good agreement between the automatically estimated \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {PPV}$$\end{document}PPV and the gold standard \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {PPV}$$\end{document}PPV obtained with manual annotations. All of the automatically estimated \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {PPV}$$\end{document}PPV index measurements (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {PPV}_{\text {auto}}$$\end{document}PPVauto) were in agreement with manual gold standard measurements (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {PPV}_{\text {manu}}$$\end{document}PPVmanu) within ±4 % accuracy. Conclusion The proposed automatic algorithm is able to give reliable estimations of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {PPV}$$\end{document}PPV given ABP signals alone during spontaneous breathing.


Methods: algorithm description
The subsequent sections explain a novel statistical signal model for ABP signals recorded from spontaneously breathing subjects and the PPV index tracking algorithm utilizing our recently developed sequential Monte Carlo estimation method.
(1) PPV (%) = 100 × PP max − PP min (PP max + PP min )/2 Notation We have adopted the notation used in [16] with minor modifications. We use boldface to denote random processes, normal face for deterministic parameters and functions, upper case letters for matrices, lower case letters for vectors and scalars, superscripts in parenthesis for particle indices, upper-case superscripts for nonlinear/linear indication, and subscripts for time indices. For example, the nonlinear portion of the state vector for the ith state trajectory (i.e., particle) is denoted as x N,(i) n where n represents the discrete time index and (i) denotes the ith particle. The unnormalized particle weights are denoted as w (i) and the normalized particle weights as w (i) . The state trajectories before resampling are denoted as x (i) n and as x (i) n after resampling.

State-space model
The proposed automatic PPV index estimator utilizes our recently developed sequential Monte Carlo estimation method which is based on the state-space model approach. The state-space model is a mathematical expression to describe the evolution of any physical system's unobservable state x n and its relation to measurement y n , where the state x n is a vector of parameters representing the system's condition. The state-space method is a technique to estimate the state x n as a function of measurement y n utilizing the statespace model. The typical state-space model can be expressed as, where (2) is a process model, (3) a measurement model, f (·) and h(·) functions of the state x n , and u n and v n uncorrelated white noises with variances q and r. A designer needs to incorporate prior domain knowledge of a system into the state-space model and define the functions f (·) and h(·). The flexibility and versatility of the state-space method are attributable to two functions, which can be either linear or nonlinear.

Measurement model
The measurement model of the ABP signal is shown in (4)(5)(6)(7), where γ n is the respiratory signal, µ n the amplitude-modulated cardiac signal, ρ k,n the amplitude modulation factor of the kth cardiac harmonic partial, κ k,n the kth cardiac harmonic partial, θ r n the instantaneous respiratory angle, θ c n the instantaneous cardiac angle, N r h the number of respiratory partials, N c h the number of cardiac partials, v n the white Gaussian measurement noise with variance r, and r ·,k,n , m ·,k,j,n , c ·,k,n the sinusoidal coefficients. This measurement model was first introduced in [17].
r 1,k,n cos kθ r n + r 2,k,n sin kθ r n Page 4 of 18 Kim et al. BioMed Eng OnLine (2016) 15:94 Process model The process model describes the evolution of each element of the state x n . In our application, x n includes the instantaneous respiratory and cardiac angles θ r n and θ c n , the instantaneous mean respiratory f r n and cardiac f c n frequency, the instantaneous respiratory f r n and cardiac f c n frequencies, and the sinusoidal coefficients {r 1,k,n , r 2,k,n , c 1,k,n , c 2,k,n , m 1,k,n , m 2,k,n }. The process model can be expressed as, where f r n is the instantaneous respiratory frequency, f c n the instantaneous cardiac frequency, T s the sampling period, f r n the instantaneous mean respiratory frequency, f c n the instantaneous mean cardiac frequency, α the autoregressive coefficient, and u r,n , u c,n , and u m,n the process noises with variances q r , q c , and q m . The clipping function g[·] limits the range of instantaneous mean frequencies, which can be written as, The range of instantaneous mean frequencies, i.e., f r n and f c n , is assumed to be known as domain knowledge.

Maximum A-posteriori marginalized PF
The proposed automated PPV index estimation method requires accurate estimates of the instantaneous respiratory frequency f r n , the instantaneous cardiac frequency f c n , (6) ρ k,n = 1 + N r h j=1 m 1,k,j,n cos jθ r n + m 2,k,j,n sin jθ r n (7) κ k,n = N c h k=1 c 1,k,n cos kθ c n + c 2,k,n sin kθ c n + u f c ,n (14) r ·,k,n+1 = r ·,k,n + u r,n (15) c ·,k,n+1 = c ·,k,n + u c,n (16) m ·,k,n+1 = m ·,k,n + u m,n (17) and the morphology of an ABP signal. In order to obtain those estimates, we utilize our recently developed particle filter technique, which is called the maximum a-posteriori adaptive marginalized particle filter (MAM-PF). The MAM-PF is a hybrid version of the marginalized particle filter (MPF) and maximum a-posteriori particle filter (MAP-PF), which leverages the advantages of the MPF and MAP-PF. In [18] we described the recursions for the MAM-PF in detail. We proposed two versions of the MAM-PF: optimal and fast MAM-PFs [18]. Within the state-space method framework, the Optimal MAM-PF computes the "optimal" trajectory of the state x n . However, its computational burden is too demanding to be practically useful. The fast MAM-PF is an approximation of the optimal MAM-PF, which requires dramatically less computational burden. However, the fast MAM-PF performs as well as the optimal MAM-PF, which we demonstrated in [8].
Recently, we proposed an automatic (PPV) estimation technique in mechanically ventilated patients by utilizing the fast MAM-PF as an ABP signal tracker [8]. Under full mechanical support, the respiratory rate of subjects is equal to the mechanical ventilation rate, which is known and constant. Therefore, the fast MAM-PF has to track only the instantaneous cardiac frequency f c n along with the signal morphology. All ABP signals included in this study were recorded from spontaneously breathing subjects. Therefore, the ABP signal tracker has to track both the instantaneous respiratory frequency f r n and the instantaneous cardiac frequency f c n along with the signal morphology. Although the fast MAM-PF based ABP signal tracker is capable of tracking multiple frequencies, there are two major issues in using the fast MAM-PF algorithm as the ABP signal tracker for ABP signals of spontaneously breathing subjects. The first issue is that the morphology of the signal, which is represented by the sinusoidal coefficients in (6, 7), does not belong to the linear state any more. Since the modulating signal ρ k,n is multiplied to the cardiac signal κ k,n , their sinusoidal coefficients c ·,k,n and m ·,k,j,n are nonlinear parameters of the measurement model in (4). The fast MAM-PF is applicable only to state-space models whose state vector can be partitioned into the linear and nonlinear portions. The second issue is that as the dimension of the state, where particle filters are used, increases the number of necessary particles to cover the state increases exponentially. As a result, the computational burden of the fast MAM-PF increases exponentially. The portion of the state space where particle filters are used is called the particle space. Since the new ABP signal tracker has to estimate both the instantaneous respiratory frequency f r n and the instantaneous cardiac frequency f c n , the dimension of the particle state becomes 2, which results in a quadruple increase of computational burden if the fast MAM-PF has to be used for the current application. In order to address these two major issues we propose a new ABP signal tracker, which is a modified version of the Fast MAM-PF. It is called, the Dual MAM-PF. The term "Dual" is borrowed from Dual Kalman filters, in which the state is divided into two portions and each portion is estimated separately assuming that the other portion is known and equal to the currently estimated value. While the fast MAM-PF treats a two-dimensional particle space as a whole, the dual MAM-PF partitions the two-dimensional particle space into two one-dimensional particle spaces assuming independence between two particle space variables, which are the instantaneous respiratory frequency f r n and the instantaneous cardiac frequency f c n . Suppose that the state vector x can be partitioned as follows, where x P n represents the particle state and x K n the Kalman state. The particle state is the portion of the state where particle filters are used as defined earlier while the Kalman state is the portion of the state where extended Kalman filters are used. The state variables whose posterior distributions are known to be multi-modal belong to the particle state while those whose posterior distributions are known to be Gaussian or uni-modal belong to the Kalman state. In [18] we demonstrated that the posterior distribution of the instantaneous frequency of a multi-harmonic signal is truly multi-modal. Given the state-space model in (4)-(7), instantaneous respiratory frequency f r n and the instantaneous cardiac frequency f c n are the particle state variables and the sinusoidal coefficients such as r ·,k,n , c ·,k,n , and m ·,k,j,n are the Kalman state variables. Assuming that the particle state variables are independent of each other the particle state x P n can be partitioned further as, where x P 1 n and x P 2 n represent the first and second particle state variables, respectively. This partitioning breaks down a two-dimensional particle space x P n into two one-dimensional particle spaces. The total posterior distribution is given by, Algorithm 1 provides a complete description of the dual MAM-PF algorithm, where N T represents the total number of signal samples, N p the number of particles for each onedimensional particle space, j the particle state variable index, and i j the particle index of the jth particle state variable. The total number of particle used in the dual MAM-PF algorithm is 2N p instead of N 2 p . At each time step n the dual MAM-PF searches for the (22) p(x 0:n |y 0:n ) = p(x K 0:n |y 0:n , x P 0:n )p(x P 0:n |y 0:n ) (23) = p(x K 0:n |y 0:n , x P 0:n )p(x P 1 0:n |y 0:n )p(x P 2 0:n |y 0:n ). best trajectory of each particle i j from the previous trajectory. This searching step can be written as, Given the best trajectory for each particle i j , corresponding Kalman state variables x P j ,(i j ) , i.e. sinusoidal coefficients, are updated utilizing the extended Kalman filter. Then, the MAP estimate of x P i n is obtained based on the value of the coefficient α j,n as follows, Since there are two groups of particles i 1 and i 2 , we need to select the best estimate of the Kalman state vector x K n among two potential estimates: . The actual estimate of the Kalman state vector x K n can be selected as follows, Then, the estimate of the entire state x n can be expressed as, In order to appreciate the algorithm of the dual MAM-PF, it is essential to understand the generic particle filter along with other variants of particle filters such as the MPF, MAP-PF, and MAM-PF. We provided detail algorithms of those particle filters in [18]. (24)

ABP signal envelope estimation
Given the estimated signal parameters in (8)(9)(10)(11)(12)(13)(14)(15)(16), it is possible to estimate the upper envelope (e µ,n ) and lower envelope (e ℓ,n ) of ABP signals by following steps below, where arg max x f (x) and arg min x f (x) are operators to obtain the value of x for which f(x) attains its maximum and minimum values, respectively. The top plot in Fig. 1 shows a five respiratory cycle period of an ABP signal y n (thick red), its estimate ŷ n (thin green), and its estimated envelopes e µ,n and e ℓ,n (blue), which are described in (32) and (33).

Pulse pressure signal envelope estimation
The pulse pressure (PP) signal is the difference between the upper envelope e µ,n and lower envelope e ℓ,n of the ABP signal. This PP signal oscillates roughly at the respiratory rate as shown in the bottom plot in Fig. 1. This oscillation is due to the respiratory effect on the variation of systemic ABP under full ventilation support [19]. Within each respiratory cycle PP reaches its maximum (PP max ) and minimum (PP min ) values, which are two critical parameters to compute the PPV index. Traditionally, the PP max and PP min values have been computed only once per each respiratory cycle. Given  Fig. 1 Top Original ABP signal (red) and its estimate (green) with automatically computed envelopes (blue). Bottom Automatically computed PP signal (red) and its envelope (blue) the estimated signal parameters in (8)(9)(10)(11)(12)(13)(14)(15)(16), however, we can compute the continuous equivalents of PP max and PP min . They are the upper ε µ,n and lower ε ℓ,n envelopes of the PP signal. The upper envelope ε µ,n is the continuous estimate of PP max and the lower envelope ε ℓ,n that of PP min . The ε µ,n and ε ℓ,n values can be estimated as described below, where 1 + ̺ k,n is equal to ρ k,n and ε µ,n and ε ℓ,n are the continuous estimates of the PP max and PP min , respectively. The blue lines in the bottom plot in Fig. 1 represent the upper ε µ,n and lower ε ℓ,n envelopes of the PP signal, which are obtained by following the method described above.

Pulse pressure variation calculation
Given the ε µ,n and lower ε ℓ,n values, it is straightforward to calculate the PPV index. It can be computed as follows, This new PPV index is different from the traditional PPV index described in (1) because the new one is continuous in time while the traditional one can be obtained only once per each respiratory cycle. Figure 2 illustrates an example of the automatically computed continuous PPV index (thick green) and the manually obtained discrete PPV index (thin red) of a real 10 min 1 + ̺ min,k,n κ max,k,n − κ min,k,n (36) PPV (%) = 100 × ε max − ε min (ε max + ε min )/2 ABP signal. Each hollow white dot represents a "discrete" PPV index, which can be obtained once per each respiratory cycle.
The subsequent sections describe how to assess the accuracy of the proposed PPV index tracking algorithm.

Assessment data
The Massachusetts General Hospital waveform database (MGHDB) on PhysioNet is a comprehensive collection of electronic recordings of hemodynamic and electrocardiographic waveforms patients in critical care units [14,15]. It consists of recordings from 250 patients representing a broad spectrum of physiologic and pathophysiologic states. The typical recording includes arterial blood pressure (ABP) signals in addition to seven other types of waveforms. By visually inspecting the spectrogram and time-series of ABP signals we identified 11 patients who breathed spontaneously. The first column in Table 2 lists the patients' identification numbers (e.g. mgh000) as in MGHDB on Physionet. Figure 3 shows the spectrogram of one of the 11 ABP signals. Each ABP signal is 10 min long and the total duration of the 11 ABP signals was 2 h. The original sample rate f s of the signals was 360 Hz, but they were downsampled by a factor of 9, so that the final sample rate f s was 40 Hz.
The number of cardiac components N c was 5 and that of respiratory components N r was 2. The total number of particles 2N p was 500. Table 1 lists the parameter values used for the PPV index estimator. Those parameter values were initialized and tuned based on the previously published work, where ABP signals were recorded during full mechanical ventilation [8]. One PPV index measurement is computed from each measurement window, which is a time period of five respiratory cycles

Manual PPV annotations (current standard)
We manually annotated the peaks and troughs of the ABP signals and calculated the PPV indices (current standard) as defined in (1). They are referred to as manual PPV indices PPV manu . PPV auto represents PPV indices obtained using the proposed PPV index tracking algorithm.

Statistical analysis
The statistical analysis used five PPV index measurements for each subject, and each measurement was separated by 2 min. Each PPV index measurement is an averaged value over 5 respiratory cycles. Figure 2 shows the 2 min apart measurement periods.  The proposed PPV index tracking algorithm was assessed by calculating the agreement (mean ± standard deviation) between PPV auto and PPV manu measurements and using Bland-Altman analysis.
A Bland-Altman plot is a statistical and visualization method that is often used in the assessment of PPV estimation algorithms in order to determine the agreement between two different PPV estimates. It has the difference PPV between PPV auto and PPV manu on the y-axis and the PPV manu on the x-axis. It visualizes the overall accuracy of estimation and estimation bias or trend versus PPV manu . We used it to compare the current standard using manual annotations with our automatic estimation algorithm. Figure 4 depicts the Bland-Altman plot of the 11 subjects. There are 5 PPV measurements available per each subject. All of PPV auto measurements were in agreement with PPV manu measurements within ±3.5 % accuracy. Table 2 lists the mean ± standard deviation of 5 PPV manu and PPV auto measurements for each subject. The second column is for PPV manu and the third column for PPV auto .

Frequency clipping function
The clipping function g[·] in (17) could be defined as follows, However, there is a major problem with the clipping function in (37) when it is incorporated into the particle filter framework. It tends to cause the boundary values, f max and f min , to be the pitfalls for the particles, x P 1 n and x P 2 n , to be trapped in when the instantaneous frequency values associated with the particles become close to f max or f min . In other words, once any particle's frequency value meets either of the boundary conditions, the particle tends to remain in the same state having the frequency value of either f max or f min . In order to address this issue, the clipping function is defined as shown in (17), which forces the frequency value to bounce back within the range, i.e., f min < f ≤ f max , once it reaches beyond f max and f min .

Algorithm's advantages
The proposed algorithm is the first automatic method described in the literature especially designed to estimate and track the PPV index in situations involving spontaneous breathing. It is important to note that the proposed algorithm is a complete new design from our previously described algorithm [20] which only worked for mechanically ventilated subjects. Our previous algorithm was made publicly available by the authors and due to its performance has been adopted by Philips Medical Systems. Currently, our previously published PPV algorithm is displayed in real-time on the Philips Intelliveu MP70 monitors (Intellivue MP70, Philips Medical Systems) and has been used in numerous studies related to PPV and fluid responsiveness. Its ability to monitor fluid responsiveness in the operating room and its accuracy against the current standard obtained by manual annotations were assessed by Cannesson [21]. Previously it was not possible to conveniently monitor the PPV index in the operating room or in the intensive care unit because it had to be manually calculated. Thus, the automatic PPV has potential clinical application for fluid management optimization in the operating room.
A limitation of our previously described algorithm adopted by Philips in their Intelliveu MP70 monitors is that it may not work adequately in regions of abrupt hemodynamic changes [20] and it is only accurate for mechanically ventilated subjects. In this paper, we provide a detailed description of a novel algorithm designed to be a robust PPV estimator during regions of abrupt hemodynamic changes and during spontaneous breathing.
The major algorithm design difference of the proposed algorithm with respect to previously published algorithms [20,22] is the fact that the proposed method is based on a statistical state-space capable of modeling spontaneous breathing, and estimation of the cardiovascular pressure signal based on this statistical model using optimal estimation methods. The state-space modeling stage results in an algorithm that is more robust to hemodynamic changes and artifacts. The statistical state-space signal model and associated model parameter estimation algorithm automatically filter out noise and artifact that cannot be captured with the model. Since the statistical signal model is based on cardiovascular physiology and pathophysiology, signal features that are not physiological in nature are automatically filtered out. Additionally, the model is general enough to accurately model both arterial blood pressure signals and plethysmogram signals. Consequently, it can also be used to calculate the pleth variability index (PVI). Figures 5 and 6 exemplify a case where signal features that are not physiological in nature are automatically filtered out resulting in more accurate PPV index estimation than manual annotation. The top plot in Fig. 5 illustrates 4 respiratory cycles of the ABP signal (red) and its estimate (green). It also shows the manually annotated signal envelopes (black) and the automated computed signal envelopes (light blue). The bottom plot in Fig. 5 depicts the PPV manu and PPV auto over the same period. Around 535 s , the PPV manu value (red) abruptly increases up to 35% while the PPV auto value (green) remains at 8 %. Around 540 s, the PPV manu value returns to 8 %. Figure 6 focuses on the time period marked with the black rectangular box in Fig. 5. The top plot in Fig. 6 shows that the heart beat between 535 and 535.5 s is contaminated by noise and has an abnormal morphology. As a result, the corresponding PP manu shown in the bottom plot reaches a large maximum value (PP min,manu : 105 mm Hg) around at 535 s. However, the automatically computed maximum PP value (PP min, auto ) at the same time is as low as 83 mm Hg. This discrepancy between the manual annotation and the proposed automatic method results from the capability of the MAM-PF algorithm, which estimates the ABP signal based on the state-space model. While the original heart beat between

Fig. 5 Top
Original ABP signal (red) with its manually annotated envelopes (black) and signal estimates (green) with its automatically computed envelopes (purple). Bottom Manual PPV indices (red) and automatic PPV indices (green) where one of the manual PPV indices has an abnormally high value 535 and 535.5 s in Fig. 6 is abnormal in a physiological sense, the estimated heart beat over the same time period shows the physiologically expected morphology and location of the heart beat.

Study limitations
The algorithm's assessment was based on only 11 subjects with pre-recorded ABP data. Additionally, for each subject five PPV estimates were used in the assessment study. This assessment was designed to be an engineering algorithm validation against current standard manual annotations, and not a clinical validation study. Consequently, a clinical validation study assessing the ability of the proposed algorithm to monitor fluid responsiveness in the operating room in situations involving spontaneously breathing subjects still needs to be conducted. This may require the proposed algorithm to be first adopted as part of a commercial system as was the case with our previous automatic PPV algorithm [20].

Conclusion
We have described the first automatic PPV tracking algorithm for spontaneously breathing subjects. This novel algorithms is based on a statistical state-space model inspired in the underlying cardiovascular and respiratory physiology. This algorithm uses our recently developed SMCM (MAM-PF) for optimal parameter estimation. The assessment results indicate good agreement against the current standard PPV. The algorithm was designed to work during regions of abrupt hemodynamic changes and spontaneous breathing. All of PPV auto measurements were in agreement with PPV manu measurements within ±4 % accuracy.  Fig. 6 Top Original ABP signal (red) with its manually annotated envelopes (black) and signal estimates (green) with its automatically computed envelopes (purple). Bottom Manual PP signal (red) and automatic PP signal (green) where the manual PP signal increases momentarily due to a contaminated heart beat in the ABP signal