2.1.

Experimental Setup and Signal Extraction

Experimental procedures were approved by the University of Pittsburgh and Carnegie Mellon University’s Institutional Animal Care and Use Committees and complied with the guidelines in the National Institute of Health’s Guide for the Care and Use of Laboratory Animals.

Details of the experimental protocol and data collection have been reported previously.^28,30 Briefly, eight anesthetized NHPs (NHPs 1 to 4 Macaca mulatta aged $7.9\pm 1.5\text{years}$, weighing $9.4\pm 0.7\mathrm{kg}$, and NHPs 5 to 8 Macaca fascicularis aged $4.2\pm 0.9\text{years}$, weighing $5.1\pm 2.1\mathrm{kg}$) were used for the experiments. Each NHP was sedated throughout experimentation. Prior to their transport into the experimentation room, the NHPs were sedated with $20\text{-}\mathrm{mg}/\mathrm{kg}$ ketamine. In some animals, an additional $0.04\text{-}\mathrm{mg}/\mathrm{kg}$ of atropine and $1\text{-}\mathrm{mg}/\mathrm{kg}$ diazepam were administered. During the experiment, the monkeys were anesthetized with a combination of 0.6% to 1.5% isoflurane (ISO) and 10- to $25.6\text{-}\mu \mathrm{g}/\mathrm{kg}/\mathrm{h}$ fentanyl. Additionally, $0.1\text{-}\mathrm{mg}/\mathrm{kg}/\mathrm{h}$ of vecuronium bromide paralytic was administered intravenously. NHPs were ventilated at 0.18 to 0.4 Hz. The ventilation frequency was kept constant for the duration of the experiment for each NHP.

A frequency-domain NIRS system (OxiplexTS, ISS Inc., United States), as shown in Fig. 1, operating at 690 and 830 nm, was used to measure cerebral hemoglobin concentration changes. The differential pathlength factor (DPF) ratio was held constant between all animals at ${\mathrm{DPF}}_{690}/{\mathrm{DPF}}_{830}=1.1$. Although DPF differences are likely between animals, the cardiac waveform shape is not expected to be influenced significantly. To overcome magnitude differences between animals, signals were normalized in steps described later. Optical probes were placed directly on the skull of the animals, above the right prefrontal cortex. For NHPs 1 to 3, the source detector distance was 2.2 cm, and for NHPs 4 to 8, it was 1.5 cm. The data for this manuscript were based on retrospective analysis. Differences in source detector distances were due to variations in overall experimental designs between animals. Light intensity changes were recorded at a sampling frequency of 50 Hz and converted to changes in oxygenated ($\mathrm{\Delta }\mathrm{HbO}$), deoxygenated ($\mathrm{\Delta }\mathrm{Hb}$), and total ($\mathrm{\Delta }\mathrm{HbT}=\mathrm{\Delta }\mathrm{HbO}+\mathrm{\Delta }\mathrm{Hb}$) hemoglobin concentrations using the modified Beer–Lambert’s law.³¹ In addition, a DCS probe was placed on the left prefrontal cortex in NHPs 1 to 4 and the right temporal cortex for NHPs 5 to 8, with source detector distances of 2.2 and 1.9 cm, respectively. DCS data from NHP 4 were excluded from further processing due to instrument instabilities. DCS sampling frequency was 100 Hz. Further details of the device are found in Ref. 27 Relative CBF was estimated using a scaled diffusion coefficient measured from the autocorrelation of light intensity changes.¹⁰ The scaled diffusion coefficient has been shown to correlate with CBF.³² CBF was then used for further processing.

Fig. 1

Experimental setup with NIRS and DCS placement on the skull above the right and left prefrontal cortex, respectively. The NIRS setup shows a defined source–detector distance. ICP was measured using a parenchymal pressure sensor. Saline was administered using a ventricular catheter. ICP was regulated by adjusting the height of the saline reservoir.

Arterial blood pressure (ABP) was recorded with an MPR1 Datalogger (Raumedic Helmbrechts, Germany) using an arterial line placed in the carotid artery. ICP was altered using a catheter (Lumbar catheter, Medtronic, Minneapolis, Minnesota, United States) placed in the lateral ventricle of the brain, with the other end connected to a saline reservoir, as shown in Fig. 1. A change in the height of the reservoir resulted in a pressure change in the head (simulating hydrocephalus). ICP was measured using a parenchymal pressure sensor recorded by the MPR1 Datalogger (Raumedic Helmbrechts, Germany). Both ICP and ABP were recorded at 100 Hz.

Each NHP experiment, lasting $\sim 22.6\mathrm{h}$ (std 2.3 h), was split between 7 and 10 separate trials. Each trial was $\sim 90\min$ long and corresponded to a particular ICP level, elevated using saline infusion. Induced ICP ranged between 5 and 60 mmHg, with natural ICP fluctuation during experimentation deviating slightly beyond these limits. The distribution of both induced and naturally fluctuating ICP levels was observed to be predominantly between 5 and 30 mmHg. Before ICP elevation, a recording at baseline ICP was also conducted.

2.2.

Signal Processing

Following the initial data collection, the signal processing procedures closely followed those outlined in Ref. 27 Analog markers in the form of voltage spikes were sent to the auxiliary ports of the DCS, NIRS, and MPR1 Datalogger devices for data alignment. In NHPs 1 to 3, an amplifier circuit was set to register a cardiac pulse when the electrocardiogram (EKG) signal exceeded an empirically determined threshold that defined the R peak of the QRS complex. When an R peak was detected, a synchronization pulse was sent to the MPR1 Datalogger. This signal was sampled at 100 Hz. For NHPs 4 to 8, EKG was measured through a separate device at 1000 Hz. A maximal overlap discrete wavelet transform was run across the EKG signal to enhance the QRS complex features of the signal (the QRS complex represents ventricular polarization). MATLAB’s (MATLAB R2020b, The MathWorks Inc., Natick, Massachusetts, United States) findpeaks function was used to index the time point of the peak of each enhanced QRS complex. The respective indices represent the R peak of the QRS complex. During data collection, occasional laser instabilities were observed. Time points in the NIRS- and DCS-based signals, for which instabilities were visually observed, were removed. For NHP 4, the DCS-based CBF signal showed consistent instability across all trials. As a result, NHP 4 was not included in the ICP estimation analysis of CBF. All signals not originally sampled at 50 Hz (ABP, ICP, CBF, and QRS complex peak indices) were downsampled to 50 Hz to match the sampling frequency of the NIRS signal.

To reduce noise in the NIRS and DCS signals, 120 consecutive cardiac pulses were averaged, and an average cardiac waveform (ACPW) was extracted. Because the heart rate of the animals varied, 120 pulses corresponded to between 39 and 78 s across all animals. The ICP values during the 120 pulses were also averaged. The 120-pulse averaging window was moved 20 pulses at a time (resulting in an 83.3% window-to-window overlap) along the entire signal of each trial. All ACPWs were normalized in time and amplitude. With the help of spline interpolation, the length of each ACPW was normalized to 66 data points, corresponding to 1.32 s. ACPW length was measured between two consecutive diastoles. The height of each ACPW (representing the amplitude of $\mathrm{\Delta }\mathrm{HbO}$, $\mathrm{\Delta }\mathrm{HbT}$, and CBF) was normalized to between 0 and 1. Normalization of the $x$- and $y$-axis removed ACPW length as a feature and thus removed heart rate as a feature of the waveform. We removed heart rate as a feature because heart rate changes can be independent of ICP.^33–35

To remove outliers, a $z$-score rejection method was applied to all ACPWs. For each trial, a $z$-score was calculated across all averaged pulses in the trial. This averaged $z$-score was then compared against each ACPW. Any $\mathrm{\Delta }\mathrm{HbO}$, $\mathrm{\Delta }\mathrm{HbT}$, or CBF ACPW with a $z$-score greater than 3 was rejected.

A Kalman filter was then used to further improve signal quality. The adaptive filter was applied to the ACPWs of each signal for each trial. Each ACPW was compared with the ideal trial pulse produced by the Kalman filter. ACPWs were corrected based on their error from the ideal pulse. In the correction method, the Kalman filter’s parameters defined the weight given to the Kalman-produced ideal pulse and the weight given to the calculated ACPW-to-ideal-pulse error. These parameters were set empirically. The output of the Kalman filter procedure was a set of high signal-to-noise ratio (SNR) ACPWs with feature morphology reflecting changes in ICP. ACPWs reflecting ICP values above 30 mmHg were removed due to their scarcity. Figure 2 shows examples of 50 Hz $\mathrm{\Delta }\mathrm{HbO}$, $\mathrm{\Delta }\mathrm{HbT}$, ICP, and CBF signals. A total of 19,000 $\mathrm{\Delta }\mathrm{HbO}$ and 19,258 $\mathrm{\Delta }\mathrm{HbT}$ ACPWs sampled from 8 NHPs, along with 17,322 CBF ACPWs sampled from 7 of the same 8 NHPs were used for feature extraction. All ACPWs represented an ICP range of 0 to 30 mmHg.

Fig. 2

Filtered time trace examples of ICP, CBF, $\mathrm{\Delta }\mathrm{HbO}$, and $\mathrm{\Delta }\mathrm{HbT}$ signals. Dashed lines indicate QRS complex peaks.

2.3.

Feature Extraction

Defining and extracting physiologically relevant features from the processed ACPWs provides the basis for estimating ICP. In our previous work, features of the CBF waveform included the percussion, tidal, and dicrotic peaks of the ACPWs.²⁷ These features were not generally discernable from the NIRS-based hemoglobin concentration signals. As such, the feature set and their extraction method were modified, as shown in Fig. 3. Figure 3 also shows an example of a $\mathrm{\Delta }\mathrm{HbO}$ ACPW.

Fig. 3

(a) Peak-based features, including peak height (${\mathrm{P}1}_{\mathrm{pk}}$) (normalized to 1), peak width (${\mathrm{P}1}_{\mathrm{w}}$), peak prominence (${\mathrm{P}1}_{\mathrm{p}}$), and peak location (${\mathrm{P}1}_{\mathrm{pos}}$), for a typical $\mathrm{\Delta }\mathrm{HbO}$ ACPW. (b) $x$- and $y$-coordinates of the center of mass of the same ACPW. MAP was used as an additional feature, not shown here.

MATLAB’s findpeaks function was used to obtain the peak height (${\mathrm{P}1}_{\mathrm{pk}}$: normalized to 1 and used as a measure for model noise), position (${\mathrm{P}1}_{\mathrm{pos}}$), prominence (${\mathrm{P}1}_{\mathrm{p}}$), and width (${\mathrm{P}1}_{\mathrm{w}}$) of the ACPWs.

The $x$- and $y$-coordinates of the center of mass (centroid) of the waveform, as shown in Fig. 3(b), defined as $\mathrm{COM}x$ and $\mathrm{COM}y$, were also extracted from individual ACPWs and used as features. This was done by turning each ACPW into a polygon and calculating its centroid using the polyshape and centroid functions in MATLAB. We hypothesized that $\mathrm{COM}x$, which describes waveform skewness, is related to blood pressure and ICP. A similar reasoning motivated $\mathrm{COM}y$. The area under the curve (AUC) of the waveform was also incorporated as a feature, as was MAP. Similar to ICP, MAP was calculated over each 120-pulse window. The feature engineering resulted in eight human-interpretable and observable features for each processed $\mathrm{\Delta }\mathrm{HbO}$, $\mathrm{\Delta }\mathrm{HbT}$, and CBF ACPW. If a feature was undetected, its value was set to 0, but it was still used.

2.4.

Machine Learning

Each of the three ACPW datasets was randomly sampled into five cross validation (CV) sets of 80% training and 20% testing. For each CV set, all animals were included. Random sampling ensured that learning became NHP and trial independent, while CV allowed us to alleviate overfitting. Python’s scikit-learn toolbox’s random forest (RF) regression algorithm was used as the ICP estimator.³⁶ The RF algorithm, explained in detail by Ref. 37, learns a set number of decision trees on a randomly sampled subset of the features using a randomly sampled subset of the training data with replacement (bootstrapping). A total of 100 trees, or estimators, were learned using this bootstrapping method. Every tree received 50% of the features and 80% of the dataset for training. All other hyperparameters were kept as default, defined in the scikit-learn documentation.³⁶ These hyperparameters were chosen empirically to maximize performance while mitigating overfitting. The hyperparameter decision-making process used a random search to gauge approximate hyperparameter ranges, after which a per-parameter and joint-parameter optimization procedure followed. Random search hyperparameter tuning is a common method used for optimizing machine learning models.³⁸ One important parameter was the number of estimators, or trees. Each of the 100 trees splits its subset of data until all leaves are pure—to max depth. Each tree receives four randomly sampled features for learning. Gini impurity was used as the measure of node split quality.³⁹ For each CV split (fold), the RF algorithm was trained on the training split and tested on the testing split. During testing, estimated ICP values (${\mathrm{ICP}}_{\mathrm{est}}$) were compared with invasively measured ICP values (${\mathrm{ICP}}_{\mathrm{inv}}$). ${\mathrm{ICP}}_{\mathrm{inv}}$ values were used as ground truth labels. The performance of the model was quantified using the coefficient of determination ($r^2$), mean squared error (MSE), and 95% confidence interval (CI).

3.1.

Waveform Features and RF Regressors

We have previously found that the CBF ACPW broadens and skews to the right when ICP increases.²⁷ We qualitatively observed similar trends in the $\mathrm{\Delta }\mathrm{HbO}$ and $\mathrm{\Delta }\mathrm{HbT}$ ACPWs.

Using an RF regressor with modified hyperparameters, we evaluated the ICP estimator’s performance on NIRS-derived ACPWs ($\mathrm{\Delta }\mathrm{HbO}$ and $\mathrm{\Delta }\mathrm{HbT}$) against its performance on DCS-based ACPWs (CBF). Although the CBF-based waveforms showed more morphological features compared with the $\mathrm{\Delta }\mathrm{HbO}$- and $\mathrm{\Delta }\mathrm{HbT}$-based waveforms, the features chosen were the same between all three modalities used. Within the 0 to 30 mmHg range of ICP values, available training and testing data were skewed toward lower ICP values, as shown in Fig. 4(a).

Fig. 4

(a) Histogram of ICP distribution for both train and test sets across all three modalities. More data were available at lower ICP values, especially between 5 and 10 mmHg. (b)–(d) Correlation plots for $\mathrm{\Delta }\mathrm{HbO}$, $\mathrm{\Delta }\mathrm{HbT}$, and CBF. Strong $r^2$ for all methods suggests that the model performs well on ICP estimation. Estimation performance drops for higher ICP values across all three modalities due to a lower availability of high ICP data for training.

To compare performance differences between the three different types of signals, the coefficient of determination ($r^2$) and MSE were used as evaluation metrics, along with a 95% CI. A fivefold CV was performed individually on the $\mathrm{\Delta }\mathrm{HbO}$, $\mathrm{\Delta }\mathrm{HbT}$, and CBF datasets. We found a mean fold $r^2=0.937$ (averaged over five folds), $r^2=0.946$, and $r^2=0.963$ for the ACPW types $\mathrm{\Delta }\mathrm{HbO}$, $\mathrm{\Delta }\mathrm{HbT}$, and CBF, respectively, with a fivefold $r^2$ standard deviation of $r_{\mathrm{std}}^2=0.003$ ($\mathrm{\Delta }\mathrm{HbO}$), $r_{\mathrm{std}}^2=0.004$ ($\mathrm{\Delta }\mathrm{HbT}$), and $r_{\mathrm{std}}^2=0.002$ (CBF). We further found a mean fold $\mathrm{MSE}=2.703$, 2.301, and $1.688{\mathrm{mmHg}}^2$ with a fivefold standard deviation of ${\mathrm{MSE}}_{\text{std}}=0.133$, 0.163, and $0.131{\mathrm{mmHg}}^2$, respectively. Our results indicate that CBF waveforms perform better than NIRS-derived parameters in terms of ICP extraction, but ICP can still be quantified within an MSE of $<3{\mathrm{mmHg}}^2$ when using hemoglobin-based waveforms.

The correlation plots between ${\mathrm{ICP}}_{\mathrm{est}}$ and ${\mathrm{ICP}}_{\mathrm{inv}}$ for all three signals are shown in Figs. 4(b)–4(d). All three ICP extractions show a good correlation between estimated and invasively measured ICP. Outliers were more common at higher ICP values for which less data were available for training, as shown in the histograms of Fig. 4(a).

A set of three Bland–Altman plots, indicating the level of agreement between ${\mathrm{ICP}}_{\mathrm{est}}$ and ${\mathrm{ICP}}_{\mathrm{inv}}$, are shown in Fig. 5 and further confirm the clear fit of our regression model to the invasively measured ICP. The plots show a 95% CI of agreement between ${\mathrm{ICP}}_{\mathrm{est}}$ and ${\mathrm{ICP}}_{\mathrm{inv}}$ of $[-3.0643.160]\mathrm{mmHg}$ with a mean of 0.048 mmHg for $\mathrm{\Delta }\mathrm{HbO}$, a 95% CI of $[-2.8412.866]\mathrm{mmHg}$ with a mean of 0.013 mmHg for $\mathrm{\Delta }\mathrm{HbT}$, and a 95% CI of $[-2.4502.397]\mathrm{mmHg}$ with a mean of $-0.026\mathrm{mmHg}$ for CBF.

Fig. 5

Bland–Altman plots for $\mathrm{\Delta }\mathrm{HbO}$, $\mathrm{\Delta }\mathrm{HbT}$, and CBF. Histograms show the distribution of datapoints across the Bland–Altman plot.

Table 1 summarizes the results across all three modalities. When considering the performance of $\mathrm{\Delta }\mathrm{HbO}$ against $\mathrm{\Delta }\mathrm{HbT}$ and CBF, we find that $\mathrm{\Delta }\mathrm{HbT}$ performs slightly better than $\mathrm{\Delta }\mathrm{HbO}$. $\mathrm{\Delta }\mathrm{HbO}$ and $\mathrm{\Delta }\mathrm{HbT}$ have a 0.009 difference in fivefold mean $r^2$ and a $0.402{\mathrm{mmHg}}^2$ difference in fivefold mean MSE. CBF performs better than the two NIRS-based methods with an increase of 0.025 fivefold mean $r^2$ and a $0.402{\mathrm{mmHg}}^2$ decrease in fivefold mean MSE against $\mathrm{\Delta }\mathrm{HbO}$. The 95% CI metric across the modalities echoes this difference in model fit performance. Using bootstrap analysis in which 50% of test samples were randomly drawn 10,000 times with replacement and a respective $r^2$ was calculated for each set of samples, we found that CBF performs statistically better than either of the NIRS methods (significance level $(\text{alpha})\le 0.05$). Despite this, NIRS performs comparatively well against DCS, along with being cheaper and easier to use. Distributions of $r^2$ scores with 5%, mean ($\mu$), and 95% thresholds are shown in Fig. 6.

Table 1

Error metrics across the three modalities.

Fig. 6

Histograms of bootstrapped $r^2$ scores for all three modalities. Bootstrapping used 50% of test samples drawn 10,000 times with replacement. The test data from the best performing CV split, per modality, are used. This results in sample mean $r^2$ slightly above the mean CV $r^2$ reported. $r^2$ is calculated for each set of samples across all three modalities ($\mathrm{\Delta }\mathrm{HbO}$, $\mathrm{\Delta }\mathrm{HbT}$, and CBF). A significance level $\le 0.05$ is used. CBF statistically outperforms both NIRS methods. NIRS methods are still comparable to CBF in overall performance.

For all three modalities, the feature relevance contributing to ICP extraction was evaluated. The relative magnitude of feature importance is shown in Fig. 7. We found that MAP and $\mathrm{COM}x$ were the most relevant features used in the training of each estimator tree. AUC, $\mathrm{COM}x$, and MAP accounted for $\sim 15\%$, $\sim 20\%$, and $\sim 28\%$ of node splits per tree, respectively. Meanwhile, $\mathrm{COM}y$ represented $\sim 11\%$ of node splits, whereas waveform peak position represented $\sim 9\%$. The remaining three features represented $<7\%$ of node splits per tree. The standard deviation of features used between trees is shown by the error bars of each feature bar. AUC importance for $\mathrm{\Delta }\mathrm{HbO}$, $\mathrm{\Delta }\mathrm{HbT}$, and CBF was approximately the same at around 15% to 18%, whereas peak-related features (peak location, width, and prominence) were more important for CBF compared with the NIRS methods. $\mathrm{COM}x$ and $\mathrm{COM}y$ were less important for $\mathrm{\Delta }\mathrm{HbT}$ and CBF compared with $\mathrm{\Delta }\mathrm{HbO}$. MAP was the most important feature across all modalities, with the level of importance being approximately the same. For $\mathrm{\Delta }\mathrm{HbT}$, MAP represented $\sim 29\%$ of splits per tree compared with $\sim 28\%$ for $\mathrm{\Delta }\mathrm{HbO}$ and $\sim 27\%$ for CBF. The importance of the height of the peak, which was normalized, is used as a proxy to measure the importance of an unrelated peak. The apparent relevance of this unrelated peak is highest for $\mathrm{\Delta }\mathrm{HbO}$ and lowest for CBF, suggesting that the model learns slightly more uncorrelated patterns on the $\mathrm{\Delta }\mathrm{HbO}$ dataset compared with the datasets of the other modalities.

Fig. 7

Comparison of feature importance across $\mathrm{\Delta }\mathrm{HbO}$, $\mathrm{\Delta }\mathrm{HbT}$, and CBF. MAP is strongest across all three modalities. AUC and $\mathrm{COM}x$ are relevant for estimation, and peak-based features (${\mathrm{P}1}_{\mathrm{pos}}$, ${\mathrm{P}1}_{\mathrm{w}}$, and ${\mathrm{P}1}_{\mathrm{p}}$) are more important in CBF-based ICP estimation compared with the hemoglobin-based methods. It is also important to note that ${\mathrm{P}1}_{\mathrm{pk}}$ is low for all methods, especially for CBF. This suggests that the model does not learn irrelevant features.

A more evenly balanced dataset was also tested for $\mathrm{\Delta }\mathrm{HbO}$. Balancing was performed by randomly removing half of all ICP values between 5 and 10 mmHg. This balanced dataset had $\sim 23\%$ less data than the skewed dataset. The balanced dataset led to a marginally worse performance compared with the skewed dataset and was not studied further for either $\mathrm{\Delta }\mathrm{HbT}$ or CBF. For $\mathrm{\Delta }\mathrm{HbO}$, $r^2$ dropped by 0.6% from 0.937 to 0.932, MSE increased by 22.1% from 2.703 to $3.299{\mathrm{mmHg}}^2$, and mean 95% CI increased by [12.3% 11.5%] from $[-3.0643.160]$ to $[-3.4413.523]\mathrm{mmHg}$ with a 14.6% change in mean from 0.048 to 0.041 mmHg. Because MAP is a dominant feature for all three signals used, we evaluated whether MAP can be removed to make DCS and NIRS standalone methods for ICP extraction. When MAP was removed as a training feature, the importance of the four peak features was increased, and the importance of AUC and $\mathrm{COM}x$ was maintained. With MAP removed, performance fell to $r^2=0.838$, $\mathrm{MSE}=6.971{\mathrm{mmHg}}^2$, and 95% CI of $[-4.9604.984]\mathrm{mmHg}$ with a mean of 0.012 mmHg for $\mathrm{\Delta }\mathrm{HbO}$; $r^2=0.865$, $\mathrm{MSE}=5.791{\mathrm{mmHg}}^2$, and 95% CI of $[-4.5564.700]\mathrm{mmHg}$ with a mean of 0.072 mmHg for $\mathrm{\Delta }\mathrm{HbT}$; and $r^2=0.916$, $\mathrm{MSE}=3.778{\mathrm{mmHg}}^2$, and 95% CI of $[-3.7853.756]\mathrm{mmHg}$ with a mean of $-0.014\mathrm{mmHg}$ for CBF. These results indicate that the performance is improved with MAP as an added feature, but ICP can still be quantified to $<\sim 5\mathrm{mmHg}$ for hemoglobin-based estimators and $<\sim 4\mathrm{mmHg}$ for CBF-based estimators.

4.1.

Dataset

In building an NHP-based dataset of ICP-dependent ACPW features, it was observed that $\mathrm{\Delta }\mathrm{HbO}$ and $\mathrm{\Delta }\mathrm{HbT}$ cardiac pulse waveforms followed a canonical cardiac pulse arch shape. Using 120-pulse averaging and Kalman filtering improved SNR to the extent that meaningful feature extraction could be performed. The number of pulses averaged over was set to 120 to include enough pulses for appropriate waveform and ICP averaging. Averaging too many pulses would likely lead to lower ICP-to-waveform-feature correlation, whereas not averaging enough pulses would not remove enough noise in the signal waveforms.

Despite selecting only ICP values up to 30 mmHg, the distribution of data was not equal for all ICP values, as shown in the histogram in Fig. 4(a). We found that balancing the dataset did not significantly improve ICP estimation. There appeared to be a more even distribution of outliers in ICP estimation between low and high ICP (0 to 30 mmHg) in the balanced dataset, compared with the skewed dataset (data not shown). Meanwhile, feature importance was not significantly different between the skewed and balanced data. It can be argued that the unbalanced dataset was more likely to cause overfitting to lower ICP values and to artificially inflate performance. Although possible, we do not have enough data points at higher ($>20\mathrm{mmHg}$) ICP values to quantify this.

4.2.

Preprocessing and Filtering

To appropriately leverage the raw signals extracted from each NHP, signal preprocessing and filtering were required. The goal was to provide the feature extraction method with low-noise signals while not overfiltering or losing valuable information. Clean signals lead to relevant feature-value-to-ICP mappings for the machine learning model to learn.

Tuning of the signal preprocessing and filtering parameters was done empirically across the cardiac waveform averaging, waveform normalization, $z$-score rejection, and Kalman filtering methods. 120 consecutive pulses were found to be a suitable averaging amount to reduce signal to noise, whereas a 20-pulse averaging window shift struck a balance between obtaining ample amounts of data for training and testing while reducing data replication. Waveform normalization, $z$-score rejection, and Kalman filtering all worked to remove data outliers while maintaining a feature distribution that supported model generalization.

4.3.

Feature Engineering

An investigation into feature importance and performance, as shown in Fig. 7, suggests that MAP on its own, though relevant, is not the only feature driving decision tree learning. This gives weight to our other chosen features, such as AUC and $\mathrm{COM}x$. We also found poor correlation between the raw values of MAP, AUC, and $\mathrm{COM}x$ with ICP. This suggests that a feature on its own is not a good indicator of ICP, giving strength to our approach of combining several features in a machine learning method to accurately estimate ICP. Interpretable feature engineering is an important part of clinically relevant and understandable machine learning tools. Explainable models ensure that conclusions made by the model are based on appropriate decisions.⁴⁰

Extracting specific features from data may improve model generalizability for small sample sizes. The FDA’s Good Machine Learning Practice for Medical Device Development: Guiding Principles⁴¹ highlights the importance of human interpretability of models and their outputs in clinical settings. As such, the features that we engineered are more interpretable when compared with machine-engineered features sometimes used in complex machine learning and neural network-based methods.²⁹ We also found that our engineered features expressed themselves relatively predictably across modalities.

4.4.

RF Regression

In our previous work on DCS-based ICP extraction, we showed a cross-validated similarity between estimated ICP and invasively measured ICP of $r^2=0.91$, an MSE of 3.3 mmHg, and a 95% CI of $[-3.73.7]\mathrm{mmHg}$. Here we compared NIRS-derived waveforms to DCS waveforms for ICP extraction and achieved a cross-validated similarity between DCS-based ICP estimation and invasively measured ICP of $r^2=0.96$, MSE of 1.7 mmHg, and 95% CI of $[-2.4502.397]\mathrm{mmHg}$. This is an improvement to our previous work. The performance increase was likely due to using more data and a revised machine learning method. Compared with Ruesch et al.’s 14,121 ACPWs across 5 NHPs,²⁷ our CBF dataset contained 17,322 ACPWs sampled from 7 NHPs. The ICP distributions were similar. Regarding the ICP estimators, Ruesch et al.²⁷ used a bagged ensemble of trees with subsampling of data. Unlike our RF model, which was chosen for its versatility, robustness, resistance to noise, and ability to manage unbalanced datasets,^42–45 the bagged tree regressor of Ruesch et al.²⁷ did not randomly subsample features for each tree in the ensemble. In an RF regression, the random sampling of features ensures that no one tree receives the full set of features for learning. This combats overfitting and improves performance.

RF also have a relatively low risk of overfitting with an increase in estimators due to their use of multiple weak and unpruned learners.⁴⁶ Each tree in the RF ensemble does overfit to the data and the features that it receives, but the averaged result of the ensembled trees produces a regression prediction that is low in variance and bias. Ultimately, this means that the more trees that are used in the ensemble, the less likely the model is to overfit. More trees, however, increase computational cost. Compared with Ruesch et al.’s 1000 trees,²⁷ we achieved comparable results with only 100 trees. This means that our approach is more computationally efficient with comparable results.

It is important to note that the signal processing, feature extraction, and training methods used on the three signal types explored in this paper ($\mathrm{\Delta }\mathrm{HbO}$, $\mathrm{\Delta }\mathrm{HbT}$, and CBF) were the same and that the features were specifically engineered with NIRS-based signals in mind. DCS-based (CBF) ICP estimation with these specific features was relatively more dependent on peak-based features than the $\mathrm{\Delta }\mathrm{HbO}$ and $\mathrm{\Delta }\mathrm{HbT}$ ICP estimation, but not to the same extent as observed in Ruesch et al.²⁷ Meanwhile, features such as AUC, $\mathrm{COM}x$, and $\mathrm{COM}y$ were important across the modalities. MAP was the defining feature for all signal types, with varying degrees of relative importance compared with other features. This suggests that the chosen features are potentially useful across signal and sensor modalities. CBF waveforms, however, do carry more waveform-related information compared with waveforms extracted using blood volume changes. Thus, waveform-based features are more relevant for CBF compared with $\mathrm{\Delta }\mathrm{HbO}$ and $\mathrm{\Delta }\mathrm{HbT}$. The features chosen for this work may therefore not be optimized for CBF, and further improvement in ICP extraction may be possible.

4.5.

Feature Differences Across NHPs

Occasionally observed differences in the morphology of the same waveform feature across subjects are likely due in part to influences on the NHPs’ natural autoregulation and their abilities to manage varying ICPs. Because ISO gas, commonly used during anesthesia, is known to suppress CA,^47,48 we used $<1\%$ to avoid anesthesia induced autoregulation impairment.³⁰ However, changes in ICP led to CPP levels outside the autoregulatory range.^8,9,30 The Windkessel effect—the elasticity of arteries that moderates blood pressure changes during a cardiac cycle—may also be a factor in observed morphological differences. Studies show that conduit arteries harden with age.⁴⁹ The stiffening of cerebral arteries results in an increased pulse pressure per stroke volume and a higher systolic pressure. This effect has primarily been studied in humans, but the effect may be similar in the NHPs used and could therefore contribute to inconsistent morphology of volume-dependent cardiac pulse waveforms between NHPs of different ages. Our youngest NHP was 3 years old, and our oldest was 10 years old. Macaques reach sexual maturity at 3 to 5 years (REFS) and can live to 25 years or more in captivity.^50,51

Small differences in probe and catheter placement between NHPs were also likely contributors to observed feature dissimilarity and variance between ${\mathrm{ICP}}_{\mathrm{inv}}$ and ${\mathrm{ICP}}_{\mathrm{est}}$. Although care was taken in the placement of the NIRS device on the skull of the NHPs, small variations did exist between subjects. In addition, ICP sensor and ventricular catheter placement varied slightly due in part to distinctive subject physiologies.

4.6.

Future Work

One avenue of future study along model selection could investigate probabilistic RF (PRF). Reis et al.⁵² showed that, compared with normal RF, PRFs improved classification by 10% on noisy features and 30% on noisy labels. This model may also be useful for regression. Natural continuations of this study would include large-scale human validation. Due to the relatively small sample size, we did not attempt a leave-one-out analysis. All animals were represented in both the training and the testing sets. A meaningful assessment of the generalizability of our ICP estimator to unseen data would be facilitated by a larger dataset.

An additional future consideration could evaluate the importance of the location of MAP estimation. In this manuscript, we measured MAP from the carotid artery. An alternative approach may the use of a blood pressure cuff, which hypothetically provides a sufficient MAP reading to be used with our approach. This, however, has not yet been evaluated. If successful, an MAP reading using a blood pressure cuff would ease clinical translation of our noninvasive ICP estimation method.