Quantitative assessment of hemodynamic and structural characteristics of in vivo brain tissue using total diffuse reflectance spectrum measured in a non-contact fashion

Here we present a new methodology that investigates the intrinsic structural and hemodynamic characteristics of in vivo brain tissue, in a non-contact fashion, and can be easily incorporated in an intra-operative environment. Within this methodology, relative total diffuse reflectance spectra (RTD(λ)) were acquired from targets using a hybrid spectroscopy imaging system. A spectral interpretation algorithm was subsequently applied to RTD(λ) to retrieve optical properties related to the compositional and structural characteristics of each target. Estimation errors of the proposed methodology were computationally evaluated using a Monte Carlo simulation model for photon migration under various conditions. It was discovered that this new methodology could handle moderate noise and achieve very high accuracy, but only if the refractive index of the target is known. The accuracy of the technique was also validated using a series of tissue phantom studies, and consistent and accurate estimates of μs’(λ)/μa(λ) were obtained from all the phantoms tested. Finally, a smallscale animal study was conducted to demonstrate the clinical utility of the reported method, wherein a forepaw stimulation model was utilized to induce transient hemodynamic responses in somatosensory cortices. With this approach, significant stimulation-related changes (p < 0.001) in cortical hemodynamic and structural characteristics were successfully measured. ©2016 Optical Society of America OCIS codes: (120.4640) Optical instruments; (170.3890) Medical optics instrumentation. References and links 1. W.-C. Lin, D. I. Sandberg, S. Bhatia, M. Johnson, S. Oh, and J. Ragheb, “Diffuse reflectance spectroscopy for in vivo pediatric brain tumor detection,” J. Biomed. Opt. 15(6), 061709 (2010). 2. B. J. Tromberg, N. Shah, R. Lanning, A. Cerussi, J. Espinoza, T. Pham, L. Svaasand, and J. Butler, “Noninvasive in vivo characterization of breast tumors using photon migration spectroscopy,” Neoplasia 2(1-2), 26– 40 (2000). 3. J. B. Fishkin, O. Coquoz, E. R. Anderson, M. Brenner, and B. J. Tromberg, “Frequency-domain photon migration measurements of normal and malignant tissue optical properties in a human subject,” Appl. Opt. 36(1), 10–20 (1997). 4. T. J. Farrell, M. S. Patterson, and B. Wilson, “A diffusion theory model of spatially resolved, steady-state diffuse reflectance for the noninvasive determination of tissue optical properties in vivo,” Med. Phys. 19(4), 879–888 (1992). 5. A. Amelink, H. J. Sterenborg, M. P. Bard, and S. A. Burgers, “In vivo measurement of the local optical properties of tissue by use of differential path-length spectroscopy,” Opt. Lett. 29(10), 1087–1089 (2004). 6. A. I. Kholodnykh, I. Y. Petrova, K. V. Larin, M. Motamedi, and R. O. Esenaliev, “Precision of measurement of tissue optical properties with optical coherence tomography,” Appl. Opt. 42(16), 3027–3037 (2003). 7. D. Levitz, L. Thrane, M. Frosz, P. Andersen, C. Andersen, S. Andersson-Engels, J. Valanciunaite, J. Swartling, and P. Hansen, “Determination of optical scattering properties of highly-scattering media in optical coherence tomography images,” Opt. Express 12(2), 249–259 (2004). 8. V. Turzhitsky, J. D. Rogers, N. N. Mutyal, H. K. Roy, and V. Backman, “Characterization of light transport in scattering media at sub-diffusion length scales with Low-coherence Enhanced Backscattering,” IEEE J. Sel. Top. Quantum Electron. 16(3), 619–626 (2010). Vol. 8, No. 1 | 1 Jan 2017 | BIOMEDICAL OPTICS EXPRESS 78


Introduction
Optical properties of biological tissues are always of great interest to researchers in the field of biomedical optics because of their intimate relationships with the intrinsic structural and compositional tissue characteristics. This is especially true for applications in investigating tissue pathology, tissue injury, and tissue functionality. Such applications in an intraoperative setting has drawn tremendous attention since it is relatively economic to set up in the operating room and efficient enough to provide real-time feedbacks [1]. Many techniques have been developed to measure optical properties from in vivo biological tissues over the past two decades [2][3][4][5][6][7][8][9][10][11][12][13]. Some of these techniques, such as optical coherence tomography [6,7] and low-coherence enhanced backscattering spectroscopy [8,9], are limited to quantify scattering-related properties. Among them, spatially and temporally resolved diffuse reflectance signals are the most commonly used to quantify both absorption and scattering properties; they can be conveniently acquired using a contact fiber optic bundle [10][11][12][13]. However, a significant drawback is associated with the utility of a contact fiber optic probe: the contact pressure introduced by the probe can significantly alter the measured diffuse reflectance spectra and hence optical properties in vivo [14,15]. This, in turn, undermines the accuracy of in vivo tissue differentiation applications based on diffuse reflectance spectroscopy, such as intraoperative brain tumor demarcation. Various probe holders have been devised to reduce hand-motion induced artifacts [15][16][17], but they are rather inconvenient to operate in an intraoperative setting.
An alternative approach to address this drawback is to acquire desired optical signals in a non-contact fashion. Over the past decade, several non-contact techniques for in vivo optical property measurements have been proposed and developed: Cuccia et al. introduced a noncontact, wide-field measurement technique of optical properties using spatially modulated illumination and achieved excellent accuracy [18]; Bish et al. used a lens system to image the illumination and collection fibers onto the tissue surface and hence achieve non-contact measurements of radially dependent diffuse reflectance [19]; Foschum et al. introduced two apparatus to carry out non-contact measurements of spatially resolved reflectance and total reflectance from that absorption and reduced scattering coefficients individually estimated [20]. However, these techniques require either a sophisticated optical system to spatially modulate the illumination or a carefully controlled collection geometry, which make them less suitable for intraoperative applications such as tumor resection guidance. Therefore, an imaging system, simple in instrumentation but accurate in estimating optical properties of biological tissue, for intraoperative applications remains to be developed.
Here we introduce a non-invasive optical technique, specifically designed for intraoperative tissue diagnosis and easily integrated into the typical surgical microscope found in most operating rooms, which allows for non-intrusive optical investigation during surgery. The methodology has both a hardware and software component. The hardware is a hybrid spectroscopy imaging system that was devised for the purpose of acquiring relative total diffuse reflectance spectra R TD (λ) from a given point within the system's field of view [21]. The system is also capable of obtaining two-dimensional relative total diffuse reflectance signals at one or two wavelengths specified by the user. By incorporating a reference measurement and a reference table established using a Monte Carlo (MC) simulation model for photon migration, R TD (λ) can be converted to absolute total diffuse reflectance signals and, subsequently, a reduced scattering coefficient to absorption coefficient ratio: μ s '(λ)/μ a (λ).
The methodology's software is a new spectral interpretation algorithm that was developed for μ s '(λ)/μ a (λ) to quantify the hemodynamic (hemoglobin concentration, [Hb], oxygen saturation level, SatO 2 ) and structural (scattering amplitude and scattering power) characteristics of biological tissues, with only a few prerequisites. In a series of experiments, the accuracy of the algorithm -in terms of measuring μ s '(λ)/μ a (λ) and quantifying the hemodynamic and structural properties -was first theoretically evaluated using total diffuse reflectance signals generated by the MC model. Next, the overall accuracy of the methodology was experimentally demonstrated using optical tissue phantoms. Finally, the clinical utility of the methodology was tested in an in vivo rat study, in which a forepaw stimulation model was utilized to induce hemodynamic responses in the somatosensory cortex of Wistar rats, and these responses were then quantified.

Instrumentation
Design details for the hybrid spectroscopy imaging system are depicted in Fig. 1. The system consisted of two signal acquisition modalities-the imaging modality and the point spectroscopic detection modality. The targeting sample was imaged through either a Nikon dSLR lens (Nikon AF 28-80 mm f/3.5-5.6 D lens with aperture ring) or a zoom imaging lens (VZM 450i, Edmund Optics Inc., Barrington, NJ). The image from the camera lens was collimated using a Hastings triplet achromatic lens (#30-229, EFL 40.3 mm, Edmund Optics Inc., Barrington, NJ) and transmitted to a beam splitter (#54-824, 50R/50T, Edmund Optics Inc., Barrington, NJ). The collimated image beam transmitting through the beam splitter would enter into the imaging modality of the system; it was further divided by a dichroic mirror with a transmission band of 400-595 nm and a reflection band of 640-750 nm (#49-471, Edmund Optics Inc., Barrington, NJ). Two narrow band-pass filters with center passing wavelengths at 500 nm (#65-149, Edmund Optics Inc., Barrington, NJ) and 700 nm (#88-012, Edmund Optics Inc., Barrington, NJ) were placed at the transmission side and the reflection side of the dichroic mirror, respectively. The collimated image beam transmitting through the band-pass filter was re-focused onto the sensor of a CCD camera (DMK 21AU04, The Imaging Source Europe GmbH) using a Hastings lens. The collimated image beam reflected by the beam splitter would enter into the point spectroscopic detection modality and be refocused using a Hastings lens. At the image plane of this lens, an optical fiber (core diameter 50 μm, NA = 0.2, GIF50C, Thorlabs Inc., Newton, NJ) was used to collect light from the center of the formed image. The distal end of the optical fiber was connected to either a red laser diode module (CPS184, 650 nm, 4.5 mW, Thorlabs Inc., Newton, NJ) assembled inhouse to track the location of investigation on the targeting sample, or a spectrometer (USB2000, Ocean Optics, Dunedin, FL) to record light reflected by the sample. The area of investigation of the point spectroscopic detection modality was determined by the magnification of the lens as well as the diameter of the optical fiber. The hybrid system was operated through a couple of LabVIEW programs developed in-house. The illumination of the sample investigated was provided by an external light source. Assuming the illumination is relatively uniform, the light received by a pixel of the camera or the point spectroscopic detection system should be considered as a fraction (f) of the total diffuse reflectance (R TD ) generated by a pencil beam illumination, in accordance with the principle of convolution [22]. The proof for this concept is provided in the Appendix A.

μ s '/μ a estimation
According to the theory proposed by Farrell et al., R TD is only a function of μ s '/μ a [4]. To easily convert R TD measured by the hybrid system to μ s '/μ a , a look-up table was constructed using a MC simulation model for photon migration in multi-layered tissue. The simulation model was developed and verified using MCML provided by Wang et al. [23]. The magnitudes of μ s '/μ a used to create the look-up table are depicted in Table 1, with a randomly selected value of μ a (2 cm −1 ). In the simulations, the tissue slab was illuminated by a pencil beam at a normal angle. The refractive index n, the anisotropy factor g, and the thickness of the simulated tissue slab were 1.395, 0.88, and 100 cm, respectively. The total number of photons used in each simulation was 10 million. The resulting look-up table is shown in Fig.  2.

Spectral interpretation
With advanced analysis on μ s '/μ a spectra, function-and structure-related information from biological tissues could be extracted. Assuming hemoglobin is the major absorber within the visible wavelength range, the absorption coefficient can be expressed as where [Hb] is the total hemoglobin concentration, ε oxy (λ) and ε deoxy (λ) are the wavelengthdependent molar extinction coefficients of oxy-and deoxy-hemoglobin, respectively, SatO 2 is the oxygen saturation level in percentage, MW is the gram molecular weight of hemoglobin (64500 g/mol) and ln(10) is the conversion factor [24]. Furthermore, a simple power law expression may be used to describe a reduced scattering spectrum μ s '(λ). That is where λ is in nm, λ 500 is 500 nm, A is a scaling factor which is equal to μ s '(500 nm) and B is the scattering power [25].
At isosbestic points of oxy-and deoxy-hemoglobin absorption spectra (i.e., 500, 529, 545, 570, and 584 nm, approximately), their molar extinction coefficients are equal, as a result of which SatO 2 will not have any influence on μ a . Hence, the least-square method could be employed to estimate A/[Hb] and B from Eq. (3) with logarithm transformation. With A/[Hb] and B, the hemoglobin absorption spectrum C(λ) at all other non-isosbestic wavelengths could be reconstructed using: 2 For any level of SatO 2 , the hemoglobin absorption spectrum C(λ) has a unique spectral profile (i.e., the number of peaks and the peak locations), which could be formulated into another look-up table specifically for SatO 2 estimation. Detailed formulations of this algorithm will be provided in Appendix B.

Validation based on Monte Carlo simulation
The Monte Carlo simulation-based validation was conducted to evaluate the effects of optical properties, other than absorption and scattering coefficients, and of random noise on the estimation accuracy of μ s '/μ a and the subsequent spectral analysis. A Matlab program was developed to generate 50 random μ s -μ a pairs for each of the six μ s '/μ a ranges depicted in Table 1. These coefficients, in conjunction with n = 1.395 and g = 0.88, were applied to the MC simulation model to generate R TD , assuming the simulated tissue phantom is infinitely thick and illuminated by a pencil beam at a normal angle. Each simulated R TD was converted to μ s '/μ a using the look-up table shown in Fig. 2. Each estimated μ s '/μ a value was subsequently compared with the corresponding input μ s '/μ a value to determine the absolute percentage error of estimation |ΔE|. That is

Look-up table for μ s '/μ a estimation
where v 1 and v 2 are the input and the estimated μ s '/μ a values, respectively. In practice, absolute R TD without noise is impossible to obtain. The effects of random noise in R TD on the μ s '/μ a estimation process, therefore, were investigated in this study. Different levels of random noise were created using a Matlab function "rand" and then added to each R TD-generated above by where the "Level" is the noise level (i.e., 1%, 5% and 10%), and the "rand" is a random number ranging from 0 to 1 and has a uniform probability density function. This is to achieve three different signal-to-noise ratios: 100, 20, and 10. The noisy R TD were converted to μ s '/μ a using the look-up table depicted in Fig. 2; the corresponding |ΔE| were calculated using Eq. .
In addition to the random noise, the effects of the refractive index n and anisotropy factor g were also evaluated with MC simulation. The parameters for setting up the MC simulation were listed in the Table 2. Each condition was randomly simulated 50 times. The simulated R TD was used as an input in the look-up table of μ s '/μ a for estimation. The discrepancies between the estimated and the original μ s '/μ a were again quantified using Eq. (5), with v 1 being the original μ s '/μ a and v 2 the estimated one. Detailed descriptions about the simulation and the following comparisons are provided in Appendix C.

μ s '/μ a spectral interpretation algorithm
In order to evaluate the efficacy of the algorithm in estimating A/[Hb], B and SatO 2 directly from a total diffuse reflectance spectrum, a validation study based on Monte Carlo simulation was designed. A Matlab program using default function "rand", which generates uniformly distributed random numbers, was developed to produce 50 sets of random values for the parameters A, B, [Hb], and SatO 2 , within the ranges depicted in Table 1, and then to compute their corresponding optical properties μ a (λ) and μ s '(λ) using Eqs. (1) and (2). These 50 sets of optical properties, along with refractive index n = 1.395 and anisotropy factor g = 0.88, were applied to a Monte Carlo (MC) simulation model for photon migration to generate total diffuse reflectance (R TD ), assuming the simulated tissue phantom was infinitely thick (~100 cm) and illuminated by a pencil beam perpendicularly. All optical properties used (Table 3) were selected in accordance to those of human brain tissue summarized by Jacques et al. [25]. The wavelength range covered in the MC simulation was also listed in Table 3, with an increment of 2 nm. The simulated R TD spectra between 450 nm and 600 nm were applied in the spectral interpretation algorithm written in Matlab to estimate the corresponding values of A/[Hb], B and SatO 2 using the interpretation methods described above. The estimated A/[Hb], B and SatO 2 subsequently compared with the original values used in the MC simulation in terms of the absolute percentage error |ΔE| as Eq. (5). In addition, for the purpose of evaluating the effects of noise on the spectra interpretation algorithm, different levels (0.1%, 0.5%, 1% and 5%) of uniformly distributed random noise were added to the simulated R TD spectra by Eq. (6). Parameters A/[Hb], B and SatO 2 were also estimated based on the noisy R TD spectra, and compared with pre-set values using Eq. (5).
In practice, the influence of random noise is usually reduced by averaging multiple measurements. To demonstrate the importance of averaging repeated measurements for parameter estimation, random noise at the same level (0.1%, 0.5%, 1% or 5%) was generated and added to the simulated R TD spectra, according to Eq. (6). The noise addition process was repeated 10 times, which yielded 10 noisy R TD (λ). The 10 noisy R TD (λ) at a given noise level were averaged to yield an averaged spectrum for spectral interpretation. The resulting parameter estimations were compared with the ones from the noisy spectra without averaging in terms of |ΔE|.

Phantom study
One of the objectives of this study was to investigate the feasibility and the accuracy of the point spectroscopic detection modality (see Fig. 1) of the hybrid system in determining μ s '/μ a from absorbing/scattering media by measuring the absolute R TD . In order to measure absolute R TD , a priori knowledge of the incident light intensity as well as the collection geometry of the hybrid system would be required, as suggested by Foschum et al. [20]. Unfortunately, these two pieces of information were difficult to obtain accurately. To circumvent this limitation, reference R TD measurements from a material with known optical properties was introduced to the process of obtaining absolute R TD and hence estimating μ s '/μ a .
Ten cylindrical polyurethane optical phantoms with different optical properties were used in this study and were made in accordance with the recipe developed by Prahl et al. [26]. The optical properties (i.e., μ s ' and μ a ) of these phantoms within the visible wavelength range were determined using the double-integrating-sphere method [27], and their ranges between 500 nm to 640 nm are shown in Table 4. Note that the ranges of μ s '/μ a in Table 4 are consistent with those in Table 3 used for simulation purpose. The refractive index of the cured polyurethane phantom is around 1.468 at 670 nm [28]. Since the optical properties of these phantoms are known, their theoretical absolute total diffuse reflectance spectra R TD (λ) could be estimated using a MC simulation. However, the challenge of using an optical phantom as a reference lies in in vivo implementation: it would be difficult to maintain the same illumination and collection geometries (for example, surface curvature) between the phantom and the in vivo tissue of investigations. One way to address the in vivo applicability of a reference phantom is to create a reflective material that is thin and flexible enough to be placed directly on top of the in vivo tissue for reference measurements. A double-layer paper-based standard was therefore created and evaluated in this study. The top layer of this standard was a piece of reflective white paper; and the bottom layer a piece of black paper absorbing all light transmitting through the top layer. The paper standard has to be pre-calibrated with a reference phantom before any new studies, following the steps depicted in Fig. 3. The reference phantom used in this study is O38. With the paper standard calibrated, it will be very easy to obtain the location-specific fraction (f) value to convert the measured reflectance spectrum into the total diffuse reflectance spectrum R TD (λ): where R paper (λ) is measured from the paper standard when it is placed on top of the studied medium (unknown phantoms or in vivo tissues), and R TD_paper (λ) is estimated during the precalibration procedure. In the phantom study, the optical properties μ s '/μ a of those phantoms (excluding O38) were assumed to be unknown and estimated using the procedures shown in Fig. 3. Each phantom was illuminated directly from above with a broad-band white light source (DC-950H, Dolan-Jenner, Edmund Optics Inc., Barrington, NJ); reflectance spectrum R unknown (λ), which was supposed to be a fraction (f) of the R TD_unknown (λ), was measured through the point spectroscopic detection modality of the hybrid system from five random locations on each "unknown" phantom at an arbitrary observation angle between 30 and 45 degrees. The side of each phantom was covered by a piece of thick black paper to ensure that only its top surface was illuminated during the study. All the measured reflectance spectra were subtracted by the dark measurement taken with the illumination light off. Therefore, the estimated total diffuse reflectance spectrum of the "unknown" phantom could be expressed as Two comparisons were made to assess the accuracy of the spectroscopic point-detection modality of the hybrid system in measuring μ s '/μ a from absorbing/scattering media using a paper-based standard. Specifically, R TD_unknown (λ) was compared with theoretical R TD (λ) derived from the look-up table (Fig. 2) using μ s '(λ)/μ a (λ) of the evaluation phantom, predetermined by the double-integrating-sphere technique, as the input. In addition, R TD_unknown (λ) was converted to μ s '(λ)/μ a (λ) using the look-up table (Fig. 2) and then compared with μ s '(λ)/μ a (λ) of the evaluation phantom (Table 4).

Forepaw stimulation in Wistar rats
Forepaw stimulation in rodents is often used to induce significant hemodynamic responses [29], including changes in cerebral blood flow, cerebral blood volume and oxygen partial pressure, within the somatosensory cortex for the forelimb (S1FL), for the purpose of investigating the neurovascular coupling mechanism. It offered a tremendous opportunity for us to verify whether the proposed new measurement scheme and accompanying spectral interpretation algorithm are capable of both measuring resting-state (baseline) hemodynamic and structural features, and detecting any stimulation-induced hemodynamic variations within the S1FL cortex.
The animal study protocols were approved by the Institutional Animal Care and Use Committee (IACUC) at Florida International University (Approval Number: 13-065) and conducted in full compliance with federal, state and local regulations and laws. All procedures were performed with anesthetized rats to minimize their suffering. Seven normal Wistar rats were used in this study ( Table 5). Three of them were used solely as the controls to measure resting-state (baseline) activity. The other four rats were used to measure both resting-state and event-related (i.e., forepaw stimulation) activity. Each rat's body temperature (~37°C), heart rate (200-300 beats per minute) and respiration rate (<50 breath per minute) were monitored continuously, using a PowerLab 8/35 data acquisition device and LabChart software (AD Instruments), to ensure that the level of anesthesia remained stable throughout the entire surgical procedure and recordings. Rats were fixed on the stereotaxic stage with their skull and dura overlying the somatosensory cortex for forelimb (S1FL) removed during the infusion of 2% isoflurane with 1L/min oxygen (14 psi). The location of the S1FL was determined based on a rat brain atlas [30]. Non-conductive paraffin oil drops (O121-1, Mineral Oil, Light, Fisher Scientific) were applied to the craniotomy to prevent the cortical surface from becoming dehydrated during the recordings. The contralateral forepaws of rats in the experimental group were then stimulated using bipolar needle-electrodes placed subcutaneously. Prior to forepaw stimulation, 0.25-mg/kg of the anesthetic/sedative Dexdomitor (active ingredient = 0.1mg/ml dexmedetomidine hydrochloride) was injected intraperitoneally to sedate the rats. Subsequently, the isoflurane was reduced to 0.25% (1L/min O 2 , 14 psi). Stimulation pulses (10-ms pulse width) were delivered by means of an isolated pulse stimulator (Model 2100, 110V, 60 Hz, A-M Systems; ~2-mA amplitude at 3 Hz). A repetitive block-design paradigm was employed, with each block consisting of an ON period of 16 seconds, during which stimulation pulses were triggered, and an OFF period (i.e., resting state) of 44 seconds.
A paper-based reference was pre-calibrated with an optical phantom prior to the animal experiments, as described in Section 2.5. During all animal experiments, the paper reference was placed on top of the craniotomy to obtain the averaged R paper (λ) spectra based on 10second continuous spectral acquisition, with an integration time of 250 ms, before and after forepaw stimulation. This secondary calibration, which will give the location-specific fraction f as per Eq. (7), is very crucial for the in vivo studies since the tissue surface is not always flat and could have significantly different illumination-collection geometry as in the precalibration step. Reflectance spectra from the S1FL cortex R cortex (λ) were acquired continuously with an integration time of 250 ms during forepaw stimulation. A dark measurement was also attained and subtracted from R cortex (λ). The time-dependent absolute total diffuse reflectance spectra of the S1FL cortex R TD_cortex was calculated as: The corresponding values for μ s '/μ a of the S1FL cortex could be derived from R TD_cortex spectra instantaneously, using the reference table for R TD versus μ s '/μ a . The event-related parameters A/[Hb], B and SatO 2 were subsequently estimated using the proposed spectral interpretation algorithm. The time-series of changes in these parameters (i.e., Δx(t)) was calculated as: were compared in old versus young rats using an unbalanced two-sample t-test.

Accuracy of μ s '/μ a estimation
The effects of varying g on the relationship between μ s '/μ a and R TD are summarized in Fig. 4. It is noticed that, within the range of investigation, an increase in g variation does elevate the magnitude of |ΔE|. Meanwhile, the maximum |ΔE| in R TD estimation using μ s '/μ a could exceed 6% when μ s '/μ a was smaller than 1 ( Fig. 4(b)), |ΔE| in μ s '/μ a estimations based on R TD is reasonable small (< 2%) over the entire range of R TD (Fig. 4(a)). This suggests that the look-up table for μ s '/μ a verse R TD established with g = 0.88 (Fig. 2) is applicable to biological tissues with different g values. The effects of varying n on the relationship between μ s '/μ a and R TD are summarized in Fig. 5. Similar to g variation, an increase in n variation also raises the magnitude of |ΔE|. In addition, |ΔE| in μ s '/μ a estimations using R TD steadily increases as R TD approaches to 1. |ΔE| in R TD estimations using μ s '/μ a , on the other hand, decreases as μ s '/μ a increases. It is noticed that |ΔE| induced by n variation is much greater than that induced by g variation; the maximum |ΔE| exceeds 40% in μ s '/μ a estimations using R TD (Fig. 5(a)) and 23% in R TD estimations using μ s '/μ a (Fig. 5(b)). This clearly indicates the look-up table depicted in Fig. 2 is not applicable to biological tissues whose n greatly deviates (i.e., more than 5%) from 1.395. The effects of varying n and g on μ s '/μ a estimations based on R TD are collectively summarized in Fig. 6. As expected, the accuracy of the μ s '/μ a estimation is insensitive to g variation when n is a constant; the plot of input and the estimated μ s '/μ a under such a condition shows a highly linear relationship (R 2 ≈1) with a slope = 1. However, the accuracy of the μ s '/μ a estimation suffers greatly when n varies.  (Table 2 and Appendix C.3), over the range of μ s '/μ a depicted in Table 1. Fig. 7. Absolute percentage error |ΔE| in μ s '/μ a estimations when n and g are fixed. Three different percentage levels of noise, 1%, 5%, and 10%, were added to the total diffuse reflectance R TD , which yielded signal-to-noise ratios of 100, 20, and 10, respectively.
The addition of random noise does have negative effects on the estimation of μ s '/μ a using R TD (Fig. 7). As expected, the higher the noise level the higher the |ΔE|. When μ s '/μ a is greater than 1, a common logarithmic (base 10) relationship between |ΔE| and μ s '/μ a is observed. The large estimation error |ΔE| in the high μ s '/μ a region is attributed to the fact that R TD is rather insensitive to μ s '/μ a variation in this region. However, based on the possible ranges of μ s '/μ a (3.70-52.16) of human brain tissues calculated in Table 3, the estimation error of μ s '/μ a should be within the acceptable ranges.

Accuracy of spectral interpretation
The simulated R TD spectra with or without noise were processed using the R TD spectral interpretation algorithm to estimate the indicative parameters A/[Hb], B and SatO 2 . The effect of the averaging process on the performance of the R TD spectral interpretation algorithm is shown in Fig. 8(a). When the averaged noisy R TD spectra were used as the input, |ΔE| of each indicative parameter was markedly reduced (Fig. 8(a) in red). This effect is especially prominent when the noise level is high. However, |ΔE| of SatO 2 estimation at the noise level of 5% was still very high (~23%) even when the averaged noisy spectra were used. It should be noted that the noise effect is not a constant for the entire range of SatO 2 ; the |ΔE| was within an acceptable range (<5%) when the theoretical value of SatO 2 was over 30%, as shown in Fig. 8(b).

Phantom study
The results of the experimental evaluation using optical phantoms demonstrate the high accuracy of the point spectroscopic detection modality of the hybrid system in terms of measuring R TD and estimating μ s '/μ a . The paper-based standard was pre-calibrated with phantom O38 by Eq. (7), since O38 has the brightest reflectance and hence highest signal-tonoise ratio in the measurements. The evaluation procedure in Fig. 3 was performed on all the remaining phantoms. The averaged |ΔE| of μ s '/μ a estimation (data size at each wavelength = 9 phantoms × 5 locations × 10 measurements = 450) within the wavelength range of 500-640 nm is 5.54 ± 4.6% (mean ± standard deviation), as shown in Fig. 9(a). When the measured μ s '(λ)/μ a (λ) is plotted against the theoretical μ s '(λ)/μ a (λ), a strong linear relationship with a slope of close to one is observed for the entire evaluated wavelength range (141 wavelengths, spectral resolution = 1 nm). A representative plot of such a relationship is provided in Fig.  9(b), where R 2 of the linear regression is 0.9924. The proposed instrumentation and look-up table for estimating μ s '/μ a have very similar results to those obtained by the doubleintegrating-sphere technique [27].
One interesting observation from this set of data is that the phantom itself is not a perfect homogenous medium. In several optical phantoms, a rather large variation is observed among R unknown (λ) acquired from five random locations on the top surface of a given phantom. Therefore, only the mean of the five acquired R unknown (λ) from a phantom was used to derive R TD_unknown (λ), which was then compared with theoretical values determined by simulation based on optical properties estimated by IAD method. Fig. 9. Representative results from the experimental evaluation study of the point spectroscopic detection modality of the hybrid system on 10 phantoms. The spectral acquisition procedure is depicted Fig. 3; a total of five investigated sites were used for each phantom. (a) Absolute percentage error |ΔE| in μ s '(λ)/μ a (λ) estimations using R TD _ unknown (λ) between 500 nm and 640 nm. The reference μ s '(λ)/μ a (λ) used here are those determined by the double-integrating-sphere technique. The middle line represents the mean |ΔE| calculated from all the phantoms studied; the error bars the standard deviation. (b) Measured μ s '(λ)/μ a (λ) were plotted against reference μ s '(λ)/μ a (λ) between 500 nm and 640 nm from all the phantoms studied (black dots). Area highlighted in yellow shows the range of μ s '/μ a (~3-60) that is relevant to human brain. Red crosses represent μ s '(λ)/μ a (λ) of the reference phantom O38. The hemodynamic baselines and event-related hemodynamic responses induced by forepaw stimulation in all rats were measured using the hybrid spectroscopy imaging system (Fig. 10), and then estimated by the spectral interpretation algorithm and Eq. (10). The area under investigation (AUI, i.e., S1FL cortex) was confirmed using the aiming function of the hybrid system (the bright spot in Fig. 10(b)). The resting-state (or baseline) values of all three indicative parameters from the S1FL cortices seem to be age-related (Fig. 11), which might indicate that an age-related difference existed in the resting-state tissue oxygenation and metabolism of the cortex of the rat's brain. Significant variations (p < 0.001) in all three indicative parameters, i.e., A/[Hb], B, and SatO 2 , within the AUI were observed during stimulation of the forepaw in all of the rats that were studied (Fig. 12(b)). It was observed that the amplitudes of the variations in two of the indicative parameters (B and SatO 2 ) were much higher in the younger rats (Rats 4 and 5) than in the older rats (Rats 1 and 2) (Fig. 12(a)).  The paired differences between the stimulated responses and baseline values of each stimulation event were significantly different from 0 (p<0.001, marked as "***"), which indicated significant responses induced by forepaw stimulation were detected by the system. It seems that there is a correlation between two indicative parameters (i.e., B and SatO 2 ) and the age/weight of the rats.

Discussion
In this paper, we have introduced a novel non-contact optical methodology that (1) acquires relative total diffuse reflectance spectra from the target, (2) converts the spectra to μ s '(λ)/μ a (λ) using a reference table, and (3) retrieves indicative physiological parameters (A/[Hb], B and SatO 2 ) from μ s '(λ)/μ a (λ). The accuracy of this technique is demonstrated using both a theoretical and experimental approach, and its clinical utility confirmed in an in vivo rat study.

R TD and μ s '/μ a
The need to accurately obtain optical properties from in vivo biological tissues is largely attributed to the fact that these properties are directly linked to tissue compositional and structural characteristics, and hence provide a diagnostic potential for diseases and injuries. In this article, a unique hybrid spectroscopy imaging system was introduced to quantify optical properties of in vivo tissue. The system possesses two detection modalities: one for point spectroscopic detection and the other for imaging. The former is used when high spectral resolution is needed, while the latter is used when investigation of an area is necessary. Both detection modalities acquire relative total diffuse reflectance signals from the target in a noncontact fashion. With assistance from a reference optical standard and a look-up table generated by the MC simulation model for photon migration, relative total diffuse reflectance measurements can be converted to absolute total diffuse reflectance and then μ s '/μ a . The results of the evaluation confirm the operating principle, as well as the high accuracy of the hybrid system.
To obtain reliable results from the hybrid spectroscopy imaging system, the following three conditions should be met: 1) the illumination intensity on the investigated subject needs to be uniform; 2) the optical properties of the investigated subject are homogeneous; 3) the anisotropy factor g and, more importantly, the refractive index n of the investigated subject should be known.
The first required condition is derived from the theoretical framework of total diffuse reflectance R TD obtained from a broad beam illumination (See Appendix A): the illumination intensity at a given point on the surface of the investigated subject should be the same as that at the point of detection in order to obtain accurate R TD . This is especially true for the illumination points that make a meaningful contribution to the diffuse reflectance signal measured at (x 0 ,y 0 ). This requirement, in turn, constitutes the minimal radius of uniform illumination of the system, which can be quantitatively assessed using the point spread function of photon migration (i.e., radially resolved diffuse reflectance from a pencil beam illumination). According to the theory of photon migration, low μ s ' and μ a generally leads to a broader point spread function. Therefore, the minimal radius of uniform illumination increases as μ s ' and μ a decrease. In order to quantify these effects, reported optical properties of pediatric brain cortex [25] were used to simulate point spread functions at 500 nm, 600 nm, and 700 nm, and the simulation results are presented in Fig. 13. It turns out that within the visible range, the radius within which 95% of total diffuse reflectance is observed (r 95 ) is 2.6 mm and 5.8 mm for the smallest and the largest point spread function at 500 nm, respectively. In other words, the minimal radius of uniform illumination at 500 nm has to exceed 5.8 mm in order to reduce the measurement error of total diffuse reflectance below 5%. This is rather easy to attain with a commercial-grade surgical light and a surgical microscope light in the operating room. As expected, r 95 increases significantly as the wavelength increases (i.e., μ s ' and μ a decrease). According to the simulation results, the largest r 95 is found at 700 nm and is about 2.9 cm. This number may not be reliable as the optical properties used may have been strongly underestimated. Nevertheless this indicates the hybrid system works better in the shorter wavelength region where the μ s ' and μ a of biological tissues are high and the corresponding r 95 , as well as the minimal radius of uniform illumination, are low. Fig. 13. Diffuse reflectance (R d ) as a function of radius (r) at three different wavelengths (a-500 nm, b-600 nm, and c-700 nm). The red curve demonstrated the narrowest diffuse reflectance profile, while the blue one is the widest one. The dash lines represent the cutoffs within which 95% of total diffuse reflectance is observed. To generate these R d (r), the Monte Carlo simulation was carried out with the optical properties of human brain [25]. Specifically, n is 1.395, g is 0.88, µ a is 1.6-4.5 cm −1 at 500 nm, 0.6-1.7 cm −1 at 600 nm and 0.07-0.19 cm −1 at 700 nm, and µ s = 83-250 cm −1 at 500 nm, 52-228 cm −1 at 600 nm and 35-211 cm −1 at 700 nm.
The second required operating condition of the hybrid system is based on the fact that the look-up table shown in Fig. 2 is derived using MC simulations with spatially homogenous media (i.e., constant optical properties). This assumption is used widely in techniques measuring tissue optical properties. For the hybrid system, the minimal radius of uniform illumination discussed in the previous section defines the minimal radius of tissue homogeneity required to ensure accurate operation of the system. To determine the minimal depth of tissue homogeneity, again MC simulation is used in conjunction with the optical properties of the pediatric brain cortex [25]. The results of the simulation indicate that the minimal depth of tissue homogeneity would be approximately 1.475 mm for 500 nm and 7.885 mm for 700 nm. Since the hybrid system reports the average μ s '/μ a within the volume of investigation defined by the minimal radius and depth of tissue homogeneity, the inhomogeneity will inevitably deviate the measured μ s '/μ a from its expected value.
The simulation results presented in this article suggest that anisotropy factor g is not as influential as refractive index n in terms of the accuracy of μ s '/μ a estimation. However, n did not introduce a significant error in estimating μ s '/μ a of the evaluation phantoms (n = 1.486) based on a look-up table designed for a different n (n = 1.395). This is attributed to the reference phantom, which has the same refractive index n as the evaluation phantoms, and thus, essentially the accuracy of the relative comparison between samples with same n would not depend on any particular look-up table. Therefore, n of the investigated tissue and the reference material are required priori knowledge for the accurate operation of the hybrid system. Refractive indices of biological tissues have been reported widely in the existing literature. For example, Jacques published a very comprehensive review on the optical properties of biological tissues [25], which could be a good reference for determining n of a specific tissue type, and as such, for creating a corresponding look-up table. If n is not readily available, it may be estimated based on its water content using the formula depicted by Jacques [25].
The hybrid system disclosed in this article is capable of obtaining absolute total diffuse reflectance at macroscopic levels and efficiently estimating a wide range of μ s '/μ a that cover the reported ranges for brain tissues. It is relatively inexpensive to construct and can be applied to various biomedical applications such as intraoperative brain tumor detection. The system may be easily used in the surgical operating field, and even could be integrated with a surgical microscope, without introducing any safety issues or disturbing normal clinical procedures. One significant limitation of the hybrid system is that it does not provide μ s ' and μ a separately. This limitation, however, may be circumvented by performing advanced analysis on μ s '/μ a spectra to extract more function-and structure-related information from biological tissues, such as the indicative parameters A/[Hb], B and SatO 2 discussed in this study.

Spectral interpretation algorithm
In this article, a new spectral interpretation algorithm was proposed to retrieve three indicative parameters related to the hemodynamic and structural characteristics of biological tissue based on total diffuse reflectance spectrum between 500 to 580 nm. Currently, the most popular spectral interpretation algorithm applied on multi-wavelength optical data to estimate the hemodynamics variations is based on modified Beer-Lambert law [31][32][33]. The scattering properties are usually assumed to be constant in order to obtain the desired hemodynamic information. A priori knowledge about the baseline hemoglobin concentration and oxygen saturation level are also needed to determine the mean optical path length in the modified Beer-Lambert law. These assumptions would introduce bias to quantify the real hemodynamic characteristics. Furthermore, these algorithms are not capable of producing absolute SatO 2 . In contrast, the spectral interpretation algorithm proposed in this study does not rely on any assumption about the baseline values of each indicative parameter. In addition, it does not assume the scattering properties are constant. The most important advantage of this algorithm is that it could objectively provide the absolute value of SatO 2 , instead of a relative estimation based on pre-set values.
Another advantage of the proposed measurement scheme and spectral interpretation algorithm is that it utilizes optical data within the short wavelength region (i.e., 500-580 nm), which leads to a small volume of investigation (VOI). This, in turn, reduces the effects of tissue inhomogeneity on the estimation accuracy of the indicative parameters. It has been demonstrated above that the VOI at 700 nm is markedly greater than the one at 500 nm, due to much lower absorption and scattering coefficients of biological tissues in the near infrared wavelength region.
In general, the proposed spectral interpretation algorithm could handle a moderate level (< 5%) of noise. By averaging the raw spectra, the estimation error could be significantly reduced. However, caution should be taken when estimated SatO 2 is below 30%, as it is more susceptible to the accuracy degradation induced by noise, as demonstrated in the validation study based on MC simulation. This susceptibility may be originated from the look-up table for SatO 2 estimation; the increase in the full width half max of the peak at the low SatO 2 levels makes the peak location detection highly inaccurate, with the presence of the strong noise.
To validate the utility and the efficacy of the new measurement scheme and the proposed algorithm in practice, an animal study with forepaw stimulation was carried out. Significant hemodynamic responses were expected to be observable in S1FL whenever the electrical current stimulation was introduced to the forepaw of the rat. The estimated indicative parameters depicted significant variations related to the stimulation, which was reproducible and clearly demonstrated the capability of the system in differentiating the variations in the hemodynamic (SatO 2 and A/[Hb]) and structural (A/ [Hb] and B) characteristics. The variations observed in the scattering power B might be a result of the increase of blood flow in the studied tissue, though it requires further investigation to gain insight to the underlying mechanism. Abookasis et al. [34] performed an open vascular occlusion to the middle cerebral artery of rats and observed changes in the scattering properties (both scattering amplitude A and scattering power B) of the cortex measured using the spatially modulated illumination methods. The pre-occlusion value of B reported by Abookasis (0.65 ± 0.1, mean ± standard error) is very close to the estimated B value at the resting-state (baseline) from the "old" rats in this study, which is 0.62 ± 0.02 (mean ± standard error). Meanwhile, this observation contradicts the assumption regarding constant scattering properties usually used in applying the modified Beer-Lambert law [32,33]. This might indicate those results obtained based on this constant scattering property assumption need to be carefully reinterpreted. One obvious limitation of the proposed spectrum interpretation algorithm is that it cannot separate [Hb] from the scattering amplitude A. Hence, it is impossible to distinguish whether a decrease in A/[Hb] (or increase in [Hb]/A, as shown in Fig. 11) during forepaw stimulation was predominantly the result of an increase in [Hb] or a decrease in A. Abookasis et al. [34] demonstrated that the scattering amplitude A in cortex decreases significantly (p < 0.05) when the middle cerebral artery is occluded. Meanwhile, the [Hb] and SatO 2 of cortex also decreased, albeit at a much higher level of statistical significance (p < 0.001). Their study results suggest that: 1) variations in [Hb] and A are directly correlated; and 2) variations in [Hb] are more contributory. Therefore, the variations in A/[Hb] observed in this study likely result more from changes in [Hb] than from changes in the scattering amplitude A. The results presented for this study are consistent with those from another previously-published study [34], in which a different imaging technique and interpretation algorithm were used.
There is one confounding factor that might render the new measurement scheme and the spectral interpretation algorithm inaccurate in quantifying the hemodynamic and structural characteristics of brain tissues, i.e., the heterogeneity of the tissue. Brain tissues have a distinct multi-layer structure and each layer has different optical properties. The thickness of the cerebral cortex usually varies from 2 to 4 mm in the human brain. Within the cortical layer, there are five or six sub-layers with unique cell structures and hence unique optical properties [35]. Our methodology was constructed based on a homogeneous medium. One may argue that such an assumption would lead significant estimation errors when the devised methodology is applied to multi-layer tissues. For this reason, we conducted a separate set of Monte Carlo simulations using the optical properties of gray and white matter as reported by multiple groups [35][36][37]. The results of these simulations (not reported here) suggested that the look-up table in Fig. 2, as well as the spectral interpretation algorithm, would remain accurate if the thickness of the cerebral cortex were greater than 2 mm. Otherwise, the total diffuse reflectance measured at cortical surface will be significantly influenced by the optical properties of the white matter, which could be the case for the present study on rat's S1FL cortex. Furthermore, the inhomogeneity in cortical layer would have a stronger impact on the estimation of A/[Hb] and B, but not on SatO 2 . A more comprehensive investigation on the heterogeneity effects on the estimation errors will be presented in the future.

"Underestimation" of SatO 2
The average value of SatO 2 measured from the S1FL cortex at the resting state ranged roughly from 36 to 44% in the anesthetized rats. These values were lower than those reported previously [34,38,39]. With the partial pressure of oxygen in the normal cortex being 35 mm Hg [40], the standard hemoglobin dissociation curve for Sprague-Dawley rats suggests that the expected SatO 2 of the brain cortex at rest should be roughly 48% [41]. This percentage was similar to, but still higher than, what we observed in the anesthetized rats.
Sakadžić et al. observed a very high heterogeneity of capillary oxygenation in the forepaw areas of the primary somatosensory cortex in mice anesthetized with isoflurane, with an SatO 2 even less than 20% at certain capillary [42]. They also found that the SatO 2 in normal pial venules could be as low as ~40% [42]. Therefore, it is possible that the "underestimation" of SatO 2 is a result of the small volume of investigation of the hybrid system ( Fig. 13(a)); the system could have picked up more information from the capillary networks as we intentionally avoid obvious pial arterials and venules. It is also highly possible that the total diffuse reflectance picked up at the cortical surface by point spectroscopic detection module had a higher contribution from the white matter, since the rat's cortex is very thin (< 2 mm). Because the information on oxy-hemoglobin was derived from total diffuse reflectance spectra using its profile characteristics between 540 and 560 nm (Fig. 15 in Appendix B), a small error in detecting the location of the peak also could lead to a significant underestimation or overestimation of SatO 2 . There are two common causes of spectral distortion, i.e., an inaccurately-calibrated spectrometer and poor resolution of the spectrometer. Therefore, it is important to control the spectral resolution of the spectrometer to 1 nm or lower and to calibrate the spectrometer carefully in order to minimize spectral distortion and, hence, achieve accurate estimates of SatO 2 .
Another major contributor to the "underestimation" of SatO 2 could be the anesthetic agents we used in this study, i.e., isoflurane and dexmedetomidine. Isoflurane has been considered to be a vasodilator, which could significantly suppress the hemodynamics, including blood flow, mean arterial pressure, and oxygenation [43][44][45]. Conversely, dexmedetomidine has been found to be a vasoconstrictor, so it could have inhibited or counteracted the cerebrovascular dilation induced by isoflurane [ Although it was not the intended scope of this study, the system also was able to detect an age-related difference in (1) the baseline levels of all three indicative parameters (A/ [Hb], B, and SatO 2 ) and (2) the magnitudes of the stimulation-induced variations in two indicative parameters (B and SatO 2 ) among the rats in this study (Rats 1, 2, and 3 vs. Rats 4, 5, 6, and 7). This phenomenon may be attributed to the aging of the brain. Age-related differences in hemodynamics have been previously identified in isoflurane-anesthetized rats [50]. The isoflurane-induced changes in hemodynamics also were found to be age-dependent [45,47]. It has been demonstrated that, even in healthy human subjects, oxygen consumption in the brain decreases with age without the effects of anesthesia [51,52], which could result from the combined effects of neuronal loss, cellular biological impairment, and functional deafferentation. Such effects were inferred to underlie or reflect the age-related cognitive changes in humans [52]. These age-related cognitive deficits also exist in rats if they are exposed to volatile anesthetic agents, such as isoflurane, which produces a sustained learning impairment in aged rats [53]. Researchers have shown that such anesthesia-induced cognitive deficits might be related to the reduced cerebral blood flow [44] and oxygenation [45], and the increased permeability of the blood-brain barrier [54] in the aged anesthetized rats.

Conclusion
In this paper, a hybrid spectroscopy imaging system and simple data interpretation algorithm were introduced to quantify the optical properties (μ s '/μ a ) of a given turbid media through the measurement of relative total diffuse reflectance. The system is capable of acquiring total diffuse reflectance (R TD ) spectra, over a broad wavelength range, from a given point on a targeted sample. With the help of reference material with known optical properties and a reference table constructed using Monte Carlo (MC) simulations, the system is able to efficiently and accurately convert the measured total diffuse reflectance to optical properties (μ s '/μ a ). The effects of other optical properties -namely the refractive index n and anisotropy factor g -on the accuracy of the estimation algorithm were also evaluated. It was found that only alterations in the refractive index could introduce significant errors to the estimates of μ s '/μ a . The proposed spectral interpretation algorithm is capable of handling moderate levels of random noise imposed on the measured data, as investigated based on MC simulation. The feasibility of using this system to estimate the hemodynamic and structural characteristics of biological tissues, through the relative R TD , was successfully demonstrated as well. In a rat forepaw stimulation study, the proposed measurement scheme and spectral interpretation algorithm were able to detect significant hemodynamic and structural variations in the somatosensory cortex for forelimb caused by electrical current stimulation. Both theoretical and experimental validations indicate that the new measurement scheme and accompanying spectral interpretation algorithm could uncover the intrinsic hemodynamic and structural characteristics of biological tissues effectively, which could then have tremendous value in both intra-operative guidance/monitoring and preclinical investigations.

Appendix A: Total diffuse reflectance
Diffuse reflectance R d (r) defines the percentage of the incident photon energy that is diffusely reflected by a tissue sample as a function of the distance between the illumination and the observation points r. Therefore, total diffuse reflectance R TD can be express as: If the surface of a tissue sample is illuminated by a large beam with a uniform intensity distribution, as illustrated in Fig. 14(a), diffuse reflectance observed from an arbitrary point (x,y) on the tissue surface will be ( ) ( ) where r is the distance between an illumination point (x',y') and the observation point (x,y).
Let (x,y) = (0,0), then 2 2 ' ' r x y = + and R(x,y) will be the same as R TD . where "xxx" could be either "paper" or "unknown" here. This relationship is applicable to both reference and evaluation measurements. It is also valid for the Monte-Carlo simulation set-up in the evaluation study, since it is always a relative comparison. As long as the systematic parameters are the same across the entire evaluation study, this f number will be the same for R TD of each μ s '-μ a pair. Therefore, the effect of f will be cancelled in the ratio operation in Eq. (7)  It should be noted that Eq. (A.6) is valid only if the angular distributions of the diffuse reflectance signals from the paper and the unknown material are identical. The changes in a look-up table of μs'/μa verse R TD resulting from g variation were investigated in a quantitative fashion. The μ s '-μ a pairs generated in the previous section were used in conjunction with n = 1.4 and g = 0.85, 0.90 or 0.95 in the MC simulation model to simulate R TD . A reference look-up table of μ s '/μ a verse R TD was established using the simulation results with g = 0.90 and, subsequently, the simulated R TD with g = 0.85 and g = 0.95 were applied to this look-up table to retrieve the estimated μ s '/μ a values. Finally, the discrepancies between the estimated and the input μ s '/μ a were quantified using Eq. (5). Additionally, the input μ s '/μ a values were applied to the look-up table of μ s '/μ a verse R TD for g = 0.90 to estimate R TD . The discrepancies between the estimated and the simulated R TD were again quantified using Eq. (5), with v 1 being the simulated R TD and v 2 the estimated one.

C.2 Varying n
To evaluate the effects of n on a look-up table of μ s '/μ a verse R TD , the input μ s '/μ a values created in the previous section were used in conjunction with g = 0.9 and n = 1.3, 1.4 or 1.5 in a MC simulation model to simulate R TD . A reference look-up table of μ s '/μ a verse R TD was established using the simulation results with n = 1.4 and, subsequently, the simulated R TD with n = 1.3 and n = 1.5 were applied to this look-up table to retrieve the estimated values of μ s '/μ a . Again, the discrepancies between the estimated and the input μ s '/μ a were quantified using Eq. (5). Additionally, the input μ s '/μ a values were applied to the look-up table of μ s '/μ a verse R TD for n = 1.4 to estimate R TD . The discrepancies between the estimated and the simulated R TD were again quantified using Eq. (5), with v 1 being the simulated R TD and v 2 the estimated R TD .

C.3 Varying g and n
A similar Matlab program was developed to generate a set of 50 random μ s ' and μ a for each of the following conditions: 1) n = 1.395 and g varying between 0.85 and 0.95; 2) n varying between 1.3 and 1.5 and g = 0.9; and 3) n varying between 1.3 and 1.5 and g varying between 0.85 and 0.95. Applying these optical property sets to the MC simulation model, simulated R TD were obtained and converted to μ s '/μ a using the look-up table shown in Fig. 2. Each estimated μ s '/μ a was compared with the corresponding input μ s '/μ a value to determine the absolute percentage error of estimation defined in Eq. (5).

Funding
Research report in this article was supported by National Cancer Institute of the National Institutes of Health under award number 1R15CA173617-01A1 and Nicklaus Children's Hospital.