Modified Beer-Lambert algorithm to measure pulsatile blood flow, critical closing pressure, and intracranial hypertension

We introduce a frequency-domain modified Beer-Lambert algorithm for diffuse correlation spectroscopy to non-invasively measure flow pulsatility and thus critical closing pressure (CrCP). Using the same optical measurements, CrCP was obtained with the new algorithm and with traditional nonlinear diffusion fitting. Results were compared to invasive determination of intracranial pressure (ICP) in piglets (n = 18). The new algorithm better predicted ICP elevations; the area under curve (AUC) from logistic regression analysis was 0.85 for ICP ≥ 20 mmHg. The corresponding AUC for traditional analysis was 0.60. Improved diagnostic performance likely results from better filtering of extra-cerebral tissue contamination and measurement noise.


Introduction
Many patients with acute brain injury experience intracranial hypertension that impairs cerebral metabolism and causes secondary hypoxic-ischemic damage [1][2][3][4].Invasive intracranial pressure (ICP) monitoring devices can be used to improve therapeutic management, but their risks often preclude or delay their use [5].Thus, the development of non-invasive ICP monitoring technologies has been investigated [6][7][8][9].One promising technology is diffuse correlation spectroscopy (DCS); DCS can non-invasively measure pulsatile heart-beat fluctuations in microvascular cerebral blood flow [10].To date, pilot studies have used DCS flow waveforms to predict ICP in non-human primates [11], adults with traumatic brain injury [12], and pediatric critical care patients [13,14].Resultant ICP predictions have been obtained from data via machine learning [11][12][13] and critical closing pressure (CrCP) models of the cerebral circulation [14].Irrespective of analysis approach, a well-known confound of DCS is its sensitivity to signal contamination from extra-cerebral tissues [15].Here we introduce a novel frequencydomain modified Beer-Lambert method to filter extra-cerebral contamination from the DCS flow waveforms, and we demonstrate its improved prediction of ICP elevations in a piglet model (equivalent to toddler-aged children).
In this work, we utilize model-based estimation of CrCP.CrCP is the sum of ICP and the compression pressure exerted by active wall tension from vasomotor tone on the arterioles [16].Since CrCP depends on ICP, the measurement of CrCP provides a proxy measure of intracranial hypertension (especially in periods of vasodilation when active wall tension is reduced).In fact, CrCP may be more clinically relevant than ICP, since it better reflects the factors which affect cerebral blood flow.While the difference between mean arterial blood pressure (MAP) and ICP is often referred to clinically as the cerebral perfusion pressure, the true pressure gradient across the arterioles also depends on vasomotor tone.Specifically, assuming that the outflow blood pressure at the distal end of the arterioles is CrCP [17], we refer to MAP -CrCP as the "actual" cerebral perfusion pressure (aCPP).Note, if pressure in the arteriole drops below CrCP, then the vessel collapses, and cerebral blood flow ceases.
By leveraging the aCPP concept, we can derive CrCP from concurrent measurements of pulsatile cerebral blood flow and arterial blood pressure during the cardiac cycle [17,18].Most previous studies have used transcranial Doppler ultrasound measurements of blood velocity in large cerebral arteries [18][19][20][21][22][23][24][25], which may not reflect arteriolar flow.DCS, however, measures blood flow in small arterioles [17], which depend more directly on the CrCP outflow pressure.DCS thus holds potential to measure CrCP with improved accuracy.To date, however, DCS measurements of CrCP have relied on nonlinear fitting of the DCS intensity autocorrelation signal to a homogeneous semi-infinite head model that ignores differences between extra-cerebral and cerebral tissues [17,[26][27][28][29][30].In addition, correlation noise in the measurements can make it difficult to resolve heart-beat fluctuations in blood flow.Herein, we report pilot results that show the modified Beer-Lambert method's mitigation of these confounds.

Vascular model of arterioles for deriving CrCP
The derivation of CrCP from pulsatile blood flow in cerebral arterioles (F) has been described elsewhere [17,27].Briefly, we employed a "resistive-only" approximation to model the cerebral arteriole vasculature compartment between the large arteries and the capillary bed.The in-flow blood pressure at the entrance to the arteriole compartment is P A , the resistance of the arterioles is R, and the outflow blood pressure at the distal end of the arteriole compartment is CrCP.CrCP is estimated based on the Ohms' law approximation relating the driving arteriolar pressure, P A −CrCP, to arteriolar blood flow (F): ( Here, t denotes time, and we assume CrCP and R remain constant over the measurement time scale.To obtain CrCP, the resistance R is first derived from pulsatile heart-beat fluctuations in P A and F. Specifically, the Fourier transform of Eq. ( 1) at the heart rate (f hr ) gives R = |P A (f hr )| / |F(f hr )|, where |P A (f hr )| and |F(f hr )| are the Fourier spectral amplitudes of P A and F at f hr (note, we assume |CrCP(f hr )| << |P A (f hr )|; throughout the manuscript, bolded variables refer to Fourier amplitudes).Next, CrCP is determined from R and the time average of Eq. ( 1): where the ⟨⟩ brackets denote the temporal means.The arteriolar flow pulsatility index, PI F ≡ |F(f hr )|/⟨F⟩, is directly measured with the optical DCS technique (see Sections 2.2, 2.3, 2.4).However, P A is not directly measured.It is instead estimated from measurement of systemic arterial blood pressure (ABP).Due to flow through the large arteries, P A is lower than ABP.In a rat model, P A was measured and found to be directly proportional to ABP, i.e., P A (t) = η ABP(t) [31].Substituting this proportionality into Eq.( 2), we obtain: Here, MAP is the systemic mean arterial blood pressure, and PI ABP ≡ |ABP(f hr )| / MAP.We assumed η = 0.6 based on the rat measurements [31].Herein, we measured PI F with 4 methods.

Measurement of flow pulsatility: nonlinear diffusion fitting (Methods 1, 2)
The traditional approach to measure PI F is via a "nonlinear diffusion fitting algorithm."Specifically, coherent near-infrared light delivered on the scalp is detected a distance ρ away.The DCS signal is the normalized temporal intensity autocorrelation function of the detected light emerging from the head (g 2 (τ)) at multiple delay-times (τ); τ typically ranges from 1 µs to ∼50 ms [32].To resolve heart-beat fluctuations in blood flow, g 2 (τ) curves are sequentially measured at high sampling rates (e.g., 20 Hz) [10]; we denote the curve measured at time t as g 2 (τ,t).Then, at each time t, F(t) is derived by nonlinear fitting of g 2 (τ,t) to the homogeneous semi-infinite medium solution of the correlation diffusion equation [32].More precisely, the correlation diffusion solution gives the normalized electric field autocorrelation function (g 1 (τ)), which is then converted to g 2 (τ) via the Siegert relation for fitting.To measure PI F , |F(f hr )| was computed from the Fourier transform of F(t), and ⟨F⟩ was derived from the nonlinear semi-infinite fit of ⟨g 2 (τ)⟩ (recall, the ⟨⟩ brackets denote the temporal means).
To explore the tradeoff between improving signal sensitivity to the brain and reducing the impact of correlation noise, we critically examined two methods that use this traditional approach to measure PI F .In Method 1, fitting is performed to the early part of the g 2 (τ,t) decay curve (improved brain sensitivity [33,34]).In Method 2, fitting is performed using essentially the whole decay curve (reduced correlation noise impact [35]).The specific delay-times used for fitting in both methods were derived from the time-average of the magnitude of g 1 (τ) across the time-scale of the PI F measurement, i.e., ⟨|g 1 (τ)|⟩.Note, ⟨|g 1 (0)|⟩ = 1 and ⟨|g 1 (τ)|⟩ decays towards zero with increasing τ [32].For Method 1, fitting was performed using τ for which ⟨|g 1 (τ)|⟩ > 0.7.For Method 2, fitting was performed using τ for which ⟨|g 1 (τ)|⟩ > 0.1.
Method 1 has improved brain sensitivity because the autocorrelation function decay times associated with long light paths from source to detector are relatively short, while the decay times associated with short light paths are relatively long [33,34].Thus, |g 1 (τ)| at shorter τ is inherently more sensitive to the deeper brain tissue.Fitting to a small number of short τ instead of to the whole curve, however, increases the impact of correlation noise [35].We chose the ⟨|g 1 (τ)|⟩ > 0.7 cutoff, instead of ⟨|g 1 (τ)|⟩ > 0.8 or ⟨|g 1 (τ)|⟩ > 0.6, to balance the tradeoff between improved depth sensitivity and reduced number of τ for fitting.This choice is supported by our simulations below (Sections 2.5, 3.1).
Method 2 is optimal for reducing errors from correlation noise because fitting is performed across the whole decay curve.Note, we chose the ⟨|g 1 (τ)|⟩ > 0.1 cutoff instead of ⟨|g 1 (τ)|⟩ > 0 for fitting because we wanted to minimize fitting to longer τ after the curve has fully decayed (the blood flow information in the signal is during the decay).
As a final aside, we are using ⟨|g 1 (τ)|⟩ to define the τ for fitting instead of using fixed τ ranges (e.g., defining the short τ range as 1µs to 10 µs) because the τ range of the |g 1 (τ)| decay depends on physiology (e.g., flow index, tissue optical properties) and source-detector distance.Using ⟨|g 1 (τ)|⟩ to define the short τ range helps ensure that the fitting is weighted to longer light paths.Further, we are using ⟨|g 1 (τ)|⟩ instead of ⟨g 2 (τ)⟩ because ⟨g 2 (τ)⟩ varies with changes in the Siegert β coefficient across experiments.

Measurement of flow pulsatility: semi-infinite modified Beer-Lambert algorithm (Method 3)
Figure 1 is a flow chart depicting the steps for an alternative method to determine PI F .This alternative method uses the DCS modified Beer-Lambert law in the semi-infinite tissue geometry [36].For each delay-time τ, and at each measurement timepoint t, we define the term "DCS optical density" as OD DCS (τ,t) ≡ −log(g 2 (τ,t) -1); log denotes the natural logarithm.In prior work, we derived the DCS modified Beer-Lambert law in the time-domain, which relates differential changes in OD DCS (τ,t) to the corresponding differential changes in F, tissue optical absorption (µ a ), and tissue optical reduced scattering (µ ′ s ) [36].The law is specifically the truncation of the Taylor series expansion of OD DCS (τ,t) to first order in F, µ a , and µ ′ s .Hence, the law is accurate for small changes in OD DCS (τ,t).  cm 2 /s) and (right) two exemplar intensity autocorrelation function curves corresponding to the temporal mean (red squares) and the flow at time t * (black circles) (there is no noise in this data).The DCS modified Beer-Lambert law relates temporal changes in OD DCS to corresponding temporal changes in flow and tissue optical absorption (µ a ) and reduced scattering (µ ′ s ) properties (the law is accurate for small OD DCS changes).(b) PI F is obtained by solving a system of equations specified by the Fourier transform of the DCS modified Beer-Lambert law at the heart rate (i.e., one Eq.for each τ; we assume that the changes in tissue optical properties during the cardiac cycle negligibly contribute to OD DCS changes).Inputs in this system of equations include ⟨F⟩ and the d F (τ) weighting factor; ⟨F⟩ is derived from a semi-infinite nonlinear fit of ⟨g 2 (τ)⟩; ⟨µ a ⟩ and ⟨µ ′ s ⟩ were inputs to this fit, and d F (τ) was calculated from ⟨F⟩ and the optical properties via Eq.(6).
We now assume that d a (τ)∆µ a and d s (τ)∆µ ′ s are negligible compared to d F (τ)∆F.The assumption that scattering factor changes during the cardiac cycle are negligible is reasonable since scattering is primarily due to particles that are extra-vascular [36].It is also reasonable to assume that d a (τ)∆µ a is small; fractional heart-beat changes in µ a are typically much smaller than fractional changes in F [10]; note also, even for the same fractional change, the size of d a (τ)∆µ a is considerably less than d F (τ)∆F [36].
Accordingly, from the frequency-domain variant of Eq. ( 4) obtained via Fourier transform, the heart-rate amplitudes of OD DCS and F are proportional, i.e., |OD DCS (τ,f hr )| = d F (τ)|F(f hr )| (note, the Fourier transforms of the constants ⟨OD DCS (τ)⟩ and ⟨F⟩ at f hr are zero).PI F , defined as |F(f hr )| / ⟨F⟩, is then readily derived from: Equation ( 7) is an over-determined system of equations; there is one equation for each delaytime τ in g 2 (τ).We used the Moore-Penrose pseudoinverse method implemented in MATLAB R2021a (pinv function, Mathworks, Natick, Massachusetts) to solve the system of equations for PI F .In principle, only one τ is needed to determine PI F , but the use of multiple τ improves signal-to-noise.For the pig experiment herein, we used the short τ for which ⟨|g 1 (τ)|⟩ > 0.7.

Measurement of flow pulsatility: two-layer modified Beer-Lambert algorithm (Method 4)
It's also straight-forward to extend the modified Beer-Lambert approach to heterogeneous tissue geometries [36].A simple useful heterogeneous geometry is the planar two-layer model of the head, which comprises a semi-infinite bottom layer (corresponding to the cortical regions of the brain) and a superficial top layer with thickness L (corresponding to extra-cerebral scalp and skull tissue).The two-layer extension of Eq. ( 4) is: where F c and F ec are the DCS flow indices in the cerebral (c) and extra-cerebral (ec) layers, respectively, µ a,c and µ a,ec are the corresponding optical absorption coefficients, and µ ′ s,c and µ ′ s,ec are the corresponding optical reduced scattering coefficients.The multiplicative weighting factors in Eq. ( 8) can be computed from the planar two-layer medium solution of the correlation diffusion equation, i.e., g 1,2L (τ) [36].We continue to assume that changes in tissue optical properties during the cardiac cycle negligibly contribute to OD DCS changes, which simplifies Eq. ( 8) to: Next, for the pig experiment herein, the pressure of the optical probe against the scalp was made large to attenuate scalp blood flow pulsatility.Further, because the brain is contained within the fixed skull, arteriolar blood flow pulsatility in the brain is expected to be larger than in the scalp and other peripheral tissues [39].For our analysis, we assumed the special case of Eq. ( 9) for which heart-beat fluctuations in extra-cerebral blood flow negligibly contribute to OD DCS changes.Then, as with the semi-infinite case, the two-layer flow pulsatility index, i.e., PI Fc ≡ |F c (f hr )| / ⟨F c ⟩, is solved from the frequency-domain variant of Eq. ( 9): Here, d F,c (τ) is computed numerically via [36] with δF c = 10 −5 × ⟨F c ⟩. ⟨F c ⟩ and ⟨F ec ⟩ were derived from nonlinear fitting of ⟨g 2 (τ)⟩ to 1+β|g 1,2L (τ)| 2 .We additionally assumed that tissue optical properties were homogeneous, i.e., ⟨µ a,c ⟩ = ⟨µ a,ec ⟩ = ⟨µ a ⟩, ⟨µ ′ s,c ⟩ = ⟨µ ′ s , ec ⟩ = ⟨µ ′ s ⟩, and we obtained them from fitting a semi-infinite head model to frequency-domain diffuse optical spectroscopy measurements.Finally, for the pig experiment herein, L was posthumously measured (see Section 2.6), and we used the short τ for which ⟨|g 1 (τ)|⟩ > 0.7 to solve Eq. (10). Figure 2 is a flow chart of the two-layer algorithm.

Impact of correlation noise on flow pulsatility measurements: simulated data
One reason we have developed this new algorithm is that correlation noise hinders determination of flow dynamics at the heart rate [10].As a result, the measured PI F underestimates the true PI F .Indeed, for severe noise levels, the observable heart rate signal can be negligible, and then the measured PI F is zero.Thus, an analysis algorithm that is more robust against correlation noise will be less prone to underestimation of PI F .As an initial test of which method can better recover the pulsatility index in the presence of noise, we generated simulated DCS data and compared the modified Beer-Lambert and nonlinear diffusion fitting estimates of PI F from the same dataset.
Briefly, we generated simulated g 2 (τ,t) with noise in the semi-infinite geometry for a 1-minutelong pulsatile F(t) waveform (2.5 cm source-detector distance, 20 Hz sampling).The generation of the data is described in detail in Appendix 1.Using a correlation noise model [40], varying levels of correlation noise were added based at different DCS intensity (photon count rate at the detector); we varied this count rate from 5 to 50 kHz in 5 kHz steps (note, the noise decreases with increasing intensity).For each intensity level, we generated 30 sets of pulsatile g 2 (τ,t) data for the cases of one detector and six detectors.For the six-detector case, which mimics our experimental data acquisition in piglets, six g 2 (τ,t) measurements from six detectors at the same source-detector distance were averaged together.
Next, we used the semi-infinite nonlinear diffusion fitting algorithm to compute PI F for each g 2 (τ,t) set (fitting performed at short τ for which ⟨|g 1 (τ)|⟩>0.7;see Section 2.2).Corresponding computations of PI F with the modified Beer-Lambert algorithm were obtained via Eq.( 7) (note, the same values of short τ were used to solve Eq. ( 7)).The means and standard deviations of PI F obtained with each algorithm were compared to the true PI F of the F(t) waveform.c ) and a superficial extracerebral layer with thickness L (flow index, F ec ; tissue optical properties, µ a,ec , µ ′ s , ec ).PI Fc is obtained by solving a system of equations specified by the Fourier transform of the 2-layer DCS modified Beer-Lambert law at the heart rate.We assume: a) small OD DCS changes solely due to F c variation, and b) homogeneous tissue optical properties.Inputs in this system of equations include ⟨F c ⟩ and the d F,c (τ) weighting factor.⟨F c ⟩ and ⟨F ec ⟩ are derived from a two-layer nonlinear fit of ⟨g 2 (τ)⟩.The measured optical properties and L were inputs in this fit; d F,c (τ) was calculated via Eq.( 11).
Finally, using the same simulated data, we tested the effects of using different ranges of τ for deriving the nonlinear diffusion and modified Beer-Lambert estimates of PI F .As mentioned in Section 2.2, restricting fitting to short τ increases the depth sensitivity [33,34], but it also increases the impact of correlation noise [35].In addition, the modified Beer-Lambert algorithm can be prone to errors from the larger OD DCS changes that occur at longer τ [36].To explore these effects, we obtained nonlinear diffusion fitting and modified Beer-Lambert estimates of PI F using other τ ranges for the fittings.Specifically, a g 1 (τ) cutoff was employed to define the τ range for fitting, i.e., the ⟨|g 1 (τ)|⟩ > cutoff, for cutoffs between 0.1 and 0.9.
Importantly, to improve the accuracy of all PI F estimates, we sought to remove the "white noise" level present at all frequencies in the |F(f)|/⟨F⟩ Fourier spectra (see Section 2.7).

In vivo piglet model of intracranial hypertension
To compare CrCP obtained with nonlinear diffusion fitting and modified Beer-Lambert algorithms to invasive ICP, we used a piglet model of intracranial hypertension [41,42].The experiment complied with the ARRIVE 2.0 guidelines, and all procedures were approved by our Institutional Animal Care and Use Committee in accordance with the National Institutes of Health Guide for the Care and Use of Laboratory Animals.
One-month-old female piglets (n = 25, weight: 10.3 ± 0.8 kg) were measured with continuous physiologic monitors (see Appendix 2).Briefly, invasive intraparenchymal ICP and microdialysis sensors were placed in the right hemisphere, a hybrid DCS and frequency-domain diffuse optical spectroscopy (FD-DOS) probe was placed on the left forehead, and arterial blood pressure was monitored via a catheter inserted in the right femoral artery (Fig. 3(a)).Details about the probe and the optical instrument are described in Appendix 3. Note, FD-DOS measures tissue blood O 2 saturation (StO 2 ), total hemoglobin concentration (THC), and the tissue absorption (µ a ) and reduced scattering (µ ′ s ) properties used for DCS fitting [37,38].Note also, the microdialysis At the end of the time-period for each ICP level, arterial blood gas and microdialysis samples were taken (black circles).(c) Exemplar arterial blood pressure (ABP, mmHg) and optical blood flow index (F, 10 −8 cm 2 /s) time-series in a piglet with invasive ICP of 20 mmHg (F was derived from semi-infinite diffusion fitting to short delay-times (Method 1), the pink circle and green square denote exemplar systolic and diastolic timepoints).(d) (Top) Individual optical intensity autocorrelation function (g 2 (τ)) measurements at the corresponding systolic (pink) and diastolic (green) timepoints labeled in panel (c) (source-detector distance, 2.5 cm; wavelength, 850 nm).(Bottom) Exemplar mean g 2 (τ) across 1 minute of data acquired at the 20 mmHg ICP level with its semi-infinite diffusion fit (solid line).(e) Normalized Fourier spectral amplitudes of 1-minute F time-series derived from semi-infinite diffusion fitting to short delay-times (Method 1, top), and the whole g 2 (τ) curve (Method 2, bottom).White noise levels (see main text) of 0.005 and 0.004, respectively, were subtracted from the heart-rate amplitudes to compute F pulsatility indices (PI F ) of 0.106 and 0.155 for the short τ and whole curve fitting methods.This results in critical closing pressure (CrCP) measurements, respectively, of −11.1 and 2.5 mmHg for the short τ and whole curve fitting methods (via Eq. ( 3); MAP = 53 mmHg, PI ABP = 0.143).(f) Normalized Fourier spectral amplitudes for blood flow obtained with the short τ semi-infinite (Eq.( 7); Method 3, top) and 2-layer (Eq.( 10); Method 4, bottom) modified Beer-Lambert algorithms, plotted against frequency.White noise levels of 0.016 and 0.027, respectively, were subtracted from the heart-rate amplitudes to compute PI F , i.e., 0.177-0.016= 0.161 and 0.306-0.027= 0.279.This results in semi-infinite and 2-layer CrCP estimates of 3.6 and 15.5 mmHg, respectively.sensor measures the concentrations of pyruvate, lactate, glycerol, and glucose concentrations in the extracellular interstitial fluid.
After a 10-minute baseline, graded elevations of ICP were induced (Fig. 3(b)) by infusing mock cerebrospinal fluid through an external ventricular drain (EVD) catheter in the right lateral ventricle (see Appendix 2).ICP was increased in ∼5 mmHg increments up to ∼60 mmHg.Each ICP level was maintained for 10 minutes.At the end of baseline and each ICP level, cerebral microdialysis and arterial blood gas samples were taken.Microdialysis samples were measured with the ISCUS Flex Analyzer (mDialysis, Stockholm, Sweden), and arterial blood gas samples were analyzed with a GEM 3000 (Instrumentation Laboratory, Bedford, MA).All animals were immediately euthanized after the final ICP level via a bolus injection of potassium chloride.Then, the scalp and skull thicknesses underneath the FD-DOS/DCS probe were posthumously measured with calipers (measurements were done at 4 distinct locations; L in Fig. 2 was the mean scalp + skull thickness across the 4 locations).
Figure 3 shows exemplar measurements in a piglet at the 20 mmHg ICP level.For Method 1, the F(t) time-series was noisy and the computed PI F gave a non-physiological negative CrCP of −11.1 mmHg.For Method 2, fitting to the whole g 2 (τ) curve resulted in a larger PI F and a CrCP of 2.5 mmHg.For Method 3, d F (τ) in Eq. ( 5) was first computed based on ⟨µ a ⟩, ⟨µ ′ s ⟩, and ⟨F⟩ measurements; note, ⟨F⟩ was obtained from short τ semi-infinite nonlinear fitting to ⟨g 2 (τ)⟩.Equation ( 7) was then solved using the short τ range (Fig. 3(f), top, shows the solution at the heart rate and other frequencies), and the CrCP derived from the PI F measurement was 3.6 mmHg.
Finally, for Method 4, we first fit ⟨g 2 (τ)⟩ to the planar two-layer correlation diffusion solution to extract ⟨F c ⟩ and ⟨F ec ⟩ (fitting was performed to the whole ⟨g 2 (τ)⟩ curve).Inputs to the fit were the measured extra-cerebral layer thickness, the Siegert β coefficient, and the homogeneous optical properties.The fit constrained ⟨F c ⟩ between 0 and 10 −6 cm 2 /s with ⟨F ec ⟩≥0.The two-layer correlation diffusion solution is shown elsewhere [36,44].Herein, we used an approximation of it to mitigate numerical errors, i.e., Eq. ( 8) in Kienle and Glanzmann [45], and we numerically evaluated the solution's integral via a Gauss-Laguerre quadrature of 1000 points.Then, d F,c (τ) was computed, and Eq. ( 10) was solved using the short τ range (Fig. 1(f)).The resultant CrCP was 15.5 mmHg.
Note, a non-zero "white noise" level is present at all frequencies in the |F(f)|/⟨F⟩ Fourier spectra (see Fig. 3).Thus, for all four methods, we sought to remove the white noise prior to computing PI F .Specifically, the white noise level was estimated as the mean |F(f)|/⟨F⟩ across the 1.2f hr to 1.8f hr frequency range to avoid contributions from harmonics.PI F was then computed from |F(f hr )|/⟨F⟩ after subtracting the white noise level.For more discussion about the white noise and its removal, see Appendix 4.

In vivo statistical analysis
To test our hypothesis that the use of novel modified Beer-Lambert methods (Method 3 and 4) improves prediction of intracranial hypertension compared to traditional nonlinear diffusion fitting methods (Methods 1 and 2), we developed univariate logistic regression models for detecting intracranial hypertension (IH) based on CrCP.We defined IH as invasive ICP ≥ 20 mmHg [2].Animals with baseline ⟨F⟩≤10 −10 cm 2 /s were excluded; this very low ⟨F⟩, which is comparable to the "biological zero" signal of DCS (e.g., the BFI observed, for example, during zero-flow circulatory arrest conditions [46,47]) is an indication that the probe was not correctly positioned over the brain.Animals were also excluded if visual inspection showed no observable Fourier spectral amplitude peak at the heart rate for F. All statistical analyses were performed on the means of data computed across each 10-minute ICP level (Fig. 3(b)).
For each CrCP analysis method, CrCP data were classified as normotension (ICP < 20 mmHg) and IH (ICP ≥ 20 mmHg) based on invasive ICP monitoring.Receiver operating characteristic (ROC) curves were computed and the area under the ROC curve (AUC) quantified to evaluate predictive performance.Specifically, a logistic regression model was fit to the normotension and IH CrCP data to derive a probability for the presence of IH, i.e., P M = 1 / (1 + e −M ), with M ≡ β 0 + β 1 ×CrCP and the weighting parameters β 0 and β 1 determined from the fit.A ROC curve was then generated by varying the P M threshold for predicting the presence of IH and plotting the sensitivity (true positive rate) versus 1−specificity (false positive rate) for each P M threshold.We report the CrCP threshold for predicting IH based on the maximal sum of sensitivity and specificity.In addition, we quantified the measurement failure rate for each CrCP analysis method as the percentage of CrCP measurements that produced negative pressures (negative CrCP values are non-physiologic).Finally, to assess the performance of CrCP in relation to that of other physiological metrics, we quantified performance for univariate IH detection based on F c , StO 2 , and MAP.
We next characterized the quantitative relationship of ICP with CrCP estimates, i.e., we performed univariate linear mixed effects regressions of invasive ICP with CrCP for each CrCP analysis method (grouping variable, animal ID #).Here, our aim was to probe the linear relationship between CrCP and ICP, so we excluded failed CrCP measurements (i.e., measurements with negative CrCP).Negative CrCP measurements were not excluded in the logistic regression analyses that assess the performance of a single CrCP measurement.

Impact of correlation noise on flow pulsatility (PI F ) measurements: simulated data
The true PI F of the 1-minute simulated F time-series without noise was 0.264 (Fig. 5(a)), and exemplar g 2 (τ) data derived from the F time-series is shown in Fig. 5(b) with and without addition of correlation noise.For the data with noise, Fig. 5(c) shows exemplar normalized Fourier spectral amplitudes obtained with the "short τ" semi-infinite modified Beer-Lambert and nonlinear diffusion fitting methods, i.e., Methods 3 and 1, respectively (see Fig. 4).The PI F recovered with the modified Beer-Lambert approach was more accurate than the PI F recovered with the nonlinear diffusion fitting approach across all DCS intensities (Fig. 5(d) and 5(e); lower intensities have higher correlation noise).A comparison of Fig. 5(d) to 5(e) also shows noticeable improvement in PI F measurement accuracy when the signals from 6 detectors are averaged in parallel to improve signal-to-noise.
In Fig. 5(f), we fix the signal intensity at 30 kHz, corresponding approximately to the regime in which the short τ modified Beer-Lambert scheme PI F approached its true value for six detectors, and we study the effects of using more delay-times in the fitting (e.g., ⟨|g 1 (τ)|⟩>0.1).In this situation, the accuracies of the Modified Beer-Lambert and nonlinear diffusion approaches worsen and improve, respectively, depending on delay-time.This is expected.The Modified Beer-Lambert approximation is worse for the larger OD DCS changes that arise at longer delaytimes; nonlinear fitting to more delay-times is better for homogeneous media.Notice also that compared to the 0.7 ⟨|g 1 (τ)|⟩ cutoff, accuracy worsens for both approaches at the 0.8 and 0.9 ⟨|g 1 (τ)|⟩ cutoffs.This supports our choice of the 0.7 cutoff for Method 3.

In vivo prediction of intracranial hypertension
Of the 25 piglets, 6 were excluded because their baseline ⟨F⟩ was ≤10 −10 cm 2 /s.One additional piglet was excluded because DCS signal-to-noise was too low to observe pulsatile heart-beat fluctuations.For the remaining 18 piglets, we acquired data across 139 ICP levels in total.On average (mean ± SD) we recorded 8 ± 2 ICP levels per piglet.The average DCS intensity and Siegert β coefficient across piglets at baseline were 36 ± 10 kHz and 0.46 ± 0.04 (mean ± SD).The mean scalp, skull, and scalp + skull thicknesses at the optical probe location were 0.28 ± 0.04 cm, 0.34 ± 0.06 cm, and 0.62 ± 0.07 cm, respectively.Finally, the average measurement of the cerebral to extra-cerebral blood flow index ratio (F c /F ec ) at baseline obtained from two-layer diffusion fitting was 14.5 (6.6, 49.0) (median (interquartile range)).
Four methods were used to compute CrCP (see Fig. 4).Logistic regression analyses showed that the two-layer and semi-infinite modified Beer-Lambert estimates of CrCP were both predictive of intracranial hypertension (Fig. 6, Table 1); the AUCs were 0.85 (95% CI: 0.78-0.90)and 0.79 (95% CI: 0.71-0.86),respectively.At their optimal CrCP decision thresholds for indicating IH detection, i.e., CrCP > 19 mmHg for the two-layer, and CrCP > 13 mmHg for semi-infinite, both methods were more specific than sensitive (Table 1).Considering the two-layer case, for example, this means that CrCP < 19 more accurately indicates intracranial normotension than CrCP > 19  measurements for the systolic and diastolic timepoints without (top) and with (bottom) noise in the semi-infinite geometry.Noise was added assuming a signal intensity of 30 kHz and averaging across 6 detectors placed at the same source-detector distance of 2.5 cm.(c) For the same noise level as in the bottom panel of (b), the normalized Fourier spectral amplitudes of F obtained with the semi-infinite modified Beer-Lambert (Method 3) and nonlinear diffusion fitting (Method 1) algorithms are plotted against frequency (short delay-times τ, corresponding to ⟨|g 1 (τ)|⟩>0.7,were used for fittings).For computation of PI F , white noise levels of 0.032 and 0.008 for modified Beer-Lambert and diffusion fit, respectively, were subtracted from the signal amplitude at the heart rate (1.45 Hz).(d) Mean PI F (±SD) obtained with the modified Beer-Lambert (red circles, Method 3) and nonlinear diffusion fitting (black squares, Method 1) algorithms on simulated data with varying noise; noise is determined based on the signal intensity and is reduced by averaging across 6 detectors.Notice that PI F approaches its true value for the modified Beer-Lambert fit at signal intensities between 20 and 30 kHz.(e) Same plot as panel (d), except simulations were done using one detector instead of six.(f) For a 30 kHz intensity signal and six detectors, the mean PI F (±SD) obtained for both algorithms are plotted using different ⟨|g 1 (τ)|⟩ fitting cutoffs (for lower cutoffs, more τ are used in the fitting).
indicates intracranial hypertension.In contrast to the modified Beer-Lambert estimates, the predictive performance of the nonlinear diffusion fitting estimates of CrCP was weak (AUC ≤ 0.60, Fig. 6 and Table 1).MAP, StO 2 , and F c metrics also exhibited weak performance for detecting IH (Table 1).
The CrCP measurement failure rates (i.e., percentage of CrCP measurements that produced negative pressures) were dramatically lower for the modified Beer-Lambert methods than the nonlinear diffusion fitting methods (Table 2).For example, 74% of the CrCP measurements obtained with short τ nonlinear diffusion fitting were negative, while the two-layer modified Beer-Lambert method produced only 1 negative measurement (0.7%).We next consider the linear relationship between CrCP and invasive ICP.Table 2 reports the slope and intercept of the linear mixed effects best-fits between CrCP and invasive ICP for the four CrCP analysis methods.For the whole curve nonlinear diffusion fitting method and both modified Beer-Lambert methods, increasing CrCP was associated with increasing ICP (p < 0.001).There was no significant association, however, for CrCP computed with the short τ nonlinear diffusion fitting method (p = 0.06).In addition, the two-layer CrCP intercept of the best-fit was positive Absorption (µa (850 nm)) (cm−1) 0.15 (0.14, 0.16) 0.15 (0.13, 0.17) 0.14 (0.13, 0.17) Lactate (mM) 1.0 (0.9, 1.2) 0.9 (0.7, 1.0) 0.9 (0.6, 1.0) a Data are reported as medians (interquartile ranges).Method 4 estimate of CrCP; CPP = MAP -ICP.
(p < 0.001), but the intercepts of the regressions for the other analysis methods were not different from zero.Finally, both the semi-infinite and two-layer modified Beer-Lambert CrCP estimates exhibited stronger correlations with ICP than the nonlinear diffusion fitting estimates.Note that the slopes of the best-fit lines for all four analysis methods were less than one, which shows that CrCP measures underestimate ICP at higher pressures.Table 3 reports average physiology during baseline, moderate IH, and severe IH.

Discussion
In this study, we introduced and demonstrated a frequency-domain modified Beer-Lambert analysis algorithm for measuring pulsatile blood flow and critical closing pressure (CrCP) with the non-invasive diffuse correlation spectroscopy (DCS) technique.CrCP measurements in piglets, derived with the new algorithm, were predictive of intracranial hypertension.The prediction was independent of mean arterial pressure and other diffuse optical metrics of blood flow and tissue oxygen saturation.

Modified Beer-Lambert versus nonlinear diffusion fitting algorithms
Traditionally, pulsatile blood flow with DCS has been obtained via nonlinear fitting of the DCS signal to the semi-infinite correlation diffusion solution.Our simulations show this approach is prone to noise-induced underestimations in blood flow pulsatility, especially when only short delay-times are used for fitting.These underestimations can in turn drive CrCP measurements negative (see Eq. ( 3)).Indeed, 74% of our CrCP measurements in piglets obtained with semiinfinite diffusion fitting to short delay-times were negative.The correlation noise confound is reduced by fitting to the whole DCS g 2 (τ) curve, but 42% of the CrCP measurements were still negative, and the correlation with ICP was weak.Since scalp/skull pulsatility is usually lower than brain pulsatility [39], extra-cerebral tissue contamination is another confound that can lead to underestimations in brain pulsatility measurements.This effect can be reduced by fitting to only short delay-times in the DCS signal to increase measurement depth sensitivity [33,34], but as discussed above, this also increases the impact of correlation noise.A key benefit of the modified Beer-Lambert approach, which we observed in our simulations, is improved signal-to-noise when using short τ ranges.This coupling of improved signal-to-noise with increased brain sensitivity likely explains the stronger correlation between ICP and the semi-infinite modified Beer-Lambert estimate of CrCP than between ICP and the nonlinear diffusion fitting estimate of CrCP.
The two-layer modified Beer-Lambert results demonstrate ease of use in heterogeneous geometries as another key benefit of the approach.Layered tissue models can also be used for nonlinear diffusion fitting to improve accuracy, but noise-induced instability in the fitting is a known problem [48][49][50][51][52][53][54][55].The modified Beer-Lambert approach, by contrast, only requires a simple linear inversion to obtain cerebral blood flow pulsatility.In future work, extending the modified Beer-Lambert approach to more realistic head models and using alternative estimations of the extra-cerebral layer thickness [51,56] may further improve accuracy.
We note that in our prior study in infants with hydrocephalus [14], we used the whole-curve non-linear diffusion fitting method to compute CrCP, and we observed a stronger correlation (R = 0.48) between CrCP and ICP than was observed in the present study with the method in piglets (i.e., R = 0.24).In this prior work, due to their thin scalp/skull, the infant measurements were carried out using a 2.0 cm source-detector distance (a 2.5 cm distance was used for the piglets).As a result of the shorter source-detector distance and thinner scalp/skull, DCS correlation noise and extra-cerebral tissue contamination might have been less confounding in the infant measurements.
We further note that promising alternative approaches to reduce the DCS correlation noise confound in the traditional nonlinear fitting analysis include averaging flow waveforms across many cardiac cycles [26,29] and combining DCS measurements with near-infrared spectroscopy photoplethysmography [30].Another proposed method to calculate CrCP uses high temporal resolution to perform linear extrapolation during diastole when the relationship between F and CrCP may be more linear [26].Our modified Beer-Lambert approach has the advantage of being less sensitive to the subtle differences between the shapes of the peripheral and cerebral ABP waveforms (the derivation of CrCP assumes that measured peripheral ABP waveforms are equivalent to cerebral ABP waveforms) [27].These approaches, however, are not mutually exclusive, and the modified Beer-Lambert law for flow could be applied in the time-domain to derive F(t) for either machine learning or diastole-only methods, which may also benefit from removal of superficial signals.

Delay-time range choices for modified Beer-Lambert algorithms
Recall, the DCS modified Beer-Lambert law (Eq.( 4)) is a good approximation for small changes in OD DCS .Importantly, if the multiplicative weighting factors in Eq. ( 4) are small, then ∆OD DCS (τ) can still be small for large tissue hemodynamic changes.For short τ, the weighting factors d F (τ), d a (τ), and d s (τ) are all close to zero, and at long τ, all 3 weighting factors deviate substantially from zero [36].Accordingly, the DCS modified Beer-Lambert law is usually accurate at short τ but breaks down at long τ.This behavior is evident in our simulations (Fig. 5(f)).When longer τ are used, i.e., τ corresponding to conditions wherein ⟨|g 1 (τ)|⟩ < 0.6, then the accuracy of the PI F estimate is diminished.
Fortunately, the short τ range is optimal for brain monitoring because the DCS signal at short τ is inherently more sensitive to the deeper brain tissue (see Section 2.2).Indeed, as we discussed above (Section 4.1), the superior in vivo performance of the semi-infinite modified Beer-Lambert method (Method 3) compared to the whole-curve nonlinear diffusion fitting method (Method 2) is likely due to this increased brain sensitivity.
Note however, we found (see Fig. 5(f)) that the use of extremely short τ ranges, i.e., corresponding to ⟨|g 1 (τ)|⟩ > 0.8 or ⟨|g 1 (τ)|⟩ > 0.9, decreased accuracy compared to ⟨|g 1 (τ)|⟩ > 0.7.This effect arises because the impact of correlation noise is substantial when only a few τ are used for linear inversion of Eq. ( 7); the smaller ∆OD DCS (τ) and higher correlation noise at the very short τ also makes it harder to distinguish signal from noise.(As an aside, when we performed simulations without adding correlation noise, the PI F errors for the 0.8 and 0.9 ⟨|g 1 (τ)|⟩ cutoffs were not significant (data not shown).)Thus, our simulations suggest the 0.7 ⟨|g 1 (τ)|⟩ cutoff is optimal for in vivo brain monitoring.Use of a 0.6 ⟨|g 1 (τ)|⟩ cutoff is another possibility, but although the noise levels for the 0.6 and 0.7 cutoffs in Fig. 5(f) are similar, there is more risk of the modified Beer-Lambert law breaking down for the 0.6 cutoff.
Finally, in principle, only one τ is needed to determine PI F with Eq. ( 7).In practice, we have found that this single τ estimate is noisier than the multiple τ estimate.For example, in our simulation (at 30 kHz intensity, average across 6 detectors), the estimate of PI F with a single τ (such that ⟨|g 1 (τ)|⟩ = 0.7) was 0.27 ± 0.03 (mean ± SD), and the corresponding estimate with multiple τ (with ⟨|g 1 (τ)|⟩ > 0.7) was 0.27 ± 0.01.
Future study is needed to confirm these findings for different source-detector separations, heterogeneous geometries, and different absolute flow indices and tissue optical properties.

CrCP predicts intracranial hypertension more robustly than other optical metrics
We have found that the CrCP prediction of IH was substantially better than the prediction achieved with the optical F c and StO 2 metrics.Many confounds, e.g., tissue optical property errors, scalp/skull thickness heterogeneity, head curvature, and hematocrit [48,50,54,[56][57][58][59], are known to impact the absolute F c metric.By contrast and importantly, the flow pulsatility index used to calculate CrCP is a fractional flow change, and fractional flow changes are more robust to these confounds, e.g., especially head curvature and baseline tissue optical property errors (scalp/skull thickness errors and temporal changes in optical properties can be more problematic) [36,50,51,54,56].Future work is needed to fully understand the effects of these DCS confounds on the flow pulsatility measurement.
In a different vein, the StO 2 metric is influenced by changes in both cerebral perfusion and oxygen metabolism [37].Thus, the minimal change in StO 2 and ∼30% decrease in cerebral blood flow observed from baseline to IH (see Table 3) suggests that the blood flow decrease was almost matched by a corresponding decrease in cerebral metabolic rate of oxygen (CMRO 2 ).The decreased CMRO 2 level, though, was not sufficient on its own to sustain brain metabolism (the increased microdialysis LPR shown in Table 3 indicates increased anaerobic brain metabolism).Interestingly, linear mixed effects regression analyses showed that F c , StO 2 , and the optical aCPP were negatively associated with increasing LPR (p < 0.001, data not shown), though the Pearson correlations were modest, i.e., R = −0.24for aCPP, R = −0.12 for F c , and R = −0.12 for StO 2 (recall, aCPP = MAP -CrCP; we used the Method 4 CrCP estimate).Despite the modest correlations, the associations motivate the future use of these optical metrics for detecting metabolic distress (as a reference comparison, the Pearson correlation between invasive CPP and LPR was R = −0.34(p < 0.001)).In a future publication, we will focus on the cerebral physiological effects of IH.

Limitations
A major limitation in the DCS measurement of CrCP is the uncertainty about the η factor in Eq. ( 3), which accounts for the blood pressure drop across the large cerebral arteries.Herein, we assumed η=0.6 based on rat measurements, but this assumption is likely in need of refinement for piglets.This factor may even be ICP-dependent.If, for example, vasodilation during severe IH resulted in a minimal pressure decrease from arteries to arterioles (i.e., η ≈ 1), the computed median CrCP during severe IH in Table 3 would be 47 mmHg, which is close to the median ICP.This suggests that assuming a constant η blunts the true CrCP change between moderate IH and severe IH.Further, CrCP might be expected to be higher than ICP at baseline because of vasomotor tone.The closeness of median CrCP and ICP at baseline (Table 3) thus suggests the true η in piglets at baseline is higher than 0.6.η may also vary across subjects.To aid clinical translation, schemes for subject-specific estimation of η are needed.
Measurement accuracy can probably be improved from better modeling of the arterial vascular tree.In the present work, arterioles were modeled as a resistor.We previously investigated an alternative two-compartment Windkessel model that incorporates compliance, i.e., the arterioles were modeled as a resistor and capacitor in parallel [17].This is a more realistic model than the "resistor-only" approach, and the phase difference (time lag) between pulsatile arteriolar blood pressure and cerebral blood flow waveforms provides information about compliance [17].Measurement of this phase difference is difficult in practice, however, as it is sensitive to location of the peripheral arterial measurement of ABP.Future work is needed to improve understanding of these transit time differences before the Windkessel or other vascular models could be routinely implemented.
A third limitation is the assumption that CrCP is constant over the time scale of the measurement.In practice, this means that the well-known heart-beat fluctuations in ICP [39], and thus in CrCP, are assumed to be negligible compared to heart-beat fluctuations in arteriolar blood pressure.Future work is needed to investigate measurement errors caused by ICP pulsatility.Finally, we assumed that scalp flow pulsatility in the piglets was negligible due to high optical probe pressure against the scalp, and we assumed that contribution of pulsatile changes in tissue optical properties to the pulsatile DCS signal were negligible.These assumptions may need to be relaxed in future work.

Conclusion
Optically measured critical closing pressure (CrCP) obtained with a frequency-domain modified Beer-Lambert algorithm was predictive of intracranial hypertension in piglets.The prediction was independent of mean arterial pressure and other diffuse optical metrics of blood flow and tissue oxygen saturation.When CrCP was obtained with traditional nonlinear diffusion fitting in piglets, measurement failure rate (i.e., fraction of negative CrCP measurements) was high, and the correlation with invasive intracranial pressure was much weaker.The key to the success of the modified Beer-Lambert approach was its robustness to correlation noise and its effective filtering of extra-cerebral tissue signals.
was placed ∼1 cm deep in the brain parenchyma through a burr hole in the right parietal bone.Finally, for continuous DCS and FD-DOS monitoring, a hybrid optical probe was sutured to the skin on the left forehead lateral to the midline and anterior to the crown [46].Throughout the study, ABP, invasive ICP, and vital signs were synchronized and recorded with a PowerLab system (1 kHz sampling, ADInstruments, Australia).The mock cerebrospinal fluid used to elevate ICP was dextrose 5% in one-fourth normal saline warmed to body temperature [41,42].The infusion rate was initially set to 4 mL/hr to increase ICP; the infusion rate needed to achieve maximal ICP was highly variable, ranging from 30 mL/hr to ∼100 mL/hr.

Appendix 3: diffuse optical instrument and probe
DCS and FD-DOS measurements were performed in parallel using a single fiberoptic probe [46].The DCS source was a 850-nm laser (DL850-100S, CrystaLaser, Reno, NV) coupled to a 90°bend borosilicate fiber bundle (1-mm diameter, 0.66 NA, 3-meter long, Fiberoptics Technology, Pomfret, CT).For DCS detection, six avalanche photodiodes (SPCM-AQ4C, Excelitas, Quebec, Canada) were connected to a fan-out bundle of single-mode fibers (780 HP, Thorlabs, Newton, NJ) placed 2.5 cm from the source; the bundle (made by Fiberoptic Systems, Simi Valley, CA) was terminated on the common end with a 3-mm right angle prism (N-SF11, Edmund Optics, Barrington, NJ) that was placed against the skin.The intensity autocorrelation function of the signal measured with each photodiode was continuously computed at 80 delay-times (between 1 µs and 4 ms) with a multiple-τ software correlator (50 ms integration time, or 20 Hz sampling) [10].The mean g 2 (τ) across the 6 detectors was then used.
FD-DOS measurements were acquired with a commercial instrument (Imagent, ISS, Champaign, IL).Specifically, the amplitude and phase of diffuse photon density waves induced by intensity-modulated sources (110 MHz) were measured at four wavelengths (690, 725, 785, 830 nm) and four source-detector distances (1.5, 2, 2.5, 3 cm).FD-DOS and DCS detection were performed through the same right-angle prism on the skin; for FD-DOS, a borosilicate fiber bundle (2 mm diameter, 0.55 NA) was placed circumferentially around the single-mode fibers at the center.Source light from each wavelength was delivered to four FD-DOS source locations via fan-out bundles of 0.66 NA borosilicate fiber (90°-bend 2.5-mm diameter common end, 1-mm diameter leg ends, Fiberoptics Technology).Note, the DCS source fiber bundle was part of the fan-out bundle at the 2.5 cm source-detector distance.An 850-nm short-pass filter (OZOptics, Ontario, Canada) in front of the FD-DOS detector blocked DCS light.FD-DOS measurements of StO 2 and THC (10 Hz sampling) were performed simultaneously with DCS.Assuming constant water volume fraction of 0.75, a semi-infinite tissue model was used to derive StO 2 and THC [46].Tissue absorption and reduced scattering at the DCS wavelength were determined from StO 2 , THC, and a Mie scattering model [57,61].
Note, since we did not have spectral filters in front of the DCS detectors to block FD-DOS light, FD-DOS light contamination in the DCS signals reduced the Siegert β coefficient from the expected value of 0.5.In practice, the amount of contamination varied with the changes in tissue optical properties within and across subjects.Thus, we decided to fit for β in the computation of CrCP (see Section 2.7; note, this entailed fitting β to the 1-minute temporal mean of the g 2 (τ,t) data used to compute CrCP, i.e., ‹g 2 (τ)›).Fitting for β, though necessary herein, does increase the risk of crosstalk between fitting parameters.Accordingly, fixing β at 0.5 (or at a value measured at baseline) is preferable for DCS signals free of contamination from other light sources.

Appendix 4: white noise in Fourier spectra
Recall, a non-zero "white noise" level, which we denote as σ, is present at all frequencies in the |F(f)|/⟨F⟩ Fourier spectra (e.g., Fig. 3(e), Fig. 5(c)).In general, a white noise level of σ arises from the Fourier transform of a random variable, X(t), with zero mean and σ 2 variance [62].Thus, the white noise in the nonlinear fitting and modified Beer-Lambert Fourier spectra originate from random noise in the F(t) and OD DCS (τ,t) time-series.Figure 3 and 5 indicate that the white noise is larger for OD DCS (τ,t) / (d F (τ)⟨F⟩) (see Eq. ( 7)) than for F(t)/⟨F⟩ obtained from traditional nonlinear diffusion fitting.Future work is needed to better understand this.
Future work is also needed to optimize the removal of the white noise level in the PI F computation.Here, we took the mean of |F(f)|/⟨F⟩ across the 1.2f hr to 1.8f hr frequency range to estimate the white noise level, which we subtracted from |F(f hr )|/⟨F⟩ to determine PIF.The selection of the 1.2f hr to 1.8f hr frequency range was somewhat arbitrary.The range is between the fundamental and first harmonics of the heart rate in the flow waveform and is also beyond higher order harmonics in lower physiologic frequencies in the flow waveform.Further, this frequency range is still plausible for white noise estimation with modestly lower DCS sampling rates than 20 Hz (e.g., 10 Hz).We did try choosing other frequency ranges for estimating the white noise, e.g., taking the average across frequencies defined by 1.2f hr to 1.8f hr and 2.2f hr to 2.8f hr .The choice of the specific frequency range did not alter our findings.

Fig. 1 .
Fig.1.Flow chart of the semi-infinite modified Beer-Lambert algorithm for measuring the cerebral blood flow pulsatility index (PI F ) with diffuse correlation spectroscopy (DCS).See main text for more details.(a) Schematics showing: (left) temporal heart-beat fluctuations in the DCS cerebral blood flow index (F, 10 −8 cm 2 /s) and (right) two exemplar intensity autocorrelation function curves corresponding to the temporal mean (red squares) and the flow at time t * (black circles) (there is no noise in this data).The DCS modified Beer-Lambert law relates temporal changes in OD DCS to corresponding temporal changes in flow and tissue optical absorption (µ a ) and reduced scattering (µ ′ s ) properties (the law is accurate for small OD DCS changes).(b) PI F is obtained by solving a system of equations specified by the Fourier transform of the DCS modified Beer-Lambert law at the heart rate (i.e., one Eq.for each τ; we assume that the changes in tissue optical properties during the cardiac cycle negligibly contribute to OD DCS changes).Inputs in this system of equations include ⟨F⟩ and the d F (τ) weighting factor; ⟨F⟩ is derived from a semi-infinite nonlinear fit of ⟨g 2 (τ)⟩; ⟨µ a ⟩ and ⟨µ ′ s ⟩ were inputs to this fit, and d F (τ) was calculated from ⟨F⟩ and the optical properties via Eq.(6).

𝑔 2 ( 1 Assump�ons•𝐿Fig. 2 .
Fig. 2. Flow chart of the two-layer modified Beer-Lambert algorithm for measuring the cerebral blood flow pulsatility index (PI Fc ≡ |F c (f hr )|/⟨F c ⟩) with diffuse correlation spectroscopy (DCS).The planar two-layer model of the head comprises a semi-infinite cerebral layer (flow index, F c ; tissue optical properties, µ a,c , µ ′ s ,c) and a superficial extracerebral layer with thickness L (flow index, F ec ; tissue optical properties, µ a,ec , µ ′ s , ec ).PI Fc is obtained by solving a system of equations specified by the Fourier transform of the 2-layer DCS modified Beer-Lambert law at the heart rate.We assume: a) small OD DCS changes solely due to F c variation, and b) homogeneous tissue optical properties.Inputs in this system of equations include ⟨F c ⟩ and the d F,c (τ) weighting factor.⟨F c ⟩ and ⟨F ec ⟩ are derived from a two-layer nonlinear fit of ⟨g 2 (τ)⟩.The measured optical properties and L were inputs in this fit; d F,c (τ) was calculated via Eq.(11).

Fig. 3 .
Fig. 3. Experimental Procedures and Exemplar Data Analysis.(a)Burr holes were drilled in piglets for the placement of an external ventricular drain (EVD) and invasive neuromonitors (intracranial pressure (ICP), microdialysis).A non-invasive diffuse optical probe was sutured to the skin of the left forehead (dashed lines indicate the midline and coronal suture).(b) Continuous optical and invasive monitoring were performed during graded ICP increases induced by EVD infusion of mock cerebrospinal fluid.At the end of the time-period for each ICP level, arterial blood gas and microdialysis samples were taken (black circles).(c) Exemplar arterial blood pressure (ABP, mmHg) and optical blood flow index (F, 10 −8 cm 2 /s) time-series in a piglet with invasive ICP of 20 mmHg (F was derived from semi-infinite diffusion fitting to short delay-times (Method 1), the pink circle and green square denote exemplar systolic and diastolic timepoints).(d) (Top) Individual optical intensity autocorrelation function (g 2 (τ)) measurements at the corresponding systolic (pink) and diastolic (green) timepoints labeled in panel (c) (source-detector distance, 2.5 cm; wavelength, 850 nm).(Bottom) Exemplar mean g 2 (τ) across 1 minute of data acquired at the 20 mmHg ICP level with its semi-infinite diffusion fit (solid line).(e) Normalized Fourier spectral amplitudes of 1-minute F time-series derived from semi-infinite diffusion fitting to short delay-times (Method 1, top), and the whole g 2 (τ) curve (Method 2, bottom).White noise levels (see main text) of 0.005 and 0.004, respectively, were subtracted from the heart-rate amplitudes to compute F pulsatility indices (PI F ) of 0.106 and 0.155 for the short τ and whole curve fitting methods.This results in critical closing pressure (CrCP) measurements, respectively, of −11.1 and 2.5 mmHg for the short τ and whole curve fitting methods (via Eq. (3); MAP = 53 mmHg, PI ABP = 0.143).(f) Normalized Fourier spectral amplitudes for blood flow obtained with the short τ semi-infinite (Eq.(7); Method 3, top) and 2-layer (Eq.(10); Method 4, bottom) modified Beer-Lambert algorithms, plotted against frequency.White noise levels of 0.016 and 0.027, respectively, were subtracted from the heart-rate amplitudes to compute PI F , i.e., 0.177-0.016= 0.161 and 0.306-0.027= 0.279.This results in semi-infinite and 2-layer CrCP estimates of 3.6 and 15.5 mmHg, respectively.

2 SelectFit 1 Fig. 4 .
Fig. 4. Schematic of the four methods used to measure the cerebral blood flow pulsatility index once every minute with diffuse correlation spectroscopy.Inputs include the intensity autocorrelation function time-series (g 2 (τ,t), 20 Hz sampling), its temporal mean (⟨g 2 (τ)⟩), the Siegert relation β coefficient, the temporal averages of the tissue optical absorption and reduced scattering coefficients (⟨µ a ⟩ and ⟨µ ′ s ⟩), and the extra-cerebral layer thickness (L).F denotes the cerebral blood flow index in the semi-infinite head model.F c and F ec denote the cerebral and extra-cerebral blood flow indices in the planar 2-layer head model.Flow indices in dark red are derived from short τ fitting/inversion.Flow indices in blue are derived from whole curve fitting.|F(f)|, |F c (f)|, and |OD DCS (τ,f)| are the Fourier spectral amplitudes of F, F c , and OD DCS (τ) at frequency f.For all 4 methods, a white noise (w noise ) correction was made to the pulsatility index measurements (see Fig. 3, Appendix 4).The white noise level is the mean |F(f)| / ⟨F⟩ (or mean |F c (f)| / ⟨F c ⟩ for Method 4) across the 1.2f hr to 1.8f hr frequency range (f hr is the heart rate).

Fig. 5 .
Fig. 5. (a)Simulated temporal blood flow index (F, 10 −8 cm 2 /s) without noise and its normalized Fourier spectral amplitudes (systolic timepoint, pink circle; diastolic timepoint, green square); ⟨F⟩ is the mean F across the 1-minute time-series (only 10 s shown), and PI F is the flow pulsatility index.(b) Simulated intensity autocorrelation function (g 2 (τ)) measurements for the systolic and diastolic timepoints without (top) and with (bottom) noise in the semi-infinite geometry.Noise was added assuming a signal intensity of 30 kHz and averaging across 6 detectors placed at the same source-detector distance of 2.5 cm.(c) For the same noise level as in the bottom panel of (b), the normalized Fourier spectral amplitudes of F obtained with the semi-infinite modified Beer-Lambert (Method 3) and nonlinear diffusion fitting (Method 1) algorithms are plotted against frequency (short delay-times τ, corresponding to ⟨|g 1 (τ)|⟩>0.7,were used for fittings).For computation of PI F , white noise levels of 0.032 and 0.008 for modified Beer-Lambert and diffusion fit, respectively, were subtracted from the signal amplitude at the heart rate (1.45 Hz).(d) Mean PI F (±SD) obtained with the modified Beer-Lambert (red circles, Method 3) and nonlinear diffusion fitting (black squares, Method 1) algorithms on simulated data with varying noise; noise is determined based on the signal intensity and is reduced by averaging across 6 detectors.Notice that PI F approaches its true value for the modified Beer-Lambert fit at signal intensities between 20 and 30 kHz.(e) Same plot as panel (d), except simulations were done using one detector instead of six.(f) For a 30 kHz intensity signal and six detectors, the mean PI F (±SD) obtained for both algorithms are plotted using different ⟨|g 1 (τ)|⟩ fitting cutoffs (for lower cutoffs, more τ are used in the fitting).

Table 1 . Univariate Thresholds, Sensitivity, and Specificity for Detection of Intracranial Hypertension (IH) a
Sensitivity is the fraction of actual IH measurements accurately detected, and specificity is the fraction of actual intracranial normotension measurements accurately detected.CrCP Methods 1−4 are described in Fig.4.
a Optimal decision thresholds (mmHg) maximized the sum of sensitivity and specificity ("> threshold", values greater than threshold indicate IH, "< threshold", values less than threshold indicate IH); IH is defined as ICP ≥ 20 mmHg.

Table 2 . Critical Closing Pressure versus Intracranial Pressure Regressions and Measurement Failure Rates a
Slope and intercept of linear mixed effects best-fit (negative measurements were excluded in the regression analyses); R, Pearson correlation; Failure Rate, percentage of negative critical closing pressure measurements. a