Open Access
Add to BibSonomy
Abstract
A multivariate method integrating time and space resolved techniques of near-infrared spectroscopy is proposed for simultaneously retrieving the absolute quantities of optical absorption and scattering properties in tissues. The time-domain feature of photon migration is advantageously constrained and regularized by its spatially-resolved amplitude patterns in the inverse procedure. Measurements of tissue-mimicking phantoms with various optical properties are analyzed with Monte-Carlo simulations to validate the method performance. The uniqueness, stability, and uncertainty of the method are discussed.
© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
Near-infrared spectroscopy (NIRS) has become a powerful tool for clinical diagnosis and medical imaging [1], owing to its non-invasiveness, high temporal resolution, and portability. It can reveal the optical absorption property of biological tissues in depth, and thereby has been widely applied for indicating the concentrations of essential biomarkers, e.g., oxy-/deoxy-hemoglobin and cytochrome-c-oxidase in human brain, kidney, and breast [2-4]. The advantages of this technology in turn also allow for more sophisticated clinical applications such as neuroimaging, diffuse optical tomography, and tissue oximetry [5–7]. By quantifying optical scattering properties, NIRS-based real-time monitoring of many in vivo mechanisms [8,9] (e.g. intracellular oxygen delivery) can also become a feasible future.
Focusing on changes in concentration of oxy- and deoxyhemoglobin in human brain, functional NIRS (fNIRS) has received an enormous amount of attention in decades [10]. This technique studies on hemodynamic and metabolic responses to human brain activities and performs on task-based or event-related experiments. It derives many fruitful ramifications in different fields, e.g. clinical cerebral metabolism [11], cognitive neuroscience [12], and brain-computer interface studies [13].
On the other hand, absolute NIRS measuring methodologies, which aim at retrieving the absolute quantitative values of optical properties in biological tissues and materials instead of changes only, are promoting progress and drawing more and more interest through the years. Absolute NIRS techniques have already been shown to be a valuable tool for non-invasive diagnosis of many pathological tissues such as breast [14], thyroid [15], traumatic brain [16], skin [17], and many others [18]. Specifically, cerebral oximeters based on absolute NIRS have found widespread use. However, most of the devices are empirically calibrated and the accuracy and comparability of results from different devices is insufficient [19]. More potentials of absolute NIRS require further exploration. By retrieving two essential optical parameters, i.e. absorption coefficient (µa) and reduced scattering coefficient (µs′), absolute NIRS methods can directly derive the absolute quantities of biomarkers’ concentrations, and in particular, oxygen saturation of hemoglobin. In order to obtain absolute values of µa and µs′, absolute NIRS methods rely on various measurands, such as spatially resolved (SR) (also termed multi-distance) continuous-wave (CW) amplitude in space-domain (SD) [20], time resolved (TR) photon distribution in time-domain (TD) [21], or demodulation and phase shift in frequency-domain (FD) [22]. The recently introduced technique of interferometric NIRS, which is based on the time-of-flight resolved measurement of coherent light, can also determine absolute values of µa and µs′ [23,24]. However, although from physics perspective µa and µs′ themselves are two totally independent quantities, in practice the external characteristics originating from the internal nature of µa and µs′ are heavily entangled in all above-mentioned measurands. Such entanglement causes that errors in one retrieved quantity can manifest as artefacts in the other and reduce the accuracy for both. Therefore, the key issue of absolute NIRS methods is to decouple absorption and scattering in detected signals.
Many analytical models (e.g. diffusion approximation, P-n approximations) and numerical tools (e.g. Monte-Carlo, finite element methods) have been developed to describe photon migration in tissues, nevertheless simultaneous and accurate determination of µa and µs′ still remains a challenging but promising task. SD approaches suffer from the non-uniqueness [25], while TD and FD approaches have difficulties regarding crosstalk among the unknowns [26] and unsatisfactory performance under low SNR circumstances like short-time acquisition or large source-detector distance [27]. Ambiguity and uncertainty of retrieved optical properties are particularly high in the presence of in vivo noise and perturbations, or in sophisticated heterogeneous structures. Some moment-based spatially-resolved TD attempts [28,29] have already been applied to handle these problems. More practical improvements aiming to increase uniqueness, stability and accuracy for clinical routine implementations of NIRS are crucial.
In this work, a new multivariate method integrating TD and SD techniques of NIRS, namely spatially-enhanced time-domain NIRS, is proposed, and tested for simultaneous quantification of µa and µs′ in homogenous turbid materials. The method is based on the principle that optical absorption and scattering effects on diffuse reflectance, despite independent, are either positively correlated in TD observations or negatively correlated in SD observations. A similar concept has been indicated by S. Arridge and M. Schweiger in theoretical basis papers [30,31]. Effective mutual complementation of TD and SD information of photon migration is realized by enhancing and optimizing TD fitting with SD constraints. Thus, a better convergence to the optimal solution in {µa, µs′} space is reachable. We first demonstrate the proof-of-principle of the method by Monte-Carlo simulations. Then the method is validated through experiments of retrieving µa, and µs′ for a set of solid homogeneous phantoms. The method could be implemented to improve the performance assessment of clinical diffuse optics instrumentation when using tissue simulating homogenous phantoms [32]. And by decreasing uncertainty, spatially-enhanced TD NIRS is capable of coping with in vivo conditions with high noise and perturbations. Especially at large source-detector separations where the low photon appearance is prevailing, as other researchers previously pointed out, it is feasible to use homogenous models to estimate the optical properties of deep regions in heterogeneous layered structures [33], which therefore may allow spatially-enhanced TD NIRS to be implemented for circumventing the influence from superficial tissues.
By measuring a series of solid tissue-mimicking homogenous phantoms, we compared the results of the proposed spatially-enhanced TD method with two conventional TD methods (single-distance TD and Global fitting TD). Considerable improvements regarding uniqueness, stability and uncertainty of the retrieved optical properties are demonstrated.
2. Theory
The contributions from absorption and scattering properties of turbid media on the diffuse reflectance are entangled in TD and SD with different ways. Effective complementation of TD and SD information can improve the retrieval performance of optical properties.
2.1 Positive correlation in time domain
In TD, time-resolved diffuse reflectance profiles are measured as histograms that represent distributions of times of flight (DTOF) curves. They are obtained by recording the time of flight from source to detector of many (typically 105 to 106) photons, by means of the time-correlated single photon counting technique. For homogeneous turbid media, the shape of the DTOFs depends on µa, µs′ and ρ (source-detector distance) [34]. An increase in µs′ or ρ broadens the DTOF, while an increase in µa narrows it. In other words, at constant ρ, the contributions from µa and µs′ are positively correlated with each other on DTOFs. The influence of µs′ and ρ on the DTOF shape is given by the product ρ2µs′ [35]. This means a nearly identical DTOF shape (disregarding amplitude) can be obtained if an increase in µs′ is balanced by a decrease in ρ. On the other hand, µa determines the DTOF slope at later time range of DTOF. The uncertainty of µa obtained by TD fit can be reduced by covering the DTOF up to late times. In the present work, this was achieved by extending the fit range to 0.1% of the DTOF’s peak amplitude. The positive correlation of µa and µs′ can also be demonstrated by considering variance (V), the second centralized statistical moment of DTOFs. In infinite media (as well as semi-infinite media with large ρ) the following relation for variance holds [36]:
where ρ is source-detector distance and c is the speed of light in the medium. According to Eq. (1), to keep V constant, i.e. the shape of DTOFs could maintain, the ratio µs′/µa3 must remain unchanged. This means for constant ρ, since µa and µs′ enter with different power, a small relative variation of µa needs to be compensated by a 3 times bigger relative variation of µs′, although in the same direction. (Note that the factor of 3 is specific for variance and does not necessarily apply to other characteristics of DTOFs).In the procedure of searching (µa, µs′) solutions and optimal fitting for TD measurements, the DTOFs are compared with corresponding theoretical curves which are convoluted with the instrument response functions (IRFs). For this comparison, the convolved theoretical curves are scaled to have the same integral as the measured DTOFs, whereby the amplitude information is lost. With the increase of µa, the curve width will become narrower whereas increasing of µs′ leads to a broadened curve. An example regarding positive correlation of contributions from µa and µs′ on normalized reflectance DTOFs will be shown in Fig. 1(a).
Fig. 1. Monte-Carlo simulations of 5 adjacent (µa, µs′) combinations. Their reflectance profiles in (a) TD: DTOFs, (b) SD: SRACs, and the χ2 distributions on the {µa, µs′} space based on: (c) TD fit, (d) SD fit, and (e) Spatially-enhanced TD method. The white X is the true value (0.127, 13.9) cm−1, and the crosses in other colors are the optimal solutions of respective methods.
2.2 Negative correlation in space domain
In SD, the contribution of scattering and absorption properties of turbid media on the diffuse reflectance profiles, i.e. spatially-resolved amplitude curves (SRACs), are negatively correlated with each other. Based on the diffusion equation under semi-infinite geometry and zero boundary condition, it can be derived that if (ρµs′) »1, i.e. the detector is far from the source, the time-resolved reflectance can be written as [21]:
In practical experiments, the light attenuation A = In(R0/R) is the relevant measurand and commonly used. By introducing A into Eq. (3) and differentiating it with ρ, one can get [37]:
To maintain the shape of a SRAC, the changes originated from µa and µs′ to SRACs will cancel each other out. Any increase of µa should be compensated by a same extent decrease of µs′, and vice versa. An example will be shown in Fig. 1(b).
2.3 Spatially-enhanced time domain NIRS
Monte-Carlo simulations (details in Section 3.3) are used here to elaborate and validate the proposed spatially-enhanced time domain method. One (µa, µs′) combination, I (0.127, 13.9) cm−1, and 4 others at vertices of its one surrounding square in {µa, µs′} space, i.e., II (0.134, 14.6), III (0.120, 13.2), IV (0.120, 14.6), and V (0.134, 13.2) cm−1, are selected to simulate the diffuse reflectance. The geometry in the simulations is a 5-cm thick slab, and source-detector distance is 2 cm. As shown in Fig. 1(a)–1(b), the normalized DTOF of I is much more similar to those of II and III than of IV and V. On the contrary, for SD the normalized SRACs of I, IV and V, which have nearly identical µeff, are much harder to split then those of II and III.
A more global view can be seen in {µa, µs′} space, as illustrated by the χ2 distributions plotted in Fig. 1. After adding Poisson noise, simulated DTOF and SRAC of I, (0.127, 13.9) cm−1, are compared with their counterparts of all potential (µa, µs′) combinations in a selected {µa, µs′} space (µa: 0.09–0.17 cm−1, µs′: 10–18 cm−1). The fit results, i.e. the distributions of χ2 values in {µa, µs′} space, are generated for TD in Fig. 1(c) and SD in 1(d). The white X is the true value of I (0.127, 13.9) cm−1, while the green (0.1305, 14.11) cm−1and blue (0.1342, 13.4) cm−1 crosses are the points where the χ2 value is minimal for TD and SD fitting, respectively. The deviations from the true value arise from the introduced noise. The red lines are the contours where the χ2 values are 2 times minimal χ2.
Clearly, the ellipse-like contour in Fig. 1(c) refers to the low-χ2 region surrounding the minimal point. Although there is a convergence in TD, i.e. a unique minimum, the presence of noises or perturbations can cause the ambiguity and move the result away from the true value, due to the small χ2 gradient within the ellipse contour. Meanwhile, the steep elongated ‘canyon’ in Fig. 1(d) refers to the low-χ2 region where all (µa, µs′) combinations can generate SRACs that are indistinguishable from I with added Poisson noise, indicating the reciprocal relation of µa and µs′ in SD and the solution’s non-uniqueness. Given these opposite characteristics, complementing TD and SD information effectively is promising to reduce ambiguity and improve uniqueness at the same time.
The above-mentioned purposes are realized by the isoperimetric inequality theorem as follows. First, the original χ2 distributions from TD and SD are normalized to their minima, then merged in an artificial spatio-temporal χST2 distribution which contains all features from high dimensional spaces.
The new χST2 distribution is the projection of all characteristics from TD and SD on the {µa, µs′} space. λ is Generalized Lagrange multiplier and will be optimized by iterations.The optimization aims at looking for the optimal λ which balances TD and SD information and eliminates the entanglement of µa and µs′, as well as at ensuring the solution’s uniqueness. To this end, the isoperimetric inequality theorem is applied. Based on analyzing the shape of any given closed contour line of χST2 distribution, the theorem states that the ratio of area and squared perimeter can reach maximum only when the shape is a circle. More specifically, the objective function for optimization is chosen to be the isoperimetric quotient Q of χST2 contours:
Fig. 2. The evolution of χ2 contour along with iterations. The solid lines are contours of 2 times minimal χ2. The dotted lines are the connected extrema of µa and µs′ at each contour. r denotes the Pearson correlation coefficient of the points on the red contours.
By using the complementarity of TD and SD information, the goal is to neutralize the correlations amongst the contributions from µa and µs′, and to improve convergence and uniqueness. As shown in Fig. 2, at the 1st iteration, the original χ2 contour presents as an oblique ellipse, when λ=0 and Pearson correlation coefficient of the points on the red contour is r = 0.8578, representing the TD positive correlation. Along with iterations, λ is updated to introduce appropriate SD information, and r is gradually reduced. The optimization evolution suppresses the deviations around the optimal convergences. To the 36th iteration, the optimal λ=0.1556 can be found to achieve r≈0, where a more well-defined convergence is reached. As already shown in Fig. 1(e), the black cross (0.1287, 13.92) cm−1 based on the spatially-enhanced TD method is much closer to the true value of I than the TD and SD solutions alone.
Essentially, applying isoperimetric inequality on the inverse procedure is based on the nature of µa and µs′. The two coefficients, from their physical definitions, are independent to each other and shouldn’t correlate. Decoupling µa and µs′ in an artificial spatio-temporal χST2 is equivalent to separating them and approaching their true values. Furthermore, the result’s uniqueness also restricts that the shape of such homocentric χST2 distribution must be a circle, where the isoperimetric inequality theorem at approaching to Q = 1 guarantees the only solution.
The orthogonality of two coordinates on optimized contours indicates that the cross-talk among µa and µs′ is suppressed. The steepness of SD contour is paved by the smoothness within TD contour, illuminating ambiguity among solutions is eliminated. Under same level of noise, the shrinking contour reveals a more confined solution.
3. Material and methods
3.1 Phantoms
The phantoms used in the experiments are five solid homogenous cylinders. The diameter and height of the cylinder are 10.5 cm and 5 cm. They were manufactured by using epoxy resin (ME500, Nils Malmgren, Sweden) as base material, aliphatic amine (H179B, Nils Malmgren, Sweden) as hardener, TiO2 powder (T-8141, Sigma-Aldrich, USA) as scatterer, and toner (4551, Infotec, UK) as absorber [38]. The amounts of absorber and scatterer added vary sequentially among the phantoms to realize gradual changes in each of the optical parameters while keeping the other constant. The details are shown in Table 1.
Table 1. The components in 5 solid homogenous phantoms measured in the experiments. Columns Scatterer and Absorber denote the added amount of scatter stock solution (1:10 mixture of TiO2 and resin) and absorber stock solution (1:100 mixture of toner and hardener).
According to the added amounts of TiO2 and toner, phantoms α2, β2, and γ2 should have nearly identical µs′ (around 9 cm−1 at 800 nm), whilst phantoms γ1, γ2 and γ3 should have nearly identical µa (around 0.125 cm−1). The final results for the optical properties of the phantoms are given in the result section below. In contrast to liquid phantoms based on intralipid and ink [39], the optical properties of epoxy based solid phantoms cannot be simply derived from the amount of the added components. In particular, the ways (e.g. mechanical, ultrasound) of mixing and dispersing the solid scatterers and absorbers and related heating and cooling processes will substantially influence the final properties. A well accepted method to determine the optical properties of solid phantoms is the application of time-resolved measurements with a single source and a single detector position. Below, we go beyond this way by applying three different TD approaches and discussing the optical properties as well as deviations between the methods.
3.2 Experimental setup
The experimental setup is a standard TD NIRS lab instrument with an SD extension. As schematically depicted in Fig. 3, a supercontinuum laser (SC500-6, Fianium, UK) equipped with an acousto-optical tunable filter (ATOF) is used as the light source, operated at 800 nm and providing picosecond light pulse with a repetition rate of 40.5 MHz. The output beam is guided by a 1.6 m multimode graded index fiber GI-Fiber 1 (core diameter 400 µm, NA 0.27, Leoni, Germany). The light is diffused in the phantoms, re-emitted from their surfaces, and then collected by the other 1.6 m multimode graded index fiber GI-Fiber 2 (core diameter 600 µm, NA 0.22, Leoni, Germany). The position of GI-Fiber 1 is fixed by a holder on the table. GI-Fiber 2 is mounted on a vertical linear stage (PT1/M, Thorlabs, USA) to adjust the height of the fiber end, and this linear stage is mounted on a horizontal long-range precision linear translation stage (length: 120 mm, OWIS, Germany) to move the detector fiber on the surfaces of phantoms to realize source-detector separations: from 15 to 25 mm with 1 mm step size. Both fibers’ ends are carefully adjusted to the same height, vertical and close to phantom surfaces. Black sponge rings and cubes are stuffed around and between the fiber ends to block stray light travelling along the surface as well as residual ambient light. The phantoms are placed on a platform adjusted horizontally with a spirit level.
Fig. 3. Schematic diagram of the experimental setup. The solid lines are optical paths and the dotted lines are electric signals paths.
The whole detection part of the setup is placed in a black metal box to avoid stray light and prevent potential damage to the detector: a hybrid photomultiplier module (HPM-100-50, Becker&Hickl, Germany) controlled by a detector control card (DCC-100, Becker&Hickl, Germany) in the computer. The collected light is guided by GI-Fiber 2 into a black box to illuminate the detector. A shutter (Melles Griot, Germany), a Labview-controlled motorized attenuator (NDC-100C-4, OD: 0 to 4.0, Thorlabs), a collimator, and a focusing lens are placed between fiber end and detector to adapt the light intensity and reshape the beam.
The detected photon pulses are recorded by a time-correlated single photon counting (TCSPC) module (SPC-150, Becker&Hickl). The Sync signal from laser system is transferred to the TCSPC to provide the timing for the CFD signal from detection system. At every source-detector separation the DTOFs are recorded with a collection time of 1 s and collected 50 times to significantly improve the SNR. The DTOF monitoring, TCSPC settings control, and real-time raw signal analysis are all performed by one LabView program. DTOFs at the various source-detector distances were recorded with the same attenuation filter settings to maintain the amplitude information. The integrals of photon counts over all time channels are used as the light amplitude information in the analysis. The dead-time compensation mode of TCSPC is applied to prevent counting loss due to dead time effect, in case of high photon exposure, particularly for short source-detector separation. The maximum count rate was limited to 800 kHz to avoid dead-time related pulse shape distortions. We should note that this was an experimental set-up for proof-of-principle studies. For in-vivo applications, only single-distance DTOF detection is needed and spatially-resolved amplitude information can be obtained by CW multiple distance detection system.
For TD measurements, the IRF measurements are essential. Both fibers are mounted into a cage system and against each other. A thin diffusor is placed in front of GI-Fiber 2 to fulfill its numerical aperture and realize a similar fiber dispersion as in the phantom measurements [40]. The distance between the two fibers’ ends is precisely measured to determine t0, the time shift between IRF and DTOF measurements. For the system shown here, at 800 nm wavelength, the time resolution of the whole experimental setup, i.e., FWHM of the IRF is about 110 ps. The IRF is always measured twice, before and after the measurements of each phantom.
3.3 Monte-Carlo simulations and forward model
The main quantitative analysis is carried out by a white Monte-Carlo (MC) [41] GPU-based CUDA program. The concept and implementation of MC in diffuse optics have been widely studied through the years. In below we give a brief overview of MC simulations in this work.
DTOF curves are simulated by assuming a homogeneous slab with the thickness of the phantoms, 5 cm. Hereby, the lateral edges of the experimental phantoms can be neglected since the DTOFs were always measured in the central part of the cylinders’ surface. For each simulated set of DTOFs, 429 billion photons are launched into the slab. Simulated photons leaving the slab in reflection are collected on concentric rings around the injection point with a step size of 0.01 cm and a time resolution of 2.44 ps. The refractive index and anisotropy factor are set as n = 1.55 and g = 0.6. Generally, the g value does not affect the simulated DTOFs for source-detector distances between 1.5 and 2.5 cm as in measurements, since under these conditions the simulations solely depend on µs′. The simulated DTOFs are then convolved with a Gaussian laser profile according to the diameter of the laser output fiber. To account for the size of the detector fiber, a circular cross section is overlaid onto the concentric rings at all detector positions as well.
As a collection of these MC simulations, a look-up table (LUT) database was created as the forward model in advance of the measurements. This database contains simulated DTOFs for µs′ range of [0.2 20] cm−1 with a step size of 0.2 cm−1. µa was set to zero according to the white MC concept. DTOFs of specific (µa, µs′) combinations are obtained from this LUT by (i) linear interpolation between the stored DTOFs with respect to µs′, and then (ii) by multiplication of this DTOF with the time dependent absorption factor exp(-µact).
3.4 χ 2 objective function and inverse models
Given the measurements and simulation curves, the universal χ2 objective function is defined as:
For analysis and comparison purposes, there are four inverse models applied in this work: (1) conventional TD model, (2) SD model, (3) global TD model, and (4) spatially-enhanced TD model. For all TD models, i ranges are always selected in the time interval from 80% height on the rising edge (left) to 0.1% height on the decay part (right), relative to the amplitude at the peak of curves. Accordingly, the simulated curves need to convolve with measured IRFs, and then be rescaled to $\mathop \sum \limits_i {m_i}$ for generating si. More specifically:
- (1) Conventional TD model calculates χ2 on DTOFs at one single source-detector distance;
- (2) SD model calculates the χ2 on the SRACs over all space channels ρ, i.e. 1.5 to 2.5 cm;
- (3) Global TD model denotes the fit procedure based on the DTOFs from all available source-detector distances to determine µa and µs′ whereby the relative amplitudes of every individual DTOFs are also considered. By expanding time channel range i to the DTOFs at all space channels ρ, the model exploits the full TD and SD information of the measured data.
- (4) Spatially-enhanced TD model is carried out as described in section 2.3.
4. Results
Figure 4 gives an example of measured DTOFs for the various source-detector distances ρ and corresponding curve fitting results from conventional TD and global TD model. The data in Fig. 4 refer to phantom β2, for its intermediate optical properties among all phantoms.
Fig. 4. (a) DTOFs measured on β2 at source-detector distances ρ: from 15 to 25 mm. (b) Conventional TD model: The optimal curve fit for DTOF at ρ=20 mm, the IRF curve and the weighted residuals. Global TD model: (c) the optimal curve fit, and (d) weighted residuals of DTOFs at various ρ (top: 15 mm, bottom: 25 mm).
In Fig. 4(a), the original DTOFs at various ρ, after background subtraction, contain all the temporal and spatial information from measurements of β2. The curve fit depicted in Fig. 4(b) is the optimal fit of the DTOF at ρ=20 mm based on the conventional TD model, and considers only the temporal information of diffuse photon migration. The flattened residuals along with time imply a good fit quality and a small (for this case the minimal) χ2 value. The curve fit in Fig. 4(c) is based on the global TD model. This approach attempts to simultaneously fit the DTOFs at all ρ and the relative amplitudes amongst DTOFs are not abandoned. Therefore the temporal and spatial information are both considered in some degree. It has to compromise between the global optimal fit and every individual optimal fit when these two contradict each other. As shown in Fig. 4(d), the fit residuals from top to bottom are from the DTOFs at ρ=15 mm to 25 mm, under the optimal global TD regime. The discrepancies at the early time channels of smaller ρ sites will be suppressed if the single distance conventional TD model were applied, but hence leads to different fit results.
Figure 5 summarizes the detailed χ2 distributions in the {µa, µs′} space under the different analysis scenarios for all phantoms. It’s worth to note that the color bars in the 1st, 3rd and 4th column (TD analysis) obey the same limit [0 300], while in 2nd column (SD analysis) it’s 10 times amplified to [0 3000], indicating the steeper χ2 valley in SD.
Fig. 5. χ2 distributions in the {µa, µs′} space. From left to right: conventional TD, SD, Global TD, and Spatially-enhanced TD; from top to bottom: Phantom α2, β2, γ1, γ2 and γ3.
For all phantoms it is obvious that µa and µs′ have opposite correlations in TD and SD profiles leading to different shapes and orientations of low value regions of χ2 distributions. Specifically, the following points can be noted from the different columns in Fig. 5: (1) In conventional TD fit, the ellipse-like low χ2 zones have small gradients along the long axis of the ellipse. Although there is an optimal solution, i.e. minimal χ2 point, the results can easily move away from the real solution or split into several regional minima in the existence of noise or various perturbations. (2) In SD fit, the extended χ2 valley indicates the non-uniqueness of an optimal solution. And the steepness of the valleys shows qualitatively the explicit reciprocal relation amongst µa and µs′, as discussed in Section 2.2. (3) The global TD reduces the ambiguity in low χ2 ellipses by shrinking their area. However, the extent is insufficient and the cross-talk and correlation among the two unknown properties are not eliminated. (4) By spatially-enhanced TD, not only the ambiguity area of low χ2 zone can be reduced, but also the cross-talk artefact can be compressed to minimum extent. The correlation coefficients among the two desired coefficients are nearly reduced to zero, which implies that the entangled contributions from scattering and absorption are decoupled by integrating SD and TD information. (5) For different phantoms, the sensitivities for µa and µs′ differ. In general, a higher µs′/µa ratio will make the inversion more sensitive to µs′, as illustrated by the contour orientations of γ1, γ2 and γ3. On the other hand, a higher µa/µs′ ratio will make the inversion less sensitive to µa, by comparing the contour orientations of α2, β2 and γ2.
Table 2 summarizes the optical properties of the 5 phantoms derived from all 3 TD methods. Suffering from the non-uniqueness, the results by SD method are arbitrary and far from those by TD methods, therefore they are excluded from the comparison. The general trends for all phantoms’ results shows a good match along with the added absorber and scatterer volume shown in Table 1. Table 2 reveals that there are small but systematic changes between the results of the three methods. The spatially-enhanced method yields the largest absorption and reduced scattering coefficients. Compared to the global TD method, absorption is larger by 0.3% to 1.5%, and scattering by 0.1% to 1.9%. The main reason for this deviation is that the center lines of the valleys in the space domain method (Fig. 5, 2nd column) do not cross the minima of the χ2 distributions from the global fits (Fig. 5, 3rd column). Instead, the center lines are slightly shifted to the upper right with respect to the ellipses from the global TD analysis. This result indicates some general deviations between the pure TD and the pure SD information. The global TD analysis exploits both types of information in the most comprehensive way. Even though, the systematic deviations of the corresponding residuals from the zero line (Fig. 4d) indicate that there is no perfect match between the global TD model and the experiment either. Hence, one cannot decide whether the global TD or the spatially-enhanced TD results are the most reliable solution. Both methods weight the measured information in different ways. Another reason for the deviations could be that both methods sample slightly different sub-volumes of the phantom, while the phantoms are not perfectly homogeneous. The absorption and reduced scattering coefficients from the conventional TD method are smaller than the results of the other TD methods. Compared to the global TD analysis, the deviation in absorption ranges from −0.2% to −2.8% and in scattering from −0.2% to -2.5%. Hence, by skipping the amplitude information of the acquired DTOFs we get a shift along the diagonal line (main axis of the ellipse in the 1st column of Fig. 5). Generally, the small deviations between the spatially-enhanced TD and the global TD method demonstrate the power of the spatially-enhanced method. Corresponding data can be acquired by a single time-domain detection channel combined with highly economical and fast CW intensity detection at selected distances ρ. In contrast, the global TD method would require either time-consuming scanning with a single TD channel or parallel detection with an expensive multichannel TD system.
Table 2. The absorption and reduced scattering coefficients derived from conventional, global, and spatially-enhanced TD methods. (units: cm−1)
5. Analysis of uncertainty, robustness, and uniqueness
The quality and performance of the 3 TD methods are quantitatively analyzed and compared from the perspectives of uncertainty, uniqueness, and stability.
5.1 Uncertainty
We propose to use the shape and size of χ2 contours as the quantitative measure to compare the uncertainty of methods. As presented in Fig. 5, the low χ2 ellipse zones from conventional TD can be squeezed to smaller zones. The uncertainty of the same level of χ2, i.e. the extremum on the contours, can be reduced. As a measure of uncertainty, the ratio of error bar on 5 times χmin2 contour lines divided by the optimal solution is defined in Eq. (8). Our analysis has shown that the results have good agreement based on different-times χmin2 contours. The choice of 5 times χmin2 contour here is used as a representative.
where µMax and µMin are the maximum and minimum values of µa or µs′ on the designated contour line. µOptima is the value of µa or µs′ at χmin2 enclosed by contours, i.e., the optimal solution.Table 3 lists the relative uncertainty ɛ of the different methods for all 5 phantoms. The ratio reveals that, under the presence of noise or perturbation, how far the maximal deviations of unknown parameters can reach. By comparison, it is apparent that Spatially-enhanced TD can significantly reduce this kind of variance.
Table 3. The uncertainty comparison about µa and µs′ of different methods.
In addition, the size of low-χ2 region is the other measure to evaluate whether a method can find a well-defined optimum and compress variance. Here we normalize all areas of 5 times χmin2 contours to that of conventional TD, as in Table 4.
Table 4. Comparison the area of 5 times χmin2 contour of different methods
Compared with the conventional TD method, spatially-enhanced TD method halves the area of the low χ2 zones. This indicates that Spatially-enhanced TD method will find more well-defined optima, since the χ2 distributions are more concentrated and the χ2 gradient in low χ2 zones becomes steeper.
5.2 Uniqueness
One of the most important attributes of TD NIRS techniques is their ability to separate effects from absorption and scattering. Although in theory all the TD methods exhibit an optimum, i.e. minimal χ2 in the {µa, µs′} space, the effects of changes in both parameters on the residuals often interfere. Error in one parameter can manifest as artefacts in the other and reduce the accuracy of both. Such correlation leads to ambiguity and cross-talk among parameters, or several regional optimal solutions which make the inversion process difficult. For quantifying the correlation strength, here Pearson correlation coefficient of points on the 5 times χmin2 contour line is used to evaluate the performance of decoupling the unknown parameters.
Table 5. Pearson correlation coefficients of 5×χmin2 contour lines
It can also be noticed from Fig. 5 that the χ2 distributions in the {µa, µs′} space for conventional TD and global TD are both positively correlated, despite global TD slightly reduces the correlation by considering the SD’s negative correlation effect. The correlation coefficients can be reduced to nearly zero by spatially-enhanced TD method for all phantoms. The essential aspect is that this method can adjust the weights between two contradicting correlation and search for an optimal balance to neutralize both. Therefore, for any µa, µs′ combinations it can always achieve a compromise between the effects from µa and µs′ and separate both parameters. The uniqueness of the optimal solution can be guaranteed by the decorrelation process.
5.3 Stability
In this work we also study the stability of the methods between successive measurements. The measurements are all recorded for 50 repetitions with 1 s per repetition. Here we split the measurements into 10 segments where each segment comprises 5 s. Then the three TD methods are applied on these segments to obtain 10 different results of µa and µs′. The mean ($\bar{\mu }$) and standard deviation (σ) are calculated. Then the coefficients of variation ($\textrm{CV} = \sigma /\bar{\mu }$) of the results from different methods are calculated and presented in Table 6.
Table 6. Comparison of coefficient of variation for different methods
The coefficient of variation can describe the degree of dispersion of the results under the same experimental conditions but at different measurement time, thus evaluate the repeatability of methods. According to Table 6, in general the dispersions are all small, owing to the fine MC database and low-noise measurements. Both global TD and spatially-enhanced TD show improved stability compared to conventional TD. For some cases CV(µa) values for spatially-enhanced TD are slightly higher than for global TD. This is mainly due to the fact that global TD processes much more information than spatially-enhanced TD. This also means, by acquiring fewer measurements, spatially-enhanced TD is able to reach the nearly a same level of stability as global TD. The stability performance of spatially-enhanced TD could be further improved, if using a reliable CW power detection system to provide amplitude information, instead of the photon counting system utilized in this experiment setup. It is essential that, in lower signal or short-time acquisition conditions, the methods which consider the SD information will compress the spread of results.
6. Conclusion
In aim to improve uniqueness, stability and uncertainty associated with the inverse problem on simultaneously retrieving the absolute quantities of optical absorption and scattering properties in turbid media, a multivariate method integrating time and space resolved techniques of near-infrared spectroscopy is proposed. The spatially-resolved distributions of times of flight contain the maximal information one can obtain from the reflectance scattered field, but also the sufficient information for extracting complete optical knowledge of an object’s interior. The essence of the problem is how numerous observations of different data types should be combined so as to give the most optimal result. We utilized the opposite correlations in time domain (positive) and space domain (negative) observations with the entangled effects that originated from two physically independent features - absorption and scattering - and dynamically merged them into an artificial spatio-temporal χST2 projection on {µa, µs′} space. By applying the isoperimetric inequality theorem on {µa, µs′} space, absorption and scattering can be separated with high accuracy. Monte-Carlo simulations are used to illustrate and validate the concept of the method. A set of tissue-mimicking phantoms have been measured in a spatially-enhanced time domain experiment. By recovering µa and µs′ of the phantoms with various inverse methods, the improvements on the result quality given in terms of uncertainty, uniqueness and stability are demonstrated through comparison. With fewer measurements and simpler system, the spatially-enhanced TD method could achieve similar level of stability as the global TD method. It is shown that, with the complementarity, Spatially-enhanced TD method can efficiently constrain the time-domain feature of photon migration by its spatially-resolved amplitude patterns. When measuring low photon counts at large source-detector separations, spatially-enhanced TD NIRS has the potential to avoid the influence from superficial tissues and allow the homogenous models to estimate the optical properties of deep regions in heterogeneous layered structures. Although the present work focuses on homogenous media so far, similar principles could be applied in more sophistical structures such as layered media, where scattering properties in deeper region are difficult to determine. Besides, the effective combination of TD and SD information may indicate the potential of integrating other kinds of data types, such as temporal moments, Mellin-Laplace moments, and various frequency-domain quantities.
Funding
H2020 Marie Skłodowska-Curie Actions (Grant No. 675332, BitMap, Innovative Training Networks).
Acknowledgment
The authors acknowledge Stefan Schieck for preparing the phantoms.
References
1. F. F. Jöbsis, “Noninvasive, infrared monitoring of cerebral and myocardial oxygen sufficiency and circulatory parameters,” Science 198(4323), 1264–1267 (1977). [CrossRef]
2. A. M. Smith, M. C. Mancini, and S. Nie, “Bioimaging: second window for in vivo imaging,” Nat. Nanotechnol. 4(11), 710–711 (2009). [CrossRef]
3. F. F. Jöbsis-VanderVliet, “Discovery of the near-infrared window into the body and the early development of near-infrared spectroscopy,” J. Biomed. Opt. 4(4), 392–397 (1999). [CrossRef]
4. G. Bale, C. E. Elwell, and I. Tachtsidis, “From Jöbsis to the present day: a review of clinical near-infrared spectroscopy measurements of cerebral cytochrome-c-oxidase,” J. Biomed. Opt. 21(9), 091307 (2016). [CrossRef]
5. T. Durduran, R. Choe, W. B. Baker, and A. G. Yodh, “Diffuse optics for tissue monitoring and tomography,” Rep. Prog. Phys. 73(7), 076701 (2010). [CrossRef]
6. M. Wolf, M. Ferrari, and V. Quaresima, “Progress of near-infrared spectroscopy and topography for brain and muscle clinical applications,” J. Biomed. Opt. 12(6), 062104 (2007). [CrossRef]
7. Y. Hoshi and Y. Yamada, “Overview of diffuse optical tomography and its clinical applications,” J. Biomed. Opt. 21(9), 091312 (2016). [CrossRef]
8. J. Mourant, J. Freyer, A. Hielscher, A. Eick, D. Shen, and T. Johnson, “Mechanisms of light scattering from biological cells relevant to noninvasive optical-tissue diagnostics,” Appl. Opt. 37(16), 3586–3593 (1998). [CrossRef]
9. A. El-Desoky, D. Delpy, B. Davidson, and A. Seifalian, “Assessment of hepatic ischaemia reperfusion injury by measuring intracellular tissue oxygenation using near infrared spectroscopy,” Liver Int. 21(1), 37–44 (2001). [CrossRef]
10. M. Ferrari and V. Quaresima, “A brief review on the history of human functional near-infrared spectroscopy (fNIRS) development and fields of application,” NeuroImage 63(2), 921–935 (2012). [CrossRef]
11. A. Torricelli, D. Contini, A. D. Mora, A. Pifferi, R. Re, L. Zucchelli, M. Caffini, A. Farina, and L. Spinelli, “Neurophotonics: non-invasive optical techniques for monitoring brain functions,” Funct. Neurol. 29(4), 223 (2014). [CrossRef]
12. M. Calderon-Arnulphi, A. Alaraj, and K. V. Slavin, “Near infrared technology in neuroscience: past, present and future,” Neurol. Res. 31(6), 605–614 (2009). [CrossRef]
13. N. Naseer and K.-S. Hong, “fNIRS-based brain-computer interfaces: a review,” Front. Hum. Neurosci. 9(3), 3 (2015). [CrossRef]
14. D. Grosenick, H. Rinneberg, R. Cubeddu, and P. Taroni, “Review of optical breast imaging and spectroscopy,” J. Biomed. Opt. 21(9), 091311 (2016). [CrossRef]
15. C. Lindner, M. Mora, P. Farzam, M. Squarcia, J. Johansson, U. M. Weigel, I. Halperin, F. A. Hanzu, and T. Durduran, “Diffuse optical characterization of the healthy human thyroid tissue and two pathological case studies,” PLoS One 11(1), e0147851–22 (2016). [CrossRef]
16. W. Weigl, D. Milej, D. Janusek, S. Wojtkiewicz, P. Sawosz, M. Kacprzak, A. Gerega, R. Maniewski, and A. Liebert, “Application of optical methods in the monitoring of traumatic brain injury: A review,” J. Cereb. Blood Flow Metab. 36(11), 1825–1843 (2016). [CrossRef]
17. N. Petitdidier, A. Koenig, H. Grateau, P. Jallon, N. Petitdidier, A. Koenig, R. Gerbelot, H. Grateau, S. Gioux, and P. Jallon, “Contact, high-resolution spatial diffuse reflectance imaging system for skin condition diagnosis,” J. Biomed. Opt. 23(11), 1 (2018). [CrossRef]
18. V. R. Kondepati, H. M. Heise, and J. Backhaus, “Recent applications of near-infrared spectroscopy in cancer diagnosis and therapy,” Anal. Bioanal. Chem. 390(1), 125–139 (2008). [CrossRef]
19. S. Kleiser, N. Nasseri, B. Andresen, G. Greisen, and M. Wolf, “Comparison of tissue oximeters on a liquid phantom with adjustable optical properties,” Biomed. Opt. Express 7(8), 2973–2992 (2016). [CrossRef]
20. F. Scholkmann, S. Kleiser, A. J. Metz, R. Zimmermann, J. M. Pavia, U. Wolf, and M. Wolf, “A review on continuous wave functional near-infrared spectroscopy and imaging instrumentation and methodology,” NeuroImage 85, 6–27 (2014). [CrossRef]
21. M. S. Patterson, B. Chance, and B. C. Wilson, “Time resolved reflectance and transmittance for the noninvasive measurement of tissue optical properties,” Appl. Opt. 28(12), 2331–2336 (1989). [CrossRef]
22. J. B. Fishkin and E. Gratton, “Propagation of photon-density waves in strongly scattering media containing an absorbing semi-infinite plane bounded by a straight edge,” J. Opt. Soc. Am. A 10(1), 127–140 (1993). [CrossRef]
23. D. Borycki, O. Kholiqov, and V. J. Srinivasan, “Reflectance-mode interferometric near-infrared spectroscopy quantifies brain absorption, scattering, and blood flow index in vivo,” Opt. Lett. 42(3), 591–594 (2017). [CrossRef]
24. D. Borycki, O. Kholiqov, and V. J. Srinivasan, “Correlation gating quantifies the optical properties of dynamic media in transmission,” Opt. Lett. 43(23), 5881–5884 (2018). [CrossRef]
25. S. R. Arridge and W. R. Lionheart, “Nonuniqueness in diffusion-based optical tomography,” Opt. Lett. 23(11), 882–884 (1998). [CrossRef]
26. J. Selb, T. M. Ogden, J. Dubb, Q. Fang, and D. A. Boas, “Comparison of a layered slab and an atlas head model for Monte Carlo fitting of time-domain near-infrared spectroscopy data of the adult head,” J. Biomed. Opt. 19(1), 016010 (2014). [CrossRef]
27. A. Pifferi, D. Contini, A. D. Mora, A. Farina, L. Spinelli, and A. Torricelli, “New frontiers in time-domain diffuse optics, a review,” J. Biomed. Opt. 21(9), 091310 (2016). [CrossRef]
28. D. Milej, A. Abdalmalak, D. Janusek, M. Diop, A. Liebert, and K. St. Lawrence, “Time-resolved subtraction method for measuring optical properties of turbid media,” Appl. Opt. 55(7), 1507–1513 (2016). [CrossRef]
29. D. Milej, A. Abdalmalak, P. McLachlan, M. Diop, A. Liebert, and K. St. Lawrence, “Subtraction-based approach for enhancing the depth sensitivity of time-resolved NIRS,” Biomed. Opt. Express 7(11), 4514–4526 (2016). [CrossRef]
30. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15(2), R41–R93 (1999). [CrossRef]
31. M. Schweiger and S. R. Arridge, “Application of temporal filters to time resolved data in optical tomography,” Phys. Med. Biol. 44(7), 1699–1717 (1999). [CrossRef]
32. A. E. Cerussi, R. Warren, B. Hill, D. Roblyer, A. Leproux, A. F. Durkin, T. D. O’Sullivan, S. Keene, H. Haghany, T. Quang, W. M. Mantulin, and B. J. Tromberg, “Tissue phantoms in multicenter clinical trials for diffuse optical technologies,” Biomed. Opt. Express 3(5), 966 (2012). [CrossRef]
33. A. Farina, A. Torricelli, I. Bargigia, L. Spinelli, R. Cubeddu, F. Foschum, M. Jäger, E. Simon, O. Fugger, A. Kienle, F. Martelli, P. Di Ninni, G. Zaccanti, D. Milej, P. Sawosz, M. Kacprzak, A. Liebert, and A. Pifferi, “In-vivo multilaboratory investigation of the optical properties of the human head,” Biomed. Opt. Express 6(7), 2609–2623 (2015). [CrossRef]
34. A. Torricelli, D. Contini, A. Pifferi, M. Caffini, R. Re, L. Zucchelli, and L. Spinelli, “Time domain functional NIRS imaging for human brain mapping,” NeuroImage 85, 28–50 (2014). [CrossRef]
35. H. Wabnitz and H. Rinneberg, “Imaging in turbid media by photon density waves: spatial resolution and scaling relations,” Appl. Opt. 36(1), 64–74 (1997). [CrossRef]
36. A. Liebert, H. Wabnitz, D. Grosenick, M. Möller, R. Macdonald, and H. Rinneberg, “Evaluation of optical properties of highly scattering media by moments of distributions of times of flight of photons,” Appl. Opt. 42(28), 5785–5792 (2003). [CrossRef]
37. S. Suzuki, S. Takasaki, T. Ozaki, and Y. Kobayashi, “Tissue oxygenation monitor using NIR spatially resolved spectroscopy,” Proc. SPIE 3597, 582–592 (1999). [CrossRef]
38. J. Swartling, J. S. Dam, and S. Andersson-Engels, “Comparison of spatially and temporally resolved diffuse-reflectance measurement systems for determination of biomedical optical properties,” Appl. Opt. 42(22), 4612–4620 (2003). [CrossRef]
39. L. Spinelli, M. Botwicz, N. Zolek, M. Kacprzak, D. Milej, P. Sawosz, A. Liebert, U. Weigel, T. Durduran, F. Foschum, A. Kienle, F. Baribeau, S. Leclair, J.-P. Bouchard, I. Noiseux, P. Gallant, O. Mermut, A. Farina, A. Pifferi, A. Torricelli, R. Cubeddu, H.-C. Ho, M. Mazurenka, H. Wabnitz, K. Klauenberg, O. Bodnar, C. Elster, M. Bénazech-Lavoué, Y. Bérubé-Lauzière, F. Lesage, D. Khoptyar, A. A. Subash, S. Andersson-Engels, P. Di Ninni, F. Martelli, and G. Zaccanti, “Determination of reference values for optical properties of liquid phantoms based on Intralipid and India ink,” Biomed. Opt. Express 5(7), 2037–2053 (2014). [CrossRef]
40. A. Liebert, H. Wabnitz, D. Grosenick, and R. Macdonald, “Fiber dispersion in time domain measurements compromising the accuracy of determination of optical properties of strongly scattering media,” J. Biomed. Opt. 8(3), 512–517 (2003). [CrossRef]
41. E. Alerstam, S. Andersson-Engels, and T. Svensson, “White Monte Carlo for time-resolved photon migration,” J. Biomed. Opt. 13(4), 041304 (2008). [CrossRef]
