Robust Peak Recognition in Intracranial Pressure Signals

Abstract

Background

The waveform morphology of intracranial pressure pulses (ICP) is an essential indicator for monitoring, and forecasting critical intracranial and cerebrovascular pathophysiological variations. While current ICP pulse analysis frameworks offer satisfying results on most of the pulses, we observed that the performance of several of them deteriorates significantly on abnormal, or simply more challenging pulses.

Methods

This paper provides two contributions to this problem. First, it introduces MOCAIP++, a generic ICP pulse processing framework that generalizes MOCAIP (Morphological Clustering and Analysis of ICP Pulse). Its strength is to integrate several peak recognition methods to describe ICP morphology, and to exploit different ICP features to improve peak recognition. Second, it investigates the effect of incorporating, automatically identified, challenging pulses into the training set of peak recognition models.

Results

Experiments on a large dataset of ICP signals, as well as on a representative collection of sampled challenging ICP pulses, demonstrate that both contributions are complementary and significantly improve peak recognition performance in clinical conditions.

Conclusion

The proposed framework allows to extract more reliable statistics about the ICP waveform morphology on challenging pulses to investigate the predictive power of these pulses on the condition of the patient.

1 Background

Traumatic Brain Injuries (TBI) affect more than 2 million people annually in the United States, and their incidence in the world keeps increasing [1]. The treatment of TBI patients in critical care units, as well as other neurological disorders, relies on the continuous measurement of intracranial pressure (ICP) (i.e. the sum of the pressures exerted within the craniospinal axis system). It is known that the management of ICP can attenuate secondary brain injuries and improve chances of recovery. Interestingly, the morphology of ICP waveform holds essential informations about the intracranial adaptive capacity (elastance), and even the outcome of head injured patients [2, 3]. For example, it has been shown that variations of the ICP pulse morphology are linked to the development of intracranial hypertension [4–6], cerebral vasospasm [7], changes in the cerebral blood carbon dioxide (CO2) levels [8, 9], decreased cerebral bloob flow (CBF) [10], and changes in the craniospinal compliance [11].

The extraction of morphological features is essential to monitor and to understand ICP in an automatic fashion with the ultimate goal of improving the treatment of pathophysiological intracranial and cerebrovascular conditions. Although ICP pulses are typically triphasic [8] (i.e. three peaks), their shape can exhibit irregular variations such that some peaks may be missing. The recognition of these top peaks is a challenging task that has recently drawn special attention from different research groups. Several algorithms have been developed to detect the first peak [12], and to recognize the three peaks of ICP pulses [13–17]. Existing methods can be divided into two categories depending if they work offline, like Morphologram [14], or online, like MOCAIP [13, 15] (M o rphological C lustering and A nalysis of I CP P ulse). These techniques offer a satisfactory accuracy to recognize the peaks in general cases. However our recent observations show that their performance deteriorates significantly when the pulses exhibit abnormalities or are simply more challenging (a pulse is considered to be challenging if any of its peaks fails to be correctly designated by the baseline MOCAIP algorithm [15], see Figure 1). Such ICP pulses are of particular interest because we suspect that they might hold essential predictive information about the patient condition.

This paper investigates how to improve peak recognition accuracy on challenging ICP pulses. The contribution is two-fold. First, to conduct this study, MOCAIP++, a robust ICP pulse processing framework that generalizes MOCAIP, is introduced. The strength of MOCAIP++ relies on its capacity to integrate different peak recognition methods, and to exploit additional features based on the derivatives of the ICP signal. Our experiments evaluate these characteristics by providing a comparative analysis of three different state-of-the-art peak recognition techniques based on Gaussian Models, Gaussian Mixture Model (GMM), and Spectral Regression Analysis (SR), and by evaluating the impact of ICP features, such as curvature, first and second derivatives on the recognition performance. Second, this paper investigates the effect of incorporating challenging pulses into the training set of peak recognition methods learned in a supervised way. A method is proposed to sample automatically a representative challenging dataset of ICP pulses from a large database of ICP signals collected from 128 neurosurgical patients. The original, and the challenging datasets allow to study how the performance of peak recognition methods, essential to extract morphological features, can be improved.

Figure 1

Illustration of two ICP pulses (the actual position of the peak is depicted in green, MOCAIP prediction in black). On the left, an ICP dominant pulse is correctly annotated with the position of the three peaks. On the right, the automatic annotation failed to correctly recognized the third peak because of the uncommon shape of the pulse. This pulse is considered as a challenging one in our study.

2 Methods

2.1 ICP Dataset

Generally, ICP signal recordings consist of several hours long segments. By reviewing those files, we observed that the majority of the recordings contain pulses whose peaks are easily recognized by automatic algorithms. A subset of ICP files, however, contains pulses that are not correctly annotated by automatic peak recognition methods. One reason for these mismatches is that the pulse morphology differs significantly from the most common ones. We consider those pulses to be challenging (an example is shown in Figure 1).

The variation in morphology of these challenging pulses might originate from a combination of external factors such as sampling rate of the ICP, noise and artifact due to the acquisition device, or coughing of the patient. It is also possible that some of these morphological variations come from the condition of the patient and that they might hold relevant predictive information. Unfortunately, the peak recognition accuracy of current techniques on the ICP recordings containing those pulses drops dramatically. It is no longer possible to extract reliable statistics about their ICP waveform morphology to perform further analysis. This observation led us to extract a challenging dataset D' (Section 2.1.2) from the dataset D (Section 2.1.1). The new dataset D' is sampled from the recordings of D that contain a large percentage of challenging pulses. Both datasets, that are described below, will be used in the experiments to evaluate the performance of our framework. In addition, we will investigate if the use of the challenging dataset as part of the training set of peak recognition methods can improve their performance.

2.1.1 Original Data

The source dataset of ICP signals originates from patients admitted to the University of California Los Angeles (UCLA) medical center. Its usage in the present study was approved by the UCLA Internal Review Board. It is a large, representative dataset that is reasonably distributed across gender, age, and type of patient (ICU or NON-ICU). A small portion of this dataset was previously used to evaluate MOCAIP [15] and its extensions based on regression analysis [18]. The ICP and ECG signals were acquired from 128 patients treated for various intracranial pressure relted conditions. ICP was monitored continuously using Codman intraparenchymal microsensors (Codman and Schurtleff, Raynaud, MA) placed in the right frontal lobe. ICP signals were recorded from bedside monitors using corporate data acquisition systems at a sampling frequency of either 240 Hz or 400 Hz. A total of 1425 recordings were extracted, each totalizing several hours. Those ICP and ECG signal recordings were subsequently pre-processed by MOCAIP so that they were first divided into 3 minutes segments. Then, a hierarchical clustering was applied on individual pulses of each segment, and the center of the dominant cluster was extracted to produce a dominant pulse. This clustering process leads to a representative set of 87,125 dominant pulses. It is referred to as the original dataset D from which a smaller, but more challenging dataset will be sampled. The actual positions of the three peaks in the ICP are obtained by manual annotation from experienced researchers following the procedure described in the next subsection.

2.1.2 Challenging Data

The selection of a challenging subset of ICP pulses D' ⊂ D is achieved using a weighted sampling procedure from the file recordings of the original dataset D. Intuitively, the sampling aims at extracting more pulses from recordings that contain a larger percentage of challenging pulses so that they are better represented in the resulting dataset. To do so, each recording is associated with a weight corresponding to its degree of difficulty which is high if MOCAIP often fails to recognize the three peaks. The procedure to weight the files is described below.

Experienced researchers establish the groundtruth by manually setting the positions of the three peaks (p ₁, p ₂, and p ₃) in each pulse. The task of the researcher is to pick the right peak candidates among those automatically detected at curve inflections (Section 2.2.2). Whenever one of the three peaks is missing, its position is labelled with the empty set. Among the set of pulses, 7173 have missing p ₁, 3699 have missing p ₂, and 4626 have missing p ₃. Researchers cross-validate their results and, if necessary, they harmonize them using the annotation of the previous and following pulses as reference. For a few difficult cases where the researchers could not agree on the position of some peaks, the pulse was removed from the dataset. This procedure ensures that the groundtruth is not biased to a specific researcher.

In parallel, MOCAIP is applied to annotate each pulse with the position of the three peaks (p ₁, p ₂, and p ₃). To find difficult files, the predictions of MOCAIP are compared with the manual groundtruth. For each ICP file f _{i = 1...F} , a weight w _{i = 1...F} is set proportional to the percentage of wrongly assigned peaks. This is done by comparing the position of each peak of the ground truth to the position obtained from the automatic files,

w_i=\frac{{ℰ}_{p1}+{ℰ}_{p2}+{ℰ}_{p3}}{{\mathcal{N}}_{p1}+{\mathcal{N}}_{p2}+{\mathcal{N}}_{p3}},

(1)

where ℰ _{p 1}, ℰ _{p 2}, ℰ _{p 3}are the number of wrongly assigned peaks and {\mathcal{N}}_p{\text{1}}_{,}{\mathcal{N}}_p{\text{2}}_{,}{\mathcal{N}}_{p\text{3}} are the total number of occurrences in the file of the peaks p ₁, p ₂, and p ₃, respectively. The distribution of the weights is illustrated in Figure 2.

Figure 2

The distribution of the weights is illustrated for the files in the original dataset of ICP pulses. The weight of each file is set proportional to the percentage of wrong peak assignations. For example, a value of 1 indicates that all the peaks in that file were not assigned correctly by the MOCAIP, this usually happens in short recordings.

Finally, the challenging dataset D' is created by extracting pulses using weighted sampling, such that a pulse has a probability v _i (Eq. 2) to be picked from file f _i . Therefore, files with large probability v _i will contribute to more pulses in the sampled dataset.

v_i=\frac{w_i}{{\displaystyle \sum _{i=1}^F{w}_i}}

(2)

To avoid redundancy from the files that contain only a few pulses, each pulse is selected at most once during sampling. The resulting dataset is made of 10638 ICP pulses among which 2816 have missing p ₁, 604 have missing p ₂, and 692 have missing p ₃. The challenging pulses are distributed among 58 patients.

2.2 MOCAIP++

This section introduces MOCAIP++, a generalization of the recently developed MOCAIP [15] which is an end-to-end framework that processes raw ICP signals to extract morphological waveform features through the recognition of the three peaks of the pulse. In its original form, MOCAIP relies on a Gaussian model to represent the prior knowledge about the position of each peak in the pulse. The Gaussian priors were replaced by a regression model in a recent extension [18].

MOCAIP++ generalizes its predecessors in two ways. First, it proposes a unifying view such that different peak recognition techniques can be integrated within the framework. Second, an additional processing step allows to exploit ICP features regardless the peak recognition method that is used. Similarly to MOCAIP, a pulse extraction technique (Section 2.2.1) first process the ICP signal to extract a reliable dominant pulse from which peak candidates are located at curves inflections (Section 2.2.2). Then, MOCAIP++ extracts different ICP features from the dominant pulse (such as curvature, first, and second derivative) (Section 2.2.3). The peak recognition module (Section 2.2.4) exploits the peak candidates and the features to recognize the peaks within the pulse. Finally, various statistics are estimated using the latency of these peaks and their ICP elevation (additional details can be found in the original papers [15, 18]). The core of the algorithm is illustrated in Figure 3 and its major components are described in the next subsections.

Figure 3

Diagram showing the different modules in MOCAIP++ framework.

2.2.1 ICP Segmentation and Dominant Pulse Extraction

The first component of the framework (ICP pulse segmentation) takes a raw, continuous ICP signal and splits it in a series of individual ICP pulses. An individual pulse is found using a pulse extraction technique [19] combined with the ECG QRS detection [20] that locates each ECG beat. Therefore, the latency of the three peaks within the ICP pulse is relative to the ECG QRS.

Because ICP recordings are subject to various noise and artifacts during the acquisition process, a robust dominant pulse S _i is extracted from a sequence of consecutive ICP pulses using hierarchical clustering [21]. It corresponds to the centroid of the largest cluster. In other words, the dominant pulse summarizes a short segment of consecutive ICP pulses.

2.2.2 Detecting Peak Candidates

Then, MOCAIP++ detects peak candidates (a ₁, a ₂, ..., a _N ) at curve inflections of the dominant ICP pulse S _i by segmenting the pulse into concave and convex regions using the second derivative of the signal. A peak is said to occur at the intersection of a convex and a concave region on a rising edge of ICP pulse, or at the intersection of a concave and a convex region on the descending edge of the pulse.

2.2.3 ICP Features

Previous MOCAIP-based studies [15, 18] exploited the dominant pulses directly as input to peak recognition techniques. In signal processing, it is common to derive features that emphasize different properties of the signal. For example, the first derivative measures the changing rate of the signal with respect to time. As illustrated in Figure 4, it is particularly interesting in our case because for a similar amplitude, a wide peak, and a narrow peak will lead to different derivative values. Therefore, features extracted from the ICP signal derivative provide additional morphological characteristics that should help to discriminate between ICP peaks. One advantage of using these features is that they are invariant to a shift of the signal elevation. Note that the framework is not restricted to these features, any other features could in principle be exploited. In our experiments, we will evaluate the impact of using the first L _x and second L _xx derivatives, as well as the curvature K extracted from the ICP signal within MOCAIP++ framework.

Figure 4

Signal L made of two Gaussian peaks with different standard deviations. Its first L _x and second L _xx derivatives are particularly usefull to discriminate peaks because their amplitude depends on the peak width but remains invariant to any global shift in elevation of the signal.

First Derivative

For more robustness, the ICP signal I(x) is first convolved with a Gaussian smoothing filter \mathcal{G}(x;\sigma ) where σ is the standard deviation of the Gaussian (σ = 3 in our experiments),

L(x,\sigma )=\mathcal{G}(x;\sigma )*I\left(x\right).

(3)

Then the derivative L _x is computed according to the smoothed version L of the ICP,

L_x=L(x,\sigma )–L(x+\text{1},\sigma ).

(4)

Second Derivative

Similarly, the computation of the second derivative L _xx relies on the first derivative L _x ,

L_{xx}=L_x(x,\sigma )–L_x(x+1,\sigma ).

(5)

Curvature

The curvature K is computed as a ratio between the first and the second derivative of the signal,

K=\frac{L_{xx}}{{(1+L_x^2)}^{3/2}}.

(6)

2.2.4 Peak Recognition

This module aims at recognizing the three peaks (p ₁, p ₂, p ₃) within an ICP pulse among the set of candidate peaks (a ₁, a ₂, ..., a _N ). Depending on the recognition technique, it can exploit the latency of the peak candidates, the raw ICP pulse, or different features extracted from the pulse. In the next, we describe three different peak recognition approaches. They are based on independent Gaussian models [15], Gaussian Mixture Models (GMM), and Spectral Regression (SR) analysis [18], respectively.

(a) Gaussian Model

The original MOCAIP algorithm exploits Gaussian priors to identify the most likely configuration of the three peaks among the set of candidates. Given P(X ₁), P(X ₂), and P(X ₃) to denote the Gaussian probability distribution of the prior position of the three peaks (p ₁, p ₂, p ₃), peak recognition amounts to searching for the maximum of the following objective function,

(7)

where P(X _i = a _k ) represents the probability of assigning a _k to the i-th peak.

In order to deal with missing peaks, an empty designation a ₀ is added to the pool of candidates. In addition, to avoid false designation, MOCAIP uses a threshold ρ such that P(X _i = a _k ) = 0, i ∈ {1, 2, 3}, k ∈ {1, 2, ..., N} if the probability of assigning a _k to p _i is less than ρ.

(b) Gaussian Mixture Models

In contrast with MOCAIP that uses a model of independent Gaussian distributions to represent the likely position of each peak, the method proposed in this paragraph exploits a multi-modal distribution to model the joint latency of the three peaks. Observed peak configurations are approximated by a Gaussian Mixture Model (GMM), where each component i represents a cluster of configurations μ _i of the three peaks. A GMM is defined as,

P(X=x|\Theta )={\displaystyle \sum _{i=1}^C{\alpha }_i}\mathcal{G}(x;{\mu }_i,{\Sigma }_i),

(8)

\mathcal{G}(x;{\mu }_i,{\Sigma }_i)=\frac{1}{\sqrt{2\pi {\Sigma }_i^2}}{\exp }^{\frac{-{(x-{\mu }_i)}^2}{2{\Sigma }_i^2}},

(9)

where α _i , μ _i , ∑ _i are the relative weight, the mean, and the variance of an individual component i, and C is the total number of components. For learning, Expectation-Maximization (EM) was used to estimate the model parameters θ = (α _1...C ; μ_1...C ; ∑_1...C ) that maximizes the likelihood of the observed peak configurations. EM was performed for a different number of components C ∈ {1, ..., 10}. The number which minimizes the Bayesian Information Criterion (BIC) [22] was selected.

The detection task amounts to find the best configuration of the three peaks among the set of peak candidates a = (a ₁, a ₂, ..., a _N ) detected in the current pulse. This can be done by finding the configuration that is the more likely on the GMM,

(10)

However, an additional difficulty is caused by missing peaks. One way to solve this problem is to use a hierarchical recognition approach where the possible configurations are first evaluated on the 3 peak model. If the largest response r ₁₂₃ = P(X = {p ₁, p ₂, p ₃}|Θ) fails to be above a given threshold τ₃, the marginals X ₁₂, X ₁₃, X ₂₃ using only two dimensions of the model are evaluated,

\{p_1,p_2\}=\underset{{\tilde{p}}_1,{\tilde{p}}_2}{\arg \max }P(X_{12}=\{{\tilde{p}}_1,{\tilde{p}}_2\}|\Theta )|{\tilde{p}}_1\in a\wedge {\tilde{p}}_2\in a,{\tilde{p}}_1\ne {\tilde{p}}_2,

(11)

\{p_1,p_3\}=\underset{{\tilde{p}}_1,{\tilde{p}}_3}{\arg \max }P(X_{13}=\{{\tilde{p}}_1,{\tilde{p}}_3\}|\Theta )|{\tilde{p}}_1\in a\wedge {\tilde{p}}_3\in a,{\tilde{p}}_1\ne {\tilde{p}}_3,

(12)

\{p_2,p_3\}=\underset{{\tilde{p}}_2,{\tilde{p}}_3}{\arg \max }P(X_{23}=\{{\tilde{p}}_2,{\tilde{p}}_3\}|\Theta )|{\tilde{p}}_2\in a\wedge {\tilde{p}}_3\in a,{\tilde{p}}_2\ne {\tilde{p}}_3,

(13)

and r ₁₂ = P(X ₁₂ = {p ₁, p ₂}|Θ), r ₁₃ = P(X ₁₃ = {p ₁, p ₃}|Θ), r ₂₃ = P(X ₂₃ = {p ₂, p ₃}|Θ). Again, if the maximum response to the GMM model of all the 2-peak configurations max(r ₁₂, r ₁₃, r ₂₃) is below a certain threshold τ ₂, 1-peak marginals X ₁, X ₂, X ₃ are evaluated, and the peak with the maximum response is marked.

p_1=\underset{{\tilde{p}}_1}{\arg \max }P(X_1=\{{\tilde{p}}_1\}|\Theta )|{\tilde{p}}_1\in a

(14)

p_2=\underset{{\tilde{p}}_2}{\arg \max }P(X_2=\{{\tilde{p}}_2\}|\Theta )|{\tilde{p}}_2\in a

(15)

p_3=\underset{{\tilde{p}}_3}{\arg \max }P(X_3=\{{\tilde{p}}_3\}|\Theta )|{\tilde{p}}_3\in a

(16)

(c) Spectral Regression

In a recent comparison of regression techniques [18], Spectral Regression (SR) [23] demonstrated excellent accuracy in peak recognition on standard ICP pulses. This motivates us to select it as the baseline regression method within MOCAIP++. The regression model y _i = f(x _i ) maps the position of the peaks as a function of the ICP dominant pulse. The model is automatically learned from training ICP pulses S = {S _i = 1...n} labeled with the latency of the peaks y _i = (p ₁, p ₂, p ₃) within the pulse. Each pulse S _i is resized to a vector x _i ∈ ℝ^sof length s = 500 ms, and normalized in amplitude between [1].

SR combines spectral graph analysis and standard linear regression to obtain a model that gives similar predictions {\widehat{y}}_i\in \widehat{Y} for data samples x _i ∈ X that are close (i.e. that are nearest neighbors in a graph representation), such that the following measure ϕ is minimized:

\varphi ={\displaystyle \sum _{i,j=1}^n{({\widehat{y}}_i-{\widehat{y}}_j)}^2}{W}_{i,j},

(17)

where W ∈ ℝ^{n × n}is the item-item similarity matrix that associates a positive value to W _i,j if the samples x _i , x _j belong to the same class. This is done by first using the eigenvectors of the matrix W,

We=\lambda De,

(18)

where D is a diagonal matrix whose entries are column sums of W, D _i,i = Σ _j W _j,i , and e ₀, e ₁, ..., e _d denote the d + 1 eigenvectors with respect to the d + 1 largest eigenvalues λ ₀≥λ ₁≥ ... ≥ λ _d.

Then SR finds d vectors {{\widehat{\beta }}_0,{\widehat{\beta }}_1,\dots ,\widehat{\beta }}that minimize the residual Sum of Square Error (SSE),

{\widehat{\beta }}_j=\arg {\min }_{\beta }{\displaystyle \sum _{i=1}^n{({\beta }^T{x}_i-y_i^j)}^2},

(19)

where y_i^j is the i-th element of e _j .

For recognition on a new pulse x _j , the regression model y _j = f(x _j ) predicts the most likely position of the three peaks y _j = (p ₁, p ₂, p ₃). A nearest neighbor search is then performed so that the nearest candidate (a ₁, a ₂, ..., a _N ) to each prediction is assigned to the peak label corresponding to the matched prediction. Additional features f _i ∈ ℝ(^s, where f _i ∈ {L _x , L _xx , K}, can be concatenated to the original input x _i ∈ ℝ(^sto create a new input vector [x _i f _i ] that combines both modalities.

Although Spectral Regression is a linear regression algorithm, it can easily be extended to become nonlinear by using a kernel projection (Radial Basis Function (RBF)) of the input vectors. We further refer to this technique as the Kernel Spectral Regression (KSR).

3 Results and Discussion

3.1 Accuracy of Peak Recognition Methods on Challenging Data

This section provides a comparative analysis of peak recognition techniques by evaluating their performance on the challenging dataset of ICP pulses. Models based on Gaussian (MOCAIP), Gaussian Mixtures (GMM), Spectral Regression (SR), and Kernel Spectral Regression (KSR) models are evaluated. A five-fold cross-validation is performed on the challenging dataset D', such that at each of the five iterations, four folds are used to train the model while the remaining one is retained for evaluation. The partitioning is randomly made with the constraint that the pulses of a given patient are grouped into the same fold. This ensures that data from the same patient are not present at the same time in the training and testing sets.

During evaluation, a predicted position {\widehat{y}}_i _i of one of the three peaks is considered to be correct if it is equal to the actual position y _i established manually. Given that peaks may be set as missing in the groundtruth, True positive (TP), false positive (FP), true negative (TN), or false negative (FN) are defined as follow,

\text{A prediticion }{\widehat{y}}_i\text{ of }{y}_i\text{ is a}\{\begin{array}{ll}TP,\hfill & \text{if}({\widehat{y}}_i=y_i\wedge {\widehat{y}}_i\ne \varnothing \wedge y_i\ne \varnothing )\hfill \\ FP,\hfill & \text{if}({\widehat{y}}_i\ne \varnothing \wedge y_i=\varnothing )\hfill \\ TN,\hfill & \text{if}({\widehat{y}}_i=\varnothing \wedge y_i=\varnothing )\hfill \\ FN,\hfill & \text{if}(({\widehat{y}}_i=\varnothing \vee ({\widehat{y}}_i\ne y_i))\wedge y_i\ne \varnothing )\hfill \end{array}

(20)

Based on these measures, the accuracy {\mathcal{A}}_p of one of the three peaks p ∈ {p ₁, p ₂, p ₃} is defined as,

{\mathcal{A}}_p\text{=}\left(TP+TN\right)/\left(TP+FP+TN+FN\right).

(21)

The accuracy {\mathcal{A}}_p of a peak p is obtained by averaging the accuracy over the five-folds. Similarly, the overall accuracy \mathcal{A} is obtained by averaging the accuracy of the three peaks, \mathcal{A}=({\mathcal{A}}_{p_1}+{\mathcal{A}}_{p_2}+{\mathcal{A}}_{p_3})/3. The learning of the recognition models is supervised in the sense that it relies on a set of manually labelled ICP pulses. As the number of training examples increases, the overall accuracy is generally expected to improve as well. We report this aspect by plotting the average prediction accuracy for each method against the number of training samples in Figure 5. To test one of the 5 folds, a model is trained by randomly extracting n pulses from the remaining 4 folds.

Figure 5

Effect of the number of training samples on the average recognition accuracy (Eq. 21) for different models (KSR, SR, GMM, MOCAIP [15]) using a five-fold crossvalidation on our challenging dataset D '. Results correspond to the average for the three peaks (p ₁, p ₂, and p ₃).

Results clearly indicates that KSR performs better by reaching a maximum accuracy of 88.78% ± 2.35, while the other techniques are less accurate; SR obtains 72.57 ± 2.6, GMM 70.47% ± 2.64, and MOCAIP 65.83% ± 2.96. It is interesting to notice that all the methods reaches their maximum accuracy before 500 training pulses. These results confirms that, besides KSR, current methods do not offer good recognition results on challenging pulses. Although KSR performs better than any other techniques, it requires all the training pulses to be kept as a part of the model to be able to compute the kernel projection. Nevertheless, KSR gives us an insight about what performance a peak recognition technique can achieve on our challenging dataset.

3.1.1 Computational Cost

One of the possible applications our framework is to be used in portable devices to monitor ICP continuously. Such an application requires real time performances of the peak recognition techniques. This section evaluates the performance of the different recognition techniques in terms of their complexity by comparing their computational time during learning and recognition.

Table 1 shows that, on a dataset of 2000 ICP pulses, MOCAIP (Gaussian), and SR are the fastest for training their model with only 60 and 90 ms, while GMM is much more slower with 33,940 ms. For recognition on a single pulse, SR ranks first with 0.19 ms, MOCAIP second (6.7 times slower), GMM is 10 times slower, and KSR is about 100 times slower than SR. Batch recognition performance is measured on a set of 2000 pulses. Under these conditions, KSR improves a lot due to the optimization of matrix operations but remains behind SR. Note that the reported durations only represents the running time of the learning, and recognition methods. Additional time is needed for MOCAIP++ to pre-process the ICP, detect peak candidates, and to compute additional features such as curvature and signal derivatives. Running time were measured using built-in MATLAB functions. These tests were performed on a DELL OPTIPLEX 760 computer equipped with INTEL DUAL-CORE E8600 cadenced with a 3.33 GHz processor and 3 GB of RAM.

Table 1 Running time for learning peak recognition models (Gaussian, SR, KSR, and GMM) from 2000 ICP pulses, and for recognition on 1 and 2000 pulses.

3.2 Feature-based Peak Recognition

This section evaluates the impact of the additional ICP features (Section 2.2.3) within MOCAIP++ on peak recognition performance. The same experimental protocol (five-fold cross-validation) of the previous section is used to evaluate the accuracy of SR, KSR, and GMM using three different features; curvature (Curv), first (L _x ) and second (L _xx ) derivatives on the challenging dataset D'.

Figure 6 shows that each feature significantly improves the overall accuracy of the SR method. While the original SR-based recognition method [18] attains an accuracy of 72.57% ± 2.6, the use of the second derivative and curvature improves it to 80.26 ± 2.29 and 80.4% ± 2.2, respectively. SR performs best when it is combined with the first derivative L _x of the ICP, reaching an accuracy of 85.81% ± 2.5. This constitutes a very significant result (+13%) in favor of our feature-based MOCAIP++ method.

Figure 6

The average recognition accuracy (Eq. 21) is reported versus the number of training samples. Results illustrate the effect of three differential features on the SR model [18]: curvature (Curv), first (L _x ) and second (L _xx ) derivatives.

When combined with derivative-based features, GMM, and KSR methods exhibit a similar ranking of improvement; first derivative offers the largest effect on accuracy, while curvature and second derivatives generally have less significant improvement. With the use of the first derivative (see Figure 7), GMM method improves from 70.47% ± 2.64 to 77.14% ± 1.85, while KSR only shows a marginal improvement from 88.78% ± 2.35 to 89.36% ± 2.51. We have also noticed in additional experiments that combining different features, such as L _x +L _xx , does not improve the performance obtained by using only the first derivative L _x of the ICP signal. These results demonstrate that the use of the first derivative within MOCAIP++ improves the recognition accuracy of the three peak recognition methods we have integrated. It can also be pointed out that the accuracy reached by SR + L _x is very close to KSR + L _x . Considering the previous remarks about the execution time and the storage of training samples for the kernel computation required for KSR, the use of SR combined with the first derivate seems to provide the right tradeoff between speed and accuracy for peak recognition on challenging ICP pulses.

Figure 7

Average recognition accuracy (Eq. 21) after a five-fold cross-validation for MOCAIP-based peak recognition methods improved with the use of the first derivative L _x of the ICP.

3.3 Impact of the Training Data Sampling Strategy

Although peak recognition models are trained in a supervised fashion such that they integrate morphological information from pulses with known peaks into models that may correctly identify peaks in new pulses, the underlying training pulses affect the estimation of the parameters and the performance of such models. Intuitively, the model should be trained on a representative range of pulses (easy, or challenging) to gain sufficient precision. This section evaluates the effect of incorporating pulses extracted from the challenging dataset into the training set of peak recognition methods.

In these experiments, peak recognition methods are estimated from two different annotated training sets ({\mathbb{T}}_1,{\mathbb{T}}_2). The first training set, named reference library {\mathbb{T}}_1, is made of 3000 randomly selected ICP pulses from the original dataset D. These pulses present a wide range of morphological variations but the majority of them are generally easily annotated. A subset of these data was used in previous works [15] to train MOCAIP. The second training set, named weighted sampling {\mathbb{T}}_2, is made of 1500 randomly selected pulses from the original dataset D, plus 1500 pulses randomly extracted from the challenging dataset D' created using a weighted sampling procedure (Section 2.1.2). Unlike previous section, where peak recognition methods were assessed against the challenging dataset D', the evaluation is now performed on the full dataset D. This allows us to see if the methods are not subject to overfitting; we verify if the methods that offer good results on challenging data also generalize well on regular pulses.

Figure 8 gives the average accuracy. It can be seen that the use of an equal number of pulses sampled from the full and challenging dataset considerably improves the performance of peak recognition methods over models exclusively learned on random pulses. The improvement is as follows: MOCAIP, from 72.31% to 87.27%, SR, from 75.96% to 82.41%, SR+Lx, from 83.67% to 92.74%, and KSR+L_x from 90.44% to 93.64%. The combination of our two contributions, the use of the first derivative and the weighted sampling for training, improves SR-based MOCAIP approach by about 17% (from 75.96% to 92.74%). This is a very significant improvement of performance that should help to extract more reliable statistics about ICP pulse morphology in real clinical conditions.

Figure 8

Effect of the source data used to train the peak recognition models on the average recognition accuracy (Eq. 21) evaluated on the large dataset D. The reference library {\mathbb{T}}_1 is made of randomly chosen ICP pulses. The second training set, weighted sampling{\mathbb{T}}_2, is made of an equal number of randomly selected pulses from the large dataset D and challenging pulses randomly extracted from the challenging dataset D'. Both datasets contain 3000 pulses

4 Conclusions

Recent works suggest that changes in the waveform morphology of ICP may provide insight to forecasting critical intracranial and cerebrovascular pathophysiological variations. However, automatic analysis of the waveform morphology of ICP acquired in clinical conditions is still beyond current ICP analysis frameworks. Their performance deteriorates significantly when the morphology of the pulse exhibits uncommon morphological changes.

This paper has described MOCAIP++, a generalization of the recently developed MOCAIP, that provides a robust framework for analyzing Intracranial Pressure signal (ICP) in terms of its waveform morphology. The proposed approach improves current methods by allowing the integration of several peak recognition methods. In addition, whereas previous MOCAIP-based studies [15, 18] exploited dominant pulses directly as input to peak recognition techniques, MOCAIP++ allows to derive additional features that capture more informative properties of the ICP signal and hence better discriminate the three peaks. The first derivative of the ICP signal has been shown to be the best among the features tested in our experiments (as shown in Figure 9). It improved all the peak recognition methods. This can be explained by its invariance to global shift in elevation from the baseline of the pulse. Performance in terms of peak recognition accuracy obtained by the proposed SR-based extension are close to the non-linear Kernel Spectral Regression (KSR). KSR can be considered really close to the best performing solution for this problem but it has the disadvantage to require to keep all the training samples, and is much slower than the SR.

Another contribution of this paper is to show that incorporating a challenging subset of ICP pulses into the training set of peak recognition methods has a positive effect on their overall accuracy.

Figure 9

Illustration of challenging ICP pulses where the original MOCAIP failed to recognize at least one of the peaks. The actual position of the peaks, correctly predicted by MOCAIP++ (SR + L _x ) are depicted by green circles, while the black diamonds correspond to the MOCAIP predictions

Experiments on a large dataset of ICP signals, as well as on a representative collection of sampled challenging ICP pulses, demonstrate that both contributions are complementary and significantly improve the recognition performance of ICP peaks in real conditions. These findings provide insight in order to potentially improve other ICP peak recognition frameworks, and will help us to collect more reliable statistics about ICP morphology to further investigate its predictive power on patient condition.