Temporal Statistical Relationship between Regional Cerebral Oxygen Saturation (rSO2) and Brain Tissue Oxygen Tension (PbtO2) in Moderate-to-Severe Traumatic Brain Injury: A Canadian High Resolution-TBI (CAHR-TBI) Cohort Study

Brain tissue oxygen tension (PbtO2) has emerged as a cerebral monitoring modality following traumatic brain injury (TBI). Near-infrared spectroscopy (NIRS)-based regional cerebral oxygen saturation (rSO2) can non-invasively examine cerebral oxygen content and has the potential for high spatial resolution. Past studies examining the relationship between PbtO2 and NIRS-based parameters have had conflicting results with varying degrees of correlation. Understanding this relationship will help guide multimodal monitoring practices and impact patient care. The aim of this study is to examine the relationship between PbtO2 and rSO2 in a cohort of TBI patients by leveraging contemporary statistical methods. A multi-institutional retrospective cohort study of prospectively collected data was performed. Moderate-to-severe adult TBI patients were included with concurrent rSO2 and PbtO2 monitoring during their stay in the intensive care unit (ICU). The high-resolution data were analyzed utilizing time series techniques to examine signal stationarity as well as the cross-correlation relationship between the change in PbtO2 and the change in rSO2 signals. Finally, modeling of the change in PbtO2 by the change in rSO2 was attempted utilizing linear methods that account for the autocorrelative nature of the data signals. A total of 20 subjects were included in the study. Cross-correlative analysis found that changes in PbtO2 were most significantly correlated with changes in rSO2 one minute earlier. Through mixed-effects and time series modeling of parameters, changes in rSO2 were found to often have a statistically significant linear relationship with changes in PbtO2 that occurred a minute later. However, changes in rSO2 were inadequate to predict changes in PbtO2. In this study, changes in PbtO2 were found to correlate most with changes in rSO2 approximately one minute earlier. While changes in rSO2 were found to contain information about future changes in PbtO2, they were not found to adequately model them. This strengthens the body of literature indicating that NIRS-based rSO2 is not an adequate substitute for PbtO2 in the management of TBI.


Introduction
Traumatic brain injury (TBI) is the preeminent form of neurotrauma globally and is a leading cause of death and disability worldwide [1,2].Contemporary management of TBI largely focuses on guideline-based management aimed at global arterial blood pressure (ABP), intracranial pressure (ICP), and cerebral perfusion pressure (CPP) treatment thresholds to try to minimize the ongoing brain injury that occurs in the acute period following the initial event [3,4].While there have been multiple iterations of this paradigm, there has been little progress in further improving outcomes following TBI [5].Given the large global burden of this disease, different approaches to the management of TBI are being explored.Attention has shifted towards a precision-medicine-based approach to the management of TBI that incorporates multi-modal cerebral monitoring [6][7][8][9].
The most widely adopted of these a"xili'ry modalities is brain tissue oxygen tension (PbtO 2 ), with recent guidelines incorporating its use into treatment algorithms [4].Additionally, at least three large phase 3 randomized control trials are underway to evaluate the ability of PbtO 2 monitoring to improve outcomes following TBI [10][11][12].PbtO 2 measures diffusible oxygen content in the extracellular space through a Clarke-type electrode.As a result, it is speculated that it has a slower response to physiologic variations than other monitoring modalities.This also necessitates its invasive placement into viable brain tissue and results in its ability to only sample a small volume of brain tissue [13].
Near-infrared spectroscopy (NIRS) as a cerebral monitoring modality, in the setting of TBI, has grown for the past three decades [14].The non-invasive nature, ease of application, and potential for high spatial resolution of NIRS monitoring are significant advantages of this modality.While there is a robust body of evidence supporting its quantitative relationship with cerebral blood flow, its precise relationship with PbtO 2 is much less clear, with some studies concluding a strong linear correlation between the modalities in retrospective observational studies [15][16][17][18][19].Other retrospective observational studies failed to identify any statistical relationship [20][21][22].
One possible etiology of this discrepancy is a failure to leverage time series analysis techniques that allow for the utilization of high-resolution data streams while accounting for the autocorrelative and hierarchical nature of this type of physiologic data.Without accounting for hierarchical and autocorrelation structures in the data, the assumptions of the regression methods utilized in these studies, mainly the independence of samples, is not entirely valid.It can be hypothesized that this may have resulted in an erroneously strong correlation between these parameters.Presented here is a study utilizing the Canadian High-Resolution Traumatic Brain Injury (CAHR-TBI) Research Collaborative database, which is a multi-institutional database of moderate-to-severe TBI patients.A highly unique feature of this dataset is its concurrently measured high-resolution PbtO 2 and rSO 2 data, the likes of which have not previously been reported in the literature.The primary objective of this study is to utilize contemporary time series analysis techniques to better characterize the relationship between PbtO 2 and NIRS-based regional cerebral oxygen saturation (rSO 2 ).The secondary objective of this study is to determine if PbtO 2 can be adequately modeled by rSO 2 in the setting of moderate-to-severe TBI.

Study Design
Using data from the Canadian High Resolution TBI (CAHR-TBI) Research Collaborative, a retrospective multicenter cohort study utilizing a prospectively collected database of critically ill TBI patients was performed, similar to recently published work from this group.Retrospective analysis was determined to be appropriate based on availability of data.Additionally, given the exploratory and observational nature of the study, an entirely prospective study would not be particularly beneficial at this time.Patients included in this study were admitted to one of the following university-affiliated hospitals: Vancouver General Hospital (University of British Columbia), Foothills Medical Centre (University of Calgary), and Health Sciences Centre Winnipeg (University of Manitoba).Institutions collaborating in this research have committed to prospective collection of high-resolution physiologic data in TBI patients.Research ethics approval at the University of Manitoba has been obtained for this database (H2017:181 and H2017:188).Additionally, ethics approval was obtained for retrospective access to the database, as well as for anonymous data transfer between centers (H2020:118, H20-03759 and REB20-0482).

Patient Population
The CAHR-TBI database is composed of moderate-to-severe TBI patients (defined as having a Glasgow Coma Scale ranging from 3 to 12) treated at an adult intensive care unit (ICU).All patients in this database had invasive ICP and ABP monitoring and were cared for using management strategies based on Brain Trauma Foundation (BTF) guidelines [3].PbtO 2 was managed based on local practice norms which varied from aggressive management to purely observational.PbtO 2 monitors were placed into viable brain tissues based on CT scans or under direct inspection in the operating room.NIRS-based rSO 2 was not actively used to guide management at any of the institutions.While granular patient-specific data are not available, similar vasopressor, sedative, and hyperosmolar/hypertonic agents were utilized at all participating institutions.Patient data were entered into the database from 2011 to 2022.
Included in this study were all patients in the CAHR-TBI database that had concurrent invasive PbtO 2 monitoring and NIRS-based rSO 2 monitoring.Those without PbtO 2 or NIRS monitoring were excluded.Age, biological sex, admission Glasgow Coma Score (GCS), admission pupil exam, and follow-up Glasgow Outcome Score (GOS) were extracted.Sample size calculations were not possible and therefore not performed due to the exploratory nature of this study.

High-Resolution Physiologic Data Collection
High-resolution physiologic data-streams included ICP, ABP, and PbtO 2 , as well as left and right rSO 2 .ABP was measured utilizing radial arterial lines.ICP was monitored using intra-parenchymal strain gauge probes (Codman ICP MicroSensor; Codman & Shurtlef Inc., Raynham, MA, USA) placed in the frontal lobe or using external ventricular drains (Medtronic, Minneapolis, MN, USA).PbtO 2 was measured using intra-parenchymal brain tissue oxygenation probes (Licox Brain Tissue Oxygen Monitoring System; Integra LifeSciences Corp., Plainsboro, NJ, USA) placed in viable frontal lobe tissue.rSO 2 was measured using NIRS regional cerebral oximetry of both the left and right frontal lobes (Covidien INVOS 5100C or 7100) when possible.
Data streams were recorded in digital high-frequency time series (≥100 Hz for ABP and ICP, 1 Hz for PbtO 2 and rSO 2 ) using analogue-to-digital signal converters (Data Translations, DT9804 or DT9826) when required.This digitized data were linked and stored in time series using Intensive Care Monitoring (ICM+) software (Version 8.5, Cambridge Enterprise Ltd., Cambridge, UK).For the purposes of this study ICP and ABP were included for the sake of cohort characterization and were not utilized in subsequent data analysis.

Physiologic Data Cleaning and Processing
High-resolution physiologic data were artifact-cleared manually by a qualified clinician utilizing ICM+ software.Artifacts were determined through the examination of waveforms for ICP and ABP.Additionally, sudden drops in PbtO 2 and rSO 2 to zero were deemed artifactual.All data were cleaned without knowledge of patient demographics or study objectives.
All high-resolution data streams were processed into both 10-seconds-by-10-seconds and minute-by-minute data utilizing ICM+ software.Data were then exported as commaseparated value (CSV) files for data analysis.The minute-by-minute data were utilized throughout the analysis as it is generally considered the standard for cerebral multimodal monitoring signal analytics as it provides a good balance between data size/computational time and temporal resolution [23,24].The exception was the cross-correlative analysis, where the 10-seconds-by-10-seconds data were also used to confirm the findings from the minute-by-minute data and in the impulse response function plots to evaluating the data at a temporal resolution and frequency in keeping with vasomotion [25].Due to the side of PbtO 2 probe placement not being available for all patients, rSO 2 values on the right were utilized for analysis, unless not available, in which case rSO 2 values from the left were utilized.

Overview
The data analysis was performed using R statistical software (Version 4.2.2,R Foundation for Statistical Computing, Vienna, Austria) with the following packages: astsa, blandr, forecast, lmtest, nlme, tidyverse, tseries, and zoo.The Intel oneAPI Math Kernel Library (Intel Corp., Santa Clara, CA, USA) was utilized for the Basic Linear Algebra Subprograms (BLAS) and the Linear Algebra Package (LAPACK) to improve computational performance.Data streams were further filtered to exclude PbtO 2 values less than 0 mmHg and greater than 60 mmHg, as these values were felt likely to be erroneous based on clinical expertise.rSO 2 values less than 25% were excluded, as this is the lower limit of output for the INVOS devices used.This, along with the previously mentioned manual artifact clearing, resulted in discontinuities in the data streams.As subsequent time series analysis required continuous data streams, discontinuities were filled through basic linear interpolation through the approx() function.For all models, alpha was set to 0.05 without correction for multiple comparisons.Additionally, no sample size or power calculations were performed.This was due to the exploratory nature of this study, which leveraged available data sets.
In this study, the relationship between PbtO 2 and rSO 2 was characterized through time series analysis and inferential statistical modeling.Linear regression assumes independent sampling.In the setting of frequent resampling from individual subjects, as is the case in all high-frequency cerebral physiologic monitoring, this assumption is invalid in two regards.Resampling from the same subject in a cohort leads to a hierarchical structure, as samples from the same subject are likely to be more similar than those taken between subjects.This is because intrasubject sampling has random unaccounted-for effects held constant that are not constant in intersubject sampling.Beyond this, when samples are taken with a high frequency, there is a tendency for samples to be correlated with previously taken samples from the same subject.This is known as autocorrelation.The simplest means by which to reinstitute the validity of these assumptions of linear regression is to average the data over large epochs of time (i.e., such as daily or over the entire recording period) for each subject and then regress over these averages.In the dynamic setting of the critically ill TBI patient, this can result in a significant loss of information.
Fortunately, there are statistical methods to account for these deviations from the assumption of independent sampling.In this study, two such methodologies were utilized.First, hierarchical linear modeling, also known as linear mixed-effects modeling, can help account for the random effects experienced by each subject.Additionally, time seriesbased autoregressive integrative moving average (ARIMA) modeling can account for the autocorrelative structure of the modeled data stream, in this case PbtO 2 .The details of this analysis are described further in the subsequent sections.However, an in-depth review of the theoretical background of these methodologies is beyond the scope of this paper.We refer the interested reader to previous works on this subject [26][27][28][29], and literature applying such methodologies to cerebral physiologic data [30].

Determination of Stationarity of Physiologic Data
Prior to performing time series modeling and analysis, it was necessary to determine if the response data streams of interest, PbtO 2 , were stationary.This is necessary, as there is an assumption of signal stationarity in the utilized time series modeling techniques.If signals are not stationary, they must be transformed to be made stationary prior to any modeling.Testing of signal stationarity was accomplished through examination of the autocorrelative function (ACF) plots for each patient's PbtO 2 data.For all ACF plots, significance levels were set to a correlation level of +/−(2/N 1/2 ), where N is the number of samples.Generally, for each patient, there was no rapid drop off to zero in lag significance, indicative of nonstationarity of the series.This was confirmed with the Augmented Dicky-Fuller (ADF) and Kwiatkowski-Phillips-Schmidt-Shin (KPSS) testing for stationarity.
To make the series stationary, the first difference was taken for each data series.Following first differencing, the stationarity of the data streams was confirmed through inspection of the ACF plots which now showed rapid decay in the significance of subsequent lags.Additionally, ADF and KPSS testing both confirmed stationarity of the first differenced signals.Physiologically, the data streams can now be thought of as the change in PbtO 2 (∆PbtO 2 ) and the change in rSO 2 (∆rSO 2 ).While this is not the same as PbtO 2 and rSO 2 , it was felt that the data obtained from examining how ∆PbtO 2 relates to ∆rSO 2 would provide insight into the relationship between these physiologic parameters.

Cross-Correlative Relationship between ∆PbtO 2 and ∆rSO 2
Given that PbtO 2 samples the extracellular fluid of brain parenchyma and NIRSbased rSO 2 reflects changes in the brain microvasculature, it is conceivable that changes in these parameters would be asynchronous [31,32].To identify and examine this potential asynchrony, a cross-correlative analysis was carried out between ∆PbtO 2 and ∆rSO 2 , utilizing the minute-by-minute data over the entire cohort.The cross-correlation function (CCF) of the minute-by-minute data shows the largest correlative magnitude at a lag of 1 (∆rSO 2 Lag1).This can be interpreted as ∆rSO 2 being most strongly correlated with ∆PbtO 2 a minute later.The CCF plot of the 10-seconds-by-10-seconds data reinforces this conclusion, as the most significant lag is seen at lag 6, corresponding to a one-minute delay between ∆rSO 2 and ∆PbtO 2 .Since ∆rSO 2 Lag1 was found to contain the most information about ∆PbtO 2 , it was used in subsequent linear modeling.

Vector Autoregressive Modeling and Impulse-Response Function Plots
To further provide insights into the relationship between ∆rSO 2 and ∆PbtO 2 , impulseresponse function (IRF) plots were created based on a multivariate vector autoregressive (VAR) model.These plots examine the modeled response of ∆rSO 2 and ∆PbtO 2 to a sudden impulse of ∆ABP.The high-frequency 10-seconds-by-10-seconds data streams of interest (∆ABP, ∆ICP, ∆PbtO 2 , and ∆rSO 2 ) were utilized as the cerebral vasoactive response was being examined in this analysis and acts on a frequency scale of approximately 0.1 Hz [23,33].Given the non-stationarity of the original data streams, the differenced data were used.These parameters were used going forward for VAR modeling and generation of the IRF plots.
To determine the appropriate autoregressive order of the VAR model, the Akaike Information Criterion (AIC) was determined for vector autoregressive models of order 1 to 15.There was diminishing marginal improvement in AIC as model order increased past 6.As such, following the principle of parsimony, a VAR model of order 6 was constructed utilizing the VAR() function in R. Finally, using the irf() function in R, this VAR model was utilized to model and plot the response in ∆PbtO 2 as well as ∆rSO 2 of an orthogonal impulse in ∆ABP over the subsequent 10 lags.

Hierarchical Linear Modeling of ∆PbtO 2 from ∆rSO 2 Lag1
Given the hierarchical nature of this dataset, a linear mixed-effects model with random slope and intercept, utilizing the lme() function, was performed with ∆PbtO 2 as the dependent variable and ∆rSO 2 Lag1 as the independent variable.Next, to evaluate whether there was any autocorrelation of the model's residuals, ACF and partial autocorrelative function (PACF) plots were made of these residuals.As with the ACF plots, for all PACF plots, significance levels were set to a correlation level of +/−(2/N 1/2 ), where N is the number of samples.There was an autocorrelative structure to the residuals, indicating that there was unaccounted for autocorrelation in the response variable, ∆PbtO 2 .
2.5.6.Modeling of ∆PbtO 2 from ∆rSO 2 Lag1 Accounting for Autocorrelative Structure Population-level linear mixed-effects models that incorporated the autocorrelative structure of ∆PbtO 2 were constructed, as has been done in previous cerebral physiology studies [34].However, given the size of this dataset, the computational times were unacceptably long, with models failing to converge even after weeks of computation.As a result, it was determined that the next best option was to construct a linear model that accounted for the autocorrelative structure of ∆PbtO 2 for each individual subject and make inferences about the relationship between PbtO 2 and rSO 2 based on the general findings of these models.
For each subject, an initial simple linear regression was performed with ∆PbtO 2 as the dependent variable and ∆rSO 2 Lag1 as the independent variable.In order to account for the autocorrelative nature of ∆PbtO 2 , the ARIMA structure of the residuals needed to be determined.This was done through the auto.arima()function for the residuals of the linear model in each subject.Next, the autoregressive and moving average orders were utilized in the arima() and sarima() functions, with ∆PbtO 2 as the response variable and ∆rSO 2 Lag1 as the external regressor, to produce a linear model that accounted for the autocorrelative structure of ∆PbtO 2 .For each subject, the significance the ∆rSO 2 Lag1 coefficient was determined.Additionally, for each model, ACF and PACF plots of the residuals were examined to determine if there was any remaining autocorrelative structure.

Evaluating Model Correlation and Agreement
To examine the correlation between the predicted values of ∆PbtO 2 and actual values of ∆PbtO 2 for each subject, a Pearson correlation coefficient was obtained.Next, a Bland-Altman plot was produced to evaluate agreement between the predicted values of ∆PbtO 2 and actual values of ∆PbtO 2 .

Cohort Demographics
A total of 20 subjects were included in the study, with a total of 114,136 min of time with concurrent PbtO 2 and rSO 2 measurements without interpolation.The full demographic data of the cohort can be found in Table 1.

Determination of the Stationarity of the Physiologic Data
The ACF plots, as well as the ADF and KPSS test results for each patient's PbtO 2 data, can be found in Supplementary File S1.ACF plots did not show a rapid drop-off of significant lags, indicating non-stationarity.An example of an ACF plot can be seen in Figure 1A.KPSS testing for each patient's PbtO 2 data also uniformly indicated nonstationarity.In most patients, ADF testing indicated no presence of a unit root.As such, a transformation to the data was required before the assumption of stationarity could be fulfilled and time series models of the data constructed.This pattern of ACF plots and KPSS and ADF testing results is consistent with a difference stationary series, and so the first-order difference of the data was taken to transform the data.
The ACF plots, as well as the ADF and KPSS test results for each patient's ∆PbtO 2 data, can be found in Supplementary File S1.Once the first difference was taken, each patient's ACF plots showed a rapid decline in the significance of lags, indicating stationarity.An example of an ACF plot can be seen in Figure 1B.Consistent with this, KPSS testing indicated stationarity of the ∆PbtO 2 data for each patient.This indicated that the firstorder-differenced data fulfilled the assumptions of stationarity.As a result, subsequent analysis was carried out using both ∆PbtO 2 and ∆rSO 2 data.
testing indicated stationarity of the ΔPbtO2 data for each patient.This indicated that the first-order-differenced data fulfilled the assumptions of stationarity.As a result, subsequent analysis was carried out using both ΔPbtO2 and ΔrSO2 data.

Cross-Correlative Relationship between ∆PbtO 2 and ∆rSO 2
Cross-correlation analysis between ∆PbtO 2 and ∆rSO 2 , over the entire cohort, indicated the ∆PbtO 2 was most strongly correlated with ∆rSO 2 one minute earlier.In other words, ∆PbtO 2 and ∆rSO 2 Lag1 shared the strongest cross-correlation.In the 10-secondsby-10-seconds data, a similar pattern was seen with ∆PbtO 2 correlating with ∆rSO 2 six lags earlier, equivalent to one minute earlier.The CCF plots of the minute-by-minute data and 10-seconds-by-10-seconds data can be seen in Figures 2A and 2B, respectively.As a result of this finding, ∆PbtO 2 and ∆rSO 2 Lag1 were utilized for modeling.

Cross-Correlative Relationship between ΔPbtO2 and ΔrSO2
Cross-correlation analysis between ΔPbtO2 and ΔrSO2, over the entire cohort, indicated the ΔPbtO2 was most strongly correlated with ΔrSO2 one minute earlier.In other words, ΔPbtO2 and ΔrSO2Lag1 shared the strongest cross-correlation.In the 10-secondsby-10-seconds data, a similar pattern was seen with ΔPbtO2 correlating with ΔrSO2 six lags earlier, equivalent to one minute earlier.The CCF plots of the minute-by-minute data and 10-seconds-by-10-seconds data can be seen in Figure 2A and Figure 2B, respectively.As a result of this finding, ΔPbtO2 and ΔrSO2Lag1 were utilized for modeling.In Panel (A), the plot of the cross-correlative function (CCF) of ΔPbtO2 vs. ΔrSO2 for the minute-by-minute data of the cohort can be seen, with the most significant lag occurring at lag 1 (i.e., 1 min).In Panel (B), the plot of the cross-correlative function of ΔPbtO2 vs. ΔrSO2 for the 10seconds-by-10-seconds data for the data of the cohort can be seen, with the most significant lag occurring at lag 6 (i.e., at 1 min).The dashed blue line represents the significance levels, which were set to a correlation level of +/−(2/N 1/2 ), where N is the number of samples.CCF = Cross-correlative In Panel (A), the plot of the cross-correlative function (CCF) of ∆PbtO 2 vs. ∆rSO 2 for the minute-by-minute data of the cohort can be seen, with the most significant lag occurring at lag 1 (i.e., 1 min).In Panel (B), the plot of the cross-correlative function of ∆PbtO 2 vs. ∆rSO 2 for the 10-seconds-by-10-seconds data for the data of the cohort can be seen, with the most significant lag occurring at lag 6 (i.e., at 1 min).The dashed blue line represents the significance levels, which were set to a correlation level of +/−(2/N 1/2 ), where N is the number of samples.CCF = Cross-correlative function, ∆PbtO 2 = Change in brain tissue oxygen tension, ∆rSO 2 = Change in region cerebral oxygen saturation.

Vector Autoregressive Modeling and Impulse-Response Function Plots
The AIC for the multivariate VAR models from order 1 to 15 can be seen in Figure 3.It is clear that there is a clear drop off in marginal improvements in AIC for VAR models with an order greater than 6.As such, a VAR model incorporating ∆ABP, ∆ICP, ∆PbtO 2 , and ∆rSO 2 with an order of 6 was utilized to construct the IRF plots.The IRF plots of the modeled response in ∆rSO 2 and ∆PbtO 2 to a sudden impulse in ∆ABP can be seen in Figures 4A and 4B, respectively.The modeled response of ∆rSO 2 indicates an almost instantaneous response to an impulse in ∆ABP followed by a sharp decrease and eventual return to equilibrium by approximately lag 6.This is in contrast to ∆PbtO 2 , where an impulse in ∆ABP results in a delayed and prolonged response that only peaks at approximately lag 6 or 7.
function, ΔPbtO2 = Change in brain tissue oxygen tension, ΔrSO2 = Change in region cerebral ox saturation.

Vector Autoregressive Modeling and Impulse-Response Function Plots
The AIC for the multivariate VAR models from order 1 to 15 can be seen in Figu It is clear that there is a clear drop off in marginal improvements in AIC for VAR mo with an order greater than 6.As such, a VAR model incorporating ΔABP, ΔICP, ΔPb and ΔrSO2 with an order of 6 was utilized to construct the IRF plots.The IRF plots o modeled response in ΔrSO2 and ΔPbtO2 to a sudden impulse in ΔABP can be seen in ure 4A and Figure 4B, respectively.The modeled response of ΔrSO2 indicates an alm instantaneous response to an impulse in ΔABP followed by a sharp decrease and even return to equilibrium by approximately lag 6.This is in contrast to ΔPbtO2, where an pulse in ΔABP results in a delayed and prolonged response that only peaks at appr mately lag 6 or 7.

Hierarchical Linear Modeling of ∆PbtO 2 from ∆rSO 2 Lag1
The linear mixed-effects model (population level model) with random slope and intercept did find ∆rSO 2 Lag1 to be a significant positive linear regressor of ∆PbtO 2 (0.35, S.E.0.10, p = 0.0002).However, the ACF and PACF plots can be seen in Figures 5A and 5B, respectively.There is a clear demonstration of autocorrelation of the residuals, indicating an unaccounted-for autocorrelative structure of ∆PbtO 2 .This indicates that the autocorrelative structure of ∆PbtO 2 needs to be accounted for, prior to valid inferences being made about the relationship between ∆PbtO 2 and ∆rSO 2 Lag1.The 95% confidence intervals are indicated by the red dashed line.Note that in both plots there is an initial rise followed by subsequent decline; however, this is prolonged in the response of ΔPbtO2.

Hierarchical Linear Modeling of ΔPbtO2 from ΔrSO2Lag1
The linear mixed-effects model (population level model) with random slope and intercept did find ΔrSO2Lag1 to be a significant positive linear regressor of ΔPbtO2 (0.35, S.E.0.10, p = 0.0002).However, the ACF and PACF plots can be seen in Figure 5A and Figure 5B, respectively.There is a clear demonstration of autocorrelation of the residuals, indicating an unaccounted-for autocorrelative structure of ΔPbtO2.This indicates that the autocorrelative structure of ΔPbtO2 needs to be accounted for, prior to valid inferences being made about the relationship between ΔPbtO2 and ΔrSO2Lag1.shows the modeled resulting response in change in regional cerebral oxygen saturation (∆PbtO 2 ) to an orthogonal impulse in change in arterial blood pressure (∆ABP); black solid line.The 95% confidence intervals are indicated by the red dashed line.Note that in both plots there is an initial rise followed by subsequent decline; however, this is prolonged in the response of ∆PbtO 2 .

Modeling of ∆PbtO 2 from ∆rSO 2 Lag1 Accounting for Autocorrelative Structure
Attempts to create a linear mixed-effects model that adequately incorporated the autocorrelative structure of ∆PbtO 2 at the population level were unsuccessful.This was due to computational complexity with failure to converge.Creating independent inferential linear models for each patient was more successful.The details of these models can be found in Table 2. Of note, ∆rSO 2 was found to be a significant regressor in 16 of the 20 patients.However, in three of these patients, the coefficient was negative, which is not consistent with the expected relationship between rSO 2 and PbtO 2 .Examination of the ACF and PACF plots, found in Supplementary File S2, showed minimal autocorrelative structure remaining in the residuals of these models.This confirms model adequacy.ΔPbtO2~ΔrSO2Lag1) with random slope and intercept are seen.In Panel (B), the partial autocorrelative function (PACF) plot is seen.In both plots, a significant correlation is seen.This is in keeping with a model that has not accounted for the autocorrelative structure of its response variable.The dashed blue line represents the significance levels which were set to a correlation level of +/−(2/N 1/2 ), where N is the number of samples.CCF = Cross-correlative function, ΔPbtO2 = Change in brain tissue oxygen tension, ΔrSO2Lag1 = The one-minute lagged change in regional cerebral oxygen saturation.

Modeling of ΔPbtO2 from ΔrSO2Lag1 Accounting for Autocorrelative Structure
Attempts to create a linear mixed-effects model that adequately incorporated the autocorrelative structure of ΔPbtO2 at the population level were unsuccessful.This was due to computational complexity with failure to converge.Creating independent inferential linear models for each patient was more successful.The details of these models can be found in Table 2. Of note, ΔrSO2 was found to be a significant regressor in 16 of the 20 patients.However, in three of these patients, the coefficient was negative, which is not ) with random slope and intercept are seen.In Panel (B), the partial autocorrelative function (PACF) plot is seen.In both plots, a significant correlation is seen.This is in keeping with a model that has not accounted for the autocorrelative structure of its response variable.The dashed blue line represents the significance levels which were set to a correlation level of +/−(2/N 1/2 ), where N is the number of samples.CCF = Cross-correlative function, ∆PbtO 2 = Change in brain tissue oxygen tension, ∆rSO 2 Lag1 = The one-minute lagged change in regional cerebral oxygen saturation.

Evaluating Model Correlation and Agreement
The results of the correlation analysis between actual and predicted values of ∆PbtO 2 are also summarized in Table 2. Notably, correlation coefficients were generally poor for each subject-based model ranging from 0.04 to 0.57.Scatter plots of actual and predicted values of ∆PbtO 2 for each subject can be found in Supplementary File S3.
An example of a Bland-Altman plot for a single model can be seen in Figure 6, with the full series available in Supplementary File S4.Uniformly, agreement was poor throughout all individual subject models of ∆PbtO 2 from ∆rSO 2 Lag1.

Discussion
A statistically rigorous exploration of the relationship between the change in PbtO2 and change in NIRS-based rSO2 was performed in this multi-institutional cohort of 20 moderate-to-severe TBI patients.There are three key insights brought about by this study.First, changes in PbtO2 are correlated with changes in rSO2 that occur one minute earlier.Second, changes in rSO2, in a linear way, contain information about changes in PbtO2 that occur one minute later.Finally, while changes in rSO2 have this delayed linear relationship

Discussion
A statistically rigorous exploration of the relationship between the change in PbtO 2 and change in NIRS-based rSO 2 was performed in this multi-institutional cohort of 20 moderate-to-severe TBI patients.There are three key insights brought about by this study.First, changes in PbtO 2 are correlated with changes in rSO 2 that occur one minute earlier.Second, changes in rSO 2 , in a linear way, contain information about changes in PbtO 2 that occur one minute later.Finally, while changes in rSO 2 have this delayed linear relationship with changes in PbtO 2 , changes in rSO 2 are not adequate for predicting changes in PbtO 2 .
In this study, changes in PbtO 2 are often best correlated with changes in rSO 2 one minute earlier.This finding is both consistent with the theoretical mechanisms of each modality and previous findings in the literature.In a study of 42 TBI patients, Budohoski and colleagues noted that NIRS reacted earlier to changes in ABP and ICP as compared to PbtO 2 [35].This was recapitulated here in the findings of the VAR-modeled IRF plot.From these plots it can be seen that the response in change in PbtO 2 is both delayed and prolonged.In the case of change in rSO 2 , it is probable that the initial step rise, subsequent overcorrection, and eventual return to equilibrium may be reflective of cerebrovascular reactivity.The delayed and prolonged nature of the change in PbtO 2 may explain why continuous cerebrovascular reactivity metrics based on PbtO 2 have been found to be so dissimilar to those based on more responsive surrogates of cerebral blood flow or volume [36].As for the mechanism of this delay, PbtO 2 is a measure of the extracellular content of oxygen in brain tissue as it only measures dissolved oxygen in the interstitial fluid of the brain.NIRS-based rSO 2 measures microvascular oxygen saturation over a volume of brain as it utilizes deoxyhemoglobin (DeOxHgB) and oxyhemoglobin (OxHgB) as chromophores to scatter and the NIR light [31,32].Oxygen is primarily delivered to the brain in the form of OxHgB through the brain's microvasculature.An increase in OxHgB in the brain's microvasculature would be detected through the absorption of near-infrared light.However, prior to observing a change in oxygen content of the extracellular space of the brain, and therefore PbtO 2 , oxygen would need to disassociate from the OxHgB and diffuse into this extracellular space.A similar delay in decreases in cerebral oxygen content might also be explained by this mechanism.This mechanism is consistent with the findings of changes in PbtO 2 being correlated with changes in rSO 2 that occurred one minute earlier.This is a significant finding that may help guide further research into the flow of oxygen through the cerebral microenvironment.
This study found that changes in rSO 2 may, in a linear way, contain information about a change in PbtO 2 approximately one minute later.The statistical significance of this relationship held true even when the autocorrelative structure of PbtO 2 was accounted for.While this was not found in every patient, there are several reasons why this may be the case, the most obvious of which is that NIRS-based rSO 2 is prone to interference from extravascular blood collections, such as those seen in subgaleal, epidural, and subdural hematomas, as well as intraparenchymal hematomas.In the setting of TBI, these forms of interference are common.
Finally, this study indicates that changes in NIRS-based rSO 2 are inadequate on their own to predict upcoming changes in PbtO 2 .Despite using a methodology that may tend towards overfitting of the models, measures of change in rSO 2 were unable to reasonably predict changes in PbtO 2 .This was primarily evident when examining the degree of agreement between actual values of change in PbtO 2 and predicted values.This adds to the body of evidence that indicates rSO 2 is not an adequate alternative to PbtO 2 [20][21][22].Once again, while rSO 2 and PbtO 2 are in some ways measures of cerebral oxygenation, they interrogate entirely different compartments.There are likely several factors that influence changes in PbtO 2 that were not utilized in the models in this study.Hemoglobin concentrations (HgB), cerebral metabolic rate of oxygen (CMRO 2 ), partial pressure of oxygen in the arterial blood (PaO 2 ), and microvascular cerebral blood flow velocity (CBFV) are all factors that may modulate how changes in rSO 2 relate to changes in PbtO 2 .It is understandable why, despite containing information about upcoming changes in PbtO 2 , changes in rSO 2 could not adequately predict changes in PbtO 2 in this cohort.
The findings of this study have significant clinical implications.The first is that rSO 2 and PbtO 2 provide related but not equivalent information about brain physiology, and as such, rSO 2 , as a raw parameter, is not a viable non-invasive alternative to PbtO 2 for the monitoring of TBI patients.Further work is needed to better elucidate the potential role of rSO 2 in the multimodal monitoring of critically ill TBI patients.The related nature of rSO 2 and PbtO 2 may be leveraged in the future to improve care, such as in the prehospital setting where invasive monitoring is not possible.Secondly, given the identifiably delayed reaction of PbtO 2 to changes in ABP as compared to rSO 2 , as demonstrated in the CCF analysis and IRF plots, it is likely that rSO 2 is better suited for continuous indices of cerebrovascular reactivity.These indices require response parameters that react quickly to changes in ABP.There is a significant degree of interest evolving in the area of continuous cerebrovascular reactivity indices in the monitoring and management of moderate-to-severe TBI.

Limitations
There are limitations to this study that need consideration when interpreting its findings.The first is the relatively small size of the cohort in this study, with only 20 patients included.This necessitates the validation of these findings in a larger cohort before the findings can be fully incorporated into patient care as generalizability might be limited.However, concurrent NIRS-based rSO 2 and PbtO 2 monitoring is relatively uncommon, with this study being the first such analysis of high-resolution concurrent recordings in the literature.This cohort represents less than 7% of the full CAHR-TBI database indicating the rarity of simultaneous recordings of these parameters.This likely reflects both the relatively recent global adoption of PbtO 2 as a means of cerebral monitoring in TBI and the paucity of evidence for the use of NIRS in monitoring moderate-to-severe TBI.
A second limitation, brought on by the type of analysis performed, was the need to interpolate data.This, obviously, injects some inherent uncertainty into the study findings.Another limitation of this study is that information about factors that may interfere with rSO 2 measurements, such as extravascular blood collections, was not available.Additionally, the side of PbtO 2 monitor placement needed to be assumed, as this was also not available for all patients.
Finally, the models utilized in this analysis assumed a linear relationship between changes in rSO 2 and changes in PbtO 2 .This was done due to the lack of evidence suggesting a more appropriate alternative structure to this relationship.It is possible that a mathematically more complex model may prove more suited to describing this relationship.Complex supervised machine learning algorithms might be an obvious method to map this relationship better.While more complex mathematical models might provide an accurate prediction of PbtO 2 values based on rSO 2 levels, they are likely to increase computational complexity.This may limit utility at the bedside if computations are not possible in real time due to this increased complexity.

Future Work
The findings of this work lay the groundwork for additional research.First, these findings need to be validated in a larger multi-institutional cohort, where information about sources of NIRS interference is also available.Ideally, additional parameters, including CMRO 2 , PaO 2 , HgB concentration, and microvascular CBFV, would also be concurrently measured to fully elucidate the relationship between these two modalities.With a larger cohort and better characterized physiology, the complexity of this relationship may be better captured.
While NIRS-based rSO 2 has not gained traction as a stand-alone parameter in the management of TBI, there is increased interest in leveraging NIRS as a means of noninvasively interrogating cerebrovascular reactivity in TBI [37][38][39][40].This may mean larger datasets with concurrent rSO 2 and PbtO 2 monitoring may become available in the future.In such complex multimodal monitoring datasets, supervised and unsupervised classification and regression machine learning algorithms may provide additional insights into how these parameters interact with one another.While complex computational models may not be deployable at the bedside, they may act as inferential models to drive our understanding further.This may ultimately lead to new treatment paradigms in the management of moderate and severe TBI in the acute phase, including specific molecular targets.

Conclusions
In this multi-institutional exploratory analysis of a cohort of 20 TBI patients with concurrent rSO 2 and PbtO 2 monitoring, changes in PbtO 2 were found to correlate most significantly with changes in rSO 2 approximately one minute earlier.Through mixedeffects and time series modeling, changes in rSO 2 were found to often have a statistically significant linear relationship with changes in PbtO 2 that occurred one minute later.This was the case even when the hierarchical and autocorrelative structure of the data was considered.However, changes in rSO 2 were inadequate on their own to predict changes in PbtO 2 based on the poor agreement between modeled and actual changes in PbtO 2 .Given the uniqueness of this dataset, only a small number of subjects were available for analysis, limiting the confidence of these findings.In the future, a larger cohort, with additional parameters that influence cerebral oxygenation, is required to validate and better explain these findings.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/bioengineering10101124/s1,Supplementary File S1: Evaluation of Stationarity of PbtO 2 Signals, Supplementary File S2: Evaluation of residuals of the linear ∆PbtO 2 ~∆rSO 2 Lag1 models to ensure autocorrelative structure was accounted for, Supplementary File S3: Scatter plots of actual vs. predicted values of ∆PbtO 2 for each subject, Supplementary File S4: Bland-Altman plots to evaluate model agreement.

Figure 1 .
Figure 1.An example autocorrelative function (ACF) plot of PbtO2 can be seen in Panel (A), indicating non-stationarity, as correlations do not reach zero.Conversely, in Panel (B), the ACF plot of ΔPbtO2 for the same subject shows a rapid decline in the significant lags, indicating stationarity.The dashed blue line represents the significance levels which were set to a correlation level of +/−(2/N 1/2 ), where N is the number of samples.ACF = Autocorrelative Function, ΔPbtO2 = Change in brain tissue oxygen tension, PbtO2 = Brain tissue oxygen tension.

Figure 1 .
Figure 1.An example autocorrelative function (ACF) plot of PbtO 2 can be seen in Panel (A), indicating non-stationarity, as correlations do not reach zero.Conversely, in Panel (B), the ACF plot of ∆PbtO 2 for the same subject shows a rapid decline in the significant lags, indicating stationarity.The dashed blue line represents the significance levels which were set to a correlation level of +/−(2/N 1/2 ), where N is the number of samples.ACF = Autocorrelative Function, ∆PbtO 2 = Change in brain tissue oxygen tension, PbtO 2 = Brain tissue oxygen tension.

Figure 2 .
Figure 2.In Panel (A), the plot of the cross-correlative function (CCF) of ΔPbtO2 vs. ΔrSO2 for the minute-by-minute data of the cohort can be seen, with the most significant lag occurring at lag 1 (i.e., 1 min).In Panel (B), the plot of the cross-correlative function of ΔPbtO2 vs. ΔrSO2 for the 10seconds-by-10-seconds data for the data of the cohort can be seen, with the most significant lag occurring at lag 6 (i.e., at 1 min).The dashed blue line represents the significance levels, which were set to a correlation level of +/−(2/N 1/2 ), where N is the number of samples.CCF = Cross-correlative

Figure 2 .
Figure 2.In Panel (A), the plot of the cross-correlative function (CCF) of ∆PbtO 2 vs. ∆rSO 2 for the minute-by-minute data of the cohort can be seen, with the most significant lag occurring at lag 1 (i.e., 1 min).In Panel (B), the plot of the cross-correlative function of ∆PbtO 2 vs. ∆rSO 2 for the 10-seconds-by-10-seconds data for the data of the cohort can be seen, with the most significant lag occurring at lag 6 (i.e., at 1 min).The dashed blue line represents the significance levels, which were set to a correlation level of +/−(2/N 1/2 ), where N is the number of samples.CCF = Cross-correlative function, ∆PbtO 2 = Change in brain tissue oxygen tension, ∆rSO 2 = Change in region cerebral oxygen saturation.

Figure 3 .
Figure 3.A plot of Akaike information criterion (AIC) versus autoregressive order of the m variate vector autoregressive (VAR) model.There is limited improvement in AIC beyond an o of 6.

Figure 3 .
Figure3.A plot of Akaike information criterion (AIC) versus autoregressive order of the multi-variate vector autoregressive (VAR) model.There is limited improvement in AIC beyond an order of 6.

Figure 4 .
Figure 4. Panel (A) shows the modeled resulting response in change in regional cerebral oxygen saturation (ΔrSO2) to an orthogonal impulse in change in arterial blood pressure (ΔABP); black solid line.Panel (B) shows the modeled resulting response in change in regional cerebral oxygen saturation (ΔPbtO2) to an orthogonal impulse in change in arterial blood pressure (ΔABP); black solid line.The 95% confidence intervals are indicated by the red dashed line.Note that in both plots there is an initial rise followed by subsequent decline; however, this is prolonged in the response of ΔPbtO2.

Figure 4 .
Figure 4. Panel (A)shows the modeled resulting response in change in regional cerebral oxygen saturation (∆rSO 2 ) to an orthogonal impulse in change in arterial blood pressure (∆ABP); black solid line.Panel (B) shows the modeled resulting response in change in regional cerebral oxygen saturation (∆PbtO 2 ) to an orthogonal impulse in change in arterial blood pressure (∆ABP); black solid line.The 95% confidence intervals are indicated by the red dashed line.Note that in both plots there is an initial rise followed by subsequent decline; however, this is prolonged in the response of ∆PbtO 2 .

Figure 5 .
Figure 5.In Panel (A), the autocorrelative function (ACF) plot for the residuals of the linear mixedeffects (LME) model (ΔPbtO2~ΔrSO2Lag1) with random slope and intercept are seen.In Panel (B), the partial autocorrelative function (PACF) plot is seen.In both plots, a significant correlation is seen.This is in keeping with a model that has not accounted for the autocorrelative structure of its response variable.The dashed blue line represents the significance levels which were set to a correlation level of +/−(2/N 1/2 ), where N is the number of samples.CCF = Cross-correlative function, ΔPbtO2 = Change in brain tissue oxygen tension, ΔrSO2Lag1 = The one-minute lagged change in regional cerebral oxygen saturation.

Figure 5 .
Figure 5.In Panel (A), the autocorrelative function (ACF) plot for the residuals of the linear mixedeffects (LME) model (∆PbtO 2 ~∆rSO 2 Lag1) with random slope and intercept are seen.In Panel (B), the partial autocorrelative function (PACF) plot is seen.In both plots, a significant correlation is seen.This is in keeping with a model that has not accounted for the autocorrelative structure of its response variable.The dashed blue line represents the significance levels which were set to a correlation level of +/−(2/N 1/2 ), where N is the number of samples.CCF = Cross-correlative function, ∆PbtO 2 = Change in brain tissue oxygen tension, ∆rSO 2 Lag1 = The one-minute lagged change in regional cerebral oxygen saturation.

Figure 6 .
Figure 6.An example from a single subject of the Bland-Altman plot comparing actual and predicted values of ΔPbtO2 from the linear ΔPbtO2~ΔRSO2Lag1 model with autocorrelative structure accounted for.Generally poor agreement can be seen.PbtO2 = Brain tissue oxygen tension, ΔPbtO2 = Change in brain tissue oxygen tension, and ΔrSO2Lag1 = The one-minute lagged change in regional cerebral oxygen saturation.

Figure 6 .
Figure 6.An example from a single subject of the Bland-Altman plot comparing actual and predicted values of ∆PbtO 2 from the linear ∆PbtO 2 ~∆RSO 2 Lag1 model with autocorrelative structure accounted for.Generally poor agreement can be seen.PbtO 2 = Brain tissue oxygen tension, ∆PbtO 2 = Change in brain tissue oxygen tension, and ∆rSO 2 Lag1 = The one-minute lagged change in regional cerebral oxygen saturation.

Table 1 .
Patient demographics for the cohort.