## 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 (

*μ′*), tissue oxygen saturation (S

_{s}_{t}O

_{2}), and total hemoglobin (HbT). Previous works have investigated the accuracy of recovered

*μ*values under this assumption. The purpose of this study was to examine the accuracy of recovered S

_{a}_{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}*μ′*, and HbT, but had minimal effects on S

_{s}_{t}O

_{2}. For the morphologically normal neonate head model, the mean absolute percent errors (MAPE) of recovered

*μ*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

_{a}_{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}*μ′*, and HbT, but had little effect on S

_{s}_{t}O

_{2}.

© 2014 Optical Society of America

## 1. 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 deoxy-hemoglobin (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

*μ′*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.

_{s}## 2. Methods

#### 2.1. 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.

#### 2.2. 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

*ψ*and Φ

_{theo}*are the theoretical intensity and phase values, and*

_{theo}*ψ*and Φ

_{slab}*are the intensity and phase measured on the calibration slab. The calibrated data is then given by where*

_{slab}*ψ*and Φ

_{cal}*are the calibrated intensity and phase values, and*

_{cal}*ψ*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].

#### 2.3. 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:

*ψ*and Φ are the amplitude and phase measurements with angular modulation frequency of

*ω*;

*ρ*is the source-detector separation;

*μ*is the absorption coefficient;

_{a}*D*is the diffusion constant equal to 1/(3

*μ*+ 3

_{a}*μ′*), in which

_{s}*μ′*is the reduced scattering coefficient;

_{s}*ν*is the speed of light in the media;

*F*is a complicated function of many variables, including

_{ψ}*ρ*,

*μ*, and

_{a}*D*that is approximately constant; Φ

_{0}is the instrument phase offset. The terms

*V*

^{+}and

*V*

^{−}are given by

*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

*μ*and

_{a}*μ′*is no longer possible.

_{s}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

*μ′*from regression of Eq. (5). The optical absorption coefficients recovered from Eq. (7a) above are related to hemoglobin by

_{s}_{2}] and [Hb] are the tissue concentrations of oxy-hemoglobin and deoxy-hemoglobin, respectively; ε

_{HbO2,λ}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

*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. The estimated values for [HbO

_{2}] and [Hb] can be used to calculate S

_{t}O

_{2}and HbT as follows:

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

*μ*values with coefficients from the matrix given by (

_{a}*E*)

^{T}E^{−1}

*E*. From Eq. (11), we can conclude that HbT is also a linear combination of recovered

^{T}*μ*values, and S

_{a}_{t}O

_{2}is a ratio with linear combinations of

*μ*values in the numerator and denominator. Thus, scaling all

_{a}*μ*values by a factor,

_{a}*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

*μ*are similar across wavelengths, the errors in HbT will be sensitive to the magnitude errors in

_{a}*μ*, but errors in S

_{a}_{t}O

_{2}will be minimal; however, in the case that the relative errors in

*μ*are dissimilar across wavelengths, the errors in S

_{a}_{t}O

_{2}can be large or unpredictable.

#### 2.4. 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.

#### 2.5. 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(a). The segmented volume had a 1 mm isotropic resolution. A single source position was simulated. Detector locations were chosen as random points on the surface of the head 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.

#### 2.6. 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.

## 3. Results

#### 3.1. Simulations on homogenous spheres

Figure 2 shows the results of simulating data on homogenous spheres and recovering *μ _{a}* (top row) and

*μ′*(bottom row) via the methods in Section 2.3. In all cases, curvature induced underestimation of

_{s}*μ*and

_{a}*μ′*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

_{s}*μ*. Figure 3 shows the results from calculating S

_{a}_{t}O

_{2}and HbT from the recovered

*μ*values in Figs. 2(a)–2(d). In general, curvature had little effect on S

_{a}_{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.

#### 3.2. Normalized partial pathlength

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 source-detector 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 results in Fig. 4 suggest that longer source-distances will be more sensitive to the brain and less sensitive to superficial layers.

#### 3.3. Simulations of varying S_{t}O_{2} and HbT in the brain

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

*μ′*(bottom row) with the slab model. There was an overall trend of increasing accuracy with increasing source-detector distance in recovering

_{s}*μ*. There was a relatively large degree of underestimation of

_{a}*μ*

_{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

*μ*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

_{a}_{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

*μ′*(bottom row) via the methods in Section 2.3. Recovered

_{s}*μ*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

_{a}_{t}O

_{2}and HbT from the recovered

*μ*values in Figs. 7(a)–7(d). There was a slight trend of increasing accuracy with source-detector distance in recovered S

_{a}_{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 9 shows selected simulations from Figs. 6 – 8 repeated at varying modulation frequency. Overall, modulation frequency had little effect on the recovery of *μ _{a}*,

*μ′*(not shown), S

_{s}_{t}O

_{2}, or HbT. There were no apparent trends in accuracy with regards to selection of modulation frequency within a range of 50–200 MHz.

#### 3.4. Effects of extra-axial fluid

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

*μ*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

_{a}_{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

*μ*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

_{a}_{t}O

_{2}, but the results were very similar and are not shown.

## 4. 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.

#### 4.1. Effects of curvature

Curvature induced underestimation of *μ _{a}*,

*μ′*, and HbT. The effects were more severe with decreasing radius of curvature. Curvature had little effect on S

_{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

*μ*and HbT. For older children and adults, curvature is expected to be less influential.

_{a}#### 4.2. 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 source-detector 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.

#### 4.3. 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.

#### 4.4. Effects of extra-axial fluid

The presence of extra-axial fluid can cause very large errors in *μ _{a}* and HbT, especially at short source-detector distances. Small amounts of extra-axial fluid do not severely affect S

_{t}O

_{2}estimates; however, when there is too much CSF, the methods described in Section 2.3 will break down resulting in negative

*μ*values, in which case estimates of S

_{a}_{t}O

_{2}are no longer feasible.

#### 4.5. Combined effects

From the results in Fig. 2, we expect curvature to influence underestimation of *μ _{a}*. Cerebrospinal fluid mostly influenced underestimation of

*μ*, although complicated geometries seemed to produce varying results in some cases. The scalp and skull can influence overestimation or underestimation depending on whether the absorption coefficient values are larger or smaller than that of the brain. This raises the possibility of opposing influences compensating for one another when the relative strengths of their influences are balanced. In general, the strengths of the influences of curvature, extra-cerebral tissues, and CSF will vary with physiological state (i.e., optical properties of the brain relative to other layers), source-detector separation distance, and thickness of extra-axial CSF. We can compare the absorption values of the extra-cerebral tissues, such as the skull or scalp, to the absorption value of the brain at a particular S

_{a}_{t}O

_{2}and HbT value to predict their influence on the errors; however, estimating the overall combined effect is not intuitive and requires systematic investigation, such as by numerical simulation.

#### 4.6. Propagation of *μ*_{a} errors to S_{t}O_{2} and HbT

Examination of Eqs. (9)–(11), as discussed in Section 2.3, gave some insight into the characteristics of errors in S_{t}O_{2} and HbT. Because S_{t}O_{2} is a ratio with *μ _{a}* values in both the numerator and dominator, it is robust to error influences that act similarly on all wavelengths, such as curvature and CSF. Conversely, the accuracy of HbT is related to the absolute accuracy of recovered

*μ*values regardless of similarity across wavelengths, since HbT is effectively a weighted sum of recovered

_{a}*μ*values. These effects are illustrated effectively in Figs. 2, 3, and 10.

_{a}#### 4.7. Typical errors in neonates

For neonates without morphological abnormalities, the mean absolute percent errors (MAPE) of recovered *μ _{a}* values were 12% and 7% for 690 nm and 830 nm, respectively. The MAPE for recovered S

_{t}O

_{2}and HbT were 6% and 9%, respectively. Larger relative errors were noted (∼20–30%), but were mostly observed at more extreme physiological values (i.e., very high/low S

_{t}O

_{2}or HbT). These results are comparable to the findings of Dehaes et. al. of 8–24% errors in recovered absorption coefficients in neonates [8]. While this level of error may be acceptable for some applications, multi-layered models may offer better accuracy when physiology significantly deviates from the normal values. For clinical populations with extra-axial fluid, advancement to more sophisticated models is necessary to obtain reliable measurements of absorption and hemoglobin.

#### 4.8. Limitations

The main limitation of this study is the reliance on a compilation of optical properties in Table 1. Furthermore, optical properties vary from person to person, such as with skin pigmentation. Since we simulated data over a wide range of physiological states (i.e., optical properties), the general results presented should represent the typical range of expected errors. Another minor limitation of this study was the use of semi-automated segmentation, which likely produced some misclassification of voxels; however, manual inspection ensured overall quality.

## References and links

**1. **F. F. Jöbsis, “Noninvasive, infrared monitoring of cerebral and myocardial oxygen sufficiency and circulatory parameters,” Science **198**, 1264–1267 (1977). [CrossRef] [PubMed]

**2. **M. A. Franceschini, S. Thaker, G. Themelis, K. K. Krishnamoorthy, H. Bortfeld, S. G. Diamond, D. A. Boas, K. Arvin, and P. E. Grant, “Assessment of infant brain development with frequency-domain near-infrared spectroscopy,” Pediatr. Res. **61**, 546–551 (2007). [CrossRef] [PubMed]

**3. **J. Wyatt, D. Delpy, M. Cope, S. Wray, and E. Reynolds, “Quantification of cerebral oxygenation and haemodynamics in sick newborn infants by near infrared spectrophotometry,” Lancet **328**, 1063–1066 (1986). [CrossRef]

**4. **J. H. Meek, C. E. Elwell, D. C. McCormick, A. D. Edwards, J. P. Townsend, A. L. Stewart, and J. S. Wyatt, “Abnormal cerebral haemodynamics in perinatally asphyxiated neonates related to outcome,” Arch. Dis. Child. Fetal. Neonatal. Ed. **81**, F110–F115 (1999). [CrossRef] [PubMed]

**5. **M. Ranucci, G. IsgrO, T. De La Torre, F. Romitti, D. Conti, and C. Carlucci, “Near-infrared spectroscopy correlates with continuous superior vena cava oxygen saturation in pediatric cardiac surgery patients,” Pediatr. Anesth. **18**, 1163–1169 (2008).

**6. **F. van Bel, P. Lemmers, and G. Naulaers, “Monitoring neonatal regional cerebral oxygen saturation in clinical practice: value and pitfalls,” Neonatology **94**, 237–244 (2008). [CrossRef] [PubMed]

**7. **S. Fantini, M. A. Franceschini, and E. Gratton, “Semi-infinite-geometry boundary problem for light migration in highly scattering media: a frequency-domain study in the diffusion approximation,” J. Opt. Soc. Am. B. **11**, 2128–2138 (1994). [CrossRef]

**8. **M. Dehaes, P. E. Grant, D. D. Sliva, N. Roche-Labarbe, R. Pienaar, D. A. Boas, M. A. Franceschini, and J. Selb, “Assessment of the frequency-domain multi-distance method to evaluate the brain optical properties: Monte carlo simulations from neonate to adult,” Biomed. Opt. Express **2**, 552–567 (2011). [CrossRef] [PubMed]

**9. **Q. Fang and D. A. Boas, “Monte carlo simulation of photon migration in 3d turbid media accelerated by graphics processing units,” Opt. Express **17**, 20178 (2009). [CrossRef] [PubMed]

**10. **G. Strangman, M. A. Franceschini, and D. A. Boas, “Factors affecting the accuracy of near-infrared spectroscopy concentration calculations for focal changes in oxygenation parameters,” Neuroimage **18**, 865–879 (2003). [CrossRef] [PubMed]

**11. **C. R. Simpson, M. Kohl, M. Essenpreis, and M. Cope, “Near-infrared optical properties of ex vivo human skin and subcutaneous tissues measured using the monte carlo inversion technique,” Phys. Med. Biol. **43**, 2465 (1998). [CrossRef] [PubMed]

**12. **S. L. Jacques, “Optical properties of biological tissues: a review,” Phys. Med. Biol. **58**, R37 (2013). [CrossRef] [PubMed]

**13. **M. Firbank, M. Hiraoka, M. Essenpreis, and D. Delpy, “Measurement of the optical properties of the skull in the wavelength range 650–950 nm,” Phys. Med. Biol. **38**, 503 (1993). [CrossRef] [PubMed]

**14. **S. A. Prahl, “Optical absorption of hemoglobin,” (1999), http://omlc.ogi.edu/spectra/hemoglobin/.

**15. **S. A. Prahl, “Optical absorption of water,” (2012), http://omlc.ogi.edu/spectra/water/.

**16. **S. Fantini, B. B. Barbieri, E. Gratton, M.-A. Franceschini, J. S. Maier, and S. A. Walker, “Frequency-domain multichannel optical detector for noninvasive tissue spectroscopy and oximetry,” Opt. Eng. **34**, 32–42 (1995). [CrossRef]

**17. **F. Shi, P.-T. Yap, G. Wu, H. Jia, J. H. Gilmore, W. Lin, and D. Shen, “Infant brain atlases from neonates to 1-and 2-year-olds,” PLOS ONE **6**, e18746 (2011). [CrossRef]

**18. **G. Nellhaus, “Head circumference from birth to eighteen years practical composite international and interracial graphs,” Pediatrics **41**, 106–114 (1968). [PubMed]

**19. **R. C. Haskell, L. O. Svaasand, T. T. Tsay, T. C. Feng, M. S. McAdams, and B. J. Tromberg, “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am. A Opt. Image Sci. Vis. **11**, 2727–2741 (1994). [CrossRef] [PubMed]