Accuracy and precision of tissue optical properties and hemodynamic parameters estimated by the BabyLux device: a hybrid time-resolved near-infrared and diffuse correlation spectroscopy neuro-monitor

: We have investigated the accuracy and precision of “the BabyLux device”, a hybrid time-resolved near-infrared (TRS) and diﬀuse correlation spectroscopy (DCS) neuro-monitor for the pre-term infant. Numerical data with realistic noise were simulated and analyzed using the BabyLux device as a reference system and diﬀerent experimental and analysis parameters. The results describe the limits for the precision and the accuracy to be expected. The dependence of these limits on diﬀerent experimental conditions and choices of the analysis method is also described. Experiments demonstrate comparable values for precision with respect to the simulation results.

Since here we have used the same method for both the simulation and the analysis, we ignore this source of error. Apart from the model used, some literature focused on the methodological errors in TRS [30,32], in FD-NIRS [31] and in DCS [33,34]. In the present work we focus not only on the accuracy but also on the precision since the latter was shown to be critical for neonatal applications [6].
Finally, the performance of the BabyLux device in terms of precision has been assessed by measurements on tissue simulating phantoms and on a piglet model. Variability over continuous measurement is calculated for both phantom and piglet measurements and compared to the findings obtained from the simulations.

Time resolved spectroscopy theory
In order to describe the propagation of pulsed light in diffusive media, such as biological tissues, we consider the analytical model based on the photon diffusion equation (PDE) [2]. The PDE is solved for a semi-infinite homogeneous medium using the so-called extrapolated boundary condition [35]. The following equation is derived for the reflectance R(ρ, t) (Wcm −2 ) , i.e. the outward flux of photons at the boundary at a distance ρ from the source and at time t [36]: where µ a (absorption coefficient) and µ s (reduced scattering coefficient) are the wavelength dependent tissue optical properties. v is the speed of light in the medium, A accounts for the index mismatch between tissue and air [15,36,37]. D = v 3µ s is the photon diffusion coefficient, ρ is the distance between the injection point and the detection point (or the source-detector separation), z + = z s and z − = −2z e − z s with z s = 1/µ s and z e = 2A/3µ s . The expression of the reflectance reported in equation (eqn.) 1 is the Green's function of the PDE for the reflectance [2,36].
TRS techniques exploit the time-correlated single-photon counting technique to measure the photon DTOFs. It measures and builds a histogram of light intensity at a certain time determined by counting the photons having the same time of flight within a time window ("bin"). This measurement is essentially primarily affected by the Poisson noise of photon statistics [38].

Diffuse correlation spectroscopy theory
DCS measures the intensity fluctuations of the diffuse laser speckles at the surface of the turbid medium. The statistics of these fluctuations are affected by the moving scatterers (i.e. the red blood cells in most tissue experiments) and DCS aims to recover information about their motion by evaluating these statistics.
In particular, the electric field autocorrelation function is calculated. Analogous to the continuous-wave photon diffusion equation, a correlation diffusion equation (CDE) can be derived for the propagation of the electric field autocorrelation functions through the tissues [13,15]. This can be solved for a semi-infinite homogeneous medium employing the extrapolated zero boundary condition as we did for TRS. The following analytic expression for the electric field autocorrelation function G 1 (τ) [15] is then obtained where τ is the delay time and the decay constant K(τ) = µ a + 2µ s k 2 0 BF Iτ v D , with k 0 the wavenumber of light in the medium, depends on the blood flow index (BFI). BFI is derived from the effective Brownian diffuse coefficient D B of moving scatterers multiplied by their fraction respect to the total scatterers α (BFI=αD B ) [39]. In eqn. 2, r + = ρ 2 + z 2 + and r − = ρ 2 + z 2 − . As mentioned at the beginning of the section, the normalized intensity autocorrelation function (g 2 (τ)) is typically measured and is related to the electric field autocorrelation function through the Siegert relation which is valid under certain assumptions [40]: where g 1 (τ) = G 1 (τ)/G 1 (0) is the normalized electric field autocorrelation function and β is a parameter that depends on the light source and the collection optics of the experiment. In order to simulated realistic data, noise must be added to the analytical curve. A model for noise in DCS curves has been previously developed [41,42] and is reported in the Appendix (eqn. 6). A schematic of the process to obtain simulated TRS and DCS curves is depicted in Fig. 1, going from left to right. The starting point of TRS data simulation process (top row of Fig. 1) was the definition of levels for HbO 2 and Hb concentrations (c HbO 2 and c Hb ) and for the percentage of water with respect to the total chromophores concentration (c H 2 O ). Thanks to the linearity between chromophore concentrations and absorption coefficients [43], the latter could then be calculated at each wavelength of interest, knowing the molar absorption coefficients [44,45] (Table 7 of the Appendix).

Curve simulation and analysis process
As explained by Fig. 1, the scattering properties were derived from Mie theory, which defines how µ s (λ) scales with the wavelength: where λ 0 is a reference wavelength set at 600 nm and a and b are parameters that depend on the tissue scatterer density and size [46,47].
Once the values for optical properties were defined, the reflectance in time resolved regime R(t) was derived exploiting eqn. 1. As was mentioned, this represents the Green's function solution of the PDE, therefore, to mimic experiments, R(t) was convoluted with the instrument response function (IRF(t)): IRF(t) is the measured distribution of the time of flight of photons when light from the injection fiber is directly coupled to the detection fiber and the one measured by the BabyLux device was used for curve generation. As suggested by Fig. 1, after the convolution with the IRF and the adjustment of the area under the curve to the desired count level, Poisson noise was added to obtain the simulated time-resolved curve R sim (t). The simulated curve could eventually be shifted by t 0 in order to simulate instabilities in the temporal TRS laser position. One curve for each wavelength of interest was generated.
The simulated curves were then analyzed to retrieve optical and hemodynamic properties by reversing the simulation procedure, going from right to left in Fig. 1. Each of the three simulated R sim (t) curves was compared to the model R model (t) (eqn. 5). Only the range between the 90% of the peak value, on the rising edge of the curve, and 1% on the falling edge was included in the fitting procedure. The absorption and the reduced scattering coefficients were considered as the fitting parameters. Optionally t 0 , the shift of the TRS laser peak, could also be considered as a fitting parameter. From the retrieved values of µ a at three wavelengths and assuming a pre-defined level of relative water concentration, concentration of HbO 2 and Hb were calculated as well as StO 2 =HbO 2 /(HbO 2 +Hb).
The simulation and the analysis process for DCS is explained in the bottom row of Fig. 1. The analytical expression for the g model 1 (τ) as expressed by eqn. 2 was calculated with the desired BFI value. Once the value for β was defined, and by employing the Siegert relation (eqn. 3), noise was then added (eqn. 6 in the appendix)) to obtain g sim 2 (τ). As was done for TRS, the simulated DCS curves were analyzed following the reverse of the simulation process. From g sim 2 (τ), g sim 1 (τ) was obtained through the reverse of the Siegert relation by estimating β from a weighted average of the first three points of the g sim 2 (τ) curves (since g 2 (0) = 1 + β) or assuming a fixed value for it. The retrieved g sim 1 (τ) was then compared to the model, with BFI as the fitting parameter, using as a fitting range the τ where g 1 (τ) was lower than 0.5. Optical properties were introduced as input parameters of the analysis process.

Parameters for simulation and analysis
We now go through the details of all the values used for the input parameters for curve simulation and analysis. In Table 1, all the experimental and tissue properties with their defined values are listed. Those are the default values used for simulations if not explicitly stated. In particular, the experimental parameters were chosen considering the BabyLux device [48] as the reference system. The tissue hemodynamics and optical properties were chosen according to what was previously measured in newborns with the BabyLux device [49].
In order to identify the contribution of each parameter in the precision and the accuracy of the recovered values, we have generated different sets of simulated curves by varying their values. Each set of curves consisted of thirty independent simulations. Those thirty curves were either generated by keeping all parameters entering the simulation process fixed or by assuming a variability for some of those parameters by defining a coefficient of variation CV = σ y /y. Here y is the mean value of an input parameter and σ y is its standard deviation which is defined a priori. Thirty values were, thus, randomly generated from a normal distribution centered at y and with a width of σ y . Those thirty values were then used for generating each of the thirty simulated curves. Only at this point noise was added to each curve, as explained above and in Fig. 1. This is photon noise and mainly depend on the number of photon detected. The whole process was done ten times to have ten sets of thirty curves for each configuration. This allowed us to have some statistics in the variability and the error of the simulated data. Figure 2 shows the different questions we wanted to be answered by our simulations (left column) and the sets of parameters that were used for addressing these questions, either in the simulation (central column) or the analysis process (right column).

• Parameters for TRS
The role of the total number of photons collected in the DTOF was investigated by simulating curves at different total photon count levels, from a minimum of 10 4 photons to a maximum of 10 7 photons. The BabyLux device performs the measurements until 10 5 photons are collected in 1 s, if reachable. Therefore we have simulated curves at lower counts, to consider measurements where that limit is not achievable. In addition we have considered higher count levels, to understand the improvement reached by an eventual change in configuration to increase the light levels. Afterwards, we have varied the absolute level of StO 2 (from a minimum of 5% to a maximum of 80%) while keeping the tHb concentration constant at 80 µM and changing accordingly the HbO 2 and Hb concentration. Then, we have proceeded by considering an undetected instability of the position of the TRS laser peak. Curves were simulated by considering a shift (t 0 ) of the peak. t 0 values were generated from a normal distribution with a mean at 0 ps and a standard deviation (σ t 0 ) of 5, 10 and 15 ps. These curves were then analyzed either by considering t 0 fixed to 0 ps or by treating it as a fitting parameters. Furthermore, we have considered the estimation of input parameters of the TRS analysis, i.e. the water content. Specifically, we have simulated curves by varying the water content from 70% to 100%, in 5% steps, of the total chromophore concentration, while the water content was always considered as 90% in the analysis. Lastly, we have simulated curves assuming a CV for HbO 2 concentration (7%), Hb concentration (3%) and for the Mie parameters a and b (3% each). Those curves were then analyzed to check how the variability increased due to noise and analysis process.

• Parameters for DCS
As suggested by Fig. 2, the questions considered for DCS were analogous to those that were mentioned above for TRS. As a first step, we have varied the count rate in the range between 20 kHz and 120 kHz, considering that 80 kHz is usually reachable for measurement on infants. While the averaging time was varied between 1 s (currently used in BabyLux) and 5 s. The following step was to vary the BFI absolute value between 0.025×10 −8 cm 2 /s and 5×10 −8 cm 2 /s. Since DCS analysis requires optical properties as input parameters, we have introduced a variability in their values. Specifically, sets of curves were simulated considering different levels of CV for both the absorption and the reduced scattering coefficients (from 0.01% to 16%), while they were considered fixed in the analysis process (at the mean of the thirty values used in the simulation). The other input analysis parameter is β, which was always kept fixed at 0.48 in the simulation. On the other hand, in the analysis, it was either estimated from the g 2 (τ) curve or was assumed fixed by varying it from 0.45 to 0.5. Last point, in accordance with what was done with TRS, focused on the role of BFI variability. Curves were generated with an assumption of variability for the BFI (from 1% to 10%) to check how the CV for BFI increased due to noise and analysis process.

The hybrid DCS and TRS device
The BabyLux device [48] was used to acquire experimental DCS and TRS curves. Briefly, it integrates both TRS and DCS. DCS uses a single long coherence length laser at 785 nm, while TRS uses three laser heads at 685 nm, 760 nm, 820 nm. A small percentage of the TRS laser power is deviated after the laser heads to be detected directly by the TRS detector in order to monitor the parameters of the injected laser pulses, e.g. the shift of the laser peak positions (t 0 ). Variable attenuators regulate TRS laser power at the probe according to the sample probed. They are adjusted by an automatic procedure in order to reach a fixed level of counts (e.g. 10 5 ) in 1 s, or the closest possible value. On the other hand, DCS light power is fixed at 20 mW at the probe level. DCS and TRS measurements are acquired simultaneously with 1 s of time resolution. The launching and detector fibers are integrated in a single probe, with a source detector separation (ρ) of 15 mm. The probe employs prisms to deflect the light towards the tissue in the perpendicular direction respect to the fiber cables. The details of the device is described elsewhere [48].

Phantom measurements
Phantoms were measured to test the precision of the device over multiple continuous measurements for both TRS and DCS. A solid phantom made by epoxy resin as solid matrix, TiO 2 particles as scatterers and ink powder as absorber was measured by TRS. It was made with nominal optical properties of µ a = 0.1 cm −1 and µ s = 10 cm −1 at 690 nm [50]. On the other hand, a liquid phantom, provided by HemoPhotonics S.L. (Castelldefels, Spain) [51], was employed for DCS, since this technique probes moving scatterers. The count rate for this measurement was 150 kHz.
Each set of phantom measurement, either solid or liquid, consisted of 30 curves.

Piglet model
The permit to perform the animal experimentation was obtained from the Danish animal experiments inspectorate (license no. 2016-15-0201-01021). Precautions were taken to reduce any animal suffering to a minimum. Anesthesia was induced shortly after arrival to our facility and sustained throughout preparation and experimentation until termination of the animal. Newborn piglets (Danish landrace) were delivered on the day of the experiment and rapidly anesthetized by intramuscular injection of Zoletil (xylazine 2 mg/kg, ketamine 10 mg/kg, methadone 1 mg/kg and butorphanol 1 mg/kg; Virbac, Denmark), intubated and put on a ventilator (Dameca, Denmark). Once the piglet was stable, the BabyLux probe was placed on the piglet head and 30 curves at 1 s time resolution were acquired. After the completion of the protocol of the main study, whose results are not presented here, the piglets were euthanized with an overdose of intravenous pentobarbital 150mg/kg. Cerebral hemodynamics was continuously monitored with the BabyLux device during the euthanasia.

Analysis of experimental curves
Experimental curves were analyzed using the same flow explained for analysis of simulated curves ( Fig. 1 from right to left, in dark-gray). As explained in the device section, the shift of the TRS lasers was monitored independently for the TRS measurement. The measured shift was used to correct the position of the R ex p (t). Afterwards, in the fitting procedure only optical properties were used as fitting parameters, while t 0 was kept fixed at 0 ps. As previously mentioned, different options are available for DCS input analysis parameters. The time average of the optical properties measured by TRS for each experimental set was used for DCS analysis. This means that they were kept constant over the whole set of measurements. An exception was made for the liquid phantom, where, in order to highlight the effect of the variability of the estimates of the optical properties, they were inserted in the analysis by defining an a priori CV for them. The other input parameter for DCS analysis is β which was estimated from each experimental g exp 2 (τ) curve.

Estimation of precision and accuracy
The precision over the thirty simulated or experimental curves was quantified considering either the standard deviation σ(X), where X = [x 1 , x 2 , ..., x 30 ] with x i as each of the values obtained from the analysis, or as the coefficient of variation (CV= σ(X)/ X ), being X the mean value of X. For the simulations, a relative error was also calculated to quantify the accuracy, defined as = ( X − x true )/x true , where x true is the parameter value used for the simulation. As mentioned before, for the simulated curves, the generation of ten data sets of thirty curves for each setting allowed us to have a standard deviation of the variability (σ(X) or CV) and of the relative error ( ). (τ) and the 30 g sim 2 (τ) and g sim 1 (τ) generated in the simulation process. The dashed vertical lines highlight the fitting range. Figure 3 shows examples of TRS and DCS data generated in the simulation process, both the model curves and the simulated ones with realistic noise are reported.

Variability and error in simulated TRS data
CV of optical and hemodynamic properties decreases drastically with increasing total photon counts (Fig. 4). The decrease with increasing counts is also noticeable for the relative error, even though it is low (below 5%) for all the count levels considered. The variability (CV) of StO 2 depends not only on the count level but also on its own absolute level. Specifically, CV and standard deviation decrease with increasing StO 2 level (Fig. 5).
The next point focuses on the role of the different available choices for t 0 determination in the fitting process. TRS data was simulated with different levels of t 0 variability while analyzed either fixing t 0 at 0 ps or treating it as a fit parameter. As expected, the uncertainty of t 0 influences the µ s estimates more drastically than the µ a . If no correction for t 0 is considered in the analysis   ( Fig. 6, left column), variability is high even for high count levels, if t 0 has a non zero standard deviation. On the other hand, considering t 0 as a fitting parameter (Fig. 6, right column) leads to a high variability for lower counts but it provides a good correction for the t 0 variation and leads to low CV for high count levels. The relative error for estimated optical properties (not shown) is around or lower than 1% for counts higher than 10 4 for both configurations of the analysis for t 0 .
On the other hand, µ a is obtained with an error of 20% and µ s of 50% at 10 4 if t 0 is treated as a fitting parameter and of about 2% if t 0 is kept fixed. The role of an error in water estimation is also investigated. StO 2 variability is found to remain  stable when the water content is either over-or under-estimated, while the error in StO 2 increase (Fig. 7 for 10 5 and 10 6 counts).
To conclude with, we have considered an a priori estimation of the physiological variability for hemodynamics and generated data at different count levels. At about 10 5 counts, the results for CV in the hemodynamic parameters match the input values (Fig. 8). The a priori 3% CV of Hb and 7% CV of HbO 2 are translated into a CV of 4% for both µ a and µ s at 10 5 counts and of 3% at 10 6 counts.

Variability and error in simulated DCS curves
As before, let us begin by looking at the pure effect of the number of detected photons. By increasing the averaging time and the count rate, CV and relative error are found to decrease (Fig. 9). If we focus on 100 kHz for count rate and 1 s for averaging time and change the absolute value of BFI, the standard deviation increases and the CV decreases at higher BFI values (Fig. 10).
The next point in the analysis is that of the effects of the input parameters. An a priori simulation of the variability for optical properties was considered in the analysis. CV of input µ s is drastically reflected in CV of the BFI, while the effect of a variability in µ a is milder (Fig. 11, right column). On the other hand the effect on the relative error is minimal (Fig. 11, left  column).
β is the other input parameter for the analysis. If β is kept fixed, an error in its assumption does not compromise CV of BFI while it does for the relative error of BFI (Fig. 12). Estimating β from each curve leads to higher CV and relative error when the count rate is not larger than, e.g., 50 kHz (Table 2).
To conclude with this section, a simulated a priori physiological variability for BFI is reflected in a comparable value, if count rate is larger than 50 kHz, while it is not reflected on the relative error, which depends only on the count rate (Fig. 13).

Variability in experimental TRS and DCS data
As described above, tissue simulating phantoms were measured by the BabyLux device to estimate the variability.
Thirty TRS measurements were acquired on a solid phantom at the three wavelengths with at approximately 10 5 total photon counts to calculate its optical properties and their CV (Table 3).
Analogously, thirty DCS measurements were acquired from a liquid phantom, with µ a of 0.17 cm −1 and a µ s of 7 cm −1 , as measured by the TRS, with a count rate of 150 kHz. They  Fig. 11. CV and relative error of BFI estimation from simulated data generated considering a certain CV for the optical properties in the data simulation, accordingly to the x and y axis. Count rate is 100 kHz and averaging time 1 s. were analyzed considering different levels of CV for input optical properties. The estimated BFI values were consistently increasing in variability with increasing variability of the input optical properties (Table 4).
In the experiments involving the piglet model, thirteen (N=13) piglets with a median weight of 3.2 kg (range 3.0 -3.4 kg) and median age of 11 days (range 9 -12 days) were measured.  All TRS measurements were acquired at 10 5 counts in a approximately 1 s integration time. While, the median of the count rate for DCS was 150 kHz, with a minimum value of 100 kHz and a maximum value of 200 kHz. The optical and hemodynamicproperties, and their CV, were retrieved and averaged over the sample population (Table 5 and Table 6).
The last experiment focused on registering the variability over a range of values changing dynamically in time by following changes after the euthanasia of the piglets. This measurement was performed for eleven (N=11) piglets. In order to quantify the variability during these dynamic changes, we have opted for selecting ten second windows of measurements taken at one second resolution, i.e. ten data points instead of thirty as was done with the simulations. We have opted  for this because choosing thirty second time windows with thirty data points would have resulted in a too variable of hemodynamics in each time window due to the rapid changes following euthanasia. The CV for BFI and StO 2 over ten seconds of measurement during the euthanasia decreased with increasing BFI and StO 2 absolute values (Fig. 14).

Discussion
The BabyLux device was developed and built in order to propose a solution for continuously and at the cot-side monitoring the cerebral hemodynamics of the preterm newborns. This work aimed, on one hand, at exploring limits for accuracy and precision of the BabyLux device in the determination of optical and hemodynamic properties and, on the other hand, at defining how different choices implemented in the analysis process could affect the results. To this purpose, TRS and DCS data with added noise were simulated and analyzed with different settings in order to investigate the influence of each parameter in the precision and the accuracy of the retrieved results. We have used reference optical and hemodynamic parameters that were in line with results of the infants studies conducted with the BabyLux device [49]. The experimental parameters were chosen accordingly to the characteristics (e.g. temporal resolution, wavelengths, noise level) of the BabyLux device. This same device was used for acquiring experimental data from phantoms and from an animal model in order to test its performance in terms of precision and compare it with the results of simulations. It must be highlighted that the limits in precision and accuracy derived from simulated data were not due to the model and to the discrepancy between the system and the diffusion approximation since simulated curves were analyzed with the same model used for their generation. In this way we could isolate the role of photon noise and of each of the effects listed in the left column of Fig. 2. We must highlight that the device has been previously tested in laboratory settings. In particular the MEDPHOT protocol [50] was used to assess its performances for optical properties estimation. Among other tests, we calculated the accuracy towards conventionally true values and precision at different count levels. Results are reported elsewhere [48].
First objective was to define the role of the count rate and the averaging time for DCS and of the total photon counts for TRS. For the latter, it is fundamental to gather enough counts in order to have accuracy and low variability in obtained optical and hemodynamic properties. The simulations showed that a substantial increase in precision could be obtained by increasing the count level from 10 5 counts (currently used in the BabyLux) to 10 6 . It has been highlighted that a standard deviation of 5% measured for current clinical, commercial devices in replacing the probe on an infant head is too high for clinical purpose [52]. If the choice of TRS aims at improving this results, it is important to reach 10 6 total counts during measurements. Differently from TRS, the accuracy and the precision of BFI estimation by DCS did not improve substantially by increasing the averaging time, especially at higher count rates. This is due to the fact that the probability to detect pairs of photons within each bin time of the correlator does not increase linearly by increasing the averaging time [42].
The lower limit to precision was found to be dependent on the absolute values of StO 2 and BFI. The simulated data suggested that the variability was not constant over the values of those two parameters. In particular, the CV increased for lower values of StO 2 and BFI. Measurements over wide ranges of physiological values have to face this limit. We decided to verify whether this behavior was confirmed by in vivo measurements. To this purpose we have monitored the variability of measurements on piglets during euthanasia when the absolute values of StO 2 and BFI decrease. As shown, the CV for StO 2 and BFI increased while lowering the absolute values of the two, confirming the results of simulation. This needs to be further explored and accurate confidence intervals must be established when reporting values in the ultimate clinical use.
The accuracy limit in the determination of optical properties by TRS was previously studied [30,32] and it was demonstrated that better performance could be achieved with samples of high reduced scattering and low absorption coefficient. In those works, the simulated curves were generated with a similar procedure as in our study. The described behavior of accuracy changing with optical properties was not then due to the limits of the diffusion approximation which is better when µ s µ a but rather due to the shape of the collected curves in this regime. It must be noted that the standard deviation of StO 2 did not vary as much as CV with the absoute StO 2 level, therefore the high CV values at low StO 2 were mainly due to the fact that the denominator in the CV calculation was small. The higher standard deviation at lower StO 2 value is due to the higher absorption coefficients compared to the one calculated for the higher values of StO 2 , which is consistent with the results in Ref. [32].
DCS being a newer technique is not so thoroughly studied. Previous studies have highlighted how the accuracy of BFI parameter depends on the accuracy of the optical properties introduced in the analysis. It has been shown that especially the error in µ s is dramatically reflected in an error in BFI [33]. In addition simulation studies have shown that it is not reliable to estimate more than one parameter from a single DCS curve and that the accuracy of BFI is improved if β parameter is estimated by more than one point at short delay times [34]. Considering the existent literature we have focused our work on studying how these and other parameters influence the precision of DCS. By looking at the BFI variability over the BFI absolute value, we have observed a higher standard deviation over multiple measurements at higher BFI values. This might be expected considering that the higher flow causes a faster decay of the DCS curve where the noise is larger due to the smaller earlier bin of the delay times [41]. This might be remedied by improving the count-rate and by optimizing the correlator design.
We have intended to simulate the eventual role of the current hardware components in the variability of results. Specifically, we have tested the influence of an unstable peak of each TRS pulse. This was summarized in the parameter t 0 which is usually determined according to how the IRF is acquired. This is often not concurrent to the TRS measurements but before or after the measurement session. In Ref. [32] it was demonstrated that an uncorrected shift of the peak of the spectrum resulted in large errors in the determination of the optical properties. For example, errors of 5% and 10% for µ a and µ s were obtained with a shift of 20 ps.
Here we have decided to not consider a shift in the t 0 since the Babylux device allows for the monitoring the position of the TRS peaks by reference branches concurrently to the measurement. This is done by directing a small amount of the power of the laser to the detector without propagating into the tissue. This works well for large and sudden change of the laser pulse but even the corrected peak position has a certain noise. A range of about 10 ps wide was measured by the BabyLux device for the corrected peak. Therefore, we have simulated variabilites of t 0 of 5, 10 and 15 ps. We have showed that this compromised the results for the variability of parameters even at high count rates if no corrections were applied. A readily available option was to consider t 0 as a fitting parameter. This option made the fitting procedure less stable, since we have increased the number of parameters to be fitted. Consequently, it gave good results only for high (> 10 5 ) count levels.
We have proceeded by analyzing the effect of input parameters of the analysis. An a priori estimation of water content is required for the analysis of the BabyLux data because the wavelengths that are used are not in a range suited for water determination. This is what is commonly done in TRS (and CW-NIRS) devices [53]. We have demonstrated that an error in the water estimation did not compromise the variability of StO 2 but of course increased the relative error.
DCS needs optical properties as input parameters. As already mentioned, it is not desirable to estimate both BFI and optical properties by fitting a single DCS curve [34] and that multi-distance or multi-wavelength measurements must be acquired for this [54,55]. One of the motivations of building and using a hybrid device was that optical properties could be simultaneously measured by TRS and used for the DCS analysis. This allowed, in principle, a better estimation of the absolute values of BFI. The drawback of using measured optical properties is that the error and the variability of TRS results is propagated into DCS results. We have shown that the CV of µ s was drastically propagated to BFI estimates while the effect of µ a was smaller. This is easily explained by eqn. 2 that shows the coupling between BFI and µ s 2 in the decay rate of the auto-correlation curve. Luckily, a change of µ s is not expected during hemodynamic changes because the bulk scattering properties of the tissue do not depend on blood flow since red blood cells account only for a tiny amount of the total scatterers. We note that the changes in µ s detected by TRS are often usually due to cross talk between µ a and µ s in TRS analysis. For this reason µ s is usually kept constant over measurements in DCS analysis and sometime in TRS analysis as well [56].
Another choice to be made in DCS analysis is the β estimation and whether to keep it fixed or to allow it to vary during the experiment. For low count rates (less than 50 kHz), it was clear that it was not convenient to estimate it from a single curve since the first part of the DCS curve is compromised by the noise and it would result in high error and variability of BFI. It is reasonable to keep β constant for an experiment as long as the laser source is known to be stable in coherence and that there is no significant leakage of external light.
The last point was to check how the physiological variability of a given hemodynamic parameter propagates into measurements and to the determination of the hemodynamic parameters by DCS and TRS. The single level of Hb and HbO 2 concentration considered in the simulation were comparable to those that were obtained on healthy term neonates by the BabyLux device [49]. At 10 5 acquired photon level, the obtained CV of Hb and HbO 2 concentration were already close to the ones used for simulation. It must be noted that this a priori variability of Hb and HbO 2 concentration (3% and 7% respectively) is translated into variability in the optical properties (4% at 10 5 counts), which can ultimately influence the blood flow index estimation, as shown in Fig. 11 and discussed. Analogously to what done for TRS data generation, different levels of blood flow variability was used for DCS. The CV of BFI results was comparable to the input values for count rates larger than 50 kHz. This provides an estimate of the realistic expectations from a device like the BabyLux.
Overall, these results define a lower value for CV and error, valid in the ideal scenario, considering the experimental parameters used for the real measurements (count rate, averaging time, temporal instabilities of TRS lasers and others). In addition, these results can be a guidance on how to choose parameters in the analysis process in order to maximize precision and accuracy. Surely, the picture here drawn is incomplete, since the error due to employing the diffusion approximation has not been considered. Nonetheless, we expect that this would minimally influence the precision of DCS and TRS measurement, while it would be more heavily reflected in accuracy. In spite of this, the conclusion on how to set the analysis process and how to select experimental parameters in order to obtain the best results are valid.
We conclude by discussing the results from experimental measurements on phantoms and piglets, performed to be compared with the simulation results. The variability of the experimental estimates of the phantom optical properties were comparable to those obtained from the simulated curves. This was an indication of the good quality of the BabyLux device in terms of precision. On the other hand, the liquid phantom results gave higher variability (4.7%) compared to the simulated curves in DCS measurements (about 2%, Fig. 9). This slight increase is due to uncontrollable and variable properties of the liquid phantom. For example, the fluctuations inside the liquid or subtle vibrations due to the placement of the probe on the Mylar window to the phantom [51] together with a non-uniform temperature could increase the noise of the measurement.
Good results were obtained for the precision of the optical property estimates from the piglet model which were equivalent to the simulation and phantom results. This was true also for HbO 2 and Hb concentration and StO 2 , especially considering the absolute level of the blood oxygen saturation, lower than what was previously measured in infants [21,29] probably because the piglets were anesthetized. Variability of BFI in piglets is higher than simulation results but comparable to phantom measurements. This could be improved by increasing the signal to noise ratio, for example adding more detectors shining in the same position.

Conclusion
We have carried out a comprehensive study of the precision and the accuracy that could be expected due to physical properties of photon counting experiments in diffuse optics. We have then compared these results to that of a newly introduce hybrid device combining near-infrared time resolved (TRS) and diffuse correlation spectroscopies (DCS) from tissue simulating phantoms and from an animal model. The methodology and the results provide quantitative and detailed means for evaluating both the existing devices and future one based on new components and designs. The BabyLux device is shown to be a promising system to achieve improved precision and accuracy in the setting of neonatal, critical care neuro-monitoring.  No financial conflicts of interest were identified. Udo Weigel is the CEO, has equity ownership in HemoPhotonics S.L. and is an employee in the company. His role in the project has been defined by the project objectives, tasks and work-packages and was reviewed by the European Commission.