### 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*
_{1}, *p*
_{2}, and *p*
_{3}) 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*
_{1}, 3699 have missing *p*
_{2}, and 4626 have missing *p*
_{3}. 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*
_{1}, *p*
_{2}, and *p*
_{3}). 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{{\mathcal{E}}_{p1}+{\mathcal{E}}_{p2}+{\mathcal{E}}_{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*
_{1}, *p*
_{2}, and *p*
_{3}, respectively. The distribution of the weights is illustrated in Figure 2.

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*
_{1}, 604 have missing *p*
_{2}, and 692 have missing *p*
_{3}. 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.

#### 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*
_{1}, *a*
_{2}, ..., *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.

##### 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 )\u2013L(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 )\u2013{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*
_{1}, *p*
_{2}, *p*
_{3}) within an ICP pulse among the set of candidate peaks (*a*
_{1}, *a*
_{2}, ..., *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*
_{1}), *P*(*X*
_{2}), and *P*(*X*
_{3}) to denote the Gaussian probability distribution of the prior position of the three peaks (*p*
_{1}, *p*
_{2}, *p*
_{3}), peak recognition amounts to searching for the maximum of the following objective function,

\begin{array}{c}J(x,y,z)=P({X}_{1}={a}_{x})+P({X}_{2}={a}_{y})+P({X}_{3}={a}_{z})|{a}_{x}\in a\wedge {a}_{y}\in a\wedge {a}_{z}\in a,\\ |x<y<z,\end{array}

(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*
_{0} 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}}}{\mathrm{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*
_{1}, *a*
_{2}, ..., *a*
_{
N
} ) detected in the current pulse. This can be done by finding the configuration that is the more likely on the GMM,

\begin{array}{c}\{{p}_{1},{p}_{2},{p}_{3}\}=\underset{{\tilde{p}}_{1},{\tilde{p}}_{2},{\tilde{p}}_{3}}{\mathrm{arg}\mathrm{max}}P(X=\{{\tilde{p}}_{1},{\tilde{p}}_{2},{\tilde{p}}_{3}\}|\Theta )|{\tilde{p}}_{1}\in a\wedge {\tilde{p}}_{2}\in a\wedge {\tilde{p}}_{3}\in a\\ |{\tilde{p}}_{1}<{\tilde{p}}_{2}<{\tilde{p}}_{3}.\end{array}

(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*
_{123} = *P*(*X* = {*p*
_{1}, *p*
_{2}, *p*
_{3}}|Θ) fails to be above a given threshold τ_{3}, the marginals *X*
_{12}, *X*
_{13}, *X*
_{23} using only two dimensions of the model are evaluated,

\{{p}_{1},{p}_{2}\}=\underset{{\tilde{p}}_{1},{\tilde{p}}_{2}}{\mathrm{arg}\mathrm{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}}{\mathrm{arg}\mathrm{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}}{\mathrm{arg}\mathrm{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*
_{12} = *P*(*X*
_{12} = {*p*
_{1}, *p*
_{2}}|Θ), *r*
_{13} = *P*(*X*
_{13} = {*p*
_{1}, *p*
_{3}}|Θ), *r*
_{23} = *P*(*X*
_{23} = {*p*
_{2}, *p*
_{3}}|Θ). Again, if the maximum response to the GMM model of all the 2-peak configurations max(*r*
_{12}, *r*
_{13}, *r*
_{23}) is below a certain threshold *τ*
_{2}, 1-peak marginals *X*
_{1}, *X*
_{2}, *X*
_{3} are evaluated, and the peak with the maximum response is marked.

{p}_{1}=\underset{{\tilde{p}}_{1}}{\mathrm{arg}\mathrm{max}}P({X}_{1}=\{{\tilde{p}}_{1}\}|\Theta )|{\tilde{p}}_{1}\in a

(14)

{p}_{2}=\underset{{\tilde{p}}_{2}}{\mathrm{arg}\mathrm{max}}P({X}_{2}=\{{\tilde{p}}_{2}\}|\Theta )|{\tilde{p}}_{2}\in a

(15)

{p}_{3}=\underset{{\tilde{p}}_{3}}{\mathrm{arg}\mathrm{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*
_{1}, *p*
_{2}, *p*
_{3}) within the pulse. Each pulse *S*
_{
i
} is resized to a vector *x*
_{
i
} ∈ ℝ^{s}of 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*,

where *D* is a diagonal matrix whose entries are column sums of *W*, *D*
_{
i,i
} = Σ_{
j
}
*W*
_{
j,i
} , and *e*
_{0}, *e*
_{1}, ..., *e*
_{
d
} denote the *d* + 1 eigenvectors with respect to the *d* + 1 largest eigenvalues *λ*
_{0}≥*λ*
_{1}≥ ... ≥ *λ*
_{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}=\mathrm{arg}{\mathrm{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*
_{1}, *p*
_{2}, *p*
_{3}). A nearest neighbor search is then performed so that the nearest candidate (*a*
_{1}, *a*
_{2}, ..., *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
} ∈ ℝ(^{s}to 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).