1. 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–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–15]). A significant improvement to NIRS tissue oximetry would be the determination of StO₂ 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–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–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–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–46]. Recent performance-assessment studies, including quantitative tests, found that moments outperform other measurands in terms of assessing Δµ_a in deep layers [47–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–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–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₂) levels from 0 to 100%. To date, all presented phantoms that used this concept had a single homogeneous liquid compartment, e.g. [63–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₂ [18–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₂ in a two-layered medium. We will utilize the phantom to demonstrate the proposed method’s capabilities of determining StO₂ in both layers.

Data and MATLAB codes used in this study were made publicly available [68].

2. Methods and materials

2.1 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₂. 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² on a skin's surface [70]. For a recent comprehensive study on IRF, see the study conducted by Pirovano et al. [73].

2.2 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.

2.2.1 Determination of µ_a and µ′_s (homogeneous model)

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.

2.2.2 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_{\textrm{tot}}} = \sum {{N_i}} $ (total number of photons), ${m_1} = \sum {{t_i}{N_i}} /{N_{\textrm{tot}}}$ (mean time of flight), and $m_2^C = V = \sum {{{({t_i} - {m_1})}^2}{N_i}} /{N_{\textrm{tot}}}$ (variance). The change in attenuation is: $\Delta A ={-} \log ({N_{\textrm{tot}}^\ast{/}{N_{\textrm{tot}}}} )$, where $N_{\textrm{tot}}^\ast $ 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_2^C$, and $m_3^C$ (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₁, 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 (χ²), 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, χ² can be calculated using the equation:

To assess the method, we will calculate and plot the distribution of $(2/{N_{\textrm{tot}}}){\chi ^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 χ² 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].

2.2.3 Calculation of concentrations of oxy- and deoxyhemoglobin

To determine the concentrations of oxyhemoglobin (HbO₂) and deoxyhemoglobin (Hb) in experiments involving blood, we used the determined µ_a values at multiple wavelengths and the known molar absorption coefficients (ε) of HbO₂ 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₂ 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.

2.3 Phantom setup

2.3.1 Construction of phantom

In recent years, various tissue-mimicking phantoms containing blood have been developed, e.g. [59,62–67]. However, to our knowledge, no presented phantom has two blood-lipid compartments with variable StO₂, 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.

 

Fig. 1. Diagrams of the phantom showing top view (a) and side view (b), and the corresponding model (c). The three containers are labeled: Deep, Superficial, and Third container, which is movable. Units are mm.

Download Full Size | PDF

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].

2.3.2 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]:

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.

2.3.3 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.

Table 1. Ingredients (in grams) and targeted µ_a values (at 750 nm) for the two experiments involving ink. SL is SMOFlipid.

View Table | View all tables in this article

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.

2.3.4 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₂) 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₂. 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.

Table 2. Protocol steps and the targeted StO₂ in two compartments for the three experiments involving blood.

View Table | View all tables in this article

Table 3. Ingredients for each of the three experiments involving blood. SL is SMOFlipid.

View Table | View all tables in this article

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.

3. Results

3.1 Method for determining Δ${\mu_{\mathrm{a}}} $ in two layers

Here we examine the proposed method for determining Δµ_a in two layers. We calculated the values of ΔA, Δm₁, 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.

 

Fig. 2. Changes in moments (ΔA, Δm₁, and ΔV) for all combinations of Δµ_a,Sup and Δµ_a,Deep (a). Error norm (χ²) for all combinations of Δµ_a,Sup and Δµ_a,Deep (white cross marks true values) when using changes in three moments (b) or using only Δm₁ and ΔV (c). Baseline parameters are µ_a,Sup = µ_a,Deep = 0.13 cm^-1, µ′_s = 10 cm^-1, L = 12 mm, and n = 1.33.

Download Full Size | PDF

We calculated the error norm (χ²) using Eq. (3) for all combinations of Δµ_a,Sup and Δµ_a,Deep. The distribution of χ² 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₁ and ΔV in Eq. (3). The presented χ² 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 χ² 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 χ² 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₁ becomes dominant, and the shape of the resulting χ² distribution resembles the line where Δm₁ = 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₁, and 10 for ΔV. Therefore, for a typical Δµ_a, the term with ΔA in Eq. (3) contributes much more to the value of χ² than the term with Δm₁, which is consistent with the one-order magnitude difference in the colorbars for χ² in Fig. 2 (b) and (c). If only a single moment is considered in Eq. (3), the shape of the χ² distribution would follow the moment’s zero contour line shown in Fig. 2 (a) and the value of χ² 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.

 

Fig. 3. Determined Δµ_a,Sup and Δµ_a,Deep for different assumed values of baseline parameters:

µ_a (a), µ′_s (b), n (c), and L (d). Green lines mark true values: Δµ_a,Sup = 0 and Δµ_a,Deep = 0.05 cm^-1, µ_a = 0.1 cm^-1, µ′_s = 10 cm^-1, L = 15 mm, and n = 1.33.

Download Full Size | PDF

3.2 Experiments involving ink

Figure 4 (a) shows the determined µ_a for two experiments involving ink, for all wavelengths, for both compartments, for three Δµ_a steps. The µ_a values for both experiments are similar and consistent with the targeted µ_a values in 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).

 

Fig. 4. (a) Determined µ_a for two experiments involving ink, for three nominal µ_a (0.05, 0.15, and 0.25 cm^-1 at 750 nm). (b) Determined µ′_s for 20 Δµ_a steps for Exp. #1. Results were obtained using measurements on each compartment (Pos. #1 or #3 in Fig. 1), analyzed using the curve-fitting method with a homogeneous model.

Download Full Size | PDF

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⁻¹ (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⁻¹) 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).

Figures 5 and 6 present results at 768 nm (8^th spectral channel) and similar results were obtained at other wavelengths. Figure 5 (a) shows the normalized IRFs and DTOFs for the two source-detector pairs located on Pos. #1 and #2, when the targeted µ_a was 0.10 cm^-1 at 750 nm in both compartments. The two source-detector pairs had similar IRFs and hence similar DTOFs. The total number of photons in the DTOF (N_tot) was between approximately 5.4 to 0.5 million photons at 768 nm for different Δµ_a steps, with lower values of N_tot at higher µ_a.

 

Fig. 5. (a) Normalized IRFs and DTOFs (when targeted µ_a = 0.10 cm^-1 at 750 nm in both compartments), for two source-detector pairs on Pos. #1 and #2, at 768 nm. (b) Changes in moments (ΔA, Δm₁, and ΔV) for 20 Δµ_a steps, measured on Pos. #2 (two-layered), at 768 nm.

Download Full Size | PDF

 

Fig. 6. Determined µ_a for 20 Δµ_a steps in first (a) and second (b) experiments involving ink, at 768 nm. Measurements on a compartment with increasing µ_a (Pos. #1 or Pos. #3) were analyzed using the curve-fitting method with a homogeneous model (black crosses). 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).

Download Full Size | PDF

Figure 5 (b) shows the changes in moments (ΔA, Δm₁, 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 (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).

3.3 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₂ and Hb (Fig. 11), and finally the determined StO₂ (Figs. 12 and 13).

 

Fig. 7. (a) DTOFs measured on Pos. #2 at different µ_a in two compartments, at 705 nm. Changes in moments, Δm₁ (b) and ΔV (c), measured on Pos. #1 and #2, for three experiments involving blood. Filled triangles indicate when yeast was added, separately in deep and superficial compartments, and hollow triangles denote the start of bubbling oxygen in the deep compartment (these triangles are also shown in Figs. 9 to 12).

Download Full Size | PDF

Figure 7 (a) shows DTOFs measured on Pos. #2, at different StO₂ values in the two compartments, at 705 nm for which Δµ_a is high. The Δµ_a due to changes in StO₂ 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₁ 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.

 

Fig. 8. Determined µ_a (a) and µ′_s (b) at maximum and minimum StO₂ 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₂ and Hb.

Download Full Size | PDF

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₁ 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₁ 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₂, i.e. at first and second baselines, respectively, which are indicated in Fig. 7 (b). We chose two baseline regions, relative to which the changes in moments will later be calculated for determining Δµ_a, so that changes in StO₂ occurred in only one compartment at a time. The first chosen baseline region is from 14.6 to 15.6 min, when StO₂ 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₂ 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.

 

Fig. 9. Determined µ_a at 705 nm in first (a), second (b), and third (c) experiments involving blood. Measurements were analyzed using the curve-fitting method with a homogeneous model (Homog.) and using the changes in moments with a two-layered model (2-L). For visualization, a moving average window (10 data points) was applied to all curves. Triangles at the top are explained in caption of Fig. 7.

Download Full Size | PDF

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.

 

Fig. 10. Reduced scattering coefficient µ′_s in three experiments involving blood, at 705 nm (a) and 805 nm (b). Measurements were analyzed using the curve-fitting method with a homogeneous model and the corresponding determined µ_a values are black and green curves in Fig. 9. Triangles at the top are explained in caption of Fig. 7.

Download Full Size | PDF

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₂ 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₂ 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₂ 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₂ and Hb that are between those of the two compartments and closer to the superficial compartment for larger L.

 

Fig. 11. Concentrations of HbO₂ and Hb in three experiments (three rows) involving blood. Curve-fitting method (Homog.) was used for measurements on Pos. #1 and on Pos. #2, and the latter was also analyzed with method based on moments (2-L. Deep and 2-L. Sup). For visualization, a moving average window (10 data points) was applied to all curves. Triangles at the top are explained in caption of Fig. 7.

Download Full Size | PDF

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₂ 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₂ are noisier than for Hb for both data analysis methods. The reason is the higher values of ε of Hb compared to ε of HbO₂ 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₂.

Figure 12 shows the time traces of StO₂, which were calculated using the concentrations in Fig. 11 and the formula StO₂ = C_HbO2 / (C_HbO2 + C_Hb), where C stands for concentration. Most StO₂ values are between 0 and 100%, except in the third cycle where negative values of C_HbO2 were obtained resulting in negative values of StO₂. The StO₂ in the deep compartment (black curve) follows five deoxygenation cycles from about 100% to 0%, as expected, but the maximum StO₂ slightly reduces with each deoxygenation cycle in all three experiments. The StO₂ values obtained from Pos. #2, analyzed with a homogeneous model (green curve), are between the StO₂ values of the two compartments and approach the StO₂ of superficial compartment for larger L.

 

Fig. 12. Determined StO₂ in first (a), second (b), and third (c) experiments involving blood for two measurements and two methods of data analysis. The concentrations corresponding to the four StO₂ curves are shown in 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.

Download Full Size | PDF

The StO₂ 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₂ 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₂ relative to the baseline and more so for larger L.

The curves of StO₂ in Fig. 12 resemble the curves of µ_a in Fig. 9, which is expected since StO₂ calculation was based on these µ_a values along with the µ_a values at other wavelengths.

Figure 13 shows a comparison of the determined StO₂ values for measurements on the deep compartment and on the two-layered medium. The StO₂ 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₂ values during five deoxygenation cycles, between the maximum StO₂ (after bubbling oxygen) and the minimum StO₂ (plateau at the end of desaturation). The third cycle was excluded from this comparison because the StO₂ 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₂ throughout all cycles, as the curves are close to the identity line. However, in the first two cycles, the results obtained at low StO₂ (i.e. furthest away from the baseline StO₂ 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₂ 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₂ (which corresponds to higher µ_a for most of the used wavelengths).

 

Fig. 13. Comparison of StO₂ time traces (from 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).

Download Full Size | PDF

4. 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.

4.1 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₂. In well-controlled experiments involving ink, we tested the method’s accuracy in determining 20 Δµ_a steps (up to ± 0.10 cm⁻¹) in the deep or the superficial compartment. The experiments involving blood demonstrated the ability to determine StO₂ in two layers during dynamic StO₂ 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–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 source-detector 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₂, e.g. as in [35].

4.2 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.

4.3 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₂ levels in each compartment followed the expected desaturation trends (StO₂ 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₂, which limits the assessment of the accuracy of the determined StO₂ values since the true values are unknown. In this study, we compared two measurements of StO₂ (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₂ values, even when they simultaneously measure on the same phantom. Furthermore, the StO₂ values also depend on the used method of data analysis. The partial pressure of oxygen (pO₂) was proposed as a reference measure in previous studies, but it has major challenges [63,67]. Obtaining reliable reference values of StO₂ would be a significant advancement in the assessment of accuracy, and for calibrating different systems and methods of data analysis.

5. 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₂ 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₂ 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₂ in the deep compartment.

The results of five experiments, two of which varied the concentration of ink and three of which varied the StO₂ 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₂ in both compartments, using time-domain NIRS data measured on two-layered medium, was confirmed experimentally.