Two-layered blood-lipid phantom and method to determine absorption and oxygenation employing changes in moments of DTOFs

Near-infrared spectroscopy (NIRS) is an established technique for measuring tissue oxygen saturation (StO2), which is of high clinical value. For tissues that have layered structures, it is challenging but clinically relevant to obtain StO2 of the different layers, e.g. brain and scalp. For this aim, we present a new method of data analysis for time-domain NIRS (TD-NIRS) and a new two-layered blood-lipid phantom. The new analysis method enables accurate determination of even large changes of the absorption coefficient (Δµa) in multiple layers. By adding Δµa to the baseline µa, this method provides absolute µa and hence StO2 in multiple layers. The method utilizes (i) changes in statistical moments of the distributions of times of flight of photons (DTOFs), (ii) an analytical solution of the diffusion equation for an N-layered medium, (iii) and the Levenberg–Marquardt algorithm (LMA) to determine Δµa in multiple layers from the changes in moments. The method is suitable for NIRS tissue oximetry (relying on µa) as well as functional NIRS (fNIRS) applications (relying on Δµa). Experiments were conducted on a new phantom, which enabled us to simulate dynamic StO2 changes in two layers for the first time. Two separate compartments, which mimic superficial and deep layers, hold blood-lipid mixtures that can be deoxygenated (using yeast) and oxygenated (by bubbling oxygen) independently. Simultaneous NIRS measurements can be performed on the two-layered medium (variable superficial layer thickness, L), the deep (homogeneous), and/or the superficial (homogeneous). In two experiments involving ink, we increased the nominal µa in one of two compartments from 0.05 to 0.25 cm−1, L set to 14.5 mm. In three experiments involving blood (L set to 12, 15, or 17 mm), we used a protocol consisting of six deoxygenation cycles. A state-of-the-art multi-wavelength TD-NIRS system measured simultaneously on the two-layered medium, as well as on the deep compartment for a reference. The new method accurately determined µa (and hence StO2) in both compartments. The method is a significant progress in overcoming the contamination from the superficial layer, which is beneficial for NIRS and fNIRS applications, and may improve the determination of StO2 in the brain from measurements on the head. The advanced phantom may assist in the ongoing effort towards more realistic standardized performance tests in NIRS tissue oximetry. Data and MATLAB codes used in this study were made publicly available.


Introduction
Near-infrared spectroscopy (NIRS) is a non-invasive and safe (even for long continuous monitoring) technique in biomedical optics for studying brain functions and its disorders. NIRS tissue oximetry enables monitoring vital parameters in real time, reporting on the balance between oxygen supply and consumption. The clinical usefulness of NIRS has been demonstrated by a wide range of clinical applications and trials (reviewed in [1][2][3][4][5][6][7][8][9]) and the popularity of NIRS continues to grow in part due to the continuous advancements in instrumentation and the resultant increasing number of NIRS systems (reviewed in [10][11][12][13][14][15]). A significant improvement to NIRS tissue oximetry would be the determination of StO 2 in multiple layers and a means of testing such method. Another field of NIRS is functional NIRS (fNIRS), which typically involves small hemodynamic changes (i.e. small ∆µ a ) related to brain activity [16,17]. The new data analysis method and the phantom introduced in this study may be useful for both fNIRS and NIRS tissue oximetry applications.
The analysis of NIRS measurements carried out on a layered tissue structure, such as the head, require a multi-layered model as the homogeneous model is inadequate [18][19][20][21][22]. Recent advances in modeling the light propagation in tissue (theories, solutions, and validations) are explained in the most recent book by Martelli et al. [23]. In this study, we used the time-domain NIRS (TD-NIRS) modality, which records the distribution of times of flight of photons (DTOFs), providing the most information about light propagation compared to continuous-wave (steady-state) and frequency-domain modalities [16]. Light in tissue can be modelled using the radiative transfer equation (RTE) and a commonly employed solution is the diffusion equation, which relies on the diffusion approximation [24][25][26]. Helton et al. [27] presented an open-access numerical routine for solving the diffusion equation for an arbitrary number of layers, for the steady-state case and in the time domain. Geiger et al. [28] modeled light propagation using spherical harmonics instead of the diffusion equation and presented the analytical solutions for the RTE.
To solve for complex geometries, the finite element method can be used, e.g. NIRFAST [29]. Another approach for modelling light propagation is based on the Monte Carlo method [30,31], which allows implementing any geometry but the main disadvantage is the long computational time. We will use the analytical solution for an N-layered medium presented by Liemert et al. [32], which applies the Laplace transform (LT) to a solution in the frequency domain (could be either numerical or analytical) to obtain the solution in the time domain. LT has important advantages over the commonly used Fourier transform, such as accuracy, efficiency (speed-up of up to several orders of magnitude), and numerical stability. Having a model for light propagation and its solution, the optical properties can be obtained from measured NIRS data by solving the inverse problem.
Various algorithms have been developed for determining the absorption (µ a ) or its changes (∆µ a ) at multiple depths of a multi-layered turbid medium using TD-NIRS measurements, e.g. more recently in Refs. [28,[33][34][35][36]. We will use the Levenberg-Marquardt algorithm (LMA), which allows solving nonlinear iterative least squares minimization problems [37,38]. The LMA is a fitting routine that searches for a solution that minimizes an error, i.e. the difference between the measured data and the theoretically-calculated data. Garcia et al. [35] proposed implementing the Kalman filter algorithm, which is another iterative algorithm that has some advantages over the LMA. Instead of using the DTOF as a measurand, Hallacoglu et al. [39] used the amplitude change and the phase shift at multiple source-detector distances as the measurands for the LMA. Steinbrink et al. [40] used changes in the number of photons in different time-channels of a DTOF for determining ∆µ a at different depths, since different time-channels have different depth sensitivity profiles. Similarly, different-order statistical moments of the DTOF have different depth sensitivity profiles and Liebert et al. [41,42] outlined the use of moments for the determination of ∆µ a at multiple depths. This has become a well-established method, e.g. [43][44][45][46]. Recent performance-assessment studies, including quantitative tests, found that moments outperform other measurands in terms of assessing ∆µ a in deep layers [47][48][49]. Using changes in moments at multiple wavelengths, which have slightly different depth sensitivity profiles, it was possible to estimate changes in the concentration of an injected indocyanine green (ICG) bolus in up to 24 layers [50]. However, up to now, the method relied on the assumption that the changes in moments are linearly proportional to ∆µ a , where the proportionality constants are the sensitivity factors [41]. This inverse problem of determining ∆µ a from changes in moments was solved using the linear least squares method or the singular value decomposition approach, i.e. non-iterative approaches.
One goal of this work is to improve the method based on the analysis of moments of DTOFs by using the analytical solution of the diffusion equation for a multi-layered medium and implementing the LMA. With this approach, we solve a nonlinear iterative least squares minimization problem, hence we remove the limiting assumption that changes in moments are linearly proportional to ∆µ a . In addition to testing the method on phantom measurements, we will investigate how errors in the values of the baseline parameters affect the determined ∆µ a values. We aim to provide a more accurate and robust method for determining ∆µ a at multiple depths in a multi-layered medium.
When introducing a new method of data analysis or a new system, phantom measurements are a necessary test before conducting in-vivo studies as phantoms allow to study the influence of the instrumental factors and uncertainties for well-defined geometries and optical properties [51]. There is an ongoing effort towards standardization of system performance tests in the field of NIRS [15,51], which has become increasingly important in part due to the increasing number of NIRS systems [10][11][12][13][14][15]. Recently, Lanka et al. [52] presented the results of the BitMap campaign, which was a multi-laboratory effort to perform a standardized performance assessment of instruments for diffuse optics, which used the BIP [53], nEUROPt [54], and MEDPHOT [55] protocols. Other standardization efforts resulted in a consensus by many groups on the guidelines for the best practices for fNIRS publications [56]. Over the past decade, many tissue-mimicking phantoms for optical spectroscopy have been developed, and they are reviewed in [51,[57][58][59], including challenges related to their construction and utilization. Konugolu Venkata Sekar et al. [60] developed solid tissue-like phantoms that mimic the anatomically correct 3D shapes of human organs for biomedical applications, e.g. [61]. For testing tissue oximeters and for other biophotonics applications, there is a need for advanced phantoms that can provide realistic simulations of tissues' optical properties. Liu et al. [62] developed a phantom that achieves simulating vascular oxygenation and perfusion, by filling with blood the tubes that mimic vessels. Afshari et al. [59] developed a phantom with an array of 148 cylindrical channels, filled with whole bovine blood, as a new approach for standardized testing of clinical oximeters. An established concept for a phantom is using whole human blood (mixed in a turbid suspension) for mimicking repeatable controlled dynamic oxygen saturation (StO 2 ) levels from 0 to 100%. To date, all presented phantoms that used this concept had a single homogeneous liquid compartment, e.g. [63][64][65][66][67]. However, when NIRS is intended to study deeper structures, e.g. the brain, it is crucial to address the influence of the superficial layer on the measured StO 2 [18][19][20][21][22].
The second goal of this work is to develop a phantom with two compartments (superficial and deep), each of which can hold a mixture containing blood that can be oxygenated or deoxygenated. The phantom is a tool that will allow quantitatively studying the contamination from the superficial layer when determining StO 2 in a two-layered medium. We will utilize the phantom to demonstrate the proposed method's capabilities of determining StO 2 in both layers.
Data and MATLAB codes used in this study were made publicly available [68].

NIRS system
The multi-wavelength TD-NIRS system was developed by Gerega et al. [69] at IBIB PAN and it was described and experimentally tested in our recent study [70], revealing its advantages compared to the current state-of-the-art systems. The system can simultaneously acquire DTOFs at 16 spectral channels (from 680 to 875 nm, with a 12.5 nm/channel width) for each of two detection modules. Each module consists of a polychromator, a multi-anode 16-channel photomultiplier tube, and a 16-channel TCSPC card, as described in [70]. Two emitting fibers and two detection fiber bundles allow simultaneous measurements on two windows of the phantom in reflectance geometry. The previous supercontinuum laser (SC480-8 WhiteLase, Fianium, UK) was replaced with a newer laser (NKT SuperK EXTREME (EXR-15), Denmark), which has similar characteristics. We used the repetition rate of 78 MHz and set the constant current power to 50% in the laser driver. The parameters of a TCSPC card were configured so that each DTOF contained 1024 time bins with a width of 12.22 ps. This width is sufficiently small as a study on accuracy suggested that a bin width of about 10 ps enables a better than 5% accuracy in quantifying optical properties using the curve fitting method [71]. In this study, the distance between the source and detector fibers was set to 30 mm. The acquisition time was 0.3 s and the sampling rate was 3 Hz. To increase the signal-to-noise ratio and reduce the computational time, 9 consecutive DTOFs were summed for measurements on phantoms with blood, resulting in a sampling time of 3 s. For measurements with ink, 128 DTOFs were summed for each ∆µ a step, resulting in a sampling time of 42.7 s. We used a shorter sampling time for analyzing measurements with blood to resolve the dynamic changes of StO 2 . We used a longer sampling time for experiments involving ink to reduce noise since the optical properties were stable.
To measure the instrumental response function (IRF), an emission fiber was aligned with a detection fiber bundle in transmission geometry at a distance of 6 cm with a neutral density filter between them and two pieces of white paper covering the source fiber and the detection fiber bundle [72]. The maximum laser power was measured as 15 mW for each emission fiber, which would equate to a power density of less than 1.5 mW/mm 2 on a skin's surface [70]. For a recent comprehensive study on IRF, see the study conducted by Pirovano et al. [73].

Data analysis
The measured NIRS data were stored and later analyzed using custom codes written in MATLAB R2020b. The LMA was implement using the function lsqcurvefit in MATLAB. To determine the values of µ a and the reduced scattering coefficient (µ ′ s ) from a measured DTOF, we used a solution of the diffusion equation for a homogeneous medium and the LMA. This is known as the curve-fitting method and it is regarded as the gold-standard approach for analyzing time-domain diffuse reflectance measurements [37,38,71], and an example of the distribution of error norm was presented in [74].
We used a homogeneous semi-infinite medium model under the extrapolated boundary conditions, which was done by using the solution of the diffusion equation for an N-layered medium [32] with N set to 1. We used equations presented in [25], which allow to account for the boundary conditions. The theoretically-obtained DTOFs are convolved with the IRF that was measured at the start of an experiment. Analyzing each measured DTOF separately, we first subtracted the background signal and then used the time-channels from 85% on the rising edge to 1% on the tail for determining µ a and µ ′ s .

Determination of ∆µ a in two layers
In this section we outline an improved method of determining ∆µ a in multiple layers. To obtain the absolute µ a with this method, we added the determined ∆µ a values to the baseline µ a values, which were obtained using the method explained in the previous section.
DTOF is represented by a histogram of photon counts (N i ) at time channels (t i ) indexed by i [41]. The first three statistical moments of DTOFs are: m 0 = N tot = ∑︁ N i (total number of photons), m 1 = ∑︁ t i N i /N tot (mean time of flight), and m C 2 = V = ∑︁ (t i − m 1 ) 2 N i /N tot (variance). The change in attenuation is: ∆A = − log(N * tot /N tot ), where N * tot is after ∆µ a . Moments have several important advantageous properties. To deconvolve the influence of the IRF, the moments of the IRF are subtracted from the moments of the DTOF; valid for m 1 , m C 2 , and m C 3 (third central moment) [46]. The resulting moments can be used to calculate the optical properties of tissue [75]. Changes in moments are independent of the IRF under the assumption that the moments of the IRF are constant. The uncertainty of moments due to photon noise, which behaves like Poisson noise (Poisson statistics), can be calculated using higher-order moments [42,46].
The method that uses the changes in moments (∆A, ∆m 1 , and/or ∆V) requires two DTOFs: during baseline and after ∆µ a . For small ∆µ a in multiple layers, the corresponding changes in moments are linearly proportional to ∆µ a in each layer, where the proportionality constants are the sensitivity factors (SF) for each layer, as explained by Steinbrink et al. [40] and Liebert et al. [41,42]. The SF for baseline optical properties can be obtained using Monte Carlo simulations [30] or an analytical solution of the diffusion equation [43], and then it is a linear inverse problem to determine ∆µ a from the measured changes in moments.
To improve the method, we overcome the assumption that the changes in moments are linearly proportional to ∆µ a by using the solution of the diffusion equation for an N-layered medium [32] and the LMA. In this study, we set N to 2, assumed the deeper layer is semi-infinite, and assumed the extrapolated boundary conditions. The use of LMA avoids the need for sensitivity factors. For a guess of ∆µ a in each layer, the analytical solution is used to generate the corresponding two DTOFs for which changes in moments are calculated. The LMA searches for ∆µ a in each layer that result in the smallest error norm (χ 2 ), which is the difference between the measured and the theoretically calculated changes in moments. The different-order moments have different units and levels of uncertainty. To account for these, we weighted each moment by its uncertainty due to photon noise [42]. When all three moments are considered, χ 2 can be calculated using the equation: where changes in moments are calculated from two DTOFs: ∆A = A − A * , ∆m 1 = m 1 − m * 1 , and ∆V = V − V * , where the star denotes the value of a moment after ∆µ a . For a guess of ∆µ a , the corresponding changes in moments (∆A sim , ∆m 1,sim , and ∆V sim ) are calculated from the theoretically obtained DTOFs. The uncertainties due to photon noise (δ∆A, δ∆m 1 , and δ∆V) can be calculated using higher-order moments [42]. We make a simplifying assumption that the noise of moments of two DTOFs (before and after ∆µ a ) is the same (true if ∆µ a = 0).
To assess the method, we will calculate and plot the distribution of (2/N tot ) χ 2 (Eq. (3)) for all combinations of ∆µ a in two layers [34,74,76,77]. N tot is a scaling factor, and it is a function of the selected laser intensity and acquisition time. Therefore, the calculated χ 2 is unitless, except for the linear dependence on N tot .
In this study, we did not use the intensity information (∆A) for determining ∆µ a because it requires careful data processing to account for possible laser drifts and changes in the optical filters. During measurements, the optical filters were adjusted (by rotating the Neutral Density filter) in order to regulate the count rate. Additionally, ∆A is more susceptible to laser drifts than other moments.
We used the developed phantom, which mimics two layers with changing µ a , for experimentally testing the new method. The method requires knowing the values of all parameters at baseline, which can be obtained for phantom measurements. The refractive index (n) was assumed to be 1.33 for both layers because the main component of the phantom mixture was distilled water. The thickness of the superficial layer (L) was measured with rulers glued to the phantom. To obtain the values of µ a and µ ′ s of each layer, we used measurements on the individual compartments (superficial and deep) and applied the curve-fitting method with the homogeneous model. For experiments involving blood, the targeted optical properties at baseline were the same for both layers, so measurements on a single layer were sufficient to determine the optical properties of both layers.
Although the changes in the first three statistical moments are independent of the IRF, their values are influenced by the chosen time-channels that are used for their computation [43]. We computed moments using the same time-channels for all theoretically obtained and measured DTOFs: from 25% on the rising edge to 3% on the tail of the DTOF that was measured at baseline. Prior to calculations of moments, the theoretically obtained DTOFs were convolved with the IRF that was measured at the start of an experiment.
The new method's computation time for each DTOF is typically a few seconds, which can be reduced with a good initial guess. Therefore, we employed the established method based on moments, which uses sensitivity factors, for obtaining the initial guess of ∆µ a in two layers. We obtained the sensitivity factors by computing the changes in moments resulting from ±2% ∆µ a in each layer, as done in [43].

Calculation of concentrations of oxy-and deoxyhemoglobin
To determine the concentrations of oxyhemoglobin (HbO 2 ) and deoxyhemoglobin (Hb) in experiments involving blood, we used the determined µ a values at multiple wavelengths and the known molar absorption coefficients (ε) of HbO 2 and Hb. We used 9 spectral channels (3 rd to 11 th out of 16) covering the spectral range from 705 to 805 nm, which is consistent with our previous study with blood-lipid phantoms [70]; the dependence on the choice of wavelengths should be further investigated.
We used the ε spectra provided by Gratzer et al. [78] and the absorption spectra of water provided by Matcher et al. [79]. We calculated the concentrations of only HbO 2 and Hb using the µ a values, which were obtained by subtracting the µ a of water (multiplied by 0.985, which was the volume fraction of water and SMOFlipid) from the determined µ a values.

Construction of phantom
In recent years, various tissue-mimicking phantoms containing blood have been developed, e.g. [59,[62][63][64][65][66][67]. However, to our knowledge, no presented phantom has two blood-lipid compartments with variable StO 2 , which could be used to mimic the extra-and intracerebral compartments of the head.
We present a new tissue-mimicking phantom that was built in-house at IBIB PAN. The phantom consists of three containers, as shown in Fig. 1. The deep container is similar to the one presented by Kleiser et al. [63] and holds the liquid for the deep compartment. The container has an irregular octagonal shape. Three walls have windows for optodes, which allow three NIRS systems to measure simultaneously, and the fourth wall has a separating window (100 mm × 100 mm). The superficial container holds the liquid for the superficial compartment and the third container separates this liquid into parts A and B. The superficial container has three walls and the side without a wall attaches to the deep container. The third container is shorter and suspended, which leaves a 25 mm gap underneath for the liquid to flow between parts A and B, as shown in Fig. 1 (b). The thickness of part A can be changed from approximately 0 to 65 mm by moving the third container. Part B of the superficial compartment allows adding and efficiently mixing ingredients. Measurements can be performed simultaneously on any of the five windows (60 mm × 40 mm) for measuring on different compartments: the deep (Pos. #1), the two-layered medium (Pos. #2), and the superficial (Pos. #3). Below we list the main features of the phantom. Materials. Five walls have windows for measurements (next to red, blue, and green arrows in Fig. 1) and these walls are made from black PVC (5 mm thick). All other walls and the covering lids are made from transparent PVC (5 mm thick), which allows seeing inside the phantom. To see and measure the height of the liquid in the deep container, one corner was cut out and replaced with transparent PVC, shown in Fig. 1 (a). To prevent ambient light interference during measurements, the entire phantom was covered with a black blanket. The deep and the superficial containers are held together by four screws, which have washers (2 mm thick solid plates cut from PVC), and in-between is a glued water-resistant rubber cord, which gets squeezed to prevent leakage. The bottom parts are made from copper (3 mm thick), which allows using hot plates for controlling the temperature. Each of the five windows (60 × 40 mm) has two latches for holding the optode holders (58 × 38 mm), which were 3D printed (pure black ABS material, Zortrax Z-ABS v2). All windows were covered with translucent Mylar film (0.050 mm thick, 785-0786 RS Pro, UK), which was permanently fixed to the walls using silicone. The effect of a thin, slightly scattering Mylar film on diffuse reflectance measurements can be considered negligible [80]. To limit ambient oxygen diffusion, the two containers were covered with lids, which have openings for fixing a thermometer, an oxygen supply, a peristaltic pump, and for adding ingredients.
Thickness of the superficial layer (L). The third container is designed to slide along one direction, which is achieved by four linear bearing guides attached to its top and two linear shafts attached to the top of the other two containers. A threaded rod is attached between the two linear shafts to enable precise movement of the container. Rulers were glued on top of the third and the superficial containers to facilitate measurement of the third container's position and hence L. To counteract buoyant forces, heavy objects can be placed inside the third container.
Separating window and bending of Mylar film. To separate mixtures contained in the deep and superficial containers, we covered the separating window (100 mm x 100 mm, as shown in Fig. 1) with Mylar film. When the Mylar film gets wet, it is able to bend easily. To maintain a constant position of the Mylar film, we applied pressure on it by pouring more liquid into the deep container, so that the liquid level was approximately 2 cm higher than in the superficial container. This created a uniform high pressure on the Mylar film, resulting in a constant and repeatable bending. During data analysis, we estimated that the Mylar film bent by approximately 0.5 mm, which we approximated as a flat boundary, given the large area of the separating window.

Stirring and heating.
Stirring and heating were accomplished by placing the deep and the superficial containers on top of two hot plates (Stuart US152D, Bibby Sterlin Ltd. and SMP-TC, SET GSM, Poland). In all experiments, the stirring was set at 500 RPM. Additionally, in the superficial container, a peristaltic pump (372.C Elpan, Poland) was used to move the liquid from the top of part B to the top of part A (shown in Fig. 1 (b)) at a flow rate of 500 ml per min through an 8 mm diameter silicone tube. In experiments involving blood, the temperature was maintained at 37 ± 1°C. To verify the efficacy of stirring, we added ink to either the deep or the superficial container and checked how quickly the measured signals reached new stable values, which typically occurred within about 10 seconds.
Controlled and monitored parameters. In experiments involving blood, the following parameters were controlled and monitored. The temperature in each compartment was recorded every 10 s using two submerged digital thermometers (RC-5+, Elitech). The pH level was kept constant at 7.4 by adding the phosphate-buffered saline powder (PBS powder, sc-24947, ChemCruz). The blood was deoxygenated by adding dry bakery yeast (Drozdze Suszone Instant, Dr. Oetker, Poland). To oxygenate the blood, medical-grade oxygen from a tank was bubbled (3 L/min flow setting on the pressure regulator) through an air stone (Mist-Air, fine bubbles, Kordon) submerged in the deep container. The air stone was held in place by a long rigid plastic tube and positioned in the middle of the deep container, approximately 3 cm above the bottom. Experiments involving ink required controlling only the stirring rate.
Simultaneous NIRS measurements can be performed on up to five windows, i.e. with up to five NIRS systems. If the source-detector distance (ρ) of the optodes fixed on the phantom is 40 mm, then the minimum distance between fibers on different windows is at least 75 mm, as in [63]. Light has a negligible chance of traversing 75 mm if the optical properties of the phantom are similar to those of biological tissue [76].
Boundary conditions. We assumed semi-infinite model with extrapolated boundary conditions. We approximated the deep container as an infinite medium (semi-infinite for measurements on the surface). For ρ = 40 mm, the distance from a fiber to the closest neighbouring wall is 60 mm, so the laser light needs to traverse at least 160 mm to reach a neighbouring wall and then a detection fiber bundle. Optodes were attached 75 mm above the ground level and the liquid was filled up to about 150 mm (or 170 mm) above the ground level in the superficial (or deep) container. We assumed that black PVC and black optode holders absorb all light that penetrates the surface. We neglected the effect of the Mylar film [80].

Characterizing ink and SMOFlipid
The use of ink and Intralipid in tissue-simulating phantoms has been studied and described in the literature, e.g. Di Ninni et al. [81,82] and Spinelli et al. [83,84]. In this section, we summarize how we determined the nominal values of µ a and µ ′ s , which we used as a guide for the experimental protocol and later for a comparison with our TD-NIRS results.
We used SMOFlipid 20% (Fresenius Kabi, Poland) as a substitute for Intralipid, which we confirmed has similar µ ′ s and µ a . We calculated the nominal µ ′ s of a phantom based on the volume fraction of SMOFlipid [81,84]: where ε ′ s,SL is the specific reduced scattering coefficient of SMOFlipid. We used the value 21.5 mm −1 for 750 nm, which was reported for Intralipid [81,83,84]. V Dilted ink , V SL , and V water are the volumes of diluted ink, SMOFlipid, and water, respectively.
We characterized the µ a of ink (Rotring, black, 23 ml, Germany) because it can have large brand-to-brand and batch-to-batch variations [82,84]. We diluted pure ink in distilled water (d ink = weight of ink / total weight), measured its absorbance (A p ) in a 1 cm thick cuvette with a spectrophotometer, and calculated its effective absorption coefficient (ε eff,ink = A p / d ink ). At 750 nm, we obtained ε eff,ink = 626 mm −1 , which is consistent with a previously reported value for this brand of ink [82]. We combined the value of ε eff,ink obtained from spectrophotometric measurements and the value of the single-scattering albedo (Λ) to predict the value of the specific absorption coefficient of ink (ε a,ink ), as proposed in [84]. By definition, Λ = (ε eff,ink -ε a,ink ) / ε eff,ink and we used the value of 0.15 for Λ. We calculated the nominal µ a of a phantom based on the volume fraction of the diluted ink [82,84]: We made a simplifying assumption that SMOFlipid has the same µ a as water, since their µ a are similar and the main constituent of the phantoms in all experiments was water.

Protocol for experiments involving ink
We conducted two experiments involving ink, in which we increased µ a in either the deep or the superficial compartment. In one compartment, the targeted µ a values ranged from 0.05 to 0.25 cm −1 (at 750 nm) in 20 steps, while in the other compartment, the targeted µ a was fixed at 0.10 cm −1 at 750 nm. The targeted µ ′ s was 10.7 cm −1 at 750 nm in both compartments. These values of optical properties are within the typical range for tissues [85,86], specifically of the human head [87]. Table 1 lists the targeted amounts of ingredients. To increase µ a without changing µ ′ s , we prepared Mixture #2 (water, diluted ink, and SMOFlipid), which has a much higher µ a (∼ 4.7 cm −1 ) but the same µ ′ s as the mixture in the phantom. We set L to 14.5 mm in both experiments, which were conducted on two consecutive days. We used the TD-NIRS system to measure on the deep compartment (Pos. #1), the two-layered medium (Pos. #2), and/or the superficial compartment (Pos.#3), according to Fig. 1. The nominal µ a values were calculated by using the weights of the added ingredients and Eq. (5). Distilled water and SMOFlipid were weighted with a precision of 0.2 g. For adding Mixture #2, we used a syringe and measured its weight before and after injecting the mixture into the phantom, with a precision of 0.1 µg. We waited at least 40 s for the ingredients to stir before taking NIRS measurements. After adding an ingredient, its measured weight was recorded and used to update the targeted weight of the next ingredient to be added, which helped achieve the targeted µ a values across the twenty ∆µ a steps.

Protocol for experiments involving blood
We conducted three experiments involving blood using the same protocol, but setting different thickness of the superficial layer L (12, 15, and 17 mm). We used TD-NIRS system to measure on the deep compartment (Pos. #1) and the two-layered medium (Pos. #2). Measurements on the deep compartment allow assessing the repeatability of the experiment, as they are independent of L.
We changed the oxygen saturation (StO 2 ) of blood separately in the superficial and the deep compartments, in one compartment at a time. The targeted µ ′ s of both compartments was 10.7 cm −1 at 750 nm. Table 2 lists the protocol steps along with the targeted StO 2 . Table 3 lists the targeted amounts of ingredients. The preparation of the phantom before starting TD-NIRS measurements comprised of six steps (P1 to P6). After step P6, the deep compartment contained 3366 g of distilled water with PBS and 180 g of SMOFlipid, as in Table 3. In step S1, we added human erythrocyte concentrate to each compartment. We used a syringe with 1.6 mm inner diameter for handling blood-containing mixtures. In step S2, we took some of the mixture (approximately 40 ml) from the deep compartment, mixed yeast in it, and poured it back into the deep compartment. We denoted this as time zero (Time = 0 min). In step S3, we bubbled oxygen in the deep compartment for 2 min. In step S4, we added yeast to the superficial compartment using a similar procedure as was used for adding yeast to the deep compartment. Then, we repeated step S3 three times (steps S5, S6, and S7). The protocol consisted of six deoxygenation cycles: one in the superficial compartment and five in the deep compartment. These cycles are indicated by six arrows at the top of figures that show time traces of the TD-NIRS results. The three experiments were conducted on different days. The first experiment (Exp. #1, L = 12 mm) and third experiment (Exp. #3, L = 17 mm) used blood from the same blood transfusion bag. The second experiment (Exp. #2, L = 15 mm) used another blood transfusion bag.

Method for determining ∆µ a in two layers
Here we examine the proposed method for determining ∆µ a in two layers. We calculated the values of ∆A, ∆m 1 , and ∆V for all combinations of ∆µ a,Sup and ∆µ a,Deep , and the results are presented in Fig. 2 (a), where the changes in each statistical moment are shown using contour lines. The lines of all pairs of two moments intersect only once, which confirms that there is a unique solution of ∆µ a,Sup and ∆µ a,Deep when changes of any two moments are analyzed. The slopes of such lines (changes in moments versus ∆µ a,Sup and ∆µ a,Deep ) are related to depth selectivity [48,88], which is a metric that reflects the sensitivity of a measurand to ∆µ a,Deep in relation to ∆µ a,Sup . The lines for ∆A are more vertical, indicating that ∆A is more sensitive to ∆µ a,Sup compared to ∆µ a,Deep , i.e. ∆A has low depth selectivity. The lines become more horizontal for higher-order moments, corresponding to increased depth selectivity, which is a well-known property of moments [48,88]. For small values of ∆µ a , the lines for all moments be approximated as linear functions. However, for larger values of ∆µ a , the lines become increasingly nonlinear, confirming the need for using the LMA, i.e. solving a nonlinear problem of determining ∆µ a in multiple layers using changes in moments. We calculated the error norm (χ 2 ) using Eq. (3) for all combinations of ∆µ a,Sup and ∆µ a,Deep . The distribution of χ 2 is shown in Fig. 2 (b) when all three moments are considered and in Fig. 2 (c) when excluding the term with ∆A, i.e. using only the terms containing ∆m 1 and ∆V in Eq. (3). The presented χ 2 distributions have a single minimum, i.e. no local minima, confirming the suitability of the LMA approach for finding this minimum. We verified that the LMA reaches the same solution (determined ∆µ a in two layers) regardless of the initial guess.
The shape of the χ 2 distribution depends on which moments are used and the weights assigned to them. When three moments are used (Fig. 2 (b)), the contribution from ∆A dominates, resulting in a χ 2 distribution that is more dependent on the superficial layer compared to the deep layer. This distribution resembles the line where ∆A = 0 in Fig. 2 (a). When the term with ∆A is excluded (Fig. 2 (c)), the contribution from ∆m 1 becomes dominant, and the shape of the resulting χ 2 distribution resembles the line where ∆m 1 = 0 in Fig. 2 (a). These results are consistent with the differences in magnitudes of the contrast-to-noise ratio (CNR), e.g. reported in [88] where a 15% increase in µ a,Sup resulted in a CNR of about 100 for ∆N tot , 30 for ∆m 1 , and 10 for ∆V. Therefore, for a typical ∆µ a , the term with ∆A in Eq. (3) contributes much more to the value of χ 2 than the term with ∆m 1 , which is consistent with the one-order magnitude difference in the colorbars for χ 2 in Fig. 2 (b) and (c). If only a single moment is considered in Eq. (3), the shape of the χ 2 distribution would follow the moment's zero contour line shown in Fig. 2 (a) and the value of χ 2 would be zero along this line, resulting in a lack of a single minimum.
We examined the impact of an error in the value of a baseline parameter on the determined ∆µ a and the results are presented in Fig. 3. We simulated two DTOFs (at baseline and after ∆µ a,Deep = 0.05 cm −1 ) and retrieved ∆µ a,Sup and ∆µ a,Deep but using different assumed values for one of the baseline parameters: µ a , µ ′ s , and n in the superficial or the deep compartment, as well as L. The findings depend on the chosen ground truth values of all parameters (for both simulated DTOFs), but the results in Fig. 3 are a representative example. Figure 3 (b) shows that µ ′ s,Deep has little effect on TD-NIRS signals, which was also found in previous studies [89]. Figure 3 (a, b) demonstrates that an error in the assumed value of µ a,Sup has an opposite effect on the determined values of ∆µ a , compared to an error in the assumed value of µ ′ s,Sup . This agrees with the general understanding that these parameters (µ a and µ ′ s ) have contrasting impacts on the shape of the DTOF [74]. Errors in µ ′ s,Sup and n Sup have similar effects (Fig. 3 (b, c)), which is expected as both impact how quickly light passes through a medium. The determined ∆µ a,Sup is minimally affected by L (Fig. 3 (d)). However, when L is underestimated (or overestimated), then the determined ∆µ a,Deep approaches (or moves away from) the true value of ∆µ a,Sup . This dependence on L remained consistent when using ground truth values of ∆µ a,Sup = 0.05 cm −1 and ∆µ a,Deep = 0.   Table 1, which demonstrates good control and repeatability of experiments. The spectral shape of µ a for the smallest nominal µ a resembles the shape of the µ a of water. For the spectral shape of µ ′ s , please refer to Fig. 8 (b). Figure 4 (b) shows the determined µ ′ s values for 20 ∆µ a steps in the deep compartment in Exp. #1, at different wavelengths. For both experiments, the µ ′ s sometimes slightly increase with µ a , especially for the first few (about four) ∆µ a steps. However, for most wavelengths, the values of µ ′ s changed by less than approximately 0.3 cm −1 (less than 3%), following a 0.20 cm −1 (500%) increase in µ a . One cause may be the coupling of ∆µ a to ∆µ ′ s , which we previously reported for our system [70] by measuring on a well-characterized array of solid homogeneous phantoms that have the same µ ′ s (either 5, 10, 15, or 20 cm −1 ) and different µ a values (ranging from 0 to 0.35 cm −1 ) [55], as part of the BitMap campaign [52]. For a comparison of the coupling with other systems, see [52]. Other causes for changes in µ ′ s could be, e.g. possible timing drifts in the laser output. Increasing the concentration of ink should not have influenced the µ ′ s because we used Mixture #2 (containing water, SMOFlipid, and ink) that had the same µ ′ s as the phantom but much higher µ a (refer to section 2.3.3).   Figure 5 (b) shows the changes in moments (∆A, ∆m 1 , and ∆V) resulting from 20 ∆µ a steps. The baseline was selected as the middle measurement (10 th ∆µ a step, when the targeted µ a = 0.15 cm −1 at 750 nm). All moments show high contrast to ∆µ a,Sup and much lower contrast to ∆µ a,Deep , which is related to the large L (14.5 mm). The contrast to ∆µ a,Deep (and hence depth selectivity) increases for higher-order moments, consistent with the results in Fig. 2 (a). Notably, ∆A appears to be close to zero even for the largest ∆µ a,Deep . However, to assess a measurand's capabilities of detecting ∆µ a , it is essential to consider also the contrast-to-noise ratio (CNR) [48,88]. Figure 6 presents the determined µ a for 20 ∆µ a steps for the two experiments involving ink. One source-detector pair was used to measure on the compartment in which µ a was changed: the deep compartment for Exp. #1 and the superficial compartment for Exp. #2. These measurements were analyzed with the curve-fitting method (homogeneous model) to determine µ a and µ ′ s . The determined µ a values (black crosses in Fig. 6) are close to the nominal µ a values (the identity line) for all 20 ∆µ a steps and for both experiments. The signal level (i.e. N tot ) decreased for higher µ a , leading to increased noise, which could be the reason for the slight discrepancies for higher µ a . For measurements on a two-layered medium (Pos. #2), the curve-fitting method Measurements on two-layered medium (Pos. #2) were analyzed using the same method (green circles) and the method based on changes in moments with a two-layered model (red and blue crosses).

Experiments involving ink
(homogeneous model) provided µ a values (green circles in Fig. 6) that were between the µ a values of the two compartments and closer to the superficial compartment's µ a . This highlights the need for a two-layered model when measuring on a two-layered medium. We applied the proposed method based on the changes in moments to determine ∆µ a,Sup and ∆µ a,Deep , and added to them the baseline µ a values. The resulting values of µ a,Sup and µ a,Deep (red and blue crosses in Fig. 6) are close to the nominal values for all 20 ∆µ a steps. These results are comparable to the results of similar experiments that were reported in e.g. [33,36,43]. Possible drifts in the laser output could affect the moments, which may account for the slight deviations from a horizontal line for the determined µ a,Deep in Exp. #2 (blue crosses in Fig. 6 (b)). Another cause could be the µ ′ s changes shown in Fig. 4 (b), which can affect only the results for the deep layer in Exp. 2, based on the findings in Fig. 3 (b).

Experiments involving blood
In this section we present the results of the three experiments involving blood. The results are presented in the order according to the data analysis: the measured DTOFs and the calculated changes in moments (Fig. 7), the determined µ a and µ ′ s (Figs. 8 to 10), the determined concentrations of HbO 2 and Hb (Fig. 11), and finally the determined StO 2 (Figs. 12 and 13).  in the two compartments, at 705 nm for which ∆µ a is high. The ∆µ a due to changes in StO 2 depends on the difference of the molar absorption coefficients (shown in Fig. 8 (a)). For example, close to the isosbestic point, e.g. at 808 nm, ∆µ a is close to zero. An increase in µ a in any layer results in a narrower DTOF, i.e. the curve after the peak decreases at a faster rate. Sato et al. [90] used the slopes of the segments along the time axis to selectively determine the µ a of different layers. The four DTOFs in Fig. 7 (a) support and demonstrate this concept. When only µ a,Sup is increased (green curve), the slope of the curve at earlier time-channels changes, while the slope at the later time-channels (from about 4 ns) appears unchanged. When only µ a,deep is increased (orange curve), the slope of the curve at later time-channels (from about 4 ns) is affected more compared to earlier time channels. The arrival time of photons is related to the mean penetration depth, and this leads to a relationship between µ a at different depths and the slope of the DTOF at different times, as explained in [90]. When µ a of both compartments is maximum (red curve), the curve has a small bump at about 4 ns (about 1 ns after the DTOF's maximum). This bump is caused by the after-pulse peak in the IRF, which occurs about 1 ns after the IRF's maximum and has a magnitude of about 1% of the IRF's maximum, as shown in Fig. 5 (a). Figure 7 (b, c) shows the time traces of changes in moments (∆m 1 and ∆V) relative to time zero, for measurements on Pos. #1 and Pos. #2, at 705 nm for which ∆µ a is high. We found that the shape of the IRF did not vary between different experiments (on different days) and the IRF shifted by less than half of a time-channel (one time-channel was 12.22 ps). Therefore, we used the same time-channels for calculating moments for all three experiments, from 25% to 3% of the maximum as shown in Fig. 7 (a). The time trace of ∆A is not shown as it was not used in data analysis and it required correction for changes in the ND filter, which we changed to obtain good signals at high µ a and avoid oversaturating detectors at low µ a .
The measurements on the deep compartment demonstrate the repeatability of the experiment, as they are independent of L. Noticeably, the rate of desaturation was slightly different in each experiment. For smaller L, the changes measured on the two-layered medium approach those measured on the deep compartment, as expected. By comparing the relative magnitudes of ∆m 1 and ∆V for the two measurement positions, we can confirm that ∆V has a higher depth selectivity [48,88]. When both compartments were desaturated (e.g. at around 43 min), we expect the same changes in moments for measurements on both windows. The observed differences may be due to different µ ′ s in the two compartments (shown in Fig. 10), which was caused by the addition of yeast as visible at 0 and 31 min in Fig. 7 (b, c) and Fig. 10. The values of m 1 and V remained stable until the fourth deoxygenation cycle, after which they gradually decreased following a drift in µ ′ s (Fig. 10). Figure 8 shows the determined µ a and µ ′ s spectra for the deep compartment at maximum and minimum StO 2 , i.e. at first and second baselines, respectively, which are indicated in Fig. 7 (b). Fig. 8. Determined µ a (a) and µ ′ s (b) at maximum and minimum StO 2 for three experiments involving blood, measured on the deep compartment. The molar absorption coefficients (ε) [78] and the µ a of water [79] were used for calculating the concentrations of HbO 2 and Hb.
We chose two baseline regions, relative to which the changes in moments will later be calculated for determining ∆µ a , so that changes in StO 2 occurred in only one compartment at a time. The first chosen baseline region is from 14.6 to 15.6 min, when StO 2 was close to 100% in both compartments. The second region is from 58.8 to 59.8 min, which is after both compartments were fully deoxygenated. We used the second baseline for analyzing measurements after 31 min, which is after yeast was added to the superficial compartment. The µ a spectra (Fig. 8 (a)) resemble the shapes of the molar absorption spectra (ε) of HbO 2 and Hb. The differences are mainly due to the contribution of water's µ a , which should be subtracted from the determined µ a values for calculating concentrations. The µ ′ s values (Fig. 8 (b)) decrease with wavelength, similar to the trend observed for Intralipid [81,84] and tissues [85,86]. The values of µ a and µ ′ s were well-repeatable for the three experiments. Figure 9 shows the time traces of the determined µ a , which were obtained for two measurement positions and two methods of data analysis. The µ a of the deep compartment (black curve) demonstrate five repeatable deoxygenation cycles, with no significant differences between the three experiments. However, after the fourth deoxygenation cycle (after 45 min), the µ a slightly drifts upwards over time, as evidenced by the increasing minimum µ a of different deoxygenation cycles, while µ ′ s drifts downwards (Fig. 10). The maximum µ a values (e.g. at around 30 min) differ for each experiment, with values of approximately 0.167, 0.173, and 0.158 cm −1 for Exp. #1, #2, and #3, respectively. The noise level expectedly increases for higher µ a values, due to the lower count rate. The comparison between the results obtained using the two methods (shown in Fig. 9) is consistent with the findings from the experiments involving ink (shown in Fig. 6). Analyzing the measurements on Pos. #2 using the curve-fitting method with a homogeneous model (green curves), produced µ a values that are between the µ a of the two compartments, and closer to the µ a of the superficial compartment for larger L. The method based on moments accurately determined ∆µ a,Sup and ∆µ a,Deep (magenta and cyan curves) in all cycles except the third cycle (between 31 and 45 min). The accuracy of determined ∆µ a,Deep (cyan curve) can be evaluated by comparing to the measurements on the deep compartment (black curve), and the accuracy of ∆µ a,Sup (magenta curve) can be assessed based on the protocol, i.e. it was constant except during the third cycle. In the third cycle, the method does not account for ∆µ ′ s,Sup caused by the yeast added at 31 min, which resulted in inaccurately determined ∆µ a,Sup and especially ∆µ a,Deep . The method is insensitive to small changes in µ ′ s,Deep , but sensitive to µ ′ s,Sup , as demonstrated in Fig. 3 (b). Figure 10 shows the time traces of the µ ′ s , which were determined together with µ a (black and green curves in Fig. 9) using the curve-fitting method with a homogeneous model. The µ ′ s measured on the deep compartment were stable up to the fourth cycle (45 min) and then began to decrease. This trend was observed for all wavelengths and in all three experiments. The trend is more similar between Exp. #1 and #3. Adding yeast at 0 min increased µ ′ s by approximately 0.3 cm −1 in all three experiments. The µ ′ s values obtained from measurements on the two-layered medium (analyzed with a homogeneous model, which is suitable only when µ a,Sup = µ a,Deep ) are overestimated in the first three deoxygenation cycles (until 45 min, when µ a,Sup ≤ µ a,Deep as shown in Fig. 9) and then underestimated in the last three cycles (after 45 min, when µ a,Sup ≥ µ a,Deep ). This is a consequence of utilizing a homogeneous model for the analysis of measurements on a two-layered medium. Adding yeast at 31 min increased µ ′ s by approximately 0.2 cm −1 . Figure 11 shows the time traces of the concentrations of HbO 2 and Hb in the deep compartment (left column) and in the superficial compartment (right column), which were calculated using the determined µ a at multiple wavelengths (shown in Fig. 9 for one wavelength). The results for measurements on Pos. #1, analyzed using the curve-fitting method with a homogeneous model (red and blue curves in the left panels), are consistent with the expected time traces of the concentrations based on the protocol used. In particular, during each deoxygenation cycle, all (or most) of HbO 2 loses oxygen and becomes Hb (due to yeast consuming oxygen), while the total hemoglobin remains constant (traces of total hemoglobin are not shown). Then, bubbling oxygen results in all (or most) of Hb gaining oxygen and becoming HbO 2 again. The total concentrations differ slightly for the three experiments, following the differences in the maximum of µ a shown in Fig. 9. In the case of Pos. #2, analyzing the measurements using the curve-fitting method with a homogeneous model (red and blue curves in the right panels), provided concentrations of HbO 2 and Hb that are between those of the two compartments and closer to the superficial compartment for larger L.
The concentrations obtained for the method based on moments with a two-layered model are consistent with the expected time traces based on the protocol. The concentrations for the deep compartment (magenta and cyan curves in the left panels) are similar to those obtained using measurements on Pos. #1 (red and blue curves). The concentrations for the superficial compartment (magenta and cyan curves in the right panels) are relatively constant during five cycles and change during the third cycle, as expected based on the protocol. The addition of yeast to the superficial compartment at 31 min had a clear influence on the determined concentrations. Negative values were obtained for HbO 2 in the deep compartment when the determined µ a,Deep had the highest error during the third cycle. Larger L lead to noisier results for the deep compartment and less noisy results for the superficial compartment. The obtained results for HbO 2 are noisier than for Hb for both data analysis methods. The reason is the higher values of ε of Hb compared to ε of HbO 2 at most of the used wavelengths, visible in Fig. 8 (a). A change in Hb has a larger effect on µ a than an equivalent change in HbO 2 . Figure 12 shows the time traces of StO 2 , which were calculated using the concentrations in The StO 2 for the superficial compartment obtained from the two-layered measurements (magenta curve) is maximum for the first two cycles and minimum after the third cycle, which follows the protocol. The StO 2 for the deep compartment (cyan curve) closely resembles the results of measurements on the deep compartment (black curve). However, the values deviate for some of the larger changes in StO 2 relative to the baseline and more so for larger L.
The curves of StO 2 in Fig. 12 resemble the curves of µ a in Fig. 9, which is expected since StO 2 calculation was based on these µ a values along with the µ a values at other wavelengths.  Fig. 11. A moving average window (10 data points) was applied to all curves. Triangles at the top are explained in caption of Fig. 7. Figure 13 shows a comparison of the determined StO 2 values for measurements on the deep compartment and on the two-layered medium. The StO 2 curves (shown in Fig. 12) were smoothed using a moving average window (10 data points), which facilitated plotting curves in Fig. 13 rather than scattered data points. We compared the determined StO 2 values during five deoxygenation cycles, between the maximum StO 2 (after bubbling oxygen) and the minimum StO 2 (plateau at the end of desaturation). The third cycle was excluded from this comparison because the StO 2 in the deep compartment was constant throughout the cycle. This method of comparing simultaneous NIRS measurements was previously applied by Kleiser et al. [63,64] and in other studies [70]. The results show that the two methods produced agreeing values of StO 2 throughout all cycles, as the curves are close to the identity line. However, in the first two cycles, the results obtained at low StO 2 (i.e. furthest away from the baseline StO 2 of approximately 100%) are noisiest and in part deviate from the identity line. In the last three cycles, the agreement is good at all StO 2 values. The results for the deep compartment for measurements on the two-layered medium are noisier for larger L (as expected) and at lower StO 2 (which corresponds to higher µ a for most of the used wavelengths).  Fig. 12) for first (a), second (b), and third (c) experiments involving blood. Measurements on deep compartment were analyzed using curve-fitting method with homogeneous model (plotted along x-axis). Measurements on two-layered medium were analyzed using method based on moments with two-layered model (along y-axis).

Discussion
This work has three main parts: (i) a new method for determining even large ∆µ a in multiple layers, (ii) a new two-layered blood-lipid phantom, and (iii) the published MATLAB codes together with the measured TD-NIRS data for two experiments involving ink and three experiments involving blood.

New method
Previous implementations of analysis based on changes in moments for determining ∆µ a were limited to a range of ∆µ a within which the linearity assumption is valid. Therefore, they were suitable only for fNIRS applications, where small ∆µ a are measured. Figures 2 (a) and 5 (b) show the changes in moments for a range of ∆µ a values, demonstrating the progressive deviation from linearity as the magnitude of ∆µ a increases. In this study, we have experimentally tested the proposed method's capability to accurately determine also large ∆µ a in two layers, i.e. in two compartments of the phantom. The method can now be applied also to NIRS applications that require obtaining absolute µ a , e.g. for determining StO 2 . In well-controlled experiments involving ink, we tested the method's accuracy in determining 20 ∆µ a steps (up to ± 0.10 cm −1 ) in the deep or the superficial compartment. The experiments involving blood demonstrated the ability to determine StO 2 in two layers during dynamic StO 2 changes. For methods that are capable of determining only small ∆µ a , the baseline optical properties had to be updated after a ∆µ a step to satisfy the linearity assumption for determining the next ∆µ a [43]. However, the new method removes the restriction on the range of ∆µ a values that can be accurately determined. Another advantageous feature of the method is its low sensitivity to the selected range of time channels used in the calculations. This is a consequence of using the same range for computing the moments of both the measured DTOF and the simulated DTOF convolved with IRF. In contrast, other methods that utilize moments [75] or changes in moments [41] suffer from errors caused by the limited integration range. These errors are sensitive to the selected range of time channels, as explained by Jelzow et al. [43].
The proposed method requires knowing the baseline optical properties, which we obtained by measuring on each compartment of the phantom separately and using the curve-fitting method with the homogeneous model. However, for analyzing in-vivo measurements, it remains a challenge to develop an accurate and robust approach for estimating the baseline optical properties of two layers as well as thickness L, e.g. [28,[33][34][35][36]. It can be particularly challenging to estimate L using a single DTOF, e.g. if the optical properties of both layers are the same, then any value of L will fit the data, i.e. the DTOF. We demonstrated in Fig. 2 (a) that for obtained changes in two moments, there is a unique solution of ∆µ a in two layers. In fNIRS applications (small ∆µ a ), it is sufficient to know approximate values of the baseline parameters (supported by Fig. 3). This is consistent with the established method based on moments, which requires estimates of the values of all baseline parameters for generating sensitivity factors [41]. However, the estimation of L and baseline optical properties is an important and needed topic for further research.
The accuracy of generating the theoretical DTOFs governs the accuracy of the recovered optical properties. One way to enhance the accuracy of theoretical DTOFs is by using Monte Carlo simulations, which is considered as the gold-standard method for simulating NIRS data [30]. The accuracy can be further improved by using a more realistic geometrical model, e.g. an anatomical head model from MRI data [91].
Sassaroli et al. [92] derived formulas for light propagation to obtain moments directly, eliminating the need to generate DTOFs and thereby reducing computational time. This approach can be particularly advantageous when incorporating changes in moments at multiple sourcedetector distances, as the computational time for the LMA may grow exponentially with the number of source-detector distances used.
The current implementation of the method based on moments uses only two DTOFs (at baseline and after ∆µ a ). To further enhance the method, more DTOFs could be incorporated in data analysis. For example, Jelzow et al. [43] fitted a third-order polynomial to the measured changes in moments before using them to determine ∆µ a , which can improve the signal-to-noise ratio. Yang et al. [74,76] proposed a new approach that exploits measurements at multiple source-detector distances, and demonstrated the benefits of incorporating spatial information in addition to time-domain measurements, which improved reliability and reduced uncertainty in determining the absolute optical properties of multiple layers. While the present study used measurements at a single source-detector distance, it is important to note that the proposed method can be extended to use changes in moments at multiple source-detector distances.
We did not use the intensity information (∆A) for determining ∆µ a because it required careful data processing to account for potential drifts and changes in the output of the laser. Incorporating ∆A information could greatly improve data analysis. We weighted changes in moments by their uncertainty due to photon noise (Poisson statistics) in Eq. (3), to account for the different units and varying levels of uncertainty. A more accurate estimate of the uncertainty could incorporate also the instrumental noise for each moment, which can be measured once for a system [43].
We did not report the uncertainties of the determined optical properties in this study, as it is dependent on several factors including the selected laser power and the acquisition time. The uncertainty of the fitted parameters (for the LMA method) due to photon statistics can be obtained using the method described in [93].
When calculating concentrations of oxy-and deoxyhemoglobin, we used the determined µ a values at multiple wavelengths, but without accounting for the differences in their uncertainties. The determined µ a at different wavelengths have varying levels of uncertainty, which relates to the µ a spectrum of the medium (higher uncertainty for higher µ a ) and to the system used (e.g. if different wavelengths have different power or detection sensitivity). Future studies should aim to incorporate information about uncertainties to improve the determination of hemoglobin concentrations and hence StO 2 , e.g. as in [35].

Measured data and MATLAB codes
The measured TD-NIRS data and the MATLAB codes used in this study were made publicly available [68]. The data is suitable for testing new data analysis methods for tissue oximeters.
The LMA is a search-based algorithm that allows entering the values of known parameters and then determining the unknown parameters. We implemented two LMA algorithms, both of which have the feature of allowing the user to choose which parameters to determine. The first algorithm uses the DTOF as an input and determines any of the eight parameters: µ a , µ ′ s , and/or L (layer's thickness), for up to three layers (the deepest layer is assumed semi-infinite). We used it to determine µ a and µ ′ s assuming a homogeneous medium. The second algorithm uses changes in moments as an input and determines any of the following: ∆µ a , ∆µ ′ s , and/or L, for up to three layers. We used it to determine ∆µ a,Sup and ∆µ a,Deep of a two-layered medium. Future studies may be interested in using this code to determine ∆µ a in two layers (superficial and deep) of a three-layered medium, with the assumption that the µ a of the middle layer is constant (mimicking the skull). Additionally, the codes can be modified to use DTOFs (in the first algorithm) or changes in moments (in the second algorithm) at multiple source-detector distances.

Phantom
The phantom allows repeatedly deoxygenating and oxygenating blood in two compartments individually (the results for oxygenating the superficial compartment were not presented), and carrying out simultaneous NIRS measurements on any of the compartments: the two-layered medium (variable superficial layer thickness), the deep (semi-infinite medium), and/or the superficial (semi-infinite medium). The measurements on individual compartments can be analyzed employing a homogeneous model, providing a reference for the optical properties within each compartment, as was done in this study (Fig. 6). The repeatability within and between experiments was demonstrated by carrying out six deoxygenation cycles in three experiments (18 cycles in total). The constructed phantom enables carrying out reproducible experiments, and has good durability and ease of use. The hemoglobin concentrations and the corresponding StO 2 levels in each compartment followed the expected desaturation trends (StO 2 decreased from 100% to 0%) during all deoxygenation cycles.
The rate of desaturation varied slightly across the three experiments, as seen from the changes in moments in Fig. 7. One potential cause could be a slight difference in the amount of yeast that was added to the phantom. After pouring the mixture containing yeast into the phantom, a small amount of yeast may have adhered to the walls of the beaker, leading to an uncertainty in the amount of yeast added. The temperature of the phantom affects the rate of desaturation, but the temperature was stable as measured with two thermometers.
We found that adding yeast to the solution resulted in an increase in µ ′ s (Fig. 10). After the start of the fourth cycle (after 45 min), the µ ′ s in the deep compartment began to steadily decrease (visible in Figs. 7 and 10). This decrease had more similar trends in the first and third experiments, in which the same blood bag was used. The decrease in µ ′ s,Deep may be related to either the yeast, which was mixing for over 45 minutes, or the blood, which was kept in a desaturated state for an extended period of time (during the third cycle). Ostojic et al. [94] studied the gradual decrease of µ ′ s over time in phantoms with blood and proposed possible explanations.
The large size of the phantom requires large amounts of ingredients, which increases the cost of conducting experiments, but more importantly, it makes the phantom easier to control. Due to the large volume, the temperature of the phantom changes at a slow rate, making it easier to maintain it at approximately 37°C. Additionally, the experiments can be carried out with higher accuracy when handling larger amounts of ingredients, improving the reproducibility of experiments. For example, we added about 3.5 g of Mixture #2 for each step of ∆µ a , whereas in another two-layered phantom, the authors needed to add only about 0.7 g to achieve similar ∆µ a [43].
The measurement windows and/or the separating window of the phantom (shown in Fig. 1) could be covered with a layer of material that has tissue-mimicking optical properties for mimicking the skin and/or the skull. This could be achieved using the silicone-based layers presented by Kleiser et al. [63,64] and/or the solid layers presented by Izzetoglu et al. [58], which can be created with variable tissue-mimicking optical properties and thicknesses.
One of the limitations of the study is the absence of an independent reference measurement of StO 2 , which limits the assessment of the accuracy of the determined StO 2 values since the true values are unknown. In this study, we compared two measurements of StO 2 (on the deep compartment and on the two-layered medium), both of which were measured by the same system. Kleiser et al. [63,64] showed that different NIRS tissue oximeters can produce different StO 2 values, even when they simultaneously measure on the same phantom. Furthermore, the StO 2 values also depend on the used method of data analysis. The partial pressure of oxygen (pO 2 ) was proposed as a reference measure in previous studies, but it has major challenges [63,67]. Obtaining reliable reference values of StO 2 would be a significant advancement in the assessment of accuracy, and for calibrating different systems and methods of data analysis.

Conclusions
We introduced and experimentally tested a new method for accurately determining even large ∆µ a in multiple layers, which utilizes (i) changes in moments (making use of two DTOFs, at baseline and after ∆µ a in two layers, and the knowledge of the optical properties at baseline), (ii) an analytical solution of the diffusion equation for a multi-layered medium, and (iii) the Levenberg-Marquardt algorithm. We applied the method for tissue oximetry for determining StO 2 in two layers (the absolute µ a was obtained by summing the baseline µ a and the determined ∆µ a ).
We used the established concept of a blood-lipid phantom and constructed a phantom with two such compartments (mimicking superficial and deep layers, with a variable superficial layer thickness), allowing for the first time, independent control of StO 2 of blood in two compartments. The phantom may be used for performance assessment, such as testing NIRS tissue oximeters or new methods of data analysis, as well as for studying the influence of the superficial compartment on the estimation of StO 2 in the deep compartment.
The results of five experiments, two of which varied the concentration of ink and three of which varied the StO 2 of blood, demonstrated a precise control of the phantom and a high reproducibility of deoxygenation cycles within and between experiments. Also, the accurate determination of µ a and/or StO 2 in both compartments, using time-domain NIRS data measured on two-layered medium, was confirmed experimentally. Disclosures. The authors declare no conflicts of interest. MW declares that he is president of the board and founder of OxyPrem AG. AL declares that he is shareholder of Brain Optics. Data availability. Data and MATLAB codes used in this study were made available under the GNU General Public License v3.0 [68].