2.1.

Hypercapnia Study in Healthy Volunteers

We conducted a human subject study with healthy volunteers to measure cerebral responses to induce hypercapnia. As part of a larger hypercapnia challenge study that had several goals beyond the scope of this paper, we enrolled a total of 27 subjects, each of which had either one or two visits, i.e. sessions. The first visit consisted of either optical-only or joint optical-TCD measurements, and the second consisted of joint optical-MRI scanning measurements; in total, 43 measurement sessions were held. Seventeen of the subjects had an MRI scanning session and 9 of the subjects had joint optical-TCD measurements, where TCD data were acquired concurrently with the optical recordings. A pressure modulation maneuver was also performed on each subject before the hypercapnia session. The study was approved by the Mass General Brigham Institutional Review Board, and all subjects provided informed consent prior to the measurements. A table of the subject information is shown in Table 1.

Table 1

Subject information. Values for age, baseline PetCO2, and increase in PetCO2 are given as mean±standard deviation.

2.2.

${\mathrm{P}}_{\mathrm{ET}}{\mathrm{CO}}_2$ Processing Methods

Raw data for ${\mathrm{P}}_{\mathrm{ET}}{\mathrm{CO}}_2$ consist of an exported ${\mathrm{CO}}_2$ partial pressure timecourse from LabChart at a sampling rate of 1 kHz. To obtain the increase in ${\mathrm{P}}_{\mathrm{ET}}{\mathrm{CO}}_2$ throughout the 5-min breathing protocol and filter out respiratory oscillations, we obtained the upper envelope of each measurement and subtracted the average value of this envelope during the baseline from the entire timecourse.

2.3.

TCD Processing Methods

Raw data from the TCD instrument consist of two timecourses of MCAV in cm/s, one each for the LMCA and RMCA at a sampling rate of 1 kHz. Heart rate and/or respiratory oscillations must be filtered out to obtain a clean MCAV timecourse usable for comparison to DCS data. We first obtained the upper and lower envelopes of each timecourse. Then, a lowpass filter of 0.01 Hz was used on all envelopes; measurements were normalized by dividing the entire timecourse by the average baseline MCAV and converting all values to percent change from baseline. Finally, the upper and lower envelopes from both (if available) L/RMCA were averaged together to produce a single TCD timecourse for each subject.

For subjects who did not have TCD measurements, we computed the expected change in MCAV by calibrating a group average of our existing TCD measurements with a group average of their corresponding ${\mathrm{P}}_{\mathrm{ET}}{\mathrm{CO}}_2$ timecourses. The travel time of the exhaled gas through the sampling tube results in a slight time lag of the ${\mathrm{P}}_{\mathrm{ET}}{\mathrm{CO}}_2$ timecourse with respect to the TCD curve. This was compounded by a physiological delay in reacting to the hypercapnia. We correct for this by calculating the $r^2$ statistic for linear regression between the TCD timecourse and multiple time-shifted ${\mathrm{P}}_{\mathrm{ET}}{\mathrm{CO}}_2$ timecourses whose delay time ranged from 1 to 10 s in 1-s increments. The offset time that resulted in the highest $r^2$ statistic between the two timecourses was chosen, and the corresponding linear regression model was then used to calibrate the TCD and shifted ${\mathrm{P}}_{\mathrm{ET}}{\mathrm{CO}}_2$ data for the 60- to 250-s time interval of the measurement. All subjects chosen with MC fitting based on systemic physiological drift and signal quality criteria (as described below in Sec. 2.4.2) did not have TCD recordings, so expected MCAV values were calculated from the ${\mathrm{P}}_{\mathrm{ET}}{\mathrm{CO}}_2$ timecourses.

2.4.

DCS Processing Methods

The DCS measurements were processed in two primary ways: first, fitting the data using the analytical solution to the semi-infinite diffusion correlation equation, which employs a single-layer, homogeneous model; and second, fitting against a multi-layered MC model of DCS measurements generated by a simulation of photon propagation through tissue. For both methods, experimentally acquired autocorrelation curves from fibers at equal s–d distances were averaged together to make a single autocorrelation at the corresponding s–d separation. Hypercapnia measurements were processed with a 10-s integration time for each data point (autocorrelation curve), and pressure modulation measurements were processed with a 3-s integration time for each data point.

2.4.3.

Monte Carlo-based fitting

MC forward simulations of light propagation through tissue were performed using the MCX software package developed by Fang et al.⁴⁴ One billion photons were launched per simulation, and relevant data saved from the software simulation included the path length and momentum transfer for each photon for each tissue layer.

Two primary volumetric geometries were used for the forward model, both with $1\mathrm{mm}\times 1\mathrm{mm}\times 1\mathrm{mm}$ resolution. The first volume was a 21-layer slab volume, with each layer 1-mm thick. After running the forward simulation, the photon path length and momentum transfer data from 21 layers were then concatenated during post-processing into either two, three, or four total layers, depending on the model tested. Two-layer models represented the extracerebral layers and the brain. Three-layer models replicated the scalp, skull, and the brain; four-layer models denoted the scalp, skull, cerebrospinal fluid (CSF), and the brain. Thicknesses of each layer were adjustable as needed for processing a given measurement. The second volume used for fitting encoded a realistic, head-like model—hereon referred to as the curved model—initially taken directly from a subject’s MRI head scan,⁴⁵ and further iteratively image-eroded to create a final 22-layer head volume, with each layer 1-mm thick (covering the full range of physiologically realistic cortical depths). As with the multi-layer slab volume, the 22 layers were concatenated into two, three, or four total layers as needed. A diagram showing how the slab and curved head models were generated is displayed in Fig. 1.

Fig. 1

Top row: the 21-layer slab is concatenated into a 2-, 3-, or 4-layer volume in post-processing (3 layer is shown). Bottom row: the segmented MRI is iteratively image-eroded to make a 22-layer head volume and is analogously post-processed into a 2-, 3-, or 4-layer volume.

Measurements were fitted using both the slab and curved models. Two-layer volumes were concatenated such that the extracerebral layer ranged from 1 to 18 mm in 1-mm increments, totaling 18 different MC model versions. Seventy different model versions were used for the three-layer volumes, with scalp thicknesses ranging from 1 to 7 mm and skull thicknesses ranging from 4 to 13 mm in 1-mm increments each; the 70 different models consist of all possible combinations of scalp and skull thicknesses throughout those ranges. Last, thickness combinations for four-layer models were taken from optimal thickness combinations chosen from the three-layer model (selection process is described in depth in Sec. 3.3), with a 1- or 2-mm-thick region removed from the skull layer to create a CSF layer.

For the vast majority of the MC models fittings, scattering was set to $0.85{\mathrm{mm}}^{-1}$ for each layer. As an exploration, the impact of using distinct scattering coefficients for each layer was also tested—values used for the scalp, skull, and brain (used with our 850-nm source and taken from Gagnon et al.¹⁸ as reduced scattering coefficients) were 0.74, 0.81, and $1.16{\mathrm{mm}}^{-1}$, respectively. For absorption coefficients, multiple literature reports using our wavelength range of interest for each of the four layers were averaged,^33,46–48 and the values 0.01, 0.033, 0.004, and $0.023{\mathrm{mm}}^{-1}$ were used for all MC models. Of note, MC simulation photon history files can simply be reprocessed to take into account changes in absorption, but the MC simulation needs to be re-run to account for changes in scattering as shown below.

Using the photon path lengths and momentum transfer obtained by the MC models, the temporal field autocorrelation function for each tissue layer was calculated as:^30,32

Eq. (5)

$$G_1(\tau )=\frac{1}{N_p}\sum _{n=1}^{N_p}\exp (-\frac{1}{3}{k}_0^2\sum _{i=1}^{N_t}{Y}_{n,i}{⟨\mathrm{\Delta }{r}^2(\tau )⟩}_i)\exp (-\sum _{i=1}^{N_t}{\mu }_{a,i}{L}_{n,i}),$$

where $k_0$ is the wavenumber of the light in the medium, $N_p$ is the number of photons detected, $N_t$ is the number of tissue layers, $Y_{n,i}$ is the total momentum transfer of photon $n$ in layer $i$, $L_{n,i}$ is the total path length of photon $n$ in layer $i$, ${\mu }_{a,i}$ is the absorption coefficient in layer $i$, and ${⟨\mathrm{\Delta }{r}^2(\tau )⟩}_i$ is the mean square displacement of the scattering particles in layer $i$, as defined in Sec. 2.4.1. Likewise, the normalized temporal autocorrelation function, $g_2(\tau )$, remains the same as in Eq. (4). However, SBF (${\mathrm{BF}}_{\mathrm{i}}$ in the topmost layer) was held to the value calculated from the 5-mm detector analytical fit for a given timepoint, and skull blood flow was held to $2\times {10}^{-8}{\mathrm{mm}}^2/\mathrm{s}$ (assumed to be 1% of typical head measurement values). As explained in Sec. 2.4.1, $\beta$ was fitted for at each timepoint, but a second pass fit was done for hypercapnia measurements in which $\beta$ was held constant for hypercapnia timecourses to the median values for the first pass at each detector.

2.5.

Estimating Probe-to-Brain Distance

As part of the larger study, some DCS measurements during hypercapnia were acquired in an MRI scanner, where both structural and concurrent functional scans were taken on the subject as validation. While we are not reporting the functional MRI data in this paper, we use the structural MRI scan to estimate extracerebral thickness and compare it to the value predicted by the optimized model in our MC processing.

Multi-echo MPRAGE T1-weighted scans were obtained on the subjects with vitamin E tabs placed on top of the DCS probe as fiducial markers to indicate the optical probe locations. From here, two methods were used to measure extracerebral thickness. The first consisted of loading the MRI DICOMs into the MicroDICOM software package and using the built-in distance tool as a ruler. In the second method, the structural DICOMs were processed into a volumetric segmentation using the MRI Freesurfer image analysis suite.⁴⁹ These volumes were then segmented into four layers—the scalp, skull, CSF, and brain—using methodology presented by Perdue and Diamond.⁵⁰ Finally, the four-layer segmented volume is converted into a mesh-based volume,⁵¹ with the fiducial marker location preserved throughout the entire process.

The node closest to the fiducial marker location was chosen on the surface of the head, and a normal vector was calculated by averaging the surface normals of all faces containing the node. From there, a line was drawn inward toward the brain along the averaged normal and was incrementally increased in length until the gray matter tissue is reached. This was done for a patch of nodes on the surface surrounding the initial chosen node, up to approximately a 1.5-mm radius, to create a final distribution of calculated probe-to-brain distances.

3.1.

${\mathrm{P}}_{\mathrm{ET}}{\mathrm{CO}}_2$ and TCD

Figure 2(a) shows the ${\mathrm{P}}_{\mathrm{ET}}{\mathrm{CO}}_2$ group average of the 9 subjects totaling 18 measurements for whom TCD measurements taken concurrently alongside optical acquisitions. Ten measurements had both left and right MCAV timecourses; four had only the left and another four had only the right. Taking the upper and lower envelopes of the total available timecourses resulted in 56 timecourses used for the group average. To remove spikes in the individual timecourses that were otherwise valid, we excluded the highest and lowest data points on a timepoint-by-timepoint basis. An increase of $\sim 11.5\pm 0.6\mathrm{mmHg}$ ($\text{mean}\pm \text{standard error}$), a value slightly above our protocol target of $10\pm 2\mathrm{mmHg}$, is reached before the sharp decrease starting at 189 s. Accordingly, Fig. 2(b) shows the group average (outliers trimmed) for the TCD measurements corresponding to the ${\mathrm{P}}_{\mathrm{ET}}{\mathrm{CO}}_2$ timecourses in Fig. 2(a). The $\mathrm{\Delta }\mathrm{MCAV}$ reaches approximately 39% $\pm 2\%$ above baseline before a sharp drop starting at about 194 s.

Fig. 2

Group averaged ${\mathrm{P}}_{\mathrm{ET}}{\mathrm{CO}}_2$ timecourse for all subjects with concurrent TCD measurements, with standard errors in light blue and the quartiles plotted in dotted lines (a), and the corresponding averaged TCD timecourse for the same subjects (b).

We regressed the individual TCD responses from their respective changes in ${\mathrm{P}}_{\mathrm{ET}}{\mathrm{CO}}_2$ for each of the nine subjects. The average slope across the nine subjects was $3.2\pm 0.9$ (given as mean and standard deviation); the average intercept was $-0.4\pm 4.0$; and the average $r^2$ value was $0.92\pm 0.07$. Figure 3 shows an example plot of the linear relationship between ${\mathrm{P}}_{\mathrm{ET}}{\mathrm{CO}}_2$ and MCAV in one subject and shows the relationship between the curves in Fig. 2 for the interval from 60 to 250 s as described in Sec. 2.3. The average ${\mathrm{P}}_{\mathrm{ET}}{\mathrm{CO}}_2$ timecourse in this particular case was offset with a time lag of 5 s. Estimation of the TCD timecourses from the ${\mathrm{etCO}}_2$ data for the remainder of this study was fitted with the following model: $\mathrm{\Delta }\mathrm{eTCD}=3.2·\mathrm{\Delta }{\mathrm{P}}_{\mathrm{ET}}{\mathrm{CO}}_2-0.4$, where eTCD is the estimated TCD MCAV derived from the ${\mathrm{P}}_{\mathrm{ET}}{\mathrm{CO}}_2$ change.

Fig. 3

Example $\mathrm{\Delta }\mathrm{TCD}$ timecourse plotted against its corresponding $\mathrm{\Delta }{\mathrm{P}}_{\mathrm{ET}}{\mathrm{CO}}_2$ timecourse. The linear regression fit is overlaid in blue.

3.2.

Observed Group-Averaged ${\mathrm{BF}}_{\mathrm{i}}$ Changes due to Hypercapnia Using Semi-Infinite Analytical Model

Figures 4(a)–4(d) show group averages of the 5-, 25-, and 30-mm s–d separation for measurements during the hypercapnia challenge from subjects for whom data were acquired at all three distances and the baseline ${\mathrm{BF}}_{\mathrm{i}}$ (consisting of the 1-min prior to ${\mathrm{CO}}_2$ administration) had a coefficient of variation (CoV) of less than 25%. Thirteen out of 59 measurements were excluded due to the CoV criteria. Outliers were excluded at the 5% level on a timepoint-by-timepoint basis again to remove spikes in individual timecourses that were otherwise clean. This resulted in two data points, the highest and the lowest, removed per timepoint for the 5-, 25-, and 30-mm distances, and no data points for the 20-mm distance. Figure 4(d) shows a group average of the subset of these measurements that contained a 20-mm separation probe attached to the subject calf (total of 28 measurements across 10 subjects). We define $\mathrm{\Delta }{\mathrm{rBF}}_i(\%)$ as $({\mathrm{rBF}}_i-1)\times 100$, where ${\mathrm{rBF}}_i$ is the ${\mathrm{BF}}_i$ timecourse divided by the average value at baseline. The timecourses peak at 220, 190, 180, and 270 s, respectively, reaching increases of $20\%\pm 3\%$, $18\%\pm 2\%$, $20\%\pm 2\%$, and $15\%\pm 4\%$, at 5-, 25-, 30-mm forehead, and 20-mm peripheral s–d separations, respectively. Recovery at 300 s reaches $11\%\pm 2\%$, $8\%\pm 1\%$, $8\%\pm 1\%$, and $10\%\pm 2\%$ of baseline.

Fig. 4

Group average of (a) 5-mm detector, (b) 25-mm detector, and (c) 30 mm detector, totaling 46 measurements. (d) A 28-measurement subset of those that contained a 20-mm probe on the subject calf. The error bars represent standard error, and the dotted timecourses above and below the group average represent the 25th and 75th quartiles, respectively.

We consider significant systemic drift to be an observed ${\mathrm{BF}}_{\mathrm{i}}$ increase in the 5-mm s–d channel that reached at least 20% during or after hypercapnia. Among the sessions containing all three s–d distances shown in Fig. 4 (where each session for a subject consisted of two hypercapnia measurements), 63% of sessions contained at least one measurement with significant systemic drift. We note that across the larger hypercapnia challenge study, which included subjects with only 5- and 30-mm s–d distances used, a 61% ratio was observed for the same phenomenon.

We extracted the peak times for each individual timecourse used in Fig. 4 and plotted the spread of the peak times for each individual detector distance (see the Supplementary Material). The spreads showed a similar result to the group averages, in that the median of the peak times was earliest for 30 mm, followed by 25, 5, and 20 mm.

3.3.

Subject 1: Exploration of Various Multi-Layer Monte Carlo DCS Models to Reduce Systemic Physiology Cross-Talk in a Sample Case

Based on the criteria outlined in Sec. 2.4.2, seven total measurements, each from a different subject, were found containing both significant systemic (scalp) ${\mathrm{BF}}_{\mathrm{i}}$ change and an observed, distinct long separation increase during hypercapnia that was discernable from the 5-mm response. Four of those measurements had long separation responses returning to within 10% of baseline at the end of the measurement, indicating limited impact of systemic physiology cross-talk. Thus, we performed multi-layer, MC-based fitting on a total of other three measurements.

Since no concurrent TCD data were collected for these measurements, hypercapnia timecourses for all three measurements are shown up to the 250-s mark, as the expected TCD timecourse synthesized from ${\mathrm{P}}_{\mathrm{ET}}{\mathrm{CO}}_2$ change using the calibration curve discussed in Sec. 3.1 does not include data past 250 s. In sections below, we first explore several different MC modeling approaches on the first subject measurement: two-layer, three-layer with uniform scattering across layers, three-layer with variable scattering, and four-layer using a range of diffusion coefficients for the CSF. Finally, based on the results from the first subject using the various approaches listed above, we report on the use of the most successful MC model approach for the other two subject measurements presented.

3.3.1.

Subject 1: reference analytical fit

Figure 5 shows the analytical fitting-based timecourse for both the pressure modulation and hypercapnia measurements on subject 1 with expected TCD changes (eTCD) overlaid.

Fig. 5

Analytical fitting-based timecourses of ΔrBFi of subject 1 during (a) pressure modulation with light pressure applied approximately between 20 and 50 s from the start, and (b) hypercapnia in terms of the percentage change in the baseline normalized $\mathrm{\Delta }{\mathrm{rBF}}_{\mathrm{i}}$. Shading indicates the pressure period for (a) and the ${\mathrm{CO}}_2$ enriched breathing gas delivery period for (b).

In the pressure modulation timecourse, we observe a 30% to 40% reduction in ${\mathrm{BF}}_{\mathrm{i}}$ in the long s–d separations and 85% to 90% decrease in the short separation measurement during the period where light pressure is applied on the probe to reduce SBF. For the hypercapnia measurement, we observe a continuous increase in ${\mathrm{BF}}_{\mathrm{i}}$ measured from the 5-mm detector, reaching 33% increase from baseline at 250 s. We note a 23% and 31% increase in the 25- and 30-mm detector timecourses, respectively, both peaking at approximately 190 s; they additionally remain elevated at slightly above 10% at 250 s. The expected TCD timecourse, estimated from the ${\mathrm{P}}_{\mathrm{ET}}{\mathrm{CO}}_2$ data, is overlaid; a 32% peak increase at approximately 180 s is reached before a return to baseline around 220 s.

3.3.2.

Subject 1: Monte Carlo two-layer (MC-2L) models

As MC methods fit for ${\mathrm{BF}}_{\mathrm{i}}$ values for the deep tissue layer as opposed to the analytical method, which fits ${\mathrm{BF}}_{\mathrm{i}}$ per detector (described in Sec. 2.4.3), we use the term ${\mathrm{CBF}}_{\mathrm{i}}$ for MC-based blood flow approximations. While minimal differences in ${\mathrm{CBF}}_{\mathrm{i}}$ timecourse were observed in using just the 5- and 30-mm separation data versus the 5-, 25-, and 30- mm separation data in the MC-based processing algorithm, the pressure modulation timecourse for subject 1 was less noisy with the former option. We therefore show MC-based fits for this subject that were processed using input from just the 5- and 30-mm separations, whereas the MC fitting for the other subjects used inputs from all three s–d distances. Thus, we proceed with MC-based two-layer (MC-2L) slab and curved models to fit the subject 1 hypercapnia timecourse, as shown in Fig. 6. Models with extracerebral thicknesses up to 18 mm were used to fit, but only a subset of them is shown. While a complete return to baseline by 250 s is reached at an unrealistic 1 mm (or 2 mm, for the curved model) extracerebral thickness, the overall response profile does not align well with eTCD. Increasing the extracerebral thickness only pulls the ${\mathrm{CBF}}_{\mathrm{i}}$ down to further unrealistic values during the recovery session (shown up to 6 mm, ${\mathrm{CBF}}_{\mathrm{i}}$ changes for thicknesses beyond 6 mm become increasingly disproportionate to what is expected). Response differences in the same extracerebral thicknesses between the slab and the curved model are observed. For instance, a 5-mm extracerebral thickness in the slab model decreases to about 80% under baseline by 250 s, whereas the same extracerebral thickness in the curved model only decreases to about 15% under baseline at the same time. In general, for the same extracerebral thickness, the curved model timecourse produces less relative drop in ${\mathrm{CBF}}_{\mathrm{i}}$ during recovery than the slab model.

Fig. 6

Subject 1 hypercapnia timecourse fitted with (a) MC-2L slab model and (b) curved model overlaid with the analytical fit at 25- and 30-mm separations and eTCD. Only fitting results using MC models with up to 6-mm extracerebral thickness are shown.

3.3.3.

Subject 1: Monte Carlo three-layer (MC-3L) models with uniform scattering across layers

Monte Carlo three-layer (MC-3L) fitting results in many more models to select from. To choose the optimal combinations of scalp and skull thickness, we refer to both the $\mathrm{\Delta }{\mathrm{rCBF}}_{\mathrm{i}}$ timecourses for the pressure modulation measurement and the absolute ${\mathrm{CBF}}_{\mathrm{i}}$ timecourses for the hypercapnia measurement. We narrow our selection range to models containing scalp and skull thickness combinations for which ${\mathrm{CBF}}_{\mathrm{i}}$ stayed the most constant before, during, and after the pressure modulation period, and also to models for which the brain-to-scalp ${\mathrm{BF}}_{\mathrm{i}}$ ratio at baseline was at physiologically relevant levels, between 3 and 8. The final model is chosen where these two criteria converge. Figure 7 shows MC-3L ${\mathrm{CBF}}_{\mathrm{i}}$ timecourses for both the pressure modulation and hypercapnia maneuvers in subject 1. For the slab model, a 1-mm scalp and a 7-mm skull thickness were used; for the curved model, a 2-mm scalp and 10-mm skull thickness were chosen. The brain-to-scalp ${\mathrm{BF}}_{\mathrm{i}}$ ratio for each model was 3.61 and 4.42, respectively.

Fig. 7

Subject 1 timecourses fitted with MC-3L slab and curved models overlaid with the analytical fit at 5-, 25-, and 30-mm separations and eTCD during (a) pressure modulation and (b) hypercapnia measurements.

In both the MC-3L slab and curved models, ${\mathrm{CBF}}_{\mathrm{i}}$ reduction during the pressure modulation reaches a maximum of 16%, significantly less than the reduction observed in the analytical long separation fit. In the hypercapnia protocol, we observe that the peak ${\mathrm{BF}}_{\mathrm{i}}$ reached in the MC-3L timecourse occurs 10 s earlier compared to the analytical fit, the MC-3L peak aligning closer to what is shown in the expected TCD timecourse. Last, similar to the eTCD timecourse, the MC-3L timecourse returns to within a few percent of baseline by 250 s, in contrast to the analytical fit, which remains at 10% to 12% above baseline.

3.3.4.

Subject 1: Monte Carlo three-layer models with variable scattering in each layer

Figure 8 displays MC-3L fits with scattering values adjusted per layer instead of using a uniform scattering value. As mentioned in Sec. 2.4.3, values for the scalp, skull, and brain used for our 850-nm source were 0.74, 0.81, and $1.16{\mathrm{mm}}^{-1}$, respectively. The same criteria for scalp and skull thickness choice were applied with this fit; the slab model used thicknesses of 2 and 5 mm, at a brain-to-scalp flow ratio of 3.75, whereas the curved model used thicknesses of 3 and 7 mm, at a flow ratio of 3.59. We observe a slight overcompensation of ${\mathrm{CBF}}_{\mathrm{i}}$ during the pressure period in the slab timecourse, whereas the curved model stays under 20% variation during applied pressure period. However, the responses of the two models are nearly the same during the hypercapnia measurement, with approximately a 10% overshoot in ${\mathrm{CBF}}_{\mathrm{i}}$ right after the ${\mathrm{CO}}_2$ period ends, but a full return to baseline by 250 s. Given the end results are very similar, this also being the case for the other subjects presented in the paper (results not shown)—we decided to use uniform scattering for the rest of the work.

Fig. 8

Monte Carlo slab and curved fits for subject 1 for (a) pressure modulation and (b) hypercapnia using variable scattering values across tissue layers.

3.3.5.

Subject 1: Monte Carlo four-layer models with variable diffusion coefficients for CSF

The final model tested for subject 1 was a four-layer MC curved geometry with uniform scattering across layers, as shown in the timecourse depicted in Fig. 9. The scalp thickness was maintained at 3 mm. One and two millimeters, for Figs. 9(a) and 9(b), respectively, were taken off the 9-mm skull to create a CSF layer between the skull and the brain. We held the diffusion coefficient for the CSF at values ranging from $1\times {10}^{-8}$ to $9\times {10}^{-8}{\mathrm{mm}}^2/\mathrm{s}$, a span of values approximately centered at the “biological-zero” of ${\mathrm{CBF}}_{\mathrm{i}}$ as measured by Busch et al.⁵²

Fig. 9

MC curved, four-layer fits with various diffusion coefficients used for the CSF layer for subject 1. A 1-mm CSF thickness was used for (a) and a 2-mm thickness was used for (b).

Both four-layer curved models using uniform scattering across layers displayed only marginal differences from the three-layer curved fit. We proceed with MC-3L fitting for the remaining subjects.

3.4.

Subjects 2 and 3: Monte Carlo-Based Three-Layer (MC-3L) Modeling Validation

Figure 10 shows subject 2 MC-3L fits for pressure modulation and hypercapnia timecourses, the latter overlaid with an expected TCD timecourse computed from the ${\mathrm{P}}_{\mathrm{ET}}{\mathrm{CO}}_2$ measurement. The slab model had a 3-mm scalp thickness and a 4-mm skull thickness, with a brain-to-scalp ${\mathrm{BF}}_{\mathrm{i}}$ ratio of 3.82; the curved model had respective thicknesses of 5 and 5 mm, with the brain/scalp ${\mathrm{BF}}_{\mathrm{i}}$ ratio of 3.35.

Fig. 10

Subject 2 timecourses fitted with MC-3L slab and curved models overlaid with the analytical fit at 5-, 25- and 30-mm separations and eTCD during (a) pressure modulation and (b) hypercapnia measurements.

In the pressure modulation maneuver, we note a 10% decrease in the 25-mm s–d ${\mathrm{BF}}_{\mathrm{i}}$ and moderate instability—possibly due to slight motion in the hand pressure—throughout the 30-mm detector timecourse fitted with the analytical model. For the hypercapnia measurement, we see a significant increase of 5-mm s–d ${\mathrm{BF}}_{\mathrm{i}}$, peaking a little over 50% of baseline at 230 s. This peak occurs about 47 s later than the eTCD peak. We observe a peak in the 25-mm ${\mathrm{BF}}_{\mathrm{i}}$ at 190 s and elevation during the recovery, failing to decrease below 25% by 250 s. We observe a similar continued increase in the post-hypercapnia period for the 30-mm detector at approximately the same percentage.

The MC-3L pressure modulation fits are more stable than what is observed in the 30 mm, although we still see mild variations, up to 10%, comparable to what is seen in the 25-mm analytical fit. In the hypercapnia measurement, both the slab and the curved models peak much closer to the expected TCD curve, albeit a little too early. A return to baseline within a few percentages for the slab and a complete return for the curved model are observed; however, a rise in ${\mathrm{CBF}}_{\mathrm{i}}$ after 220 s, similar to what is shown in the 30-mm s–d ${\mathrm{BF}}_{\mathrm{i}}$, is also noted for both models. This is likely related to subject motion due to hypercapnia discomfort.

Since MR scans were available for subject 2, we compare the total extracerebral thicknesses used in our MC fitting to what is observed in the MR scan. Figure 11 shows the estimated probe to brain distance as measured using the two methods described above in Sec. 2.5. The optical and MR data for this subject were acquired in one measurement session—that is, the fiducial marker denoting the probe location was preserved between the pressure modulation, hypercapnia, and MR structural scan. We observe a $\sim 12\text{-}\mathrm{mm}$ distance in Fig. 11(a); this falls in approximately the center of the distribution in Fig. 11(b). The thickness combination used for the MC-based slab model is about 5 mm less than what is measured, whereas the MC curved model underestimates the measured extracerebral thickness by 2 mm.

Fig. 11

Total extracerebral thickness for subject 2 as measured using a distance ruler on image software (a) and a mesh-based MATLAB function after segmentation (b).

In this subject, we were also able to run MC fitting that used a head segmented volume derived directly from the MRI scan taken as an input to the MC forward simulation. The head volume was segmented into four layers (scalp, skull, CSF, and brain tissue). However, the ${\mathrm{CBF}}_{\mathrm{i}}$ results from these fits were very unrealistic (high noise and exaggerated variation), suggesting there may be a systematic bias in the MC modeling of actual anatomy, a potential subject for future investigation.

Last, Fig. 12 displays the MC-3L slab and curved fits for subject 3. A significant drop in the 25- and 30-mm ${\mathrm{BF}}_{\mathrm{i}}$ is observed during the pressure period ($\sim 70\%$ and 60%, respectively) in the first plot. We observe $\sim 30\%$ less of a decrease in ${\mathrm{CBF}}_{\mathrm{i}}$ with MC fitting, though it is unable to keep estimated ${\mathrm{CBF}}_{\mathrm{i}}$ flat throughout the timecourse. Similarly, while the MC-3L plots do not show ${\mathrm{CBF}}_{\mathrm{i}}$ returning to baseline completely at the end of the hypercapnia measurement, they are still able to reduce the ${\mathrm{BF}}_{\mathrm{i}}$ elevation by 50% at 250 s compared to the analytical fit. The slab model used a 2-mm scalp and 7-mm skull thickness, whereas the curved model used 3 and 9 mm, respectively. The brain-to-scalp ${\mathrm{BF}}_{\mathrm{i}}$ ratio was 3.58 for the slab and 3.12 for the curved model.

Fig. 12

Subject 3 timecourses fitted with MC-3L slab and curved models overlaid with the analytical fit at 5-, 25- and 30-mm separations and eTCD during (a) pressure modulation and (b) hypercapnia measurements.

4.1.

Group Response to Hypercapnia

A gradual, significant increase in systemic blood flow due to ${\mathrm{CO}}_2$ administration can be observed in both the group average responses shown in 5-mm scalp probe and the 20-mm calf probe. A number of previous studies of hypercapnia-induced CBF changes with DCS have either used only a long separation detector in the DCS probe or an alternate modality for measuring SBF, such as the laser Doppler flowmetry used by Durduran et al.^15,43,53,54 Our study, as well as recent work by Milej et al.,⁵⁵ shows that measurements without a short separation detector may fail to reveal a possible SBF response that could be driving a significant portion of the observed long separation ${\mathrm{BF}}_{\mathrm{i}}$ timecourse. While Durduran et al. reported negligible SBF during the hypercapnia measurement, we note that the data shown were from a single subject, whereas the 5- and 20-mm data shown in our group average span 46 measurements across 18 subjects. As also mentioned in Sec. 3.2, over half of the sessions across our larger hypercapnia study contained at least one measurement with significant ($>20\%$ increase) 5-mm response. Selb et al.³³ have posited the possibility of SBF increase in their DCS measurements during hypercapnia, and several other studies indicate a positive correlation between systemic flow or flow-related parameters and percentage of ${\mathrm{CO}}_2$ inhaled.^56–58

On the other hand, our cerebrovascular reactivity as measured by TCD is consistent with what Coverdale et al.⁵⁹ showed in a hypercapnia-based validation study between TCD and phase contrast MRI (PC-MRI); they report an approximately $2.5\%\pm 2\%$ increase in CBFV per mmHg with TCD and a $3.6\%\pm 2\%\mathrm{CBFV}/\mathrm{mmHg}$ increase with PC-MRI. As shown in Sec. 3.1, we observe a 3.3% increase in CBV/mmHg with TCD. However, we note that Coverdale et al. also reported an $18\%\pm 8\%$ higher increase in their calculated arterial CBF versus MCAV during hypercapnia. In our study, the analytical model estimates of ${\mathrm{CBF}}_{\mathrm{i}}$ changes from DCS are below TCD changes, likely due to partial volume effects and the fact that a non-negligible minority of our DCS measurements showed little to no ${\mathrm{BF}}_{\mathrm{i}}$ increase in the long separation (probably due to low/no brain sensitivity,). The correlates of this apparent lack of sensitivity are still under investigation and a larger sample is likely needed to draw definite conclusions. Moreover, while studies such as Coverdale et al.⁵⁹ and Poulin et al.⁶⁰ are able to derive blood flow from the measured MCAV in TCD, the blood vessels measured by TCD and DCS are not the same; TCD measures the middle cerebral artery, and the estimated velocity is largely determined by the angle between the ultrasound beam and the direction of the blood flow.⁶¹ DCS, on the other hand, probes microvasculature in the frontal cortex. Apart from the size of the vessels themselves being sources of potentially different CBF values, and Liu et al.⁶² and Duffin et al.⁶³ have shown that the cerebrovascular reactivity to increased ${\mathrm{CO}}_2$ may not be homogeneous throughout parts of the brain.

Thus, we restrict the comparison of our observed relative ${\mathrm{CBF}}_{\mathrm{i}}$ increase during hypercapnia to previous studies that primarily used DCS. Durduran et al.¹⁵ observed an average of 2.4% increase of ${\mathrm{CBF}}_{\mathrm{i}}/\mathrm{mmHg}$ in five adult subjects, similar to our long separation group average shown in Sec. 3.2. Buckley et al.⁴³ reported a 49% increase to a 30-min long hypercapnia challenge containing $\sim 3\%$ inhaled ${\mathrm{CO}}_2$, a change significantly higher than what we observe. However, we note that their study was performed on children and not adults, along with a much longer hypercapnia period. Selb et al.³³ observed an unusually high average ${\mathrm{CBF}}_{\mathrm{i}}$ increase of 17% per mmHg; as noted above, this may be a result of a moderate systemic response as shown in their 8-mm separation data. Last, Milej et al.⁵⁵ performed hypercapnia on a group of five subjects and observed approximately a $0.35{\mathrm{cm}}^2/\mathrm{s}$ increase in long separation ${\mathrm{BF}}_{\mathrm{i}}$ for a 12-mmHg increase in ${\mathrm{etCO}}_2$. While the statistics given in their study cannot be directly compared to ours, the similarity of their observed extracerebral contamination will be noted shortly.

We abstain from directly comparing relative increases observed between our TCD and DCS data. The interpretations we make instead rest on the assumption that because the frontal cortex receives its blood supply from the middle cerebral arteries,⁶⁴ the temporal profile of blood flow in this cerebral region is likely to follow MCAV changes fairly closely. We first note that systemic responses shown in our DCS group averages (increases in 5-mm s–d separation on the forehead and 20-mm s–d separation on the calf, as discussed earlier) display distinct temporal curves and peaks than what is seen in our TCD group average. In particular, the systemic response peaks begin a gradual decrease approximately 30 and 90 s in the 5- and 20-mm s–d separations, respectively, after the end of the ${\mathrm{CO}}_2$ administration; the MCAV response as shown in TCD group average peaks drops sharply within 10 s after ${\mathrm{CO}}_2$ administration is stopped. Additionally, MCAV returns to baseline within 90 to 120 s at the end of hypercapnia, whereas the 5-mm and 20-mm response remains elevated 10% to 15% above baseline at 300 s. This elevated ${\mathrm{BF}}_{\mathrm{i}}$ response during recovery is also observed in 25- and 30-mm data. The 25-mm timecourse peaks earlier than the 5-mm (190 versus 210 s, closer to what we see in the TCD data) and recovers several percentages nearer to baseline at the end of the measurement, but we do not observe the full recovery as shown in TCD. The 30-mm s–d separation, probing deeper and thus containing more signal from the brain, peaks even sooner than the 25-mm curve at approximately 180 s, at the end of hypercapnia. Even so, it still recovers by about the same amount as the 25-mm distance at 300 s. As such, we interpret the group average ${\mathrm{BF}}_{\mathrm{i}}$ results as evidence of a systemic response to hypercapnia contaminating results shown in the long-separation detector timecourses. Milej et al.⁵⁵ demonstrate a similar phenomenon, where a lasting elevated response in the SBF is shown to contaminate the long separation ${\mathrm{BF}}_{\mathrm{i}}$, and was confirmed via repeating the hypercapnia challenge with a tourniquet tightened over the probe to suppress SBF. They estimated the extracerebral contamination level to be approximately $48\%\pm 18\%$ at 3 cm.

As mentioned above, significant scalp perfusion increases in our study were observed in a little over half the measurements. These measurements showed a wide range of scalp reaction magnitudes. The degree of scalp contamination in the long separation ${\mathrm{BF}}_{\mathrm{i}}$ will likely vary from subject-to-subject depending on factors, such as extracerebral thickness, amount of scalp reaction, and signal-to-noise ratio in the long separation signal. This phenomenon should be taken into account for future studies using hypercapnic validation and most importantly for studies using DCS to measure brain blood flow changes in adults clinically.

4.2.

Comparison of Monte Carlo Modeling to Previous Scalp Reduction Approaches

In the NIRS field, short-separation regression approaches have been extensively studied;^{21–26,65,66} however, for DCS, the non-linear nature of how the blood flow affects the different parts of the autocorrelation curve makes it less suitable for short separation regression. Thus, post-processing scalp physiology contamination reduction methods for DCS usually rely on a layered version of the analytical model to fit for ${\mathrm{CBF}}_{\mathrm{i}}$.^18,27,28 Our technique can be considered as an MC-based analog of Baker et al.’s analytically derived pressure modulation algorithm,²⁹ applied alongside hypercapnia measurements as further model validation. The primary advantages of MC-based modeling are the numerical stability and, as Li et al.³² and Shang et al.⁶⁷ have shown, the ability to accommodate realistic geometries as inputs into the forward simulation. As the development of MC photon propagation simulations becomes more sophisticated, MC modeling would allow us to be able to incorporate more nuanced physiological parameters into DCS fitting. This could include investigating the effect of probe location with respect to cortical folds or mapping DCS sensitivity across the entire forehead. Overall, data from MC simulations need not be restricted to in-silico validation of analytical models but can be used as inputs to fit in-vivo measurements themselves, particularly as MC methods are becoming increasingly faster and more accurate.^44,68

4.3.

Qualifications of Successful Monte Carlo Enhanced ${\mathrm{CBF}}_{\mathrm{i}}$ Quantification; Limitations of the Technique

We observe that pressure modulation, when used in conjunction with reasonable physiological estimates of brain-to-scalp ${\mathrm{BF}}_{\mathrm{i}}$ ratio and scalp/skull thicknesses, has the potential to help guide optimal models for MC fitting on a case-by-case basis. However, the maneuver may also serve as an indicator of brain sensitivity—measurements in which the long separation ${\mathrm{BF}}_{\mathrm{i}}$ shows equal or comparable decrease to the short separation ${\mathrm{BF}}_{\mathrm{i}}$ may indicate little or no cerebral sensitivity. Ideally, this should be checked in real time when setting up the measurement, and the DCS probe moved as needed to minimize the relative drop in the long separation measurement compared to the short separation. We interpret subject 3 as an example of this; the long separation ${\mathrm{BF}}_{\mathrm{i}}$ decreases significantly more during the pressure period as compared to the other two subjects. Correspondingly, the MC fits for this subject are unable to extract a full recovery to baseline at the end of the measurement. Like subject 3, baseline levels of sensitivity for future subjects may be inferred through protocols, such as pressure modulation or a hypercapnia challenge.

The criteria previously discussed in Sec. 2.4.2 serve as qualifications outlining the extent to which our study of multi-layered MC-based fitting can remove extracerebral influence from ${\mathrm{CBF}}_{\mathrm{i}}$. As a brief investigation (results not shown), we performed MC fits on subjects from the protocol with little to no evidence of brain sensitivity in the measurements, i.e., the long separation response did not differ much from the short in either the pressure modulation or hypercapnia measurements. We observed that the MC results were at best the same as the analytical fits; model geometries that showed reasonable brain-to-SBF ratios or increased ${\mathrm{CBF}}_{\mathrm{i}}$ during hypercapnia sometimes contained much noisier timecourses (presumably due to model mismatch). Analogously, we conducted MC modeling on measurements that already showed high levels of brain sensitivity in the analytical fit, i.e., significantly higher and timelier increase during hypercapnia and nearly full return to baseline in the long separation. The MC timecourses for these measurements matched the analytical fits. With the current resolution and accuracy in MC photon simulation, fitting DCS data using MC modeling seem to prove most advantageous under specific conditions, where the long separation analytical fit is showing evidence of both brain sensitivity and cross-talk from systemic physiology (verified with short-separation measurements).

Last, we acknowledge the computational limitations imposed by the use of the MC-based models. The two computationally intensive steps are the MC simulation itself, and the post-processing of MC history files to generate autocorrelation curves. With GPU acceleration, the MC simulation for 1 billion photons can completed within a few minutes (and is getting faster every year). Post-processing MC history files to obtain $g_2$ curves for given ${\mathrm{BF}}_{\mathrm{i}}$ values for each layer take on order of 50 to 100 ms in MATLAB if double-precision GPU acceleration is available or about 5 times longer using CPU only. Once the MC geometry is selected data can be fit in real time on a high end workstation computer as long as scattering is assumed to remain constant (this allows the reuse of the MC history file without re-running the simulation). However, there is a preparatory phase, in which data from one or more pressure modulation maneuvers are processed using a variety of geometries (varying scalp and skull thickness as explained in Sec. 2.4.3)—this could take half an hour or so on today’s hardware using MATLAB code. With highly optimized dedicated code, we believe the preparatory phase can be reduced to a few minutes, after which real-time MC-based monitoring will be available, e.g. if the goal is to offer brain perfusion monitoring in a clinical setting.

4.4.

Monte Carlo Geometry Choices

As seen in subject 1, two-layered MC-based models appear to fail—in particular, increasing extracerebral thickness closer to what is physiologically accurate leads to overcompensation and a sharp decrease in ${\mathrm{CBF}}_{\mathrm{i}}$ during recovery period post-hypercapnia, and no thickness setting follows the TCD measurement temporal profile effectively. The only two-layer models that appear to have reasonable post-hypercapnia recovery responses have 1-2 mm extracerebral thicknesses, which are too thin to be considered realistic. The three-layered models containing a constant diffusion coefficient value for a 4- to 10-mm slice—in other words, with a layer representing the solid slice of the skull—between the probe and the brain succeed in revealing more accurate ${\mathrm{CBF}}_{\mathrm{i}}$ responses despite contrasting superficial influence at much more reasonable extracerebral thicknesses. Disregarding the presence of low blood flow in the skull renders the two-layer models as highly unrealistic, contributing to their failure. Using varying scattering coefficient values across the three-layers in the MC-3L model seems to give a similar result as to their respective homogeneous scattering models, though at different scalp and skull thickness combinations. This is expected, as adjusting the scattering and adjusting the superficial layer thickness can both be interpreted as modifying sensitivity to the brain; changes in the former should result in changes to the latter, and vice versa. Given the substantial uncertainty regarding the actual scattering properties of inner tissue layers, the use of uniform scattering values appears warranted in DCS measurement analysis. Last, adding a 1- to 2-mm extra layer for CSF and using a four-layer model seems to have little impact on the final results.

We observe similar responses in the slab and the curved models in both pressure modulation and hypercapnia timecourses after tuning. In all subjects, the difference in the timecourses between the two geometries is close to negligible; however, the optimized slab extracerebral thickness combinations across subjects fall a few millimeters short of their respective curved thickness combinations. Comparing these values to the probe-to-brain distances measured on the structural scan for subject 2 reveals a closer match to the curved extracerebral distances than to the slab. In situations where either the scalp, skull, or total extracerebral thickness can be acquired and subsequently used as an input for model tuning, a curved geometry may be a more accurate and useful model. In other applications where this type of anatomical information is unknown and may never be acquired, the slab and curved models can be seen as functionally equal. With either case, we note the importance of model tuning in MC-based fitting or similar methods where the superficial layer thicknesses are adjusted over a range of values. Some a priori or other objective criteria—in this case pressure modulation and the analytical brain-to-scalp ${\mathrm{BF}}_{\mathrm{i}}$ ratio—must be used to guide the model to the correct scalp and skull thicknesses. Future research or clinical applications of DCS will involve situations where the expected CBF response cannot be predicted or is unknown; these circumstances necessitate some objective methodology to choose the most accurate model out of all possible scalp and skull thicknesses.

While minor differences are observed in the slab and the curved geometry fits, both can do several things. First, as seen in all subjects, they can provide a more stable ${\mathrm{CBF}}_{\mathrm{i}}$ during the pressure modulation maneuver than what is seen in the analytical long-separation fit. Second, the MC models can bring ${\mathrm{CBF}}_{\mathrm{i}}$ closer to baseline in post-hypercapnic recovery when superficial contamination is observed in the analytical long-separation fit. Third, as seen especially in subject 1, MC models can substantially remove of scalp influence to the temporal curve of the fit, i.e., they can remove a time lag caused by the scalp flow increase in the DCS ${\mathrm{CBF}}_{\mathrm{i}}$ estimate peaking and then decreasing.