Accuracy of oxygen saturation and total hemoglobin estimates in the neonatal brain using the semi-inﬁnite slab model for FD-NIRS data analysis

: Frequency domain near-infrared spectroscopy (FD-NIRS) is a non-invasive method for measuring optical absorption in the brain. Common data analysis procedures for FD-NIRS data assume the head is a semi-inﬁnite, homogenous medium. This assumption introduces bias in estimates of absorption ( µ a ), scattering ( µ ′ s ), tissue oxygen saturation (S t O 2 ), and total hemoglobin (HbT). Previous works have investigated the accuracy of recovered µ a values under this assumption. The purpose of this study was to examine the accuracy of recovered S t O 2 and HbT values in FD-NIRS measurements of the neonatal brain. We used Monte Carlo methods to compute light propagation through a neonate head model in order to simulate FD-NIRS measurements at 690 nm and 830 nm. We recovered µ a , µ ′ s , S t O 2 , and HbT using common analysis procedures that assume a semi-inﬁnite, homogenous medium and compared the recovered values to simulated values. Additionally, we characterized the effects of curvature via simulations on homogenous spheres of varying radius. Lastly, we investigated the effects of varying amounts of extra-axial ﬂuid. Curvature induced underestimation of µ a , µ ′ s , and HbT, but had minimal effects on S t O 2 . For the morphologically normal neonate head model, the mean absolute percent errors (MAPE) of recovered µ a values were 12% and 7% for 690 nm and 830 nm, respectively, when source-detector separation was at least 20 mm. The MAPE for recovered S t


Abstract:
Frequency domain near-infrared spectroscopy (FD-NIRS) is a non-invasive method for measuring optical absorption in the brain. Common data analysis procedures for FD-NIRS data assume the head is a semi-infinite, homogenous medium. This assumption introduces bias in estimates of absorption (µ a ), scattering (µ ′ s ), tissue oxygen saturation (S t O 2 ), and total hemoglobin (HbT). Previous works have investigated the accuracy of recovered µ a values under this assumption. The purpose of this study was to examine the accuracy of recovered S t O 2 and HbT values in FD-NIRS measurements of the neonatal brain. We used Monte Carlo methods to compute light propagation through a neonate head model in order to simulate FD-NIRS measurements at 690 nm and 830 nm. We recovered µ a , µ ′ s , S t O 2 , and HbT using common analysis procedures that assume a semi-infinite, homogenous medium and compared the recovered values to simulated values. Additionally, we characterized the effects of curvature via simulations on homogenous spheres of varying radius. Lastly, we investigated the effects of varying amounts of extra-axial fluid. Curvature induced underestimation of µ a , µ ′ s , and HbT, but had minimal effects on S t O 2 . For the morphologically normal neonate head model, the mean absolute percent errors (MAPE) of recovered µ a values were 12% and 7% for 690 nm and 830 nm, respectively, when source-detector separation was at least 20 mm. The MAPE for recovered S t O 2 and HbT were 6% and 9%, respectively. Larger relative errors were observed (∼20-30%), especially as S t O 2 and HbT deviated from normal values. Excess CSF around the brain caused very large errors in µ a , µ ′ s , and HbT, but had little effect on S t O 2 .

Introduction
Near-infrared spectroscopy (NIRS) is a non-invasive method for measuring the optical absorption of cerebral blood due to hemoglobin [1]. Measurements acquired at multiple wavelengths of light can be used for spectral estimation of both oxy-hemoglobin (HbO 2 ) and deoxyhemoglobin (Hb), which subsequently provide estimates of tissue oxygen saturation (S t O 2 ) and total hemoglobin (HbT). Because NIRS is portable, non-ionizing, and non-restraining, it is an attractive technique for bedside measurements of cerebral physiology in neonates with applications in assessing brain development [2], detecting hypoxic-ischemic brain injuries [3,4], and monitoring during pediatric cardiac surgery [5].
Current clinical NIRS systems use continuous wave (CW-) light sources, in which light intensity is constant. These systems have shown promising evidence in monitoring trends in cerebral physiology, but lack quantitative accuracy in absolute estimates of baseline physiology [6]. Frequency domain (FD-) NIRS systems use a light source that is modulated in intensity by a sinusoidal function, which provides measurements of both amplitude and phase of the modulated light. The phase measurement contains information about the optical pathlength, which improves quantification.
A common method for analysis of FD-NIRS data assumes a homogenous, semi-infinite medium [7]. The existence of an analytical solution to light propagating in such a medium provides a convenient and computationally efficient method for data analysis, but introduces bias in recovered absorption values due to curvature/finite volume and the heterogeneous structure of tissues in the head. Dehaes et. al. investigated the accuracy of recovered optical absorption coefficients using this model by simulating FD-NIRS measurements on one-, two-, and three-layered segmentations of structural magnetic resonance images (MRI) and comparing recovered absorption coefficients to simulated values [8]. They found typical errors of 8-24% in recovered absorption coefficients in neonates. Based on the errors in absorption, they predicted errors of 8% and 11% in S t O 2 and HbT, respectively. In this work, we used simulations of light propagation through tissue to directly investigate the accuracy of S t O 2 and HbT estimates in neonates. We used a four-layer segmentation of a structural MRI from a full term, 24 day old neonate to simulate FD-NIRS measurements. We used literature values for the scalp, skull, and cerebrospinal fluid (CSF) and parameterized the brain properties in terms of S t O 2 and HbT. Simulations were performed for varying S t O 2 (30-100%) and HbT (15-120 µM) in the brain layer. Recovered S t O 2 , HbT, µ a , and µ ′ s values were compared to the simulated values. Additionally, we characterized the effects of curvature on recovered values via simulations on homogenous spheres of varying radius (30-120 mm). Lastly, we examined the effects of extra-axial fluid increasing the thickness of the CSF layer.

Simulation of measurements
Simulated NIRS measurements were obtained by computing light propagation through volumetric images via the Monte Carlo Extreme software (MCXLAB 0.7.9 [9]). The volumetric images contained tissue labels, specifying the optical properties compiled in Table 1. Simulations were performed using reduced scattering coefficients with anisotropy factor set to zero (isotropic scattering). For all simulations, wavelengths of 690 nm and 830 nm were used according to the findings in [10]. The optical properties of the brain were calculated based on tissue oxygen saturation (S t O 2 ) and total hemoglobin (HbT) content. Normal baseline values of 70% S t O 2 and 60 µM HbT were assumed for the brain. Additionally, background absorption from 70% water content was assumed for the brain voxels. The time-domain signal output from MCXLAB was converted to a frequency-domain measurement via discrete Fourier transform at a modulation frequencies of 50 MHz, 100 MHz, 150 MHz, and 200 MHz.

Calibration
The simulated FD-NIRS data were calibrated in the same manner as the calibration procedure for experimental data (e.g., [16]) by simulating an additional calibration data set on a homogenous slab of 180 mm × 180 mm × 100 mm with optical properties matching the brain at 70% S t O 2 and 60 µM HbT. The respective calibration terms for intensity and phase are given by where ψ theo and Φ theo are the theoretical intensity and phase values, and ψ slab and Φ slab are the intensity and phase measured on the calibration slab. The calibrated data is then given by where ψ cal and Φ cal are the calibrated intensity and phase values, and ψ and Φ are the measured intensity and phase values. The calibration procedure reduces errors due to the diffusion approximation, especially at short source-detector separation distances [8].

Recovery of optical properties and physiological parameters
The solution to light propagation through a semi-infinite, homogenous medium is described in detail by Fantini et. al. [7]. The exact equations used in this work were as follows: where ψ and Φ are the amplitude and phase measurements with angular modulation frequency of ω; ρ is the source-detector separation; µ a is the absorption coefficient; D is the diffusion constant equal to 1/(3µ a + 3µ ′ s ), in which µ ′ s is the reduced scattering coefficient; ν is the speed of light in the media; F ψ is a complicated function of many variables, including ρ, µ a , and D that is approximately constant; Φ 0 is the instrument phase offset. The terms V + and V − are given by Note that Eq. (3) is an approximation of the full solution in [7]. This is often further approximated to the following equations: where F Φ replaces the additional terms in Eq. (3b), and is assumed to be constant. We chose to use Eq. (3) for estimation of optical properties for this work, since we found Eq. (3) to produce more accurate results with uncalibrated data from simulations on a homogenous slab than Eq. (5). With either set of equations, the left-hand side quantities form a linear relationship with ρ with the following slope values: for the equations for ψ and Φ, respectively. The optical absorption and scattering coefficients are solved for algebraically in terms of the slope values: Note that in the continuous wave case (ω = 0), Φ is a constant and solving for µ a and µ ′ s is no longer possible.
The regression of Eqs. (5a) and (5b) is straightforward. The regression of Eq. (3a) and (3b) is performed iteratively, starting with an initial estimate of µ a and µ ′ s from regression of Eq. (5). The optical absorption coefficients recovered from Eq. (7a) above are related to hemoglobin by where the [HbO 2 ] and [Hb] are the tissue concentrations of oxy-hemoglobin and deoxyhemoglobin, respectively; ε HbO 2 ,λ and ε Hb,λ are the respective extinction coefficients for HbO 2 and Hb at a wavelength of λ ; B(λ ) is the background absorption at the wavelength λ . In this study, background absorption was assumed to be from 70% water content, in which B(λ ) was set to 70% of the absorption coefficient of water (absorption spectrum from [15]). With two or more wavelengths, Eq. (8) can be used to form a system of linear equations, which can be organized in matrix notation as where u is the vector containing the absorption coefficients, corrected for background absorption; E is the matrix gathering the extinction coefficients; h is the vector containing oxy-and deoxy-hemoglobin concentrations. The least-squares solution to the matrix expression is given by where the superscript T is the transpose operation.
Equations (10) and (11) give some insight into the propagation of errors from recovered µ a values to S t O 2 and HbT. In the limit that background absorption goes to zero (B(λ ) = 0), we find that [HbO 2 ] and [Hb] are linear combinations of the recovered µ a values with coefficients from the matrix given by (E T E) −1 E T . From Eq. (11), we can conclude that HbT is also a linear combination of recovered µ a values, and S t O 2 is a ratio with linear combinations of µ a values in the numerator and denominator. Thus, scaling all µ a values by a factor, c, will scale HbT by c, but have no effect on S t O 2 (when B(λ ) = 0). We can use this fact to conclude that if relative errors in µ a are similar across wavelengths, the errors in HbT will be sensitive to the magnitude errors in µ a , but errors in S t O 2 will be minimal; however, in the case that the relative errors in µ a are dissimilar across wavelengths, the errors in S t O 2 can be large or unpredictable.

Effects of curvature
In order to examine the effects of head curvature on estimates of S t O 2 and HbT, data were simulated on homogenous spheres of varying radius from 30-120 mm, as well as a flat slab (infinite radius) with dimensions of 180 mm × 180 mm × 100 mm. All voxels within the sphere or slab were labeled with optical properties matching the brain at 70% S t O 2 and 60 µM HbT. The volumetric images defining the optical properties were all created with 1 mm isotropic resolution. A single source position was simulated with detectors arranged in a linear fashion along the surface at distances from 10-40 mm in 1 mm increments. In order to decrease the effects of discretization of the spherical geometry, the position and orientation of the probe were randomized for each iteration. A total of 128 repetitions of 10 7 photon packets were simulated per wavelength for each simulation. Recovery of S t O 2 and HbT was performed for four subsets of the measurements based on source-detector separation: 10-25 mm, 15-30 mm, 20-35 mm, and 25-40 mm. Recovery of µ a , µ ′ s , S t O 2 , and HbT was performed via the methods in Section 2.3 and compared to the simulated values.

Simulations of varying S t O 2 and HbT in the brain
In order to investigate the magnitude of errors in a realistic model of an infant head with normal morphology, simulations were performed on a segmented structural MRI of a neonate (full term, 24 days old). Segmentation of gray matter, white matter, and CSF was performed via SPM8 (Wellcome Trust Centre for Neuroimaging, London, UK) using the UNC neonate atlas [17]. Segmentation of the scalp and skull layers was performed via manually thresholding the structural image after masking out the brain. A sample of the segmentation results is shown in Fig. 1 with uniformly distributed distances between 10-40 mm. The source and detector locations are shown in Fig. 1(b). Optical properties of the skin/scalp, skull, and CSF layers were fixed according to Table 1 for all simulations. The white and gray matter regions were combined to a single brain region with optical properties based on S t O 2 and HbT. Simulations were performed with varying S t O 2 (30-100%) and fixed HbT (60 µM). A second set of simulations were performed with varying HbT (15-120 µM) and fixed S t O 2 (70%). A total of 128 repetitions of 10 7 photon packets were simulated per wavelength for each simulation. Recovery of µ a , µ ′ s , S t O 2 , and HbT was performed via the methods in Section 2.3 for four subsets of the measurements according to source-detector separation: 10-25 mm, 15-30 mm, 20-35 mm, and 25-40 mm. Recovered values were compared with simulated values for the brain.

Effects of increased extra-axial fluid
In order to simulate the effects of extra-axial fluid, we increased the thickness of the CSF layer by performing a binary erosion of the brain mask with varying sized kernels and replacing the outer brain voxels with CSF in the resulting segmented volumes. This is in-effect similar to the actual morphological changes associated with extra-axial fluid, in which the build up of excess CSF in the subarachnoid space pushes the brain away from the the skull. The resulting segmentations had CSF layers that varied in thickness from approximately 1 mm to 4 mm with some spatial heterogeneity due to anatomy. The probe shown in Fig. 1(b) was used for simulation as in the previous section. Simulations were performed for varying S t O 2 (30-90%) and fixed HbT (60 µM). Recovery of µ a , µ ′ s , S t O 2 , and HbT was performed via the methods in Section 2.3 for four subsets of the measurements according to source-detector separation: 10-25 mm, 15-30 mm, 20-35 mm, and 25-40 mm. Figure 2 shows the results of simulating data on homogenous spheres and recovering µ a (top row) and µ ′ s (bottom row) via the methods in Section 2.3. In all cases, curvature induced underestimation of µ a and µ ′ s with more severe effects as radius of curvature decreased. There were no strong trends between the magnitude of errors and source-detector separation distance. Errors in µ ′ s were more severe than errors in µ a . Figure 3 shows the results from calculating S t O 2 and HbT from the recovered µ a values in Figs. 2(a)-2(d). In general, curvature had little effect on S t O 2 , but caused significant underestimation of HbT. There were no strong trends across source-detector distance for errors in HbT. The average circumference of the neonatal head was reported to be 343 mm [18], giving an approximate radius of curvature of 55 mm using a circular approximation. At a radius of curvature of 60 mm, the percent errors were 5-7%, 6-7%, 7-8%, and less than 1.5% for µ a,690 , µ a,830 , HbT, and S t O 2 , respectively. Figure 4 shows the normalized partial pathlength (i.e., fraction of the total pathlength) of the simulated photons through each tissue type in the neonate model as a function of sourcedetector distance. The data shown was for 70% S t O 2 and 60 µM HbT in the brain. Each point on the graph represents one of the source-detector measurements. In general, the fraction through the brain increased as a function of source-detector distance, while the fraction through the scalp and skull decreased with source-detector distance. There appeared to be a split between increasing and decreasing trends in the partial pathlength through CSF with source-detector separation, potentially due to anatomical variations in CSF across the head. Overall, the re-sults in Fig. 4 suggest that longer source-distances will be more sensitive to the brain and less sensitive to superficial layers. Figure 5 shows the results of simulating data with varying brain S t O 2 in the neonate model and recovering µ a (top row) and µ ′ s (bottom row) with the slab model. There was an overall trend of increasing accuracy with increasing source-detector distance in recovering µ a . There was a relatively large degree of underestimation of µ a,690 at low S t O 2 where µ a,690 is relatively large. Scattering coefficients were severely underestimated. Figure 6 shows the results from calculating S t O 2 and HbT from the recovered µ a values in Figs. 5(a)-5(d). There was an overall trend of increasing accuracy with increasing source-detector distance. A summary of the percent errors for simulations with varying S t O 2 is shown in Table 2.  Figure 7 shows the results of simulating data with varying HbT in the neonate model and recovering µ a (top row) and µ ′ s (bottom row) via the methods in Section 2.3. Recovered µ a values showed an overall trend of increasing accuracy with increasing source-detector distance, especially for high HbT values. Again, scattering coefficients were severely underestimated. Figure 8 shows the results from calculating S t O 2 and HbT from the recovered µ a values in Figs. 7(a)-7(d). There was a slight trend of increasing accuracy with source-detector distance in recovered S t O 2 values. There was an overall trend of increasing accuracy in recovered HbT values with increasing source-detector distance, especially at high HbT values. The percent errors for simulations with varying HbT is shown in Table 3.        Figure 10 shows the results of simulations with increased extra-axial fluid. Excess CSF caused very large errors in recovered µ a values, eventually breaking the analysis methods when CSF was too thick (dependent on source-detector distance). The errors in µ a propagated large errors in recovered HbT values. The CSF caused similar effects on µ a for 690 nm and 830 nm, which had minimal effects on recovered S t O 2 . This is consistent with expectations based on the discussion in Section 2.3. The first three source-detector distances (10-25 mm, 15-30 mm, 20-35 mm) showed a consistent pattern of underestimation of µ a due to excess CSF, with shorter source-detector distance being more sensitive. The longest source-detector distance (25-40 mm) showed an unusual pattern, in which increasing CSF initially caused recovered µ a values to increase before rapidly declining. We found this behavior be related to the specific anatomy of the head model, since this behavior was not reproducible on a regular geometry consisting of concentric spheres (data not shown). We also performed the 30%, 50%, and 90% S t O 2 , but the results were very similar and are not shown.

Discussion
In this study, we simulated FD-NIRS measurements to assess the magnitude of errors in recovered S t O 2 and HbT values using a semi-infinite, homogenous medium model of the head. We characterized the influence of curvature on errors via simulations on spheres of varying radius. We also simulated data on a segmented structural MRI of a full term neonate in order to estimate the magnitude of errors in a realistic model of a neonate with no morphological abnormalities. Lastly, we investigated the effects of extra-axial fluid by increasing the thickness of CSF layer surrounding the brain and repeating the simulations. The results presented here help to understand the influence of model violations on the accuracy of recovered S t O 2 , HbT, and µ a values and assess the validity of common analysis procedures for FD-NIRS data.

Effects of curvature
Curvature induced underestimation of µ a , µ ′ s , and HbT. The effects were more severe with decreasing radius of curvature. Curvature had little effect on S t O 2 estimates. There did not appear to be a strong trend with source-detector separation distance. For the neonate, head curvature can be a potentially significant source of error in recovering µ a and HbT. For older children and adults, curvature is expected to be less influential.

Effects of source-detector distance
Normalized partial pathlength showed an increasing pathlength through the brain with increasing source-detector separation. The pathlength through superficial layers showed a decreasing trend with source-detector separation. Together these results support the use of longer sourcedetector distances if signal to noise ratio (SNR) is sufficient. Simulations in the neonate model generally showed an overall trend of increasing accuracy with increasing source-detector distance.

Modulation frequency
The analysis methods described in Section 2.3 are based on the diffusion approximation, which is expected to break down at modulation frequencies on the order of ∼1 GHz [19]. In practice, phase wrapping occurs at ∼250 MHz at a distance of 40 mm for the optical properties used in this study. Our results showed only minor differences in the recovered values for modulation frequencies ranging from 50-200 MHz, suggesting that selection of source-detector distance is more important than selection of modulation frequency.