1.

Introduction

Cerebral perfusion is a key parameter for the monitoring of cerebral ischemia. Position emission tomography (PET), magnetic resonance imaging (MRI), and lately also computer tomography (CT) can approximate this parameter with different tracer and bolus techniques. However, they rely on complex equipment and cannot be used at the bedside. Potentially, near infrared spectroscopy (NIRS) offers an inexpensive bedside monitoring tool when used for tracking an optical dye bolus. The reasoning is that the arrival time, i.e., the delay between the bolus signal and the time of intraveneous injection and the shape (width) of the bolus allow perfusion parameters to be reliably approximated.

For the last few years the optical dye indocyanine green (ICG) has been used in NIRS bolus measurements. When ICG signals are measured in babies and children, the bolus peak is less than 10 s in width.^{1 2 3} In this case Fick’s principle is applied to quantify the blood flow. Based on an experimental study in pigs, Kuebler et al.;⁴ derived a so-called blood flow index from the initial slope of the bolus signal. Alternatively, measurements of the ICG bolus signals at the inlet (aorta) and the outlet (jugular bulb) of the brain were exploited for perfusion quantification.^{5 6} However, this provides a measure of the global cerebral blood flow only.

The key obstacle to wide application of NIRS bolus methods in adults is the limited photon penetration depth into the tissue. On the one hand, this will always restrict the detection of perfusion impairments to the upper cortical layers when applied to adults. Also, the optical data will contain a signal not only from the cortex but one from the skin and skull as well. Therefore a method by which to separate the intracranial signal from extracranial contamination is required. This problem is obviously less severe in flow monitoring of babies or young children where the skull is considerably thinner.

The problem of extracerebral signals is highlighted by analyzing standard MR-bolus signals using gadolinium-DTPA as the contrast agent. In Figure 1 extra- and intracerebral bolus signals measured in a healthy volunteer are compared. The intracerebral time course [Figure 1(b)] at the six measurement positions (cf. the inset of Figure 1) has a minimum at t=15 s with a width [full width at half maximum (FWHM)] <10 s. In contrast, the extracerebral signal [Figure 1(a)] has a different shape with a peak delayed by a few seconds and a much slower recovery to the baseline. These extracerebral signals strongly depend on the measurement position. We found similar shapes in all volunteers analyzed. It is obvious that these different time courses impose a problem for perfusion measurements performed by NIRS.

Figure 1

Time courses of MRI perfusion signals using Gd-DTPA as a tracer normalized at the time of injection (t=0 s). Bolus signals were (a) analyzed at six extracerebral sites and (b) compared to those of adjacent intracerebral sites (marked by the arrows in the inset).

To solve this problem and to establish NIRS perfusion measurements, it has been suggested that one should subtract the signal/attenuation changes received from different source–detector positions. This is based on the assumption that light traveling through the extracranial tissue is almost identical for all source–detector distances.⁷ Moreover, McCormick et al.; injected ICG boli into both the exposed internal carotid artery (ICA) and the external carotid artery (ECA).⁸ When detecting at detector spacings of 1 and 2.7 cm, a bolus injected into the ECA resulted in a similar time course and magnitude in both detectors, while injection into the ICA was followed by a signal at larger detector spacing only.

The major approach we follow in this paper is based on the observation that different measurement parameters like intensity and photon time of flight have different sensitivity to different depths of tissue.^{9 10} Therefore such measurements do provide information on the separation of extra- and intracerebral signal. The main objective of this paper is to find out whether frequency-domain or multidistance measurements can help to better separate extra- from intracerebral absorption changes. The concept is that the photon mean time-of-flight parameter (phase) obtained from the frequency-domain technique is more sensitive to absorption changes in deeper tissue layers than intensity measurements. This sensitivity is estimated from a physical model for light transport in tissue based on Monte Carlo simulations. In principle time-domain measurements contain the maximum amount of available information for a depth resolved measurement.^{11 12} Here we limit ourselves to frequency domain and the following questions. (1) Is the information contained in the two measurement parameters, intensity and mean time of flight (phase) sufficient for signal separation? (2) How reliable is this separation when the basic assumption for the underlying model of photon transport is varied? (3) For its lower demand on instrumentation, multidistance intensity-only measurements seem preferable. Its potential use for the separation of cortical ICG bolus measurements is therefore examined as well.

2.

Methods

2.2.

Depth Resolved Analysis

NIRS intensities (I) were converted into attenuation changes ΔA=log(I₀/I), where I₀ is the intensity at t=0 s (injection time). Changes in phase Δϕ (in degrees) are equivalent to changes in the mean time of flight Δ〈t〉, ¹⁴

Eq. (1)

$$\Delta 〈t〉=\Delta \text{Φ}/({\mathit{\nu }}_{\mathit{M}}\cdot 360°)\text{,}$$

where ν_M is the modulation frequency (110 MHz) of the spectrometer. The conversion of the recorded changes of ΔA and Δ〈t〉 into changes in absorption coefficient Δμ_a (i.e., the ICG concentration) is based on the following assumption. First, the tissue is separated into an upper and a lower compartment separated at depth z_s with (unknown) absorption changes Δμ_a,low and Δμ_a,up . Second, the experimental measurement parameters ΔA and Δ〈t〉 are written as

Eq. (2)

$$\Delta \mathit{A}={\mathit{l}}_{\text{low}}\cdot \Delta {\mu }_{a,\text{low}}+{\mathit{l}}_{\text{up}}\cdot \Delta {\mu }_{a,\text{up}}\text{,}$$

$$\Delta 〈t〉={\mathit{m}}_{\text{low}}\cdot \Delta {\mu }_{a,\text{low}}+{\mathit{m}}_{\text{up}}\cdot \Delta {\mu }_{a,\text{up}}\text{.}$$

Here, l _low , l _up , m _low , and m _up are the sensitivities of A and 〈t〉 with respect to changes in μ_a in the lower and upper compartments. Based on these equations, Δμ_a can readily be obtained for both compartments once l_j and m_j are known and ΔA and Δ〈t〉 are measured.

Monte Carlo (MC) simulations were performed to estimate the sensitivity factors l_j and m_j for layered media. Details of the procedure are found in papers by Steinbrink et al.; and by Uludag et al.;^{11 12 15} l_j is equivalent to the (mean) optical pathlength of the photons in tissue layer j (j for low or up) and is calculated from

Eq. (3)

$${\mathit{l}}_{\mathit{j}}={\displaystyle \sum _i}{\mathit{l}}_{\mathit{ij}}{\mathit{W}}_{\mathit{i}}/{\displaystyle \sum _{\mathit{i}}}{\mathit{W}}_{\mathit{i}}\text{,}$$

where W_i is the survival weight of the ith photon bundle that reaches the detector and l_ij is the pathlength of the ith photon bundle in layer j. The summation is over all photon bundles (i). In a homogeneous medium, m_j corresponds to the variance in photon temporal distribution m_j=(〈t〉²−〈t²〉)⋅v_j, where v_j is the velocity of light in layer j. In an inhomogeneous medium it is¹¹

Eq. (4)

$${\mathit{m}}_{\mathit{j}}={\mathit{l}}_{\mathit{j}}\cdot 〈t〉-\sum _k\frac{\mathit{l}}{{\mathit{v}}_{\mathit{k}}}〈{\mathit{l}}_{\mathit{k}}{\mathit{l}}_{\mathit{j}}〉\text{,}$$

where the cross term is

$$〈{\mathit{l}}_{\mathit{k}}{\mathit{l}}_{\mathit{j}}〉=\frac{{\mathit{\sum }}_{\mathit{i}}{\mathit{l}}_{\mathit{ij}}{\mathit{l}}_{\mathit{ik}}{\mathit{W}}_{\mathit{i}}}{{\mathit{\sum }}_{\mathit{i}}{\mathit{W}}_{\mathit{i}}}\text{.}$$

The medium was assumed to consist of nine layers 2 mm thick followed by a semi-infinite half space. The medium was homogeneous, i.e., it had the same tissue optical properties (μ_a,μ_s ^′,n) independent of the depth. Anatomical structures such as the skull or the cerebro-spinal fluid (CSF) were disregarded. Photons (typically 4.2×10⁷ ) were launched at one position and the photons reflected were detected in concentric rings 3 mm in width at various distances from the source. Scattering anisotropies were disregarded. The refractive index of the medium was n=1.4. Fresnel reflection at the boundary of the medium was not included. The survival weight of the detected photons was calculated as well as the pathlength distribution (equivalent to the time distribution) for each layer.

In Figure 2 l_j and m_j are plotted for the layered medium assuming μ_s ^′=1 mm ⁻¹ and μ_a=0.01, 0.02, and 0.03 mm⁻¹. These μ_a values sufficiently cover the range expected for biological tissue in the NIR wavelength range. l_j is given for source–detector distances of r=10 and 30 mm, while m_j is plotted for r=30 mm only. It is apparent that measurements of 〈t〉 are more sensitive to variations of Δμ_a in deeper layers. While l_j has the highest sensitivity in the layer with depth between 4 and 6 mm with regard to r=30 mm, m_j peaks for the layer at 8–10 mm depth. By choosing a short source–detector distance of 10 mm, the sensitivity of l_j to changes at depth >8 mm is very low. For the analysis outlined in Eq. (2), the sensitivities of the lower and upper compartments separated at a depth z_s were obtained by adding the l_j and m_j values for z<z_s and z>z_s, respectively. Similar to in Eq. (2), a separation of the upper and lower absorption changes based on attenuation measurements (ΔA_r1,ΔA_r2) at different source–detector separations, r₁ and r₂, is theoretically feasible when the pathlength in the layers ( l _low,r1, l _low,r2, l _up,r1, and l _up,r2 ) are included. The whole set of estimated values of l_j and m_j as a function of the depth for 0.01 mm ⁻¹<μ_a<0.07 mm ⁻¹, 1 mm ⁻¹<μ_s ^′<3 mm, and 6 mm<r<60 mm can be obtained from the authors.

Figure 2

Sensitivity of attenuation ( l_j; closed symbols) and mean time of flight ( m_j; open symbols) with respect to changes in μ_a in a layer of thickness 2 mm at various depths based on Monte Carlo simulations. μ_s ^′=1 mm ⁻¹ and μ_a values of 0.01, 0.02, and 0.03 mm⁻¹ were assumed. Source–detector distances of r=10 and 30 mm were chosen for l_j and r=30 mm for m_j.

Equations (2 3 4) are based on the linear approximation which is accurate only when the fractional change in μ_a is small. Estimation of higher order terms with sufficient signal to noise ratio would require considerably more computational time in the MC simulation so it was not considered here. Therefore we have no means by which to estimate the error when the fractional change in μ_a is of the order of 20–30 as calculated for the experiments discussed below.

3.

Results

Models for the separation of ICG concentrations into different tissue compartments based on a single wavelength analysis require that there are no concurrent changes in hemoglobin. This was tested by spectral analysis of attenuation changes within the range of 740–880 nm. An example is given in Figure 3(a) where the changes in hemoglobin concentration (deoxy-Hb, oxy-Hb) and ICG are plotted for the time period following three bolus injections. There are fluctuations in the hemoglobin traces, however, within the noise level they are not correlated with the ICG bolus. This is apparent from the bolus attenuation spectrum [see Figure 3(b)] that closely fits the in vitro ICG absorption spectrum. The small discrepancies in the peak absorption wavelength are due to uncertainties in the ICG absorption spectrum itself or to the wavelength dependence of the optical pathlength. The in vitro ICG absorption spectrum is not corrected for the wavelength dependence of the optical pathlength since the major ICG signal contribution originates from the upper tissue layers. In these compartments the effect of the wavelength dependent absorption spectrum on the pathlength is small and was therefore neglected.

Figure 3

(a) Calculated concentration changes in ICG, oxy-Hb, and deoxy-Hb following three ICG bolus injections. For ICG the values are arbitrarily scaled. (b) Attenuation spectrum averaged over the period of the first ICG signal (t=490–570 s). For comparison, the arbitrarily scaled in vitro absorption spectrum of ICG is shown.

An example of the bolus time course measured with the frequency-domain spectrometer in a volunteer is given in Figure 4. Averaged data of six consecutive ICG injections are shown for a source–detector separation of 30 mm at λ=780 nm. The maximum change in attenuation [Figure 4(a)] is reached about 27 s after injection. As already seen in Figure 3 for the ICG trace, the recovery to baseline values is slow with an initially faster decay. The 〈t〉 time course rises a few seconds earlier and peaks about 22 s after injection. Note that the mean time of flight signal decreases (at the peak, Δ〈t〉=−23 ps ), consistent with a rise in the absorption coefficient. For t>40 s, the recovery to the baseline of ΔA and Δ〈t〉 is very similar. The bolus signals (Δμ_a) for an upper compartment (0–12 mm) and a lower compartment (>12 mm) were calculated from these data based on Eq. (2) and are shown in Figure 4(b). A transport scattering coefficient of μ_s ^′=1 mm ⁻¹ and an absorption coefficient of μ_a=0.02 mm ⁻¹ were assumed for estimation of l_j and m_j. To separate the lines in Figure 4(b), the baseline of Δμ_a,low was shifted by −0.002 mm⁻¹. The signal of the lower compartment peaks at about 20 s with a FWHM of ≈10 s, followed by a slower recovery to the baseline (t≈30 s). In contrast, the Δμ_a time course of the upper compartment is delayed with a slower recovery to the baseline. These differences closely resemble those of the Gd bolus for the extra- and intracerebral compartments shown in Figure 1.

Figure 4

(a) Measured change in attenuation and mean time of flight averaged over six consecutive injections of ICG. (b) Conversion of the data in (a) into changes in absorption coefficient of an upper and a lower compartment based on the assumption μ_a=0.02 mm ⁻¹, μ_s ^′=1 mm ⁻¹ and that the upper compartment has a thickness of 12 mm and the lower compartment is semi-infinite.

The obvious question is whether the calculated time course of μ_a is affected by basic assumptions of the tissue optical properties and the layer geometry. In Figure 5 the calculated μ_a changes are shown for the two compartments when the thickness of the upper compartment (z_s) is varied between 6 and 18 mm. The calculations are based on the experimental data shown in Figure 4(a). In all cases μ_s ^′=1 mm ⁻¹ and μ_a=0.02 mm ⁻¹ were assumed. The results indicate that variation of z_s has only a minor effect on the shape of the Δμ_a time course of the upper compartment whereas the magnitude increases with an increase of z_s. The Δμ_a time course of the lower compartment has a clear bolus peak in all cases. However, for z_s⩽8 mm the recovery to the baseline is incomplete for times >25 s. For z_s⩾14 mm the calculated Δμ_a changes (i.e., ICG concentration changes) are negative for t>25 s. This is clearly unreasonable for an increase in dye concentration, and puts a theoretical limit on depth penetration.

Figure 5

Calculated changes in μ_a for (a) the upper compartment and (b) the corresponding lower compartment based on the experimental data shown in Fig. 4(a). The two compartments were separated at depth z_s (between 6 and 18 mm). In all cases μ_s ^′=1 mm ⁻¹ and μ_a=0.02 mm ⁻¹ were assumed.

Next, the influence of the assumed background tissue absorption coefficient was examined. Figure 6 gives the results for μ_a values of 0.01, 0.02, and 0.03 mm⁻¹, i.e., a variation that sufficiently covers the range expected for NIR wave lengths. In all cases a thickness of the upper compartment of z_s=12 mm was assumed and a scattering coefficient of μ_s ^′=1 mm ⁻¹. Again, the Δμ_a time course of the upper compartment is hardly affected by the choice of parameters. For the lower compartment, however, there is a strong variation in magnitude that can be explained by the strong decline in the m_j sensitivity for deeper layers when the absorption coefficient is increased (see Figure 2). Again, the recovery to the baseline directly after the bolus peak is influenced by the variations in background absorption. In a similar way, variations in the scattering coefficient of the model tissue were considered. By changing μ_s ^′ from 1 to 2 mm⁻¹ the effect is by and large on the magnitude of the calculated Δμ_a changes only (data not shown).

Figure 6

Calculated changes in μ_a for an upper layer and the corresponding lower layer for the data of Fig. 4(a). For the calculation μ_a=0.01, 0.02, and 0.03 mm⁻¹ were assumed. In all three cases the two layers were separated at a depth z_s=12 mm; μ_s ^′ was 1 mm⁻¹.

The shape of the bolus signal strongly varied among different volunteers. In some experiments there was no apparent difference between the time course of attenuation and the mean time of flight with a response similar to the attenuation data shown above. In other data sets the short bolus was already apparent in the raw data of the Δ〈t〉 signal. An example is given in Figure 7, where Δ〈t〉 has a clear ICG peak of about 10 s width that is earlier by more than 5 s compared to ΔA.

Figure 7

Raw data of Δ〈t〉 and ΔA averaged over two ICG boli injected at t=0 s measured in an adult volunteer. The data are consistent with a faster and shorter bolus signal in the brain that is primarily detected by the mean time-of-flight signal compared to a slower bolus signal in the skin and skull primarily measured by intensity changes.

To further describe the bolus signals the rising edge was analyzed as outlined in Figure 8(a). The raw data of ΔA or Δ〈t〉 were first approximated by a polynomial of fifth order between the time of injection and time of maximal ICG signal. From these smoothed lines, the rise time t_r defined as the difference between the 10 and 90 points in ΔA or Δ〈t〉 were obtained. Furthermore, based on the polynomial fit, the time of the maximal slope t _ms was calculated. For each subject t _ms and t_r were determined and gave mean values of 〈t _ms (ΔA)〉=20.5±5.5 s and 〈t_r(ΔA)〉=10.7±3.5 s for attenuation and 〈t _ms (Δ〈t〉)〉=18.4±3.5 s and 〈t_r(Δ〈t〉)〉=10.4±4.1 s for the mean time of flight. In Figure 8(b) the differences in t _ms and t_r for attenuation and mean time of flight data are shown for all subjects. While there is no significant difference in the mean rise time [Δ〈t_r(ΔA,Δ〈t〉)〉=0.3±2.5 s], evaluation of the time of maximal slope does indicate a more pronounced difference [Δ〈t _ms (ΔA,Δ〈t〉)〉=2.1±3.1 s] for the two measurement parameters.

Figure 8

(a) Evaluation of rise time t_r defined as the difference between 10 and 90 points of the ICG induced attenuation increase and time of maximal slope t _ms . (b) t _ms (ΔA)−t _ms (Δ〈t〉) vs t_r(ΔA)−t_r(Δ〈t〉) for measurements in adult volunteers.

The finding that analysis of ΔA and Δ〈t〉 allows different tissue compartments to be separated is further supported by data previously measured in babies during a cardiac pulmonary by-pass (CPB) operation.¹⁶ In these babies ICG was injected into the CPB circuit at concentrations of 0.1 mg/kg body weight. It was found that the shape of the bolus time course of ΔA and Δ〈t〉 is very similar (see the example for a 13 week old baby in Figure 9) with no apparent differences in the ICG arrival time or bolus width. The likely explanation is that in babies the skin and skull are thin and therefore both optical measurement parameters sample predominantly cere- bral tissue with only small contributions from extracerebral compartments (cf. Figure 2).

Figure 9

Changes in attenuation and mean time of flight following two ICG injections measured on the head of a 13 week old baby during cardiac by-pass operation (see Ref. 16). λ=806 nm, r=30 mm. For both injections ( t=2166 and 2298 s) the second passage of the dye is visible.

The separation of extra- and intracranial bolus signals was also attempted by analyzing multidistance measurements ( r=10, 20, 30, and 40 mm). An example of bolus time courses is shown in Figure 10 for source–detector separations of 10 and 30 mm. The signals are similar with a slightly higher amplitude for the larger separation. The assumption is that at larger distances a larger signal is sampled from deeper layers and at the same time a moderately larger signal from layers close to the surface (cf. Figure 2). The concentration changes (Δc) of both separations were first arbitrarily scaled to give similar overall signal magnitudes and then subtracted. When a scaling factor of 1.3 is used, the time course of Δc(30 mm)−Δc(10 mm) 1.3 peaks about 25 s after bolus injection with values close to zero for larger times. This is consistent with a greater depth sensitivity for the larger interoptode distance (IOD) (see Figure 2). However, the choice of scaling factor is somewhat arbitrary. It is consistent with the physical model (Figure 2) only when assuming z_s=2 mm [ l_j(z=0–2 mm)=17.94 and 14.65 mm for r=30 and 10 mm, respectively]. Slight variations in z_s, i.e., in the thickness of the skin, result in large changes in the scaling factor derived from the Monte Carlo model. While for separations based on ΔA and Δ〈t〉 variations in z_s (Figure 5), μ_a (Figure 6) or μ_s ^′ have a minor influence on the time course of the bolus signals, similar variations strongly affect the multidistance approach. Introducing a model derived rather than arbitrary scaling is therefore not justified as long as there are not precise data on the tissue’s optical properties and anatomy. In most of the data we found the analysis based on multidistance attenuation data to be unreliable for the recovery of peak bolus time courses when based on pathlength estimates as shown in Figure 2.

Figure 10

ICG concentration changes Δc at interoptode distances of r=30 and 10 mm for two consecutive bolus injections (injection times set to t=0 s ). In the lower section the traces Δc (30–10 mm) 1.3 are shown. For graphical reasons an offset of −0.002 was chosen.

4.

Discussion

The data presented here demonstrate the capability of NIRS to noninvasively monitor cerebral perfusion in the human adult based on an ICG-bolus tracking technique similar to the Gd-DTPA bolus technique established for MRI. While in simple intensity measurements the signal from the cortex is masked by that from the skin and skull, the application of frequency-domain measurements in combination with a model based analysis allows the signals from an upper layer corresponding to the skin and skull and a lower layer that reaches the cerebral cortex to be separated. In volunteers it was regularly found that the arrival of the bolus signal is delayed in the upper layer (skin and skull) compared to in the lower layer. This is consistent with a lower perfusion rate and a longer retention time in the skin compared to the brain which is one of the most perfused organs. This finding is supported by analysis of MRI-perfusion data, which show a similar difference between the MTT in extra versus intracerebral tissue (Figure 1). Such a measure for extra- and intracerebral perfusion with a bedside method is of great potential relevance for the monitoring of vascularly diseased or traumatized patients. It should be noted that the approach presented here is not an alternative but may serve as an extremely helpful extension to the MTT assessment performed by Gd-DTPA enhanced functional MRI. Since the two techniques allow simultaneous application there is the perspective of an initial combined approach that is followed up on the intensive care or stroke unit at short intervals. The initial delay between the extra- and intracerebral bolus will serve as a baseline. To offset the effect of skin perfusion changes, a multisite measurement may be required. Changes of this delay may readily be explained by worsening or normalization of cerebral perfusion in the tissue illuminated. We are currently exploring the technique described here in patients undergoing cardial by-pass surgery.

Our data suggest that more demanding time- or frequency-domain equipment is a prerequisite for reliable monitoring. However, preliminary measurements with a multisource multidetector topography system based on pure intensity data indicate that the spatial resolution might offset a lack of depth resolution (these data will be discussed in a separate paper). For reliable monitoring the combination of multisite and frequency- (time-) domain measurements seems even more desirable.

The core argument for the potential to differentiate between the extracerebral and the truly intracerebral compartments is the differences between the bolus shape and their relative latency. However, we regularly observed slow oscillations, best seen in the unaveraged data (see, e.g., Figures 6 and 7) with frequencies close to 0.1 Hz. They are much more prominent in the Δ〈t〉 traces compared to in the ΔA traces. The oscillations have recently been examined by NIRS¹⁷ and are likely a cortical phenomenon, thus providing an additional argument for the stronger contribution of cerebral hemodynamics to the Δ〈t〉 measurements. Another argument can be given by a comparison between the data from infants and adults. The lack of a delay between the two measurement parameters in infants corresponds to their much thinner skull and skin, and provides an anatomical explanation for this difference from adults.