Enhancing diffuse correlation spectroscopy pulsatile cerebral blood flow signal with near-infrared spectroscopy photoplethysmography

Abstract. Significance Combining near-infrared spectroscopy (NIRS) and diffuse correlation spectroscopy (DCS) allows for quantifying cerebral blood volume, flow, and oxygenation changes continuously and non-invasively. As recently shown, the DCS pulsatile cerebral blood flow index (pCBFi) can be used to quantify critical closing pressure (CrCP) and cerebrovascular resistance (CVRi). Aim Although current DCS technology allows for reliable monitoring of the slow hemodynamic changes, resolving pulsatile blood flow at large source–detector separations, which is needed to ensure cerebral sensitivity, is challenging because of its low signal-to-noise ratio (SNR). Cardiac-gated averaging of several arterial pulse cycles is required to obtain a meaningful waveform. Approach Taking advantage of the high SNR of NIRS, we demonstrate a method that uses the NIRS photoplethysmography (NIRS-PPG) pulsatile signal to model DCS pCBFi, reducing the coefficient of variation of the recovered pulsatile waveform (pCBFi-fit) and allowing for an unprecedented temporal resolution (266 Hz) at a large source-detector separation (>3  cm). Results In 10 healthy subjects, we verified the quality of the NIRS-PPG pCBFi-fit during common tasks, showing high fidelity against pCBFi (R2 0.98±0.01). We recovered CrCP and CVRi at 0.25 Hz, >10 times faster than previously achieved with DCS. Conclusions NIRS-PPG improves DCS pCBFi SNR, reducing the number of gate-averaged heartbeats required to recover CrCP and CVRi.


Introduction
Near-infrared spectroscopy (NIRS) is an established non-invasive diffuse optical method used to quantify hemoglobin oxygen saturation (SO 2 ) and hemoglobin concentration (Hb) changes from light attenuation at different wavelengths. 1,2Diffuse correlation spectroscopy (DCS) is a more recent optical technology that uses speckle fluctuations to measure an index for blood flow (BF i ).][5] Combining DCS with NIRS allows for calculating the cerebral metabolic rate of oxygen 6 and furthering our understanding of the relationship between hemoglobin concentration and cerebral blood flow (CBF) changes in healthy 7,8 and pathological conditions. 9,102][13][14][15][16] Despite the encouraging results, the low signal-to-noise ratio (SNR) of current DCS devices limits pCBF i to sourcedetector separations (SDsep) of up to 2.5 cm, which reduces brain sensitivity in adults, 17 and to achieve sufficient time-points within a pulsatile waveform, it requires cardiac-gated averaging of 50 arterial pulses, 11 which dampens the pulsatile peaks and provides CrCP and CVR i estimates at 0.02 to 0.07 Hz, rates that are too low to investigate the dynamic pressure-flow relationship of the cerebral vasculature. 18To overcome the DCS noise, increase SDsep, and recover pCBF i with much less averaging, we propose a new method based on the combination of NIRS and DCS pulsatile signals.
Because NIRS measurements at the same sampling rate typically detect multiple orders of magnitude more photons, the SNR of NIRS is much better than the SNR of DCS, 19,20 allowing for the measurement of pulsatile blood volume waveforms with high temporal resolution at long SDsep (≥3 cm). 21In particular, we recently developed an open-source, wearable, and wireless NIRS device, called FlexNIRS, with a low noise equivalent power (NEP < 70 fW∕ p Hz), able to acquire 10 channels at up to 266 Hz sampling rate. 22The high SNR performance of this device allows us to resolve the pulsatile blood volume and its time derivative to perform cerebral pulse waveform analysis (PWA) on the pulsatile light absorbance of NIRS photoplethysmography (PPG) at a 3.3 cm SDsep (NIRS-PPG) with zero to a few beats averaging.PWA generally refers to the study of the morphology of the PPG waveform measured at short SDsep with pulse oximetry devices. 23Morphological features extracted from the superficial PPG and its time derivative have been investigated in the literature and often include amplitude, latency, and width of the PPG wave contour.These features are generally located with algorithms that find the local maxima and minima in the signal and its first to third time derivatives. 239][30] Further, when simultaneously measuring with DCS and NIRS, by exploiting the pulsatile blood volume and blood flow relationship, [31][32][33] we can separate the pulsatile inflow and outflow and model pCBF i as a linear contribution of NIRS-PPG and its first time derivative [d(NIRS-PPG)/dt].The resulting fitted pCBF i-fit displays an exceeding SNR over DCS, while accurately matching DCS pulsatile flow, allowing us to estimate PI, CrCP, and CVR i at the cardiac frequency.
To validate this model, we measured 12 healthy subjects simultaneously with the FlexNIRS and a state-of-the-art DCS prototype available in our lab, which operates at 1064 nm and employs superconducting nanowire single-photon detectors (SNSPD).The SNSPD-DCS system provides a >16-fold SNR increase versus the standard DCS technology, 17 allowing us to resolve pCBF i at larger separations and to use lower cardiac-gated averaging.We performed the NIRS and DCS measurements on the subjects while performing standard tasks that change cerebral and systemic physiology and recovered both pulsatile and slow varying signals under various conditions.

Participants and Measurement Protocol
For the study, we enrolled and measured 12 healthy subjects.Two subjects were excluded due to low DCS data quality: low DCS photon counts (<20 kcps∕s) and a low coherence factor (Siegert relationship β < 0.3). 34Data from these two subjects were not included in this paper and were not numbered.The remaining 10 subjects included 6 males and 4 females of diverse skin tones, with an age range of 21 to 68 years old (Table S1 in the Supplementary Material).The study protocol was reviewed and approved by the Mass General Brigham Institutional Review Board (IRB No. 2020P002463) in accordance with Ethical Principles and Guidelines for the Protection of Human Subjects (Belmont Report).All subjects signed the informed consent before participating in the study.
The protocol consisted of cold pressor tests (CPTs), paced hyperventilation (HV), and brief repeated breath-holding (BH) tasks.For the CPT, subjects submerged one hand in warm water (30°C to 36°C) for 2 minutes; then they moved the hand into icy cold water (3°C to 5°C) for 1 min; finally, they moved the hand back into the warm water for 2 more min.CPT is known for raising blood pressure, 35 and the protocol barred any subject with high blood pressure from performing this task.For the paced HV, after a 1-minute baseline of sale-paced breathing, subjects breathed at a rate of 70 breaths per minute following a metronome cue for a minute and then resumed their self-paced breathing for two more minutes.For the BH task, after a 1-minute baseline, subjects repeated a set of 20 s BH and 40 s normal breathing for four times, ending with an extended recovery period of 1 min.

Devices
The wearable NIRS system used in the study is a continuous wave, open-source, low-cost, and wearable wireless cerebral oximeter called FlexNIRS, recently developed in our group. 22It consists of a wearable NIRS headband with an embedded flexible optical probe that communicates with a tablet through Bluetooth low energy (BLE) and is powered by a rechargeable lithium-ion battery (Fig. 1).The probe includes two dual-wavelength LEDs (735 and 850 nm, SMT735D/ 850D, Marubeni, Japan) and three PIN-photodiodes (VEMD5060X01, Vishay Intertechnology, Inc., United States), forming two short SDsep channels of 0.8 cm and two pairs of symmetrically arranged long SDsep channels of 2.8 and 3.3 cm.This geometry allows for a self-calibrating multi-distance scheme that uses two values at each separation to eliminate the effect of different component efficiencies and therefore provides a robust measurement of cerebral hemoglobin oxygenation without the need for external calibration. 36,37The use of surface-mount source and detector components in direct contact with the skin minimizes signal attenuation, and the use of a dedicated analog front-end chip that amplifies and digitizes the signals on one integrated circuit results in a device with low NEP (<70 fW∕ p Hz).The optimized custom software allows for a sampling rate of 266.66 Hz for all channels.This was a >2.5 times improvement in the sampling rate over the original description of the FlexNIRS device (100 Hz). 22This was achieved through shortening the sampling durations (LED-on time for all was reduced from 366 to 117 μs) and eliminating repeated ambient and short SDsep channels.
The DCS system used in this study operates at 1064 nm using a long coherence length laser (CrystaLaser, United States) and four high photon detection efficiency (88%) SNSPD detectors (Opus One, Quantum Opus, United States) 17 (Fig. 1).To further increase SDsep while maintaining SNR, we used a dual-source illumination scheme, splitting the laser into two 3.5 mm diameter illumination spots, each emitting up to the maximal permissible exposure limit of 103 mW at 1064 nm (ANSI Z136.1).Source and detector fibers were terminated with prisms into a custom optical probe made of flexible 3D-printing material (NinjaTek, United States) and arranged to form SDsep of 2.0 cm (one detector fiber) and 3.0 cm (three detector fibers) (Fig. 1).In addition, to acquire scalp signals, we used a silicon single-photon avalanche diode (Si SPAD) detector at 0.5 cm from the sources.A custom 20-channel field-programmable gate array board documented the arrival time of each photon at 150 MHz resolution, and a custom data acquisition software calculated and displayed the g 2 and the blood flow index for each channel in real time.
To avoid possible crosstalk between FlexNIRS and the DCS, we positioned the two probes symmetrically on opposite sides of the forehead with the source placed laterally to maximize the distance between them.

Data Analysis
FlexNIRS and DCS data were analyzed using standard routines to recover slow varying hemodynamic changes and with cardiac-gated averaging algorithms to extract pulsatile waveforms.

Slow varying signals
To quantify slow changes in SO 2 and hemoglobin concentration, we used a mixed approach based on the multi-distance self-calibrating method described in Ref. 22 to calculate initial values and a single-distance modified Beer-Lambert law method to calculate hemoglobin changes.After ambient light subtraction, followed by a 1-s moving averaging, light intensity attenuation at 2.8 and 3.3 cm was used to recover absorption coefficients (μ a ), assuming reduced scattering coefficients (μ s 0 ) of 6.8 cm −1 at 735 nm and 5.9 cm −1 at 850 nm. 22,38A 75% water fraction was used to estimate hemoglobin concentration and oxygenation. 39Changes in SO 2 and hemoglobin concentration with respect to the initial values were calculated for each subject by considering each source-detector pair independently.This hybrid method was used to eliminate spurious results of the multi-distance method due to tissue inhomogeneity (seeFig.S1 in the Supplementary Material).
To estimate slow varying blood flow changes, we computed the g 2 curves at 1 Hz smoothed with a moving average of 3 s for each SDsep.The g 2 obtained by three colocalized detectors at 3.0 cm were averaged together to increase SNR.We calculated BF i by fitting the g 2 to the semi-infinite correlation diffusion equation, 40 using optical properties μ s 0 of 4.7 cm −1 and μ a of 0.18 cm −1 .The scattering was extrapolated from the values used for FlexNIRS at shorter wavelengths to 1064 nm, and the absorption was obtained assuming typical hemoglobin and water concentrations. 41n addition, we calculated heart rate (HR) as the inverse of the R-R interval in the ECG signals and extracted the time series of ABP from calibrated Finapres measurements and SpO 2 from the pulse oximeter.For group averaging, optical and systemic measured parameters were down-sampled to 1 Hz, subjected to a 3-s window moving average, and normalized with respect to the initial baseline.The four repetitions of BH were averaged together after detrending and normalizing the data to the three-second periods before the BH.Before averaging across subjects, the missing points that were rejected due to motion artifacts were interpolated.

Pulsatile waveform analysis
The pulsatile waveform of NIRS-PPG was extracted from the 850 nm channels where the pulsatile signal amplitude is larger.A wavelength at the hemoglobin isosbestic point would be ideal to measure pulsatile blood volume.Because the arterial oxygenation is relatively constant, the use of non-isosbestic wavelengths is acceptable.In fact, assuming an arterial oxygenation drop of 10% from 100% to 90%, larger than the arterial oxygenation changes measured in this study, at 850 nm the amplitude of the pulsation changes by <3%.Moreover, in our analysis we only consider the relative amplitude of features of the pulsatile waveform, for which the wavelength dependence is further minimized as arterial oxygenation changes within a heartbeat are negligible.
To extract pulsatile waveforms from the FlexNIRS data, we first subtracted ambient light from the raw intensity measured at 266 Hz and removed data segments with motion artifacts as detected by the embedded accelerometer and gyroscope.We used a high pass filter (f c ¼ 0.4 Hz) to remove breathing and Mayer waves and a low pass filter (f c ¼ 14 Hz) to remove highfrequency noise.We then calculated the delta absorbance, or delta optical density (ΔOD ¼ ln½I 0 ∕IðtÞ, I 0 : baseline light intensity) 42 at 850 nm as NIRS-PPG.In addition to motion artifacts, we rejected distorted heartbeats that had outliers, which were defined as data points deviating from the median waveform by more than a 2.5 interquartile range (IQR).The median and IQR were calculated with cardiac-gated heartbeats within a short moving window of 15 s with a 50% overlap.To establish condition-representing NIRS-PPG waveforms, we gate-averaged heartbeats measured during the various tasks (conditions): each 1-minute baseline was considered to be one condition; the same was done for each recovery segment, excluding the first 20-s transition period; during the HV and CPT, we excluded the initial transitions and considered the last 40 s with more stable hemodynamic responses; and for the BH condition, we considered a period of 20 s centered around the end of the BH and averaged the four repetitions together.The FlexNIRS optical data were aligned with the ECG signals by matching the beat-to-beat HR variability.For the gate averaging, the cardiac cycles were identified using the diastolic onsets in the NIRS-PPG signals, and each heartbeat was normalized by stretching, through interpolation, to a uniform period equivalent to the subject average HR during that condition.We performed these calculations for both NIRS wavelengths and all source-detector pairs; for each subject NIRS-PPG PWA analysis, we report the 3.3 cm pair at 850 nm with the highest SNR (SNR ¼ σ pulsation ∕σ noise ; σ: standard deviation).The SNR was calculated based on high-pass filtered light intensity signals (f c ¼ 0.4 Hz): the pulsation is defined as the low-pass filtered (f c ¼ 14 Hz) intensity, and the noise is the difference between pre and post low-pass filtering.In two subjects (Nos.3 and 8), a 2.8 cm channel was used because the 3.3 cm channels had a low SNR (<15).
To obtain DCS pulsatile blood flow waveforms, we computed g 2 curves at 100 Hz using 40 ms of integration time, i.e., a 20 ms photon-inclusion radius.The data periods excluded in the FlexNIRS due to motion artifacts were also excluded for the DCS time series.Representative DCS pulsatile blood flow waveforms during the various conditions were calculated through cardiac-gated averaging of g 2 curves instead of the CBF i .Each heartbeat onset was determined by the co-acquired ECG R-peaks; across heartbeats, g 2 curves of the same location in the cardiac cycle were averaged together as done in a previous study. 11As for the slow CBF i , the gateaveraged g 2 curves were fitted using the semi-infinite correlation diffusion equation, 40 assuming the optical properties listed previously.During the fitting iterations for the pCBF i , as systolic g 2 curves of high flow lack a clear plateau, β was kept constant.For each measurement, this value was set to be the median β from the traditional free-β fit.
Pulse waveform analysis was performed to quantify the peak amplitudes and their corresponding latencies of NIRS-PPG, d(NIRS-PPG)/dt, and pCBF i .We used a time derivative method 23 modified for brain PPG.When comparing pulsatile NIRS and DCS, to match the NIRS sampling rate, pCBF i was up-sampled to 266 Hz through cubic spline interpolation; to compensate for the longer integration time of DCS (40 ms), NIRS-PPG and d(NIRS-PPG)/dt were subjected to a moving average of 11 points (i.e., 11∕266.66Hz ≈ 40 ms).

Modeling of pCBF i as a Function of NIRS-PPG
Herein we examine the relationship between pCBF i and NIRS-PPG when the two are simultaneously measured.The tissue is assumed to be a semi-infinite homogenous medium, and a similar depth sensitivity is assumed for the two methods.The pulsatile light absorbance of NIRS-PPG is proportional to the pulsatile blood volume changes.As a result, its first time derivative represents the changing rate of blood volume, which is the difference between the blood inflow and the blood outflow during a heartbeat and is given as [31][32][33] E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 4 ; 6 1 6 As suggested by a prior study, if we assume the outflow blood drains out passively, BF out is modeled to be proportional to the pulsatile volume of blood within the area interrogated by diffusely propagating photons in the tissue bed, which is given as 31 ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 4 ; 5 4 7 DCS pCBF i represents the total flow through the vasculature, which gives E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 4 ; 5 1 0 with k being the scaling factor that converts the DCS pCBF i units (cm∕s 2 ) to the blood flow units ½mL∕ð100 g Ã minÞ for perfusion per tissue volume.DCS BF i measures the motion of red blood cells, which is proportional to blood flow, as shown in previous validation studies 43 and recent modeling studies. 44,45y combining these equations, we derive pCBF i-fit , the estimated pCBF i based on the linear combination of the NIRS-PPG signal and its derivative, as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 4 ; 4 1 4 An additional parameter, the alignment lag C 0 was applied to temporally align the signals acquired with independent devices.This is due to the uncertainty of absolute time bases of the devices introduced by the BLE communication and the clocks of the tablet and the computer.A non-linear fit, Levenberg-Marquardt, was performed to minimize the least-squares of residuals and determine the four parameters (C 0 , C 1 , C 2 , and C 3 ).We found that this approach results in many local minima, to which, based on the initial guess, non-linear fitting algorithms are known to converge, as shown in Ref. 46.Hence, we used an initial guess grid search for Levenberg-Marquardt to find the proper initial guess, which leads to the convergence to the global minimum.Finally, by combining Eqs. ( 1), (3), and (4), we were able to separate the inflow (BF in ) and outflow (BF out ) contributions to the pulsatile flow, given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 4 ; 2 3 4 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 4 ; 1 8 5 We calculated C 1 and C 2 for each task condition by fitting the averaged pCBF i and NIRS-PPG signals and then projected these parameters back to NIRS-PPG to estimate pCBF i-fit with no cardiac-gated averaging.C 1 and C 2 were interpolated during task transitions.C 3 was replaced by the slow varying mean CBF i time traces calculated above.We evaluated the SNR improvement with respect to SNSPD-DCS pCBF i , by calculating the coefficient of variation (CV) as σ∕μ (σ: standard deviation and μ: mean) of each point along the cardiac cycle across heartbeats during the 60-s baselines.
2.5 Critical Closing Pressure, Pulsatility Index, and Cerebrovascular Resistance We obtained the PI, CrCP, and CVR i using the pCBF i and the alternative NIRS-PPGreconstructed pCBF i-fit with much less cardiac-gated averaging.With the SNSPD-DCS-only pCBF i , due to its higher noise, we used 20 heartbeats; with NIRS-PPG-reconstructed pCBF i-fit , we used 4 heartbeats.The signals were aligned with the pulsatile ABP (pABP) obtained from the Finapres using the same cardiac-gated averaging.CrCP was obtained by linearly extrapolating the pCBF i versus pABP relationship to the pABP-axis intercept.Because of the nonlinearities during the systolic phase, only the diastolic run-off phase of the pulse cycle was considered, as done previously. 11In the past, the dicrotic notch (DN) was manually found for each subject; however, now due to the high SNR of the pCBF i-fit , we can automatically find the diastolic run-off using a PWA algorithm, 23 which we customized for brain NIRS-PPG.The CVR i was calculated as the inverse of the slope of the pulsatile pressure-flow relationship, which is given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 7 ; 5 6 8 We rejected outliers with loose criteria that only include data points that have pressure-flow regression R 2 of higher than 0.8, CrCP within −30 to 150 mmHg, and positive CVR i .The PI, defined as the difference between the systolic and diastolic flow velocity divided by the mean flow velocity as measured by transcranial Doppler (TCD) ultrasound, is used to assess vascular resistance, 47 intracranial pressure, and cerebral perfusion pressure. 48Herein we define PI of DCS as the difference between the systolic and diastolic pCBF i divided by the mean.The same was calculated for pCBF i-fit .
When comparing results between conditions, after excluding outliers (exceeding 1.5 IQR above Q3 or below Q1), we examined normality with the Kolmogorov-Smirnov test (α ¼ 0.05).If normality was not found, instead of the pairwise t-test, we applied the paired-sample sign test.and cold pressor conditions.We report the pulse waveforms in normalized z-scores, which is defined as the number of standard deviations away from the mean (calculated as z ¼ ðx i − μÞ∕σ.x i : data points, μ: mean, and σ: standard deviation).We observed a strong similarity between pCBF i and d(NIRS-PPG)/dt, with comparable morphology changes of the main peaks with different conditions.Across all subjects and conditions, the R 2 between pCBF i and d(NIRS-PPG)/dt was 0.89 AE 0.10.
The difference between pCBF i and d(NIRS-PPG)/dt can be accounted for by incorporating the NIRS-PPG waveform contribution.Using Eqs. ( 4)-( 6), we were able to reconstruct pCBF i as pCBF i-fit using NIRS-PPG and separate inflow and outflow [Fig.2(b), results for the same representative subject].We found pCBF i-fit to be nearly identical to pCBF i , with R 2 of 0.98 AE 0.01 across all subjects and conditions [Fig.3(a)].The fitting parameters C 1 and C 2 are reported in Fig. S2 in the Supplementary Material, and the ratio C 2 ∕C 1 is shown in Fig. 3(b).The parameters C 1 and C 2 differ between subjects, but weakly depend on different conditions, except for HV, where C 1 , C 2 , and C 2 ∕C 1 dropped (for C 2 ∕C 1 relative to baseline, p < 0.05, paired-sample sign test).
With the higher SNR of the pCBF i-fit , we could identify waveform features automatically.We quantified the systolic peak (P 1 ), secondary systolic peaks (P 2 , tidal wave), DN and diastolic peak amplitudes, and time delay.As an example, the values of P 1 , P 2 , and P 2 ∕P 1 across subjects and conditions are reported in Fig. S3 in the Supplementary Material.We found that P 2 ∕P 1 dropped during HV (relative to recovery, p < 0.01, paired-sample sign test) and increased during CPT (relative to baseline, p < 0.01, paired-sample sign test).

Minimal Cardiac-Gated Averaging Noise Characterization
To quantify the improvement in SNR of pCBF i-fit with respect to SNSPD-DCS pCBF i , we evaluated the CV of the pulsatile waveforms during a 60-s baseline.Figure 4(a) shows the average and standard deviation of pCBF i and pCBF i-fit for each subject.The shapes of the pulsatile waveforms were diverse among subjects, yet the NIRS-PPG fits showed a strong agreement with the DCS-measured pCBF i .pCBF i-fit provided an improvement on the CV, as the CV of pCBF i-fit was consistently lower than the CV of pCBF i , and on average, CV pCBFi ∕CV pCBFi-fit ¼ 3.1 AE 1.7 [Fig.4(b)].The paired-sample sign test rejected the null hypothesis with p < 0.01.

Slow Changes in Cerebral and Systemic Measured Parameters with Tasks
In addition to the NIRS hemoglobin concentration and oxygenation changes and DCS relative changes of the CBF index during the three tasks, we calculated the PI, CrCP, and CVR i changes with time using pCBF i (one point every 20 s) and pCBF i-fit (one point every 4 s).Across subjects, for the two methods, baseline CrCP was 32 AE 13 mm-Hg and 36 AE 14 mm-Hg, respectively, and baseline PI was 1.14 AE 0.35 and 1.37 AE 0.41, respectively.In Fig. 5, we show the slow cerebral responses to each task averaged across all subjects along with the systemic physiology responses, including ΔHb, ΔSO 2 , scalp rBF i , CBF i , PI, CVR i , CrCP, mean ABP (MAP), HR, and SpO 2 .Subject No. 2 was excluded because of poor attachment of the FlexNIRS probe, which caused large light intensity jumps during task transitions.Subject No. 5 was excluded due to technical problems with the Finapres pABP measurement, which prevented us from calculating CrCP and CVR i .For the cold pressor task, subject No. 8 did not perform the task due to high blood pressure.

Hyperventilation
HV causes hypocapnia due to over-breathing and is expected to cause vasoconstriction and a decrease in blood flow. 49,50Because cerebral metabolism is maintained during this period, the reduction in blood flow results in a decrease in SO 2 , which induces a vasoactive reaction to return the blood flow to the baseline level. 51This biphasic response was observed in HbT, SO 2 , and CBF i , whereas scalp blood flow (0.5 cm rBF i ) showed minimal changes.The pulse amplitude of both NIRS-PPG and its time derivative increased steadily.A similar but smaller biphasic response was observed in CVR i (opposite direction, i.e., initial increase and subsequent decrease), whereas PI and CrCP increased during HV and slowly returned to baseline during recovery.The observed systemic physiologic response to the HV was a significant increase in HR (13.5 AE 8.6 BPM) and an initial increase (8 AE 4 mmHg) and subsequent decrease (−8 AE 11 mmHg) in MAP, as seen in Fig. 5(a).

Breath-holding
A series of repeated brief breath holds caused an increase in blood pressure 52,53 and a hypercapnic state. 54,55As shown in Fig. 5(b), the average response to a 20s BH task show increases in HbT (0.4 AE 0.8 μM), SO 2 (0.5% AE 0.8%), CBF i (14% AE 11%), and scalp BF i (23% AE 28%).The pulse amplitude of NIRS-PPG and its time derivative dropped.From the pulsatile parameters analysis, we observed a decrease in CVR i and PI during BH.MAP increased (9 AE 5 mmHg), whereas the HR remained relatively constant.SpO 2 measured on the finger showed a drop (to ≤95%) in 4 out of 8 subjects.This drop was observed to lag behind the head and brain responses, with a nadir of about 22 s after the resumption of breathing.This delay of finger responses was expected and is consistent with literature findings. 56

Cold pressor test
The CPT induces both pain and a response of the autonomous nervous system, including an elevation in blood pressure. 35In our experiments, MAP increased substantially (23 AE 8 mmHg), whereas HR and SpO 2 remained relatively constant.After an initial small drop in HbT, SO 2 , and CBF i , the parameters increased throughout the remainder of the task.CBF i increased by 41% AE 17%, whereas scalp BF i reached a 143% AE 90% increase.The large increase in CBF i and blood pressure caused a small increase in CrCP toward the end of the stimulus and a decrease in CVR i .PI showed a biphasic response: a decrease during CPT and a step increase when the hand was moved back into the warm water.The step increase was also observed in the pulse amplitude of both NIRS-PPG and its time derivative.

Discussion
We demonstrated that the simultaneous use of NIRS and DCS allows for the recovery of pulsatile blood flow waveforms with a higher SNR and a higher temporal resolution than with DCS alone, even when using a high-end DCS system such as the SNSPD-DCS at 1064 nm here.The use of a low noise, high sampling rate NIRS device, such as FlexNIRS (266 Hz), allowed for a three times improvement in pCBF i-fit CV with respect to SNSPD DCS pCBFi (Fig. 4) and an increase in the sampling rate from 100 to 266 Hz at >3 cm SD separation.Moreover, the low noise of the recovered pCBF i-fit allowed us to use automatic pulsatile feature extraction algorithms for PWA.This method was demonstrated with SNSPD DCS in 10 healthy subjects, showing that a three-parameter linear fitting model using the NIRS-PPG and its time derivative signals is sufficient to recover pCBFi with high accuracy (R 2 ¼ 0.98) (Fig. 2).We showed that the fitting parameters are quite stable within the same subject, despite induced physiological changes, except for HV (Fig. 3).The changes in the fitting parameter C 2 during HV imply a change in the vasculature capacitance and resistance.Ideally, we should update the fitted parameters continuously during a monitoring session.This may be a potential issue when using a common DCS system, for which a much larger cardiac gated averaging is needed to provide a stable pCBF i waveform, especially at 3 cm source-detector separations.
Because of the similarity of pCBF i and d(NIRS-PPG)/dt (R 2 ¼ 0.89), when performing PWA, if necessary, d(NIRS-PPG)/dt alone can be used to approximate pCBF i .Instead, to calculate CrCP, CVR i , and PI, the pCBF i mean amplitude is needed to calibrate d(NIRS-PPG)/dt, so a combined NIRS-DCS approach is necessary.The use of pCBF i-fit allows for calculating CrCP, CVR i , and PI with minimal averaging and detecting changes during transients.The 4 s average was done to reduce pulse-to-pulse variability and to obtain a more stable fitting when estimating CrCP.In the literature, CVR i is often calculated using the mean CBF and the mean blood pressure values, which requires the assumption of a zero intercept on the pressure-flow relationship.This is not always the case as shown by the high CrCP values (the pressure-intercept) and changes with tasks.It is interesting to see that, although the PI is often used as a measure of cerebrovascular resistance, CVR i and PI do not always match during task conditions, likely because the PI is also affected by the compliance of the cerebral arterial bed, ABP, and HR, 57 which differentially vary during the selected tasks.
As shown in Fig. 2, this method allows for the separation of inflow and outflow contributions to the pulsatile flow.The separation of the individual contributions may be useful when trying to estimate intracranial pressure using pulsatile features because increases in ICP would likely affect the lower-pressure venous outflow more than, and in a different way than, the higherpressure arterial inflow.
With respect to the slow hemodynamic signal's changes measured during the various tasks, our findings are in agreement with what was expected based on literature reports.Hemoglobin concentration changes during HV 58,59 and BH 60 are consistent with previous NIRS findings.The CBF i response during HV and BH is consistent with blood flow results previously reported by our group. 172][63] During the cold pressor tasks, the measured increases in CBF, hemoglobin oxygenation, and total hemoglobin concentration are due to systemic blood pressure increases above cerebral autoregulation thresholds and, in the frontal area, are also due to functional activity caused by pain stimulation. 647][68] Similarly, during cold pressor tasks, decreases in PI have been observed with TCD, 69,70 and the increases that we observed in CVR i are consistent with TCD-derived resistance area product (RAP) increases reported by Panerai et al. 71 This study has several limitations.Because of optical crosstalk, we could not measure NIRS and DCS in the same location, but instead measured them contralaterally across the forehead.This may account for the small differences that we observed between the pCBF i and pCBF i-fit waveforms.This problem can be resolved in the future by unifying the NIRS and DCS systems using a single probe and temporal multiplexing by pulsing the DCS laser source.
We assumed a homogenous medium.To further separate the scalp and brain signals, the addition of multi-layer models for pulsatile NIRS and DCS is required.
The NIRS-PPG fitting procedure is affected by having several local minima, especially for the NIRS-PPG contribution term C 2 .This is because the NIRS-PPG waveform has fewer distinct features and contributes less than the first derivative.We overcame this problem by grid-searching for the initial guess.
To demonstrate the method, we used a high-end DCS device with a more than 16 times higher SNR than conventional DCS.When using conventional DCS at traditional NIRS wavelengths, the recovered pCBF i will be much noisier.With current systems, we can use larger pulse averaging and shorter source-detector separations to reduce the errors in determining the C 1 and C 2 fitting parameters.Novel DCS devices are emerging with improved SNR, comparable to the one of the SNSPD-DCS system used here, for which shorter averaging should be sufficient.
Cardiac-gated averaging introduces error, especially when the HR varies.Using the combined NIRS-DCS approach, we were able to substantially reduce the number of cardiacgated averaging needed to obtain clean pulsatile blood flow signals.
In conclusion, the use of NIRS-PPG and d(NIRS-PPG)/dt allows for achieving pulsatile blood flow with high temporal resolution and minimal cardiac-gated averaging.Combined with DCS, to obtain calibration, it allows us to derive CrCP and CVR i at a rate of 4 s.The combination of FlexNIRS and SNSPD-DCS at 1064 nm allows for the comprehensive probing of hemodynamic response during breathing and other tasks.

Disclosures
At the time of the study, MAF had a financial interest in 149 Medical, Inc., a company developing DCS technology for assessing and monitoring CBF in newborn infants, which is now dissolved.MAF's interests were reviewed and managed by Massachusetts General Hospital and Mass General Brigham in accordance with their conflict-of-interest policies.

Code, Data, and Materials
De-identified data will be provided upon request to the authors following MGB regulations.
FlexNIRS is an open-source NIRS cerebral oximeter that can be used for research purposes.Information can be found on our website.

Fig. 1
Fig. 1 Measurement setup.Right: FlexNIRS for hemoglobin concentration and SO 2 .A piece of white electrical tape was applied to the short SDsep PD to attenuate the light.Left: SNSPD-based DCS for longer separations (2.0 and 3.0 cm) and Si SPAD-based DCS for the short separation (0.5 cm).Middle: pictures of the FlexNIRS and DCS probes.Bottom: the physiological signals from BIOPAC and Finapres recorded with PowerLab.
Fig. 2 (a) Average pulsatile waveforms of a representative subject (No. 5) during the measured conditions.The features for PWA are listed on the baseline panel.(b) pCBF i and pCBF i-fit comparison and inflow and outflow contribution (in BF i unit cm 2 ∕s) to the pulsatile flow on the same subject.

Fig. 5
Fig.5Group average of all measured parameters during (a) hyperventilation task (n ¼ 8, HV period in yellow), (b) BH task (n ¼ 8, block average of four repetitions, BH period in green), and (c) cold pressor test (n ¼ 7, CPT period in cyan).r-prefix indicates relative to, and Δ-prefix indicates delta changes; both are with respect to the average of 60-s initial baseline.The standard deviation across subjects is shown in lighter colors.rPI, ΔCrCP, and rCVR i were calculated every 4 s using pCBF i-fit (continuous trace) and every 20 s (circles) using pCBF i .All other parameters were calculated at 1 Hz with a 3 s moving average.For ΔHb graphs: red ΔHbO, blue ΔHbR, and black ΔHbT.The pulse amplitude of NIRS-PPG and its derivative are shown in Fig.S4in the Supplementary Material.