Influence of the Signal-To-Noise Ratio on Variance of Chromophore Concentration Quantification in Broadband Near-Infrared Spectroscopy

This study presented a theoretical or analytical approach to quantify how the signal-to-noise ratio (SNR) of a near infrared spectroscopy (NIRS) device influences the accuracy on calculated changes of oxy-hemoglobin (Δ[HbO]), deoxy-hemoglobin (Δ[HHb]), and oxidized cytochrome c oxidase (Δ[oxCCO]). In theory, all NIRS experimental measurements include variations due to thermal or electrical noise, drifts, and disturbance of the device. Since the computed concentration results are highly associated with device-driven variations, in this study, we applied the error propagation analysis to compute the variability or variance of Δ[HbO], Δ[HHb], and Δ[oxCCO] depending on the system SNR. The quantitative expressions of variance or standard deviations of changes in chromophore concentrations were derived based on the error propagation analysis and the modified Beer-Lambert law. In order to compare and confirm the derived variances versus those from the actual measurements, we conducted two sets of broadband NIRS (bbNIRS) measurements using a solid tissue phantom and the human forearm. A Monte Carlo framework was also executed to simulate the bbNIRS data under two physiological conditions for further confirmation of the theoretical analysis. Finally, the confirmed expression for error propagation was utilized for quantitative analyses to guide optimal selections of wavelength ranges and different wavelength combinations for minimal variances of Δ[HbO], Δ[HHb], and Δ[oxCCO] in actual experiments.


INTRODUCTION
Near-infrared spectroscopy (NIRS) (Quaresima et al., 2012;Boas et al., 2014;Scholkmann et al., 2014;Bigio and Fantini, 2016;Shi and Alfano, 2017) has been considered an effective and non-invasive tool for functional imaging and diagnostics in medical applications. Besides the conventional NIR window from 650 to 950 nm (Villringer et al., 1993), the advancement of new technology made other NIR windows within the 1,000-2,500 nm wavelength range become possible Wang W et al., 2016;Alfano et al., 2018). These NIR windows have gained much attention recently because of their capacity to obtain greater imaging and sensing depth, and thus capable of investigating both soft tissue constituents (such as with high lipid content or cancers) and hard tissue constituents (such as bones) Sordillo et al., 2017).
The modified Beer-Lambert law (Liu et al., 2000;Sassaroli and Fantini, 2004) is mostly employed as the mathematical fundamental of a NIRS system to quantify changes in chromophore concentrations, namely, changes of oxy-hemoglobin (Δ[HbO]), deoxy-hemoglobin (Δ[HHb]), oxidized cytochrome-c-oxidase (Δ[oxCCO]). Theoretically, the minimum number of wavelengths required for the calculation of chromophore concentration changes must be minimally equal to the number of chromophores. When Δ[oxCCO] is the parameter of interest, however, such an approach is extremely sensitive to the system noise and may lead to inaccurate quantification because of its low concentration in tissue with respect to Δ[HbO] and Δ [HHb]. Thus, an adequate signal-to-noise ratio (SNR) of a NIRS measurement becomes a key determinant to judge or predict accuracy of experimental results of chromophore concentration changes. Broadband NIRS (bbNIRS) (Heekeren et al., 1999;Kolyva et al., 2012;Bale et al., 2014;Kolyva et al., 2014), which provides a full or broad wavelength range of measurement, is expected to alleviate the influence of system noise and offer more reliable quantification results. Several studies have also been carried out to identify optimal wavelength combinations in order to minimize the measurement redundancy but still ensure the accuracy of the quantification results (Wobst et al., 2001;Arifler et al., 2015).
Since Δ[HbO], Δ [HHb], and Δ[oxCCO] must be inferred from the optical densities measured at multiple wavelengths, measurement variations caused by the system's thermal or electrical noise, drifts, and disturbance will lead to variability in the respective quantities. To our knowledge, very few studies have addressed such a problem. Most existing works have focused on the effect of extinction coefficients, optical pathlengths, or wavelength combinations on the accuracy of estimated changes of chromophore concentrations. Kim and Liu (Kim and Liu, 2007) proved that small variations in hemoglobin extinction coefficients would lead to large variations in quantifications of hemoglobin concentration changes. Funane et al. (Funane et al., 2009) Sudakou et al. (2019) presented an error propagation analysis method to estimate the deviations of the recovered tissue constituent concentrations using a time-resolved NIRS.
In this study, we sought to investigate the influence of SNR of a bbNIRS system on quantifications of Δ[HbO], Δ[HHb], and Δ [oxCCO]. Practically, all measured data are not noise-free; they must consist of natural variability or uncertainty from the measurement system. The accuracy of the results derived from these measurements will depend on the measurement errors or SNR (Goodman, 1960;Bevington et al., 1993). In other words, measurement variances or errors will propagate to the quantified results. Note that the objective of this study was to quantify the SNR-derived variance of Δ[HbO], Δ[HHb], and Δ[oxCCO] caused only by the noise of the bbNIRS instruments/devices. Other parameters, such as the extinction coefficients and the optical pathlengths, are not the concern of variables in this study. Being aware of how much system noise propagates into the calculated chromophore concentration changes is essential when considering a bbNIRS system. Since the system SNR can be calculated easily throughout multiple baseline measurements, we can estimate or predict the SNR-derived variance of Δ[HbO], Δ[HHb], and Δ[oxCCO] based on the results of error propagation analysis. Consequently, one may want to improve the SNR of the bbNIRS system to lessen the uncertainties of calculated chromophore concentration changes by considering warming up the system to reach a stable state, optimizing the light exposure time, or selecting an optimal wavelength range.
Specifically, we performed the error propagation analysis for bbNIRS to calculate variances of Δ[HbO], Δ[HHb], and Δ[oxCCO] induced by the measurement system's noise. We considered only the case where the number of wavelengths is larger than the number of chromophores. Thus, the chromophore concentration changes were estimated by fitting the model to the measured data using the least-squares method (Heekeren et al., 1999;Kolyva et al., 2012;Bale et al., 2014;Kolyva et al., 2014). As a consequence, the analytical expressions of the error propagation in estimating Δ[HbO], Δ[HHb], and Δ[oxCCO] were derived directly from the best-fit model. In order to compare the analytical results with experimental results, we first carried out two real bbNIRS experiments, namely, 1) one taken from a solid tissue phantom using two spectrometers concurrently and 2) the other from the human forearm at resting state. Then, a Monte Carlo (MC) simulation framework that permits to replicate or simulate the measurement of a bbNIRS system was conducted. To demonstrate good consistency of influence of SNR on chromophore concentrations among all three cases, the theoretical quantitation of the SNR-derived variances of Δ[HbO], Δ[HHb], and Δ[oxCCO] were compared with the respective variances calculated directly from the measured and simulated data. Finally, we also performed analyses to assess optimal selections and ranges of wavelengths for bbNIRS to minimize variances of Δ[HbO], Δ[HHb], and Δ[oxCCO].

Modified Beer-Lambert Law
It is known that changes in chromophore concentrations are calculated based on the modified Beer-Lambert law (Liu et al., 2000;Sassaroli and Fantini, 2004;Kocsis et al., 2006): where ΔOD(λ i ) is the change in optical density at wavelength λ i , I 0 and I are the detected light intensities of the baseline and transient conditions, respectively. ϵ C k (λ i ) is the extinction coefficient of the k-th chromophore at wavelength λ i , Δ[C k ] is the change of the k-th chromophore concentration, and L(λ i ) is the optical pathlength of wavelength λ i . Such optical pathlength is estimated as L(λ i ) = r ×DPF(λ i ), where r is the source-detector separation distance, and DPF(λ i ) is the differential path-length factor, which is introduced to consider light scattering effects in Beer-Lambert's law. If the number of chromophores and wavelengths are equal, the chromophore concentration changes are calculated by simultaneously solving the number of equations (i.e., Eq. 1 with respective parameters). If the wavelength number is larger than the number of chromophores, the chromophore concentration changes are estimated by fitting the model (i.e., Eq. 1) to measured data using the least-squares method (Heekeren et al., 1999;Kolyva et al., 2012;Bale et al., 2014;Kolyva et al., 2014). In this study, we consider the general case using multiple wavelengths to estimate Δ[HbO], Δ[HHb], and Δ[oxCCO]. In such a case, the modified Beer-Lambert law, namely, Eq. 1, can be rewritten as follows: where m is the total number of wavelengths, ϵ HbO (λ i ), ϵ HHb (λ i ), and ϵ diffCCO (λ i ) are the extinction coefficients of HbO, HHb, and oxidized-reduced difference of CCO (Bale et al., 2016) at wavelength λ i , respectively.

Error Propagation From
Eq. 2 is equivalent to:

Δ[HbO], Δ[HHb]
and Δ[oxCCO] in Eq. 2 are typically estimated using the least-squares fitting method. The principle is to fit the results from the linear modelŷ Ca to the measured data y [ΔOD(λ 1 ) / ΔOD(λ m )] T as closely as possible. Such a best-fit minimizes the chi-squared optimization function defined as follows: where V y σ 2 is the diagonal matrix of variances of y and σ 2 yi is the variance of the measurement y i . The χ 2 optimization function is minimized with respect to the parameters a by solving the equation of zχ 2 za | a â 0, which leads to the linear least-square estimate of the parameters a as follows:â Note that the linear least square estimate of Δ[HbO], Δ[HHb] and Δ[oxCCO] defined in Eq. 5 is a function of the measured ΔOD, the ΔOD covariance matrix, the extinction coefficients, and the optical pathlength. As mentioned above, we consider the error propagation caused only by the noise from instruments. Thus, matrix C, which includes the extinction coefficients and the optical pathlengths at all wavelengths, is considered a constant matrix. The covariance matrix of the estimates of Δ[HbO], Δ[HHb] and Δ[oxCCO] is calculated as follows: By substituting the matrices C and V y in Eq. 6 with the original definitions, the covariance matrix of the estimated Δ[HbO], Δ[HHb] and Δ[oxCCO] can be rewritten as: For simplification, let ξ HbO (λ i ) = L(λ i )ϵ HbO (λ i ), ξ HHb (λ i ) = L(λ i )ϵ HHb (λ i ), and ξ diffCCO (λ i ) = L(λ i )ϵ diffCCO (λ i ), the covariance matrix Vâ can be simplified as follows: This covariance matrix is known as the error propagation matrix, which indicates how measurement errors in ΔOD, described by If multiple measurements are repeated for both baseline and transient instant, the detected light intensities I 0 and I can be replaced by the mean values of all measurements' detected light intensities. Mathematically, the fraction I 2 /σ 2 I is equivalent to the SNR of the measurement system. Furthermore, assuming that the measurement system has consistent SNR for both baseline and transient instant measurements, Eq. 9 can be rewritten as: Substituting Eq. 10 into Eq. 8, we obtain the error propagation matrix calculated from the SNR of the measurement system.
The square root of the diagonal of this covariance matrix corresponds to the standard deviation of Δ[HbO], Δ[HHb] and Δ[oxCCO], which indicates the uncertainties of the estimated chromophore concentration changes caused by the measurement system noise. We note the standard deviation of the chromophore concentration changes derived by the proposed error propagation analysis as σ T

Actual Measurements From a Solid Phantom and the Human Arm Using bbNIRS
In order to evaluate how well the derived error propagation calculations matched with the real bbNIRS data, we first conducted actual measurements from a solid phantom using two separate bbNIRS spectrometers combined with one light source, as shown in Figure 1A. The two spectrometers utilized 1) a twodimensional CCD spectrograph (Teledyne Princeton Instrument, 3660 Quakerbridge Road Trenton, NJ 08619 United States) and 2) a back-thinned cool-down CCD spectrometer (QE-Pro, Ocean Optics Inc.). A tungsten halogen lamp (Model 3900, Illumination Technologies Inc., East Syracuse, NY) covering 400-1,500 nm light was used as the light source. We used the solid phantom provided by ISS (ISS Inc., Champaign, IL, United States) with the absorption coefficient μ a of 0.155 cm −1 at 690 nm and 0.15 cm −1 at 830 nm. Three fiber bundles were used, one for the light delivery and two for light collection from the solid phantom. The distance between each detection bundle to the source bundle was 3 cm. The bbNIRS data were collected concurrently by these two spectrometers to ensure the identical environmental condition. The exposure time of both spectrometers was set to 1 s, and the data were collected continuously over 30 min, resulting in a total of 1800 data points for each spectrometer. We repeated the measurements seven times on different days to ensure the arbitrary nature of the data collection.
The second experiment was conducted on the human forearm under the resting state with no stimulation. The experimental protocol was approved by the Institutional Review Board (IRB) of the University of Texas at Arlington. Two healthy adults participated in this experiment. Each subject attended eight measurement sessions on different days, which led to n = 16 measurements. Informed consent was acquired prior to all measurements.
The experimental setup was illustrated in Figure 1B. The bbNIRS system consisted of a tungsten halogen lamp (Model 3900, Illumination Technologies Inc., East Syracuse, NY) and a CCD spectrometer (QE-Pro, Ocean Optics Inc.), both of which were used in the phantom measurements. A flexible probe holder was used to firmly hold the two optical fiber bundles on the subject's forearm. An optical shutter was also employed to switch on the light only during the exposure time to minimize the tissue heating effect from the broadband light source. An experiment lasted 15 min with a 5 s single light-exposure/data-acquisition time per minute, giving rise to 16 spectra per measurement.

Monte Carlo Simulations
We also conducted a Monte Carlo (MC) framework that simulated bbNIRS data of the hemodynamic and metabolic states of the exposed tissue. Two physiological conditions were defined within the MC simulation: 1) a baseline state corresponding to the normal resting state of the tissue and 2) a stimulated state where oxygenation and metabolism improved notably. We employed MCmatlab (Marti et al., 2018), an open-source MC program for light propagation in three-dimensional (3-D) media. Specifically, a (600 × 600 × 300) voxel model corresponding to a (4 × 4 × 2 cm) 3-D volume was first defined as the simulation geometry. The medium that replicates general tissue was considered to be made of water, fat, blood with different oxygen saturation levels, and concentrations of oxCCO and reduced CCO (redCCO). The composition and the optical properties of the medium were defined based on the properties of biological tissues reported in (Jacques, 2013;Jacques et al., 2013). The scattering coefficient μ s (λ) was considered to be dependent only on the wavelength λ, while the absorption coefficient μ a (λ) was estimated as the sum of absorption coefficients of all major chromophores composing the medium (Jacques, 2013;Bale et al., 2016). Table 1 summarizes the simulation medium composition of two physiological conditions, while Figure 2 depicts the corresponding absorption coefficient μ a (λ) (Figure 2A) and the scattering coefficient μ s (λ) ( Figure 2B). The Henyey-Greenstein scattering anisotropy factor g was set to 0.9.
Note that the concentrations of HbO and HHb were estimated from the predefined percentages of blood and oxygen saturation of hemoglobin, with an average concentration of hemoglobin in blood equal to 150 g/L.
A light source and a detector of 2 cm distance were also defined and placed on the top surface of the 3D voxel model. MC simulation was repeated for different wavelengths within the NIRS range from 780 to 900 nm with 1-nm wavelength resolution, leading to a total of 121 wavelengths, for two physiological conditions. For each MC simulation execution, 50 million (i.e., 5 × 10 7 ) photon packets were launched from the light source. The number of photons reaching the detector and their partial pathlengths were recorded for later use when calculating the changes in chromophore concentrations. We repeated the MC simulation 20 times for each wavelength and for each physiological condition. The recorded data were further permuted within each wavelength and each condition 50 times. Gaussian noise was also added to enrich the simulated data. Thus, the whole simulation framework can be considered to generate bbNIRS data of 50 measurements. Each consisted of 20 spectra for each condition, namely, I 0 for the baseline condition and I for the stimulated state.   and σ T Δ [oxCCO] were finally compared to verify/confirm the accuracy of the quantitative analysis. from the theoretical analysis (Eq. 11). Specifically, the mean and standard deviation spectra of SNR estimated from the detected optical spectra of multiple measurements are depicted in Figures 4A-C, taken from (A) spectrometer 1 in the phantom measurement, (B) spectrometer 2 also in the concurrent phantom experiment, and (C) the human forearm measurement, respectively. Consequently, Figures 4D-F  by theoretically derived EP (blue bars). As seen in Figures 4A-C, the spectral shapes of SNR appeared to be different from one another. This observation is understandable since each of them was obtained from a specific bbNIRS system. Specifically, SNR spectra in Figures  4A,B were derived from the data collected concurrently by two spectrometers under the same experimental conditions, and the differences in the spectral sensitivity of these two spectrometers led to distinctive SNR spectral shapes. The bbNIRS data from the human arm experiment ( Figure 4C) were acquired by the same QE-Pro spectrometer as in the solid phantom/spectrometer 2 case ( Figure 4B). However, different exposure times for data acquisition (5 vs. 1 s) resulted in un-identical SNR spectra.

Standard Deviation of Δ[HbO], Δ[HHb] and Δ[oxCCO] Derived From Error Propagation Analysis and Actual Measurement Data
Despite the variations of SNR spectra, no significant differences were found between measured σ M Δ [HbO] , σ M

Δ[oxCCO]
calculated directly from the EP analysis (Eq. 11) in all the cases. Specifically, the human forearm measurement led to a high SNR range (i.e., 57-58 dB) ( Figure 4C) (Figure 4F) because of a prolonged exposure time for data acquisition (5 s). The same spectrometer was used in the solid phantom measurement, but a shorter exposure time of 1 s resulted in a reduced SNR range of the measured data ( Figure 4B)

Standard Deviation of Δ[HbO], Δ[HHb] and Δ[oxCCO] Derived From Error Propagation Analysis and MC Simulation Data
We also compared the uncertainties of Δ[HbO], Δ[HHb] and Δ[oxCCO] obtained from the MC simulation and EP results derived from the corresponding SNR. Figure 5A depicts the mean and standard deviation spectra of SNR of the MC computed from the MC simulation data (red bars). One more time, no significant difference was found in the standard deviation of chromophore concentration changes calculated from the MC simulation data and from the EP analysis. The SNR of the MC simulation ( Figure 5A) was lower than those of the real measurements ( Figures 4A-C) due to the MC simulation settings. Nevertheless, it could be improved by having a larger tissue volume and a larger number of photons. Generally, the higher the average SNR of the measurement system, the lower the error propagates to the calculated changes in chromophore concentration.
These results confirmed that the analytically derived expressions of σ given an SNR of a bbNIRS system. In the following two sub-sections, we applied the EP analysis to examine the influence of wavelength selection and spectral range of bbNIRS on the variance of chromophore concentration changes. Note that the values of the optical pathlength L(λ i ) = r ×DPF(λ i ) can be estimated by diffusion theory, as given in refs. Patterson et al. (1989), Sevick et al. (1991).   . To address this problem, broadband spectroscopy with a full wavelength range of measurement has been suggested and utilized in recent studies, such as 740-900 nm in studies of (Heekeren et al., 1999;Kolyva et al., 2012;Kolyva et al., 2014;Wang X. et al., 2016), 770-905 nm in publications of (Bale et al., 2014(Bale et al., , 2019, and 780-900 nm in reports of (Heekeren et al., 1999;Kolyva et al., 2012;Kolyva et al., 2014;Giannoni et al., 2020). Since a large number of wavelengths required in bbNIRS can be an obstacle for development of a compact or portable system (Chitnis et al., 2016), various studies made efforts to find an optimal set of wavelengths to minimize the number of required wavelengths while keeping satisfactory accuracy for the quantifications (Arifler et al., 2015). Bale et al. (2016) reviewed existing clinical NIRS systems and summarized different broadband ranges or wavelength combinations used in these systems. We followed this review paper and selected/updated some representative wavelength ranges or wavelength combinations to illustrate how wavelength selections affect the accuracy of the calculated chromophore concentration changes. uncertainty using 16 consecutive wavelengths separated by an equal interval of 12.5 nm with a moving spectral window covering from 650 to 950 nm. However, such a study focused on the timeresolved NIRS and utilized a 3-layer model (scalp, skull, and brain). Consequently, the error propagation analysis included the moment method for analyzing time-resolved statistical uncertainty and multilayer considerations. Thus, the variance of Δ[CCO] and respective analysis by Sudakou et al. were distinct from this work and thus excluded in Table 2 for comparison. Figure 6 depicts σ T Δ [HbO] , σ T Δ[HHb] and σ T Δ[oxCCO] calculated using the EP analytical expression (i.e., Eq. 11) for different broadband wavelength ranges ( Figure 6A) or different wavelength combinations ( Figure 6B). Selections of the ranges or combinations were based on the exact values given in Table 2. In TABLE 2 | Summary of different ranges of wavelengths or optimal wavelength combination sets.
Frontiers in Photonics | www.frontiersin.org June 2022 | Volume 3 | Article 908931 8 these calculations, we temporally assumed a constant SNR value across all wavelengths and performed the EP calculation for different SNR values varying from 20 to 60 dB (equivalently, noise level (NL) being from 10% to 0.1%). As shown in Figure 6A, the broadest bandwidth of 720-920 nm (with 201 wavelengths) led to the smallest σ  Figures 6A,B. For instance, at an SNR of 40 dB, EP varied from 0.08-0.12 μM for Δ[HbO] and from 0.05-0.1 μM for Δ[oxCCO] when using all wavelengths of bbNIRS, while the EP increased by at least four times when using a limited number of wavelengths; namely, EP varied from 0.25 to 0.64 μM for Δ[HbO] and from 0.2 to 0.5 μM for Δ [oxCCO]. In most cases, the more wavelengths used to calculate concentration changes of each chromophore, the smaller error propagation in the results. For the case that the same number of wavelengths was used (e.g., 8 wavelengths used in (Arifler et al., 2015;Giannoni et al., 2020) or in (Chitnis et al., 2016)), the distinct wavelength combinations led to different error propagation results. For instance, these two sets of 8 wavelengths led to approximate σ T Δ[HbO] value but distinct σ T

Δ[HHb]
and σ T Δ[oxCCO] (green and pink curves in Figure 6B). Note that the results presented here should be considered as approximate references to roughly estimate σ T Δ[HbO] , σ T Δ[HHb] and σ T Δ[oxCCO] measured from a bbNIRS system with an approximate SNR varying within the given range. For a specific bbNIRS system with a known SNR spectrum, more precise estimations of error propagation in calculating chromophore concentration changes can be carried out by substituting the system's SNR spectrum into Eq. 11 according to the respective wavelength ranges or wavelength combinations. HHb] and σ T

Δ[oxCCO]
As revealed in Figure 6, using a broader wavelength range in bbNIRS led to significantly smaller σ T Δ [HbO] , σ T Δ[HHb] and σ T Δ [oxCCO] . However, to our knowledge, no study has reported the spectral-bandwidth dependence of variance in calculating chromophore concentration changes. In the following, we demonstrate how selections of spectral ranges would impact the variances of Δ[HbO], Δ[HHb] and Δ[oxCCO] based on the derived Eq. 11.
As an example, three panels of Figure 7 depict σ T Δ[HbO] , σ T Δ[HHb] and σ T Δ[oxCCO] with different spectral ranges at the SNR of 40 dB (for simplicity). Each curve in each panel plots EP dependence on spectral ranges with a starting wavelength of λ i and an ending wavelength of λ j , where λ i varies from 700 to 800 nm in a step of 10 nm (as noted on top of the figures), and λ j corresponds to the x-axis from 820 to 920 nm in each panel, also in a step of 10 nm. For instance, the red curves represent σ concentration changes within wavelength ranges of 710-820 nm, 710-830 nm, . . . , 710-920 nm.
The three zoomed panels in Figure 7 provide more detailed information. Specifically, σ T Δ[HbO] varied around or below 0.15 μM if the starting λ i was shorter than 780 nm and the ending λ j equal to or longer than 880 nm. Meanwhile, both σ T Δ [HHb] and σ T

Δ[oxCCO]
had a large reduction if the start λ i was shorter than 760 nm. Overall, the spectral range of 760-890 nm can be considered as an optimal range to minimize variances of calculated  (Kim and Liu, 2007) or different wavelength combinations (Funane et al., 2009;Sudakou et al., 2019;Caredda et al., 2020) Δ[oxCCO] caused by uncertainties of the bbNIRS measurement system. The derivation of EP of respective concentration changes was achieved by fitting the measured data using the least-squares method (York, 1966). The analytical expressions and experimental/MC-simulation results were compared and confirmed statistically nonsignificant difference between the theoretical and measurement-based variances of chromophore concentration changes. Experimental results also indicated that the larger the SNR of a bbNIRS system, the smaller the quantified σ T Δ [HbO] , σ T Δ[HHb] and σ T Δ[oxCCO] ( Figure 4D-F). Note that in both actual experiments presented in Section 2.3 and Section 3.1, we tried to enhance the SNR of the bbNIRS system by warming up the system to reach a stable state before acquiring the data. For a new or self-designed bbNIRS system, the thermal and electrical fluctuations of the system over time should be first investigated to ensure the required warm-up time to reach the desired system stability.
In the case of bbNIRS, a large number of wavelengths with different spectral ranges or bandwidths have been considered when calculating the chromophore concentration changes (Heekeren et al., 1999;Kolyva et al., 2012;Bale et al., 2014;Kolyva et al., 2014;Wang X. et al., 2016;Bale et al., 2019;Giannoni et al., 2020), as summarized in . Accordingly, Figure 7 may serve as a useful guide or predictor to select an optimal spectral range/bandwidth for accurate determination of changes in chromophore concentrations by a bbNIRS measurement device. It is worth noting that although the results presented in Figures 6, 7 were calculated using a constant SNR value, this does not imply that the derived EP calculation cannot be applied to a real bbNIRS system where SNR has a spectral feature depending on the wavelength. We presented Figures 6,  7 for the purpose of giving the reader an approximated amount/order of σ T Δ [HbO] , σ T Δ[HHb] and σ T Δ[oxCCO] caused by a certain SNR level. For instance, in our previous studies on the effects of transcranial photobiomodulation (tPBM) (Wang et al., 2017;Pruitt et al., 2020), we employed the same bbNIRS system as that used in this study for the human arm measurement (Figures 4C, F) Some efforts have also been carried out to select a limited number of wavelengths to lessen the complexity of the spectral hardware or system but still ensure the calculation or quantification accuracy (Wobst et al., 2001;Arifler et al., 2015;de Roever et al., 2018;Caredda et al., 2020). One objective of such studies was to design wearable NIRS systems (Wyser et al., 2017;Chitnis et al., 2016 Figure 6B), compared to 0.07-0.1 μM for the full bandwidth of 720-920 nm ( Figure 6A, red curve). Moreover, distinct wavelength combinations led to different error propagation results. Thus, one may want to estimate the SNR-derived variance of Δ[HbO], Δ[HHb] and Δ[oxCCO] using different available wavelength combinations before deciding the optimal set that leads to minimal error propagation. Additional data processing or fitting steps, such as the genetic algorithm (Arifler et al., 2015), may be added to help improve the quantification accuracy. If such a step is taken, further error analysis should be carried out to consider the contribution of additional steps in the EP calculation.
Last, to better understand the scientific reasoning why the variance of chromophore concentration changes depends on the spectral range and wavelength selection, we performed a mathematical expansion of the EP matrix inverse (Eq. 11), as given in Appendix. The final results demonstrate that the variance for each of Δ[HbO], Δ[HHb], and Δ[oxCCO] is computed from the extinction coefficients of all three chromophores at given multiple wavelengths. Since these expressions are highly wavelength dependent and extinction coefficient nested, it is impossible to directly infer an optimal spectral bandwidth or wavelength combination that can lead to minimal error propagation. However, the knowledge learned from these equations is that the device-driven errors for Δ[HbO], Δ[HHb], and Δ[oxCCO] in a bbNIRS system depend closely on the extinction coefficient spectra of the three chromophores. Also, it might be theoretically possible to optimally select wavelengths by maximizing Eq. A2 and minimizing Eqs A4-A6 simultaneously for a given bandwidth. However, the latter point is beyond the scope of this work and needs to be investigated in future studies.

CONCLUSION
This study investigated the influence of SNR on variance of Δ[HbO], Δ[HHb] and Δ[oxCCO] measured by a bbNIRS device or system. Since all measured data contain inevitable uncertainties caused by thermal or electrical fluctuations or disturbance of the devices, such uncertainties or errors must impact the accuracy of the calculated Δ[HbO], Δ[HHb] and Δ [oxCCO]. Based on error propagation analysis, we derived analytical expressions of EP for all three chromophore concentrations depending on the SNR spectral curve of the bbNIRS measurement system. To compare the quantitative results and those obtained from actual bbNIRS measurements, we performed two sets of experiments on a solid tissue phantom and the human forearm using two bbNIRS systems. We also introduced an MC framework mimicking a set of bbNIRS measurements with two predefined physiological states of tissue. Both experimental and MC simulation results statistically confirmed and supported the analytical expression of the variance or EP of Δ

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by The Institutional Review Board of the University of Texas at Arlington. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
NT derived theoretical expressions, performed both phantom experiments and MC simulations, analyzed the data, interpreted the results, and prepared the manuscript. SS assisted and discussed theoretical formulation. SK performed human forearm measurements. XW discussed data analysis and the results. HL initiated and supervised the study, discussed and interpreted the results, as well as reviewed and revised the manuscript. All authors reviewed and approved the manuscript.

FUNDING
This work was supported in part by the National Institute of Mental Health/National Institutes of Health under the BRAIN Initiative (RF1MH114285).