Time-domain NIRS system based on supercontinuum light source and multi-wavelength detection: validation for tissue oxygenation studies

We present and validate a multi-wavelength time-domain near-infrared spectroscopy (TD-NIRS) system that avoids switching wavelengths and instead exploits the full capability of a supercontinuum light source by emitting and acquiring signals for the whole chosen range of wavelengths. The system was designed for muscle and brain oxygenation monitoring in a clinical environment. A pulsed supercontinuum laser emits broadband light and each of two detection modules acquires the distributions of times of flight of photons (DTOFs) for 16 spectral channels (used width 12.5 nm / channel), providing a total of 32 DTOFs at up to 3 Hz. Two emitting fibers and two detection fiber bundles allow simultaneous measurements at two positions on the tissue or at two source-detector separations. Three established protocols (BIP, MEDPHOT, and nEUROPt) were used to quantitatively assess the system’s performance, including linearity, coupling, accuracy, and depth sensitivity. Measurements were performed on 32 homogeneous phantoms and two inhomogeneous phantoms (solid and liquid). Furthermore, measurements on two blood-lipid phantoms with a varied amount of blood and Intralipid provide the strongest validation for accurate tissue oximetry. The retrieved hemoglobin concentrations and oxygen saturation match well with the reference values that were obtained using a commercially available NIRS system (OxiplexTS) and a blood gas analyzer (ABL90 FLEX), except a discrepancy occurs for the lowest amount of Intralipid. In-vivo measurements on the forearm of three healthy volunteers during arterial (250 mmHg) and venous (60 mmHg) cuff occlusions provide an example of tissue monitoring during the expected hemodynamic changes that follow previously well-described physiologies. All results, including quantitative parameters, can be compared to other systems that report similar tests. Overall, the presented TD-NIRS system has an exemplary performance evaluated with state-of-the-art performance assessment methods.


Introduction
Time-domain near-infrared spectroscopy (TD-NIRS) uses short pulses of light to non-invasively and in real-time monitor tissue optical properties, i.e. absorption (µ a ) and reduced scattering (µ ′ s ) coefficients, which are related to the concentrations of chromophores in tissue, e.g. oxyhemoglobin (HbO 2 ) and deoxyhemoglobin (Hb). TD systems acquire the highest amount of information, for a defined number of source-detector pairs, compared to continuous-wave and frequency-domain acquisition modes [1]. The monitored parameters, e.g. cerebral metabolism [2] (related to changes of chromophores concentrations) or tissue oxygenation [3] (related to absolute values), help diagnose physiological conditions and monitor results of clinical interventions. NIRS is safe and non-ionizing, which makes it suitable for studies on neonates [4][5][6] and on adults [7].
A number of recent publications review instrumental developments for various applications of NIRS, including the past and the current status of TD-NIRS [8], clinical brain monitoring with TD-NIRS [1], metabolic brain measurements on neonates [9], neuromonitoring for neonatal encephalopathy [10], functional NIRS [7], instrumentation for functional NIRS [11], and applications related to stroke patients [12], studies on muscles [3], and evaluating the level of consciousness typically after a brain injury [13]. Continuous advancements in optical technologies and electronics for TD systems are leading to better measurements and a wider range of applications. The advancements have enabled TD-NIRS measurements at larger source-detector separations, e.g. 9 cm [14], at longer wavelengths, e.g. in the range from 600 to 1350 nm [15], with more sources and detectors, e.g. 16 sources and 8 detectors with fibers [16], or 1032 detectors [17], or 36288 detectors [18], with fast acquisitions, e.g. wavelength switching at 160 Hz [19], more compact, e.g. 200 × 160 × 50 mm 3 [20], wearable [22,21] and also wireless [23], or with an improved light harvesting capability that allows collecting more photons than restricted by the pile-up limit [24]. Another direction of advancements involves integrating various imaging techniques, e.g. TD-NIRS and diffuse correlation spectroscopy [25,26].
Measurements at more wavelengths allow estimating absolute concentrations of more chromophores, e.g. lipid, water, and collagen [15], in addition to HbO 2 and Hb. Furthermore, to follow dynamic changes associated with physiology requires fast acquisitions on a time scale of typically 1 s. Of particular interest is monitoring concentration changes of cytochrome-c-oxidase (CCO) enzyme, which requires many wavelengths as well as fast acquisitions [2]. The importance of monitoring changes in CCO, which relates to tissue metabolism, and the relevant advances in optical technologies were reviewed in [9]. Recently, Lange et al. [19] presented an assessment study of a multi-wavelength TD-NIRS system that was designed for monitoring changes in CCO in addition to HbO 2 and Hb. The system uses a broadband light source coupled with two acousto-optic tunable filters that allow switching between many different wavelengths. Other common approaches for multi-wavelength TD-NIRS are time multiplexing [27] and space multiplexing [16]. In the first approach, multiple pulsed lasers with different wavelengths are shined in parallel and the optical pulses of different lasers are delayed relative to each other such that they arrive with different time offsets to a single TCSPC card and get recorded in the same TAC, which allows a limited number of wavelengths, typically two. In the second approach, multiple pulsed lasers are sequentially shined using an optical switch, which allows employing many wavelengths at the cost of acquisition time. Renna et al. [27] presented such system with 8 pulsed diode lasers, but recommend sequential scanning of no more than 2 or 3 wavelengths for monitoring dynamic processes.
In the present work, the multi-wavelength TD-NIRS system avoids switching wavelengths. The principle of this kind of instrument was introduced in [28]; the presented system is an advanced version and with modern electronics. We report the most relevant performance tests from three established protocols for photon migration instruments (MEDPHOT [29], BIP [30], and nEUROPt [31]), which were intended for quantitatively assessing the performance of timedomain instruments that use pulsed laser sources, single-photon detectors, and time-correlated single-photon counting electronics. Measurements on 32 homogeneous phantoms and a solid inhomogeneous phantom [32] were carried out as part of the BitMap campaign [33], which is a multi-laboratory effort to standardize the performance assessment of instruments for diffuse optics. Furthermore, the ability to estimate absolute concentrations of HbO 2 and Hb, and hence oxygen saturation, is assessed on a series of blood-lipid phantom measurements using the phantom and the protocol developed by Kleiser et al. [34,35]. Lastly, in-vivo monitoring during hemodynamic changes is demonstrated using measurements on the forearm of three healthy volunteers during arterial and venous cuff occlusions. This study is intended to quantify the capabilities of the system and support its clinical use for diverse applications.

Overview
The system measures the distributions of times of flight of photons (DTOFs) in parallel for 16 wavelengths with a sampling time of 0.3 s. Two emission fibers allow delivering light at two locations. Two detection fibre bundles allow measuring at up to two source-detector distances (ρ). Therefore, the system provides 32 DTOFs (16 wavelengths and 2 detection modules) at a sampling rate of up to 3 Hz. A previous version of the system was used for hemodynamic studies involving the injection of an optical contrast agent Indocyanine Green [36,37]. The system was designed for the clinical environment with the possibility to continuously monitor patients before, during, and after surgery. The current fiber holder design is suitable for adults and can be potentially adjusted for neonates.

Housing setup
A photo of the arrangement of the system is shown in Fig. 1 (a) and a schematic diagram of the inner components is shown in Fig. 1 (b). All components were mounted inside a 19-inch metal rack case with four wheels, which is a commonly used housing setup for TD-NIRS systems because it is compact, portable, and robust [19,38]. An uninterruptible power supply unit (APC Smart UPS SMC1500I-2U) protects against interruption of external power supply and also allows the system to operate without an external power supply for a brief period.

Emission module
The results in this study were obtained with fiber laser SC480-8 WhiteLase (Fianium, UK), which operates at 80 MHz repetition frequency. The fundamental pulse width is 6 ps, the maximum output power is up to 8 W, and the spectral bandwidth range is from 460 to 2000nm.
(1) The output beam of the laser first passes through a system of two polarizers with a half-wave plate in-between, which serves as a variable attenuator for unpolarized light. For a chosen output power of the laser, the half-wave plate is manually rotated and then secured with a screw, such that the power of light exiting any of the two emission fibers is always lower than 20 mW, which equates to a power density of less than 2 mW/mm 2 on the surface of the skin due to angular spread and distance, which is in line with the safety regulations for measurements on patients. Next, the beam passes through two edgepass filters, which remove the wavelengths that are outside of a desired spectral range. The current chosen spectral range is from 650 to 850 nm, which is achieved with a longpass filter (FEL0650 with a cut-on wavelength of 650 nm, Thorlabs, Sweden) and a shortpass filter (FES0850 with a cut-off wavelength of 850 nm, Thorlabs, Sweden). The two filters allow some unwanted light above 1200 nm to be emitted, but this light does not pass through the polychromator on the detection side.
(2) A beam splitter directs a fraction of the laser light to a reference photodiode (PHD 400, Becker & Hickl GmbH, Germany), which feeds a signal into the SYNC input of each TCSPC card. The details of using reference photodiodes can be found in the TCSPC handbook, pages 76 and 202 [39].
(3) The rest of the light passes through one of four openings in a rotatable wheel: i) an unimpeded opening that lets all light pass through, ii) a solid wall that blocks all light, iii) a band-pass filter centered around 780 nm (FL780-10, Thorlabs), which is used to calibrate diffraction gratings of the two polychromators, and iv) a bandpass interference filter (750 nm with FWHM of 50 nm, Edmund Optics), which is used for the excitation of the ICG in the fluorescence mode measurements.
(4) Final filtering is achieved with the use of a continuously variable neutral density filter (NDC-50C-2-B, Thorlabs), which removes a variable amount of excessive light power.
(5) A beam splitter directs the filtered light into two fiber couplers that guide the light into two optical fibers, which deliver the light to the surface of a medium, e.g. tissue or phantom, at two chosen locations. The parameters of each emission fiber are: step index, 0.39 NA, 2 m length, and 400 µm diameter (FT400EMT, Thorlabs).

Detection module
After travelling through a diffusive medium, the remitted light is collected by two separate detection fiber bundles (2 m length, 0.22 NA, and a 90°bent tip on the patient side with a 30 mm long-term bending radius, custom-made, Ceram Optec), which are terminated with a circular tip (diameter 3.6 mm) on the medium side and a line-shaped tip (1.4 × 7.3 mm) that fits into the slit of a polychromator; illustrated in the TCSPC handbook, page 171 [39]. Each fiber bundle guides light to one of the two separate detection modules (PML Spec, Becker & Hickl GmbH, Germany). Each detection module is comprised of two parts: a polychromator (MS125, Oriel Instruments, USA), which has 0.135 NA and uses 77414 diffraction grating (Grating Groove Density is 600 lines / mm), and a multi-anode photomultiplier tube detector (PML-16-1-C, Becker & Hickl GmbH, Germany). The power to both photomultiplier tubes is provided by a detector control card (DCC-100, Becker & Hickl GmbH, Germany), which also allows controlling gain parameters and overload shutdown. The width of each spectral channel is about 12.5 nm, which can be changed by using another diffraction grating, e.g. 35, 18, 9, or 4.5 nm as in [28]. The spectral range, which is a parameter of a polychromator, was set from 674 to 874 nm. Consequently, the central wavelengths of spectral channels were: 680, 692.5, 705, 717.5, 730, 742.5, 755, 767.5, 780.5, 792.5, 805, 817.5, 830, 842.5, 855, and 867.5 nm. For a further description of the detection module, we refer to the instrumentation section in the previous version of the system [36].
The output electrical signals of two PML-16-1-C are passed to two 16-channel TCSPC cards (SPC 150, Becker & Hickl GmbH, Germany). The synchronization between the acquisitions of the two TCSPC cards was achieved using the PCI Express I/O card (PCIe-6321, National Instruments), which was configured with LabView software (National Instruments). The card generates a trigger signal at 3 Hz for the start of the recording of the DTOF. The two TCSPC cards were used in continuous flow mode with a collection time of 0.3 s. During data analyses, the collection time can be increased by summing consecutive measured DTOFs. The parameters of a TCSPC card were configured such that each DTOF contained 1024 time bins with about 9.77 ps width each. For typical measurements, a temporal resolution of a recorded DTOF must be on the order of 10 ps to quantify the optical properties of diffuse media with accuracy better than 5% using a curve fitting method [40].

Software
All components were operated using the software provided by the manufacturers and run on a PC with Windows 7 Professional. The laser functions were controlled using the Fianium Laser Controller software. The DCC card and the two TCSPC cards were operated using the TCSPC Package 6.5.0 software (Becker & Hickl GmbH, Germany), which allows to adjust all TCSPC-related parameters, record the data, and assess the quality of the data before, during, and after a measurement. All data were saved and later processed offline using custom codes written in MatLab R2020b.

Optode holder
The two emission fibers and the two detection fiber bundles can be attached to the surface using custom 3D printed optode holders with holes at desired ρ (typically a few cm apart). The chosen printing materials were pure black ABS filament (Zortrax) for phantom measurements and black Fiberflex 40D filament (Fiberlogy) for in-vivo measurements, which is a more flexible material and achieves better contact with the skin. For phantom measurements, the optode holder was held in place using a support stand and covered with tape. For in-vivo measurements, the optode holder was fixed in place using Velcro strips and covered with a cohesive bandage.

BIP protocol
The basic instrument protocol (BIP) [30] allows the assessment of parameters that are related to the hardware components of a time-domain instrument. In particular, we report the instrument response function (IRF), the differential nonlinearity (DNL), the warm-up, and the temporal fluctuations.
To measure the IRF [41], an emission fiber was aligned with a detection fiber bundle at a distance of 6 cm with a neutral density filter between them and two pieces of white paper covering the fiber and the fiber bundle. The DNL was assessed by recording the response to a battery-operated light-emitting diode, which produced a continuous light signal. Repeated measurements were summed to obtain more than 10 5 photon counts per time bin for a good signal-to-noise ratio. The deviation in the number of photons detected in different time bins was quantified by calculating the peak-to-peak difference normalized to the mean value (ε DNL ) [30]. To assess the warm-up, the IRF was repeatedly acquired for up to 6 hours. Afterward, to assess the temporal fluctuations, the DTOF was repeatedly acquired on a homogeneous solid phantom (label B2 in the MEDPHOT protocol) at ρ = 3 cm. The first three statistical moments, i.e. the total number of photons (N), the mean time of flight (m 1 ), and the variance (V) [31], were calculated after background subtraction and using time bins that have photon counts higher than 1% of the maximum. Each DTOF contained about 10 6 photons, which was achieved by integrating repeated acquisitions. The temporal fluctuations of the moments were assessed using the standard deviation of 100 values.

MEDPHOT protocol
The MEDPHOT protocol [29] allows the assessment of a system's capability to recover µ a and µ ′ s in terms of linearity, coupling, and accuracy. The protocol includes measuring on 32 solid homogeneous phantoms that are labeled with a letter and a number. The phantoms make up all combinations of four nominal µ ′ s values and eight nominal µ a values (see the caption of Fig. 3). The base component of these phantoms is epoxy resin and we assumed the refractive index as 1.55. A measurement on each phantom was acquired for 30 s.
Four scatter plots were obtained: the measured µ a (or µ ′ s ) versus the nominal µ a (or µ ′ s ) and the measured µ a (or µ ′ s ) versus the nominal µ ′ s (or µ a ). Using the first two plots, the mean deviation in µ a (or µ ′ s ) from linearity was quantified by fitting a straight line and calculating the mean percentage difference between the measured data and the fitted line. This calculation does not use the values of nominal optical properties. Using the other two plots, the coupling of ∆µ a to ∆µ ′ s (or ∆µ ′ s to ∆µ a ) was quantified by fitting a straight line and calculating the slope, which requires using the changes of nominal values (assumed 0.05 cm −1 for ∆µ a and 5 cm −1 for ∆µ ′ s ).
The nominal values of µ a and µ ′ s of the 32 phantoms do not match the true values, and there is no established consensus yet on the agreed true values, as pointed out in [24]. Also, µ a and µ ′ s depend on wavelength, especially µ ′ s [27]. Therefore, to assess the accuracy, we present the obtained spectra of µ a and µ ′ s for three well-characterized phantoms: phantom B2, the solid inhomogeneous phantom, and the 1% aqueous solution of Intralipid (the two inhomogeneous phantoms are presented in the next section). The obtained values can be compared to the reported results of other systems that measured on the same phantoms at multiple wavelengths.

nEUROPt protocol
The nEUROPt protocol [31] allows the assessment of a system's sensitivity to a small localized µ a perturbation positioned at various depths and later positions. Measurements were repeated on two inhomogeneous phantoms (solid and liquid) that had similar µ ′ s but different µ a . Monte Carlo simulations were performed for a qualitative comparison.
The mechanically switchable solid inhomogeneous phantom was described in detail in [32]. The phantom has a cavity (diameter 14 mm) for a movable cylindrical rod that has a localized cylindrical µ a perturbation (we used the rod with the following perturbation: 5 mm in diameter and length, i.e. 98 mm 3 volume, and ∆µ a equivalent to 0.017 cm −1 ). The rod can be precisely moved by a step motor.
The liquid phantom was a 1% aqueous solution of Intralipid (Fresenius Kabi AG, Bad Homburg, Germany) poured inside a dedicated cell [42]. A localized µ a perturbation was introduced using a submerged small black PVC cylinder (4 mm in diameter and length, i.e. 50 mm 3 volume), which was fixed to an end of a 0.5 mm thin, rigid, white metallic wire, which was verified to have a negligible effect on TD-NIRS measurements [42]. The depth and the lateral position of the metallic wire can be manually changed with a precision of about 1 mm.
The reported optical properties of the two phantoms without a perturbation are: µ a ∼0.1 cm −1 (∼0.026 cm −1 ) and µ ′ s ∼8 cm −1 (∼10.5 cm −1 ) at 800 nm for the solid [32] (liquid [19,43]). The solid phantom mimics the typical µ a and µ ′ s of biological tissue. The liquid phantom, which contained only water and Intralipid and hence can be easily replicated, is used for a comparison with the results on the solid phantom. The refractive index was assumed as 1.55 for the solid and 1.33 for the liquid.
One emission fiber and one detection fiber bundle were fixed at ρ = 3 cm. Their positions on the solid phantom were illustrated in [32,38]. When varying the depth, a perturbation was positioned in the midplane between the emitter and the detector, and it was moved in steps of 1 mm (2.5 mm) in the solid (liquid). When varying the lateral position, a perturbation was 1.5 cm (1.2 cm) below the surface and it was moved in steps of 2 mm (2 mm) in the solid (liquid).
The position of a perturbation was defined as the center of a perturbation, as in [31,32]. A measurement was acquired at each position of a perturbation for 60 s (18 s) on the solid (liquid).
The contrasts that result from introducing a perturbation were calculated for five measurands: N, m 1 , V, and the number of photons in early and late time windows (from 0.5 to 1 ns for N Eearly and from 3 to 4 ns for N Late ). The zero time (0 ps) was defined at the maximum of the IRF. Measurands were calculated from DTOFs after background subtraction and using time bins that have photon counts ≥ 0.1% (2%) of the maximum for measurements on the solid (liquid) phantom. The contrasts were calculated as: where ∆A TW (x) is the change in the attenuation for a time window TW (N, N Early , and N Late ) due to a perturbation located at x, which is either the depth or the lateral position. The zero subscript represents the reference measure when a perturbation was furthest away from the emitter and the detector.
Monte Carlo (MC) simulations were performed to obtain the sensitivity factors for the first three statistical moments: the mean partial pathlength (MPP), the mean time of flight sensitivity factor (MTSF), and the variance sensitivity factor (VSF). The MC code was developed by Wojtkiewicz et al. [44]. The sensitivity factors relate changes in moments to localized changes in µ a (MPP = ∆A / ∆µ a , MTSF = ∆m 1 / ∆µ a , and VSF = ∆V / ∆µ a ) at different depths and lateral positions, and hence they can be used for a comparison with the measured contrasts. For a detailed description of the sensitivity factors and their use, we refer to [44,45]. The MC parameters were chosen as in [46]. The model parameters were set to mimic the solid and the liquid inhomogeneous phantoms, using their optical properties obtained at 805 nm (values listed in Table 2).

Blood-lipid phantoms
Blood-lipid phantoms allow the assessment of a system's capability to retrieve absolute concentrations of HbO 2 and Hb, and hence total hemoglobin (HbT = HbO 2 + Hb) and oxygen saturation (StO 2 = HbO 2 / HbT), during controlled dynamic StO 2 changes.
The presented system was transported to University Hospital Zurich, University of Zurich, in Switzerland, to perform measurements on their blood-lipid phantom in parallel with a commercially available oximeter (OxiplexTSTM, frequency-domain, 692 and 834 nm, ISS, Champaign, Illinois, USA). OxiplexTS provides the absolute values of µ a and µ ′ s at two wavelengths, the concentrations of HbO 2 , Hb, and HbT, and StO 2 . A rigid sensor with four source-detector separations (2.5, 3.0, 3.5, and 4.0 cm) was used and the sampling rate was set to 1 s. Co-oximetry was used as another reference measure of HbT, by analyzing blood samples from the erythrocyte bag using a blood gas analyzer (ABL90 FLEX, Radiometer Medical ApS, Brønshøj, Denmark) and calculating the degree of dilution in the phantom, which is the same method as used in [47].
A detailed description of the phantom can be found in [34,35] and the following is a brief overview. The phantom and a similar protocol were used in previous studies that also used OxiplexTS [48][49][50]. The phantom consisted of a black 3D printed cubic container with windows on four sides, which allows simultaneously measuring with up to four systems. The windows were made from silicone and had tailored properties to mimic a layer of fat: µ a = 0.063 cm −1 , µ ′ s = 6 cm −1 at 690 nm, and 1 mm thickness. The three main components that made up the liquid were: 2.5 L of phosphate-buffered saline (PBS, after Kreis, pH = 7.4, Kantonsapotheke Zurich, Zurich, Switzerland), a varying amount of Intralipid (Fresenius Kabi AG, Bad Homburg, Germany), and a varying amount of human blood from a human erythrocyte concentrate bag.
Fresh baker's yeast (3.0 g) was added to deoxygenate the phantom (StO 2 = 0%) and the rate was increased by adding glucose 50% (AlleMan Pharma GmbH, Reutlingen, Germany). Blood was fully oxygenated (StO 2 = 100%) by bubbling oxygen gas from an oxygen tank (typically for 1 min). Sodium bicarbonate buffer was added to keep the pH level close to 7.4. The temperature (37°C) and the mixing speed (500 rpm) were controlled with a hotplate that has a magnetic stirrer.
Measurements were repeated on two blood-lipid phantoms in which either blood or Intralipid was added at the end of each deoxygenation cycle before bubbling oxygen, i.e. when StO 2 = 0%, similar as in [50]. The first phantom contained 74 ml of Intralipid and four amounts of blood: 20, 35, 55, and 70 ml. The second phantom contained 45 ml of blood and five amounts of Intralipid: 25, 50, 75, 100, and 125 ml. The refractive index was assumed as 1.33. The optical properties resulting from 74 ml of Intralipid and 45 ml of blood mixed with 2.5 L of saline are close to the typical optical properties of a neonatal brain (µ ′ s ≈ 5.5 cm −1 ) [34].

In-vivo measurements
In-vivo measurements during cuff occlusions were performed to provide an example of a system's ability to monitor hemodynamic responses in tissue. It is an easy-to-repeat test that allows checking the performance of oximeters in-vivo [51]. Measurements were repeated on three healthy volunteers: skinfold thickness (2.1, 2.4, and 2.5 mm), age (25,26, and 29 years), one male. The measured blood pressure before the start of the experiment was around 120/80 mmHg. Written informed consent was filled in by participants.
A subject was relaxed in a Fowler's position on a medical bed. A probe holder was attached on the posterior side of the left forearm, over the extensor muscles about midpoint between elbow and wrist. An arm-intended pneumatic cuff was loosely placed around the upper arm and rapid inflation was achieved with an electro-pneumatic regulator (SMC ITV2010-31F2N3). The inflation time was a few seconds; the recommended upper limit for a calf muscle was is 6 s [52]. 250 mmHg pressure was applied to occlude arterial and venous flow [25,38,53,54], and 60 mmHg pressure was applied to partially restrict venous flow without obstructing arterial flow [52,55]. The pressure was applied for 2 min. The refractive index of tissue was assumed as 1.4 [56].

TD-NIRS data analysis
Most measurements were performed with two detection modules at two source-detector distances (ρ = 2 and 3 cm, or 3 and 4 cm) away from the same emission fiber. The acquisition time was 0.3 s and the sampling rate was 3 Hz. During post-processing, 3 consecutive DTOFs were summed to increase the signal-to-noise ratio, resulting in a sampling rate of 1 s. µ a and µ ′ s were estimated using the curve fitting method [57,58] as implemented in NIRFAST software [59]. It is considered a gold standard method to determine the absolute optical properties. A homogeneous semi-infinite medium model under extrapolated boundary conditions was used. The fitting region of the DTOFs was limited based on the percentage of its maximum, i.e. from 85% on the rising edge to 1% on the tail for phantom measurements (Section 4.2), and from 75% to 3% for blood-lipid phantoms and in-vivo measurements (Sections 4.4 and 4.5).
The concentrations of HbO 2 and Hb were calculated using the estimated µ a at multiple wavelengths and the Beer-Lambert law. For in-vivo measurements, another method was also used to calculate the changes in concentrations relative to a baseline, and the results were compared to the first method, as in [60]. The second method uses changes in light attenuation and the modified Beer-Lambert law (MBLL). The mean optical pathlength was calculated for each spectral channel using the method demonstrated by Delpy et al. [61], which relies on the measured mean time of flight and the speed of light in tissue. Zhao et al. [62] compared the hemoglobin absorption spectra from different sources and found that the most consistent performance was achieved with the values provided by Moaveni et al. [63], which were used in this study. We used the absorption spectra of water provided by Matcher et al. [64]. The concentration of water was assumed constant and the contribution of water's µ a for the estimation of absolute concentrations was addressed as follows. The concentrations of three chromophores (HbO 2 , Hb, and water) were estimated in the first step. A mean value of water concentration over the whole measurement was calculated and its contribution to the measured µ a was subtracted. Finally, the concentrations of only two chromophores (HbO 2 and Hb) were estimated using the water-subtracted µ a . Resolving for concentrations of only two chromophores leads to less noisy results. The concentrations were calculated using 11 spectral channels that covered the spectral range from 705 to 830 nm.

Results and discussions
Results were similar for both detection modules and hence some results are shown for only one.  Figure 2 (c) shows the IRF at 805 nm and the DTOFs measured on two inhomogeneous phantoms that were used for the nEUROPt protocol (before introducing a perturbation). Similar DTOFs for the same solid phantom were reported in [25]. The IRFs match the reported IRF for the same PML detector [39] since the FWHM is close to 150 ps and there is a small afterpulse peak (1 to 2% of the maximum) about 1 ns after the main peak. It is common for conventional PMTs to have such afterpulse peaks, which are probably due to photoelectrons that are reflected at the first dynode but then forced back to the dynode by the electric field, page 169 in [39]. The reported FWHM of multi-wavelength TD-NIRS systems are typically broader, e.g. 465 ps [19], 370 [38], 230 ps [53], and 83 to 263 ps depending on the detector [15], at around 800 nm. The narrowest reported IRF (83 ps) was obtained with a detector that has a much smaller active area (0.1 mm diameter) and hence a much lower responsivity. The ε DNL was found to be 3% in the interval between 1 and 10 ns; the maximum deviation of the photon count was 1.5%. This is consistent with the previously reported measures of ε DNL , e.g. 6.77% [16] and 3.5% [24]. The warmup period was found to be between 30 to 60 min, during which time the parameters of the IRF can change drastically. Afterward, long-term drifts can be observed during the next few hours (up to about 5 hours), which can be accounted for by periodically measuring the IRF (e.g. every 30 min). The duration of most clinical measurements is typically less than 1 hour, in which case, after a 30 to 60 min warmup period, it is sufficient to measure the IRF at the start and at the end, i.e. before and after attaching probes to a patient. The two IRFs can be used in data analysis to remove drift in N, m 1 , and V, by assuming that the drift was constant over time. Also, the first and second IRFs can be used for analyzing the first and second half of the measurement, respectively. A more rigorous approach to account for drifts involves directing a small portion of the laser beam to the detector along with the signal from the tissue [15], which can be adapted in the next version of the system.

BIP protocol
The calculated temporal fluctuations of moments, which were measured on phantom B2, were about: ±0.2% for N, ±1.0 ps for m 1 , and ±700 ps 2 for V, at 805 nm. The values at other wavelengths were similar. These reported standard deviations contain the photon noise, which depends on the count rate and shape of the DTOF, and the instrumental noise, which depends on the system [65]. Figure 3 shows the four scatter plots that illustrate the linearity and the coupling of µ a and µ ′ s (for ρ = 2 cm and 805 nm). Qualitatively, the scatter plots illustrate an exemplary level of linearity and coupling, for both µ a and µ ′ s . The same was observed for ρ = 3 cm and for other wavelengths (not shown). Important to note, for measurements at ρ = 3 cm, the signals were too low for some of the more absorbing phantoms and hence phantoms with nominal µ a ≥ 0.25 cm −1 (labels 6 to 8) were excluded. For measurements at ρ = 2 cm, 31 phantoms were used, excluding only phantom C8. The power of the laser was within the safe limit for measurements on patients (2 mW/mm 2 ), although it could be increased for phantom measurements.

MEDPHOT protocol
The mean deviations from linearity and the coupling are reported in Table 1 for two ρ (2 and 3 cm) and two wavelengths (705 and 805 nm). The results are presented for each subset of the 32 phantoms, which provides more information about the system's performance and also allows comparing with other systems that measured on only a subset of phantoms. Measurements on phantoms with the lowest nominal µ ′ s (5 cm −1 , label A) and lowest nominal µ a (0 cm −1 , label #1) have the highest deviations from linearity for ρ = 2 cm but not so much for ρ = 3 cm. When estimating low values of optical properties, even small deviations have a high magnitude in units of percentage. Increasing ρ reduces the uncertainty due to the finite width of time bins [40]. Furthermore, the deviation from linearity will be lower if using fewer phantoms since a straight line is easier to fit with fewer data points, while for the same reason the coupling of ∆µ a to ∆µ ′ s will be worse if using fewer phantoms. This is observed when comparing the results for ρ = 2 cm (31 phantoms were used) and ρ = 3 cm (20 phantoms were used). The coupling is similar for different phantom sets, except for those that had a low signal level (at ρ = 3 cm on phantoms with nominal µ a = 0.20 cm −1 , label 5). The coupling could be in part due to the theoretical model that was used for the fitting method.
The overall measure of linearity and coupling can be taken as the mean of different phantom sets. For ρ = 2 cm and 805 nm, the mean deviation from linearity is 4.3% for µ a and 1.4% for µ ′ s , and the mean coupling is 2.8 for ∆µ a to ∆µ ′ s and 2.2 × 10 −4 for ∆µ ′ s to ∆µ a . In other words, if µ ′ s increases by ∼1 cm −1 , the system and the used method will show an increase also in µ a by ∼2.2 × 10 −4 cm −1 . The obtained linearity is consistent with previously reported systems, e.g. mean deviation of 4.2% for µ a and 4.9% for µ ′ s [29]. The linearity and coupling reported in [33] were calculated using the median instead of the mean of different phantom sets, which lead to smaller values of linearity (as low as 1% for µ a and 1.5% for µ ′ s ) although similar values of coupling (between 1 and 10 for ∆µ a to ∆µ ′ s and between 2 × 10 −4 and 15 × 10 −4 for ∆µ ′ s to ∆µ a for the 29 instruments that were enrolled in the BitMap exercise). The obtained coupling is consistent with the best performance of the reported TD-NIRS systems.
The calculation of linearity requires an array of phantoms with linear increments of µ a and µ ′ s . The calculation of coupling requires also a rough estimate of the step size of the increments. Therefore, similar results could be obtained using e.g. a liquid phantom as in [66], where µ a and µ ′ s were linearly increased by adding ink and Intralipid, respectively. The assessment of linearity and coupling can be enhanced if the true optical properties of the 32 phantoms get established.   Here, we report the system's accuracy of retrieving the absolute µ a and µ ′ s . The values shown in Fig. 3 match the values in similar scatter plots in [23,25,27,67] for a similar wavelength. The obtained values for ρ = 3 cm are similar and one example is shown in Fig. 4 (a). However, due to a lack of true values for the 32 phantoms, we focus on the three well-characterized phantoms that are commonly used for assessing a system's performance. Figure 4 shows the spectra of the obtained µ a and µ ′ s for the three phantoms, and Table 2 reports these values for two ρ along with the previously reported values around the same four wavelengths. Note, most referenced values were roughly estimated from the reported spectral figures for the chosen four wavelengths. The results for phantom B2 match with the previously reported spectra [27][28][29] since µ ′ s is close to 10 cm −1 and decreases with wavelength and µ a is close to 0.07 cm −1 at all wavelengths. The obtained values of µ a , and especially µ ′ s , are slightly lower than the previously reported, but consistent with reproducibility of modern TD-NIRS systems, e.g. maximum day-to-day variations with respect to the mean: -5% to +6% for µ a and -4% to +3% for µ ′ s [25]. Similarly, the results for the solid inhomogeneous phantom match the spectra reported in [32], although now the obtained values of µ ′ s , and especially µ a , are slightly higher than the previously reported.   Figure 4 (c) shows the results for a 1% aqueous solution of Intralipid alongside the reported µ a of water. The shape of the obtained µ a spectra matches that of water, which supports that the system's spectral channels correspond to the specified wavelengths. The µ a and µ ′ s values are close to those reported in [15,19,43] for four wavelengths, and to those reported in [68] for 760 nm (µ a = 0.029 cm −1 and µ ′ s = 10 cm −1 ). The system and the used method can retrieve spectra of µ a and µ ′ s with high accuracy, as compared to the previously reported spectra. Figure 5 shows the obtained contrasts of five measurands and the sensitivity factors (SF) of three statistical moments, for the two inhomogeneous phantoms. The quantitative results are  Solid inhomogeneous phantom, data in Fig. 5 (a, b) ∆A  Table 3. The findings for the two phantoms are similar, except for the liquid phantom the contrasts reach deeper and have a higher magnitude. The contrasts for moments agree with SF and with their known behavior [45] since the position of the maximum contrast is deeper for higher order moments and the FWHM of the lateral scan increases following the broadening of 'banana' shapes. The moments m 1 and V increase when a perturbation is positioned near the surface and decrease when it is positioned deeper; this effect has been well described [45]. Oddly, the contrasts for ∆A can be negative at the smallest depths (the same was reported in [23,24,53]), which corresponds to an increase in the photon count rate after introducing a perturbation, although at these depths a perturbation is outside of the phantom (negative depth). The phantom inherently has a weak heterogeneity in the structure due to a movable part, and importantly, it was found that a homogeneous rod causes a noticeable perturbation [32]. Also, the measurement geometry is not semi-infinite because the rod extends beyond the phantom's surface. These reasons could have contributed to the difference between the obtained ∆m 1 (37.3 and 28.7 ps, shown in Table 3) during two scans when a perturbation was at the same depth and lateral position. These reasons are not present for the liquid phantom and the obtained ∆m 1 was similar during two scans (35.1 and 34.2 ps). The measure of ∆A is more prone to a slight drift, as seen in Fig. 5 (e) by comparing the values at -4 and 4 cm. Also, the sensitivity of ∆A is focused at a shallow depth, which is a known challenge for continuous-wave NIRS [69]. The results for ∆V are consistent with previous studies on similar phantoms [31,42].

nEUROPt protocol
The contrasts of the early ∆A Early and late ∆A Late time windows confirm that the later photons on average travel deeper and hence have a higher sensitivity to deeper perturbations. For an ideal system, the latest time windows could probe well beyond a depth of 3 cm [70]. However, an IRF with a substantial late afterpeak can have a devastating effect for time windows [31], although the changes in moments are independent of the IRF if the time range encompasses the whole DTOF [45]. The presented system's IRF allows to confidently use time windows and relevant methods, as confirmed by the IRF parameters (Section 4.1) and ∆A Late in Fig. 5. Figure 6 shows time traces of the obtained concentrations and StO 2 for two blood-lipid phantoms during dynamic StO 2 changes. The concentrations obtained with the presented system (TD-NIRS) at 3 and 4 cm are similar, and the values for ρ = 3 cm are shown because they contain lower noise. The signal level, i.e. photon count rate, greatly reduced after adding blood, which is a consequence of an increase in µ a .

Blood-lipid phantoms
The bubbling of oxygen can form bubbles inside the phantom, which then dissipate over time. In Fig. 6 (a, b), HbT always drops by a few µM during the bubbling of oxygen and afterward gradually returns to the initial value (prior to bubbling). The magnitude of the drop is higher for a higher blood concentration and reaches up to ∼5 µM for TD-NIRS at 3 cm and OxiplexTS, and up to ∼10 µM for TD-NIRS at 4 cm. The time to return to the initial value is about 12 min for TD-NIRS and always longer for OxiplexTS. The chosen zero-time is after oxygen was bubbled, and hence the first cycles show an increasing HbT as it returns to the value that it had. Apart from this dip, HbT remains constant during all desaturation cycles as intended.
The concentration of oxyhemoglobin HbO 2 obtained with TD-NIRS at 3 cm approaches 0 µM at the end of all desaturation events, except at the lowest scattering (25 ml Intralipid) where it is always above 5 µM. Consequently, the corresponding StO 2 is in the range from 100% to 0%, except at the lowest scattering where it is always above 10%. The amount of blood and Intralipid have no apparent effect on the estimation of StO 2 with TD-NIRS at 3 cm, except for the lowest amount of Intralipid. TD-NIRS at 3 and 4 cm obtain similar values of Hb, but the values of HbO 2 for 4 cm are slightly lower at the maximum and higher at the minimum compared with the values for 3 cm, which is amplified for the following cycles. Consequently, StO 2 for 4 cm is in the range from 95% to 10% during the first cycle, and the maximum decreases reducing the range with each addition of blood or Intralipid. Also, StO 2 for 4 cm is significantly noisier, which is visible in Fig. 6 despite the higher temporal averaging. Under most conditions, the system reliably retrieves concentrations and StO 2 with good precision.
For an assessment of the accuracy, we compare TD-NIRS with the results of OxiplexTS and co-oximetry. The closest match between the TD-NIRS and the OxiplexTS results is observed for 55 ml of blood ( Fig. 6 (a, c)) and 75, 100, and 125 ml of Intralipid ( Fig. 6 (b, d)). The largest difference is observed for the lowest concentrations of blood (20 and 35 ml) and Intralipid (25 and 50 ml), where OxiplexTS obtained values slightly below zero. When blood or Intralipid is increased, for the results of OxiplexTS, the maximum of HbO 2 remains the same while the minimum increases, such that HbO 2 and hence StO 2 become non-negative and match closer with TD-NIRS at 3 cm. The previously reported measurements with OxiplexTS on a similar phantom do not have negative values [50], which is probably because the authors used windows that Fig. 6. Measurements on two blood-lipid phantoms that contained a varying amount of blood (left panels) or Intralipid (right panels); the amount is expressed as the background color and written in the bottom panels. Concentrations of HbO 2 , Hb, and HbT obtained with TD-NIRS at 3 cm and OxiplexTS (a, b). StO 2 obtained with TD-NIRS at 3 and 4 cm and OxiplexTS (c, d). For displaying purposes, 10 data points were averaged for TD-NIRS at 3 cm and OxiplexTS, and 30 data points were averaged for TD-NIRS at 4 cm due to higher noise.   mimick a neonatal head (4 mm thickness, µ a = 0.10 cm −1 , and µ ′ s = 9.6 cm −1 ) and the windows in the present study mimick a layer of fat (1 mm thickness, µ a = 0.063 cm −1 , and µ ′ s = 6 cm −1 ). Table 4 reports the mean values of HbT during the plateau regions at the end of desaturation events, which lasted a few minutes, and the results of co-oximetry. The standard deviations of these values are negligible, except for TD-NIRS at 4 cm (∼1 µM). OxiplexTS values are close to the results of co-oximetry; the biggest discrepancies are for the highest concentrations of blood and Intralipid. The values of TD-NIRS at 4 cm match with OxiplexTS and co-oximetry within a  Fig. 6. TD-NIRS at 3 cm versus OxiplexTS (Left panels). TD-NIRS at 4 cm versus OxiplexTS (Middle panels). TD-NIRS at 4 cm versus at 3 cm (Right panels). The legends correspond to the amount of blood or Intralipid inside the phantom. Linear fits were obtained for data points 30% < StO 2 < 95% and reported in Table 5. Table 5. Coefficients for linear transformation of data in Fig. 7:

TD-NIRS (4 cm) OxiplexTS
TD-NIRS (4 cm) TH-NIRS (3 cm) The values of TD-NIRS at 3 cm are only slightly lower than the values at 4 cm. In summary, the system can retrieve HbO 2 , Hb, HbT, and StO 2 with an accuracy that is consistent with OxiplexTS and co-oximetry, except for the lowest concentration of Intralipid. The StO 2 data in Fig. 6 were resampled to 10 s on a common time base for comparing the values obtained by different systems during deoxygenation events. The StO 2 values obtained with TD-NIRS are almost always higher than obtained with OxiplexTS ( Fig. 7(a, b, d, e)), which could be observed also in Fig. 6. The values obtained with the two systems better match at higher StO 2 , which was also found for most systems that were similarly compared with OxiplexTS in [34,35]. Increasing the amount of blood brings the data points closer to the unity line, and varying the amount of Intralipid (excluding 25 ml) does not significantly affect the results. For the TD-NIRS measurements at 3 cm versus 4 cm (Fig. 7(c, f)), most data points are slightly above the unity line at low StO 2 and below the unity line at high StO 2 .
The differences between the traces of the obtained StO 2 can be quantified using the approach proposed by Kleiser et al. [34,35]: apply a first degree polynomial fit (lowest least square error) to the data in Fig. 7 for values of StO 2 from 30 to 90% (Kleiser et al. used 15 to 95%). The obtained parameters of the fit are presented in Table 5 for each deoxygenation cycle. There is no noticeable dependence on the concentration of blood and Intralipid. The slope is always lower than 1 and the y-intercept is always positive. The fits are better than for most systems that were used in [34,35].
The duration of deoxygenation cycles ranged from 23 to 38 min (mean 28 min) and it depended mostly on the added amounts of glucose and blood. A previous study showed a significant decrease in µ ′ s (around 2 cm −1 ) following deoxygenations inside a similar phantom [48]. The µ ′ s (obtained with TD-NIRS and OxiplexTS, not shown) was stable in both phantoms and decreased by an insignificant amount. The greatest reduction (about 0.2 cm −1 in 30 min) was during the last deoxygenation cycle that had the highest amount of blood. Similar results for Fig. 6 and 7 and Tables 4 and 5 were obtained using a subset of spectral channels of the TD-NIRS system, as long as the first and last channels were included, apart from an increased level of noise when using fewer channels. Figure 8 shows the retrieved dynamic changes in concentrations, which were calculated using the fitting method (time-domain information) and the MBLL method (intensity information), during arterial and venous cuff occlusions on the forearm of one subject. For the other two subjects, the following findings were also observed. During the venous occlusion ( Fig. 8(a, b)), both HbO 2 and Hb increase at two rates, first rapidly and then much slower (close to 0 for the shown subject, and even negative for HbO 2 for another subject). The observed changes are consistent with the physiology of a venous occlusion [52,55]. During the arterial occlusion ( Fig. 8(c, d)), the obtained increase in Hb is greater than the decrease in HbO 2 , which leads to an increase in HbT similar as found in [19,38,53,54,56]. After releasing the pressure, the concentrations 'overshoot' and then slowly return to baseline. The expected hemodynamic responses for both occlusions were retrieved using measurements at either ρ = 3 or 4 cm.

In-vivo measurements
The measurements at two ρ lead to similar changes in concentrations if using the fitting method, but not if using the MBLL method. The magnitudes are higher for ρ = 4 cm and closer to the results of the fitting method. A difference in the trend is noticeable in Fig. 8 (b), where HbO 2 for 3 cm decreases between 2.5 and 4 min but HbO 2 for 4 cm remains constant. An arm is a heterogeneous medium and greater hemodynamic changes are expected to occur in the deeper tissues, i.e. in muscles. The sensitivity is more focused inside the superficial layer for the MBLL method compared to the fitting method [68] and depth sensitivity can be improved for both methods by increasing ρ. The results were obtained with the fitting method (a, c) and the MBLL method (b, d). The red background represents the start and end of a 2-min occlusion period. The results of the fitting method were smoothed for displaying purposes: a third order median filter (13.8 s) for ρ = 4 cm, and a smoothing filter (1.5 s) for both ρ.

Discussion and future perspectives
This work demonstrated the performance of the system, which is intended for muscle and brain monitoring in a clinical environment. The IRF, the temporal fluctuations, and the warmup time confirm the good specifications of modern electronics (BIP protocol). The warmup time must be considered when planning clinical measurements, and a good IRF is necessary for using time windows and relevant methods [31]. The linearity for µ a and µ ′ s , the coupling between ∆µ a and ∆µ ′ s , and the accuracy of retrieving µ a and µ ′ s were quantitatively assessed (MEDPHOT protocol). The results are consistent with the best performances of previously reported multi-wavelength TD-NIRS systems for a big range of µ a and µ ′ s ; the nominal optical properties of the tested phantoms ranged from 0 to 0.04 cm −1 for µ a and from 5 to 20 cm −1 for µ ′ s . The contrasts measured on the solid inhomogeneous phantom, which has µ a and µ ′ s similar to tissue, confirm the system's sensitivity at great depths for measurements on tissues (if tissue has similar µ a and µ ′ s ) since the obtained ∆V is non-zero at up to ∼2.5 cm depth. The contrasts measured on the liquid inhomogeneous phantom, which had a much lower µ a but similar µ ′ s , showed a different behavior of measurands but which was consistent with the theory as illustrated with Monte Carlo simulations. Although using later time windows improves depth sensitivity, it also increases noise which hinders the overall performance of detecting deeper changes in µ a [46].
Blood-lipid phantom measurements demonstrated the system's ability to reliably and accurately retrieve absolute concentrations of hemoglobin, and hence StO 2 , for phantoms with a varying amount of blood and Intralipid. The results obtained using the reported system (TD-NIRS) are close to the reference measurements, which were obtained using OxiplexTS and co-oximetry. The biggest discrepancy is for the lowest amount of Intralipid, which corresponds to µ ′ s lower than a typical biological tissue. Improvements in the method of retrieving optical properties could improve the results. The used blood-lipid phantom is well suited for assessing the performance of oximeters: the ingredients can be precisely varied, the controlled deoxygenation cycles are highly repeatable, and up to four systems can measure in parallel without experiencing cross-talk. The used windows were mimicking a layer of fat, but they could be replaced to mimic other layers, e.g. neonatal head [48]. Previous studies have reported that occasionally the scattering properties of blood-lipid phantoms change over time, which is addressed in [48] and where possible reasons are discussed. In the present work, µ ′ s was obtained by both NIRS systems and it did not change significantly over time. The challenges of validating oximeters were described in detail in a review paper [71], including the pros and cons of in-vitro and in-vivo methods. Recently, Kovacsova et al. [72] presented a novel algorithm for measuring StO 2 and compared its performance to other published algorithms using data from a blood-lipid phantom that was developed in their group, as well as in-vivo data on a neonate, which demonstrates an example of another blood-lipid phantom and the latest improvements in the assessment of StO 2 .
The description of the system explained how 32 DTOFs (16 spectral channels and 2 detection modules) get acquired at up to 3 Hz. A drawback of measuring at 16 wavelengths with a single TCSPC card is the limited count rate, e.g. because of the pile-up effect, which inherently limits the count rate per wavelength compared to measuring at a single wavelength. The optical power is spread over a chosen wavelength range and the current version of the system cannot redistribute it to specific, desired wavelengths. For accurately estimating concentrations of chromophores, the extinction coefficients can be modified to account for the finite spectral bandwidth [73]. An advantage of the system is the excessive number of spectral channels and the fast sampling rate, which combined, may potentially allow monitoring the oxidation state of the CCO enzyme. The time-domain information allows applying depth-resolved methods for separating signals that come from different depths, e.g. brain and extracerebral tissue. Additionally, different spectral channels have slightly different depth profiles, which was exploited by Gerega et al. [37] to estimate changes in concentrations of ICG in many layers. The depth-resolved assessment will be investigated in the future and it has the potential to overcome the problems related to the extracerebral layer [69], which will make the measurements more reliable as well as more informative of the cerebral tissue.
The system is well suited for detecting fluorescence during an injection of Indocyanine Green [37], which allows assessing brain perfusion. The current study supports using the system also for other applications that involve monitoring constant and changing optical properties of tissue, e.g. tissue oximetry, tissue spectroscopy, and functional imaging.

Conclusions
We presented a description of the multi-wavelength time-domain NIRS system that obtains 32 DTOFs at up to 3 Hz. The quantitative results for the three established protocols are consistent with the best performances of the reported systems. The results for measurements on blood-lipid phantoms are close to the reference values obtained with OxiplexTS and co-oximetry, except a discrepancy was found for the lowest concentration of Intralipid. Venous and arterial cuff occlusions on the forearm resulted in the expected hemoglobin changes that are consistent with the physiology of such hemodynamic challenges. The presented quantitative parameters can be compared to other systems that measured similar tests.