1.

Introduction

A major challenge in tissue optics is the quantification of the optical properties of the different types of biological tissues, their spatial inhomogeneity, and their evolution in time. This task is inherent to different applications such as functional brain imaging,^1,2 cancer diagnostics,^3–5 therapy monitoring,^6,7 and, in general, tissue spectroscopy.^8–18 Biological tissue is a complex random medium that, in the red and in the near-infrared (NIR) therapeutic window ($600\mathrm{nm}<\lambda <1000\mathrm{nm}$) and for dimensions in the centimeter range, typically acts as a diffusive medium. It can be optically characterized by an absorption coefficient, ${\mu }_a$, carrying information about tissue chromophores, such as oxy- and deoxy-hemoglobin, by a reduced scattering coefficient, ${\mu }_s^{\prime }$ and by a refractive index $n$.

In the simplest case, the medium can be modeled as being homogeneous. The determination of ${\mu }_a$ and ${\mu }_s^{\prime }$ of homogeneous turbid media has been approached through the years with several techniques that can be divided into three main categories:¹⁹ continuous wave (cw) techniques,^20,21 frequency-domain techniques,¹⁵ and time-domain techniques.^11,13,22,23 Most of these methods use indirect procedures to obtain the optical properties in which ${\mu }_a$ and ${\mu }_s^{\prime }$ are retrieved by fitting an appropriate model to the measured data. The three acquisition techniques have an intrinsically different level of information content so that the accuracy of the retrieved optical properties differs. For instance, the use of cw data with a measurement at single source-detector separation does not permit uncoupling of ${\mu }_a$ from ${\mu }_s^{\prime }$ and thus the ability to determine them separately. In contrast, from time-domain measurements for homogeneous media, ${\mu }_a$ and ${\mu }_s^{\prime }$ can be obtained from a single-distance measurement by analyzing the shape of the complete distribution of time of flight (DTOF) of photons. Theoretically, frequency-domain and time-domain approaches are expected to be equivalent. However, in their practical implementation, mainly due to the limited number of modulation frequencies, the information content of frequency-domain measurements is usually lower than that of time-domain approaches.

The different relative information content of cw, frequency, and time-domain approaches is also relevant when investigating more complex, heterogeneous media. In fact, biological tissues are complex media with heterogeneous structures. A number of investigations utilized the homogeneous model as an approximation to obtain an estimate of ${\mu }_a$ of heterogeneous tissues.^14,16,24,25 However, the neglected presence of the superficial layers (e.g., fat, scalp, and skull) may alter the retrieved ${\mu }_a$ values of the deep target tissues (e.g., muscle and brain). In other works, trying to take into account the actual structure of tissues, the optical properties were obtained in the time domain by using layered models.^10,26 Layered model geometries were employed, motivated by the fact that some tissues have a more or less layered architecture. In particular, this is the case for muscle underneath a superficial fat layer¹⁰ or for the head with compartments like scalp, skull, and brain.²⁶ Thus, the introduction of layered models in the inversion procedures allowed optical properties of the deep layers to be distinguished from those of the superficial layers. This point has been addressed in a number of papers dedicated to tissue spectroscopy^26–32 with investigations on simulated and experimental data. The layered models refer to two-layered^26–29 or multilayered^30–32 geometries with plane interfaces. Also, real anatomical structures were considered.²⁹ The experimental data of these investigations were obtained from the following: (1) layered solutions of the diffusion equation (DE)^{27,28,30–32} or Monte Carlo (MC) simulations based on real anatomical structures²⁹; (2) phantom measurements;^26,28,31,32 and (3) in vivo measurements on head and muscle tissues.^10,26 It should be noted that phantom measurements are a necessary test before the applications to in vivo data. They allow one to study the influence of the instrumental factors and uncertainties while the geometry and the optical properties are well defined and known.

In the following, we discuss the most relevant approaches in more detail. In a pioneering work, Kienle et al.²⁸ showed the reconstruction of ${\mu }_a$ and ${\mu }_s^{\prime }$ of both layers of a two-layered medium if the thickness of the first layer is known. In the inversion procedure, they exploited a solution of the DE for this geometry. To validate their approach, they performed time-resolved (TR) reflectance measurements at two source-detector separations on two-layered solid tissue phantoms. Sato et al.³¹ developed a different time-domain technique to determine the absorption coefficients of a two-layered medium, based on a single-distance measurement. This method does not employ traditional nonlinear least-square fitting procedures. It rather acts as a differential method, requiring two sets of measurements, one on the object under investigation and the other on a reference object. The same authors recently generalized their approach to a four-layered human head model.³² In both these works, the methods were validated making use of TR measurements on layered epoxy-resin-based phantoms. The advantage of this method is that it does not imply the complexity of inversion algorithms. The disadvantages are that it requires two sets of measurements (object and reference data) and that the approach does not enable a full reconstruction of all the parameters of the layered model. Shimada et al.³⁰ exploited an inversion procedure that is based on scaling properties of the DE. It combines TR reflectance measurements at short and long source-detector separations to retrieve the absorption coefficients of a two-layered diffusive medium. Again, TR reflectance measurements at different source-detector separations are needed. Finally, we mention the fitting routine for multidistance TR data that were implemented by Selb et al.²⁹ and validated on data generated by MC simulations on segmented anatomical images. This approach has two novelties: the application of an atlas head model with two tissue types (intra- and extra-cerebral) as model geometry, in comparison with a simple two-layered model, and the use of a MC-based fitting procedure (exploiting a graphics processing unit).

It is worth noting that in some cases, such as functional NIR spectroscopy, it can be sufficient to determine absorption changes ($\mathrm{\Delta }{\mu }_a$) occurring in different tissue compartments instead of the absolute values of ${\mu }_a$. This information can be obtained more easily than ${\mu }_a$, in particular, when small values of $\mathrm{\Delta }{\mu }_a$ are addressed and a linear approximation is applicable. Still, even in this case, accuracy and cross-talk between tissue absorption changes in the various compartments are not trivial issues. Limiting to the time-domain approach, different strategies have been developed. The complete DTOF is exploited by reconstruction methods utilizing TR mean partial path lengths in the compartments.^33,34 Other methods are based on moments^35–37 or the Mellin-Laplace transform^38,39 of the measured DTOFs. Also, simpler heuristic approaches based on the selection of specific time windows of DTOFs have been developed.^40,41

In our previous works, dedicated to the two-layer medium^26,27 and to TR measurements, we performed a reconstruction of the optical properties exploiting multidistance TR measurements on a two-layer phantom. The retrieval was obtained with a nonlinear fitting procedure employing the Levenberg–Marquardt (LM) algorithm⁴² based on a closed-form solution of the DE for two layers. We showed that it is actually possible to retrieve ${\mu }_a$ of both layers, ${\mu }_s^{\prime }$ of the first layer and, less accurately, also the first layer thickness. With the aim to further developing this successful approach in the time domain, we have introduced these substantial improvements: (1) we have replaced the LM algorithm with an optimal estimation (OE) algorithm^43,44 that can also account for a priori information on the range of variability of the optical and geometrical characteristics of the medium, and (2) we have replaced the multidistance measurements with a single-distance measurement. The purpose of this strategy was twofold. With the OE, we aimed at improving the performances of the inversion procedure exploiting a priori information while using single-distance measurements, we wanted to investigate the information content encoded in such a basic dataset. Moreover, since in reality, the tissues layers are not exactly homogeneous, the advantage of single-distance over multidistance in-vivo measurements is that the latter can be more compromised by the heterogeneity of the first tissue layer (for instance, nonuniform distribution of skin pigmentation and blood vessels). However, we also stress that when dealing with a single-distance measurement, the simultaneous reconstruction of all optical properties of a layered medium is inherently challenging since the inverse problem is ill posed. Therefore, this study addresses the limiting performances achievable with a single TR measurement and a two-layered geometry. The choice of the two-layered model is a reasonable compromise between robustness and fineness in the description of the medium, while the use of a multilayered model could result in a low reliability and stability of the inversion procedure. Finally, we note that for applications related to the oxygen saturation of brain and muscle tissues, the only relevant parameter is often the ${\mu }_a$ of the deeper layer. With this work, we show that with the OE, the ${\mu }_a$ of the lower layer can be retrieved with a high level of reliability and accuracy.

In the following, we first describe (Sec. 2) the algorithm for the reconstruction, the two-layer tissue phantom used, the experimental setup employed, and the kinds of measurements and retrievals carried out. The results of the various kinds of retrievals carried out with the OE and LM methods are presented and compared in Sec. 3. In Sec. 4, the results are discussed and conclusions are provided in Sec. 5.

2.

Methods

2.2.

Layered Liquid Phantom

The layered liquid phantom was implemented by means of a modular scattering cell made of black polyvinyl chloride (PVC) with small plexiglass windows that was previously developed.^46–48 The inner dimensions of the cell were 120 mm width, 145 mm height, and 70 mm thickness. The 2-mm-thick front wall contained transparent circular windows of 7 mm diameter. The two-layered structure was implemented with the help of a Mylar foil of $30\mu \mathrm{m}$ thickness that was inserted with a spacer of 10 mm thickness to separate the two compartments of the cell with minimal perturbation of light propagation.⁴⁶ Difficulties in flattening the inelastic Mylar foil led to variations in the thickness of the first layer. The thickness was estimated to be 9 mm in the particular experiment reported here, see Ref. 48. The effects of the finite lateral dimensions of the cell were studied in detail by MC simulations.⁴⁷ Under the given conditions, these effects were negligible so that the layers can be regarded as infinitely extended. Moreover, the second layer of about 60 mm thickness can be assumed to be infinitely thick.

The scattering and absorbing liquid was prepared from Intralipid-20%, India ink, and water. The scattering and absorption properties of Intralipid and ink had been accurately characterized in a multicenter study by Spinelli et al.⁴⁹ using different methods.^50–57 The amounts of components to be mixed are solely determined by accurate weighing, which ensures that the concentrations are known with high accuracy.

Based on the known scattering and absorption coefficients of the components, a suspension S0 with the initial optical properties ${\mu }_s^{\prime }=1{\mathrm{mm}}^{-1}$ and ${\mu }_a=0.01{\mathrm{mm}}^{-1}$ was prepared. In the first part of the phantom experiment, the ${\mu }_a$ of the upper layer was increased from $0.01{\mathrm{mm}}^{-1}$ to about $0.02{\mathrm{mm}}^{-1}$ in 11 steps, while the optical properties of the lower layer remained unchanged. Then, the upper layer was emptied and refilled with the suspension S0 before the absorption in the lower layer was increased in the same steps approximately. Care was taken to add absorber without inducing a change in scattering. In both parts of the experiment, the first two steps of absorption changes were about $0.0005{\mathrm{mm}}^{-1}$ and the remaining steps were about $0.001{\mathrm{mm}}^{-1}$. The exact amount of liquid added in each step was determined by a precision balance so that the nominal values of ${\mu }_a$ in each step were accurately quantified. The procedure of phantom preparation was described in detail in Ref. 37. From now onward, we denote the absorption coefficients of the first and second layers by ${\mu }_{a1}$ and ${\mu }_{a2}$, respectively, and the reduced scattering coefficients by ${\mu }_{s1}^{\prime }$ and ${\mu }_{s2}^{\prime }$. The thicknesses of the first and second layers were $s_1=9\mathrm{mm}$ and $s_2=60\mathrm{mm}$. The refractive index of the medium was assumed equal to that of water, $n=1.33$, while the refractive index of the external medium (PVC wall of the scattering cell) was $n=1.54$. For all measurements, we fixed ${\mu }_{s1}^{\prime }={\mu }_{s2}^{\prime }=1{\mathrm{mm}}^{-1}$, while ${\mu }_{a1}$ and ${\mu }_{a2}$ were changed as described earlier. In Fig. 1, their nominal values are displayed versus the measurement number. The uncertainty of the nominal values can be estimated with 3% for the optical properties and roughly 10% for $s_1$.

Fig. 1

Nominal values of the absorption coefficients, (a) ${\mu }_{a1}$ and (b) ${\mu }_{a2}$, inside the two layers of the phantom for a set of 24 independent measurements.

3.

Retrievals with Optimal Estimation and Levenberg–Marquardt

To test the OE method versus the LM method, we have implemented six types of retrievals (cases 1 through 6) on the set of 24 measurements identified by the nominal values illustrated in Fig. 1. For the retrievals carried out here, the relevant parameters to be considered are ${\mu }_{a1}$, ${\mu }_{a2}$, ${\mu }_{s1}^{\prime }$, ${\mu }_{s2}^{\prime }$, $s_1$, and the origin of time, $t_0$. Actually, the reconstruction of the optical properties from TR measurements may be strongly influenced by the value estimated for $t_0$. The fit parameters of the retrievals were chosen among these parameters. The characteristics of each retrieval in terms of the fit parameters and the fixed forward model parameters can be found in Table 1. The table specifies the a priori information used in the OE retrievals for each fit parameter or fixed forward model parameter providing its expectation a priori value together with the related standard deviation. For the LM retrievals, we considered the same initial values of the fit parameters used for the OE retrievals, while the fixed forward model parameters were assumed to be known and equal to their expectation values.

Table 1

Kind of retrievals carried out and priors used in the retrievals. For each “Fit” parameter and “Fixed” forward model parameter, the expectation value and the standard deviation σ of the a priori information used in the retrieval are shown.

The aim of this investigation was to find out whether the use of a priori information can really improve stability and accuracy of the inversion procedures when dealing with real experimental data, while the ranges in which the parameters can vary are still reasonably large. The sample of six kinds of retrievals selected was aimed at testing the effect of a priori information on the reconstruction of ${\mu }_{a1}$, ${\mu }_{a2}$, ${\mu }_{s1}^{\prime }$, and ${\mu }_{s2}^{\prime }$ that are in most cases the physically relevant unknowns of the inversion procedure. The number of fit parameters used in each retrieval was varied from two up to six while trying to figure out possible scenarios of reconstructions carried out in real applications. The curve fit was usually performed selecting a time range in percent of the curve peak starting at 85% (0.85 multiplied by the maximum intensity) on the rising edge up to 0.01% (0.0001 multiplied by the maximum intensity) in the tail of the measured DTOFs. The SNR of the measured DTOFs was sufficient to extend the curve fitting up to ${10}^{-4}$ of the maximum. It should be noted that the inclusion of photons with long flight times is important for the retrieval of the parameters of the lower layer. The measured and theoretical DTOFs were both normalized to the whole area subtended by the curves within the fit limits. This normalization of experimental and theoretical DTOFs superseded the use of an amplitude parameter as extra fit parameter. In this work, we have also considered the origin of time $t_0$ as fit parameter or fixed forward model parameter since $t_0$ is important and critical for TR measurements and analysis while its determination is usually subjected to several error sources.

We have performed OE and LM retrievals. For the OE retrieval, the a priori information on the fit parameters is given by their initial values (estimated expected value) and their range of variation (estimated standard deviation) as shown in Table 1. In case 1, we have addressed a four fit parameters retrieval aimed at reconstructing ${\mu }_{a1}$, ${\mu }_{a2}$, ${\mu }_{s1}^{\prime }$, and ${\mu }_{s2}^{\prime }$. The thickness $s_1$ and the origin of time $t_0$ were taken as fixed forward model parameters exactly known and equal to their nominal values determined from the phantom and the experimental setup. The a priori information of the fit parameters ($3\sigma$) roughly spans a range of $\pm 30\%$ to 40% around the initial values (a priori estimate) of the parameters. In Fig. 2, the retrieved values and the nominal values of the fit parameters are plotted versus the measurement number. The fixed forward model parameters $s_1$ and $t_0$ are also illustrated. The results show an excellent reconstruction of ${\mu }_{a2}$, but less accurate reconstructions for the other parameters. In particular, the worst reconstruction was obtained for ${\mu }_{s2}^{\prime }$ and ${\mu }_{a1}$, which is in agreement with the results reported in previous papers.^26,27 In particular, the LM values exhibited a remarkably larger scatter around the nominal values.

Fig. 2

Case 1: retrieval with four fit parameters, i.e., ${\mu }_{a1}$, ${\mu }_{a2}$, ${\mu }_{s1}^{\prime }$, and ${\mu }_{s2}^{\prime }$. The retrieved values are given by the markers, the nominal values by the continuous line.

In case 2, $s_1$ was included as additional fit parameter, apart from ${\mu }_{a1}$, ${\mu }_{a2}$, ${\mu }_{s1}^{\prime }$, and ${\mu }_{s2}^{\prime }$. The a priori information was roughly spanning a range of $\pm 30\%$ to 40% around the initial value (a priori estimate) of the parameters. For this case, the convergence achieved with the LM method was quite poor and instable with an evident sensitiveness to the initial values of the fit parameters. For some measurements, the convergence was not even obtained; and for this reason, the LM results are not shown in the figure. In Fig. 3, the values retrieved with the OE, their error bars (three times the standard deviation of the retrieved parameters), and the nominal values of the fit parameters and the fixed forward model parameters are plotted versus the measurement number. The results show a quality of the OE retrieval that is even better than the previous case with four fit parameters. It is important to note the stability and the accuracy of the retrieval of the absorption coefficient of the second layer ${\mu }_{a2}$ since for applications such as brain oxygenation measurements, this parameter is the most relevant one. It is also interesting to note the good stability obtained for the reconstruction of $s_1$.

Fig. 3

Case 2: retrieval with five fit parameters, i.e., ${\mu }_{a1}$, ${\mu }_{a2}$, ${\mu }_{s1}^{\prime }$, ${\mu }_{s2}^{\prime }$, and $s_1$. The retrieved values are given by the markers, the nominal values by the continuous line.

Case 3 pertains to a retrieval with six fit parameters, i.e., ${\mu }_{a1}$, ${\mu }_{a2}$, ${\mu }_{s1}^{\prime }$, ${\mu }_{s2}^{\prime }$, $s_1$, and $t_0$. The a priori information was roughly spanning a range of $\pm 30\%$ to 40% around the initial value (a priori estimate) for ${\mu }_{a1}$, ${\mu }_{a2}$, ${\mu }_{s1}^{\prime }$, ${\mu }_{s2}^{\prime }$, and $s_1$. For $t_0$, we selected a standard deviation of 10 ps around the initial a priori estimate. As for the previous case, we do not report the LM retrievals due to the poor convergence obtained. In Fig. 4, the values retrieved by the OE, their error bars, and the nominal values of the fit parameters are plotted versus the measurement number. The results show a quality of the OE retrieval similar to the case of five fit parameters. We note the stability of the reconstruction for the origin of time $t_0$ that shows a spread of its values within 20 ps for all the measurements. A good stability of the reconstruction is also observed for $s_1$. Also in this case, we note the rather good robustness and accuracy of the reconstructed values of the absorption of the second layer ${\mu }_{a2}$.

Fig. 4

Case 3: retrieval with six fit parameters, i.e., ${\mu }_{a1}$, ${\mu }_{a2}$, ${\mu }_{s1}^{\prime }$, ${\mu }_{s2}^{\prime }$, $s_1$, and $t_0$. The retrieved values are given by the markers, the nominal values by the continuous line.

Case 4 pertains to the same six fit parameters as in case 3, but some of the a priori information used for case 3 was released. We used a range of about $\pm 60\%$ to 80% ($\pm 3\sigma$) around the initial estimate of ${\mu }_{a1}$, ${\mu }_{a2}$, ${\mu }_{s1}^{\prime }$, ${\mu }_{s2}^{\prime }$, and $s_1$ (for details see Table 1). The initial values were $0.8{\mathrm{mm}}^{-1}$ for ${\mu }_{s1}^{\prime }$ and ${\mu }_{s2}^{\prime }$ and 11 mm for $s_1$. For $t_0$, the same a priori information as in case 3 was applied. As for the previous cases, the results of the LM retrievals are not shown due to the poor convergence obtained. In Fig. 5, the retrieved values with the OE, their error bars, and the nominal values of the fit parameters are plotted versus the measurement number. The figure reveals a larger spread of the retrieved values around the nominal values compared to case 3, and larger error bars. This change can be ascribed to the less strict a priori information. For ${\mu }_{a2}$, although larger error bars are observed, the retrieved values are still very close to the nominal values. We note that, for all the cases 1 to 4 presented here, the retrieved values obtained with the OE are distributed within the whole range of three standard deviations centered at the initial values of the fit parameters (see Table 1 and the above figures). This confirms the consistency and reliability of the a priori information used. Even if the nominal values are close to the limits of the $3\sigma$ intervals around the expected values, the fit does not seem to be misled by the a priori information.

Fig. 5

Case 4: same fit parameters as in case 3, but with a released a priori information according to Table 1. The retrieved values are given by the markers, the nominal values by the continuous line.

In case 5, we have addressed a two-parameter retrieval where the fit parameters were ${\mu }_{a1}$ and ${\mu }_{a2}$ only. The a priori information used for these parameters is the same as for the cases 1 to 3. Further, in this retrieval, we have simulated a systematic error on the origin of time $t_0$ equal to $+40\mathrm{ps}$. An error of this magnitude could easily occur in real experimental conditions. The origin of time $t_0$ was therefore considered in the OE as a fixed parameter of the forward model with an expectation value equal to 40 ps, and with a standard deviation ${\sigma }_{t_0}=20\mathrm{ps}$. The aim of this retrieval was to test the OE and LM algorithms in the presence of a systematic deviation as well as of an uncertainty of a forward model parameter. An uncertainty of a forward model parameter can only be taken into account with the OE method, while with the LM method any forward model parameter must be assumed to be exactly known. The other fixed forward model parameters ${\mu }_{s1}^{\prime }$, ${\mu }_{s2}^{\prime }$, and $s_1$ were set to their nominal values. In case 6, we have repeated the analysis of case 5 for a $t_0$ equal to $-40\mathrm{ps}$.

The results in Figs. 6 and 7 show that the OE retrieval is able to correctly reconstruct ${\mu }_{a2}$ even in the presence of systematic errors, while the LM retrieval leads to less accurate results. For what concerns ${\mu }_{a1}$, both algorithms have some difficulties in the retrieval that is largely due to the systematic error on $t_0$. Additional retrievals with fit parameters ${\mu }_{a1}$ and ${\mu }_{a2}$ were carried out introducing systematic errors on the other forward model parameters ${\mu }_{s1}^{\prime }$, ${\mu }_{s2}^{\prime }$, and $s_1$ (data not shown). The results obtained were essentially similar to those of Figs. 6 and 7 obtained with a systematic error on $t_0$. The absorption coefficient of the second layer ${\mu }_{a2}$ was always retrieved with the lowest uncertainty and remained rather close to its nominal values. Furthermore, we always observed a better estimate obtained with the OE compared to the LM.

Fig. 6

Case 5: retrieval with two fit parameters, ${\mu }_{a1}$ and ${\mu }_{a2}$. A systematic error on the origin of time $t_0$ equal to $+40\mathrm{ps}$ was assumed. The retrieved values are given by the markers, the nominal values by the continuous line.

Fig. 7

Case 6: retrieval with two fit parameters, ${\mu }_{a1}$ and ${\mu }_{a2}$. A systematic error on the origin of time $t_0$ equal to $-40\mathrm{ps}$ was assumed. The retrieved values are given by the markers, the nominal values by the continuous line.

The results shown in the preceding figures are a sample of some thousands of retrievals done on the same set of 24 measurements. In addition to the retrievals presented here, we have carried out other sets of retrievals where we reduced the number of decades of the trailing edge of the DTOF signal taken into account in the inversion procedure. When reducing the fit range to include the data down to 1% of the curve maximum only, we did not observe significant variations in the retrieved values. This fact indicates the robustness of the results obtained that are only slightly depending on the temporal range considered in the retrieval as long as at least two decades are covered. This statement is also valid for the reconstruction of ${\mu }_{a2}$, which was expected to be particularly sensitive to the coverage of long times of flight.

The results obtained would be meaningless if they were not supported by a high quality of the fit between model and experimental data. To this end, we investigated the reduced chi-squared, ${\widehat{\chi }}^2={[\mathbf{y}-\mathbf{F}(\widehat{\mathbf{x}},\widehat{\mathbf{b}})]}^T{\mathbf{S}}_T^{-1}[\mathbf{y}-\mathbf{F}(\widehat{\mathbf{x}},\widehat{\mathbf{b}})]/\nu$ ($\nu$ degrees of freedom; for $\mathbf{y}$, $\mathbf{F}$, and ${\mathbf{S}}_T$, see Sec. 2.1) as a test of the goodness of the fit. Figure 8 depicts ${\widehat{\chi }}^2$ of each OE retrieval versus the measurement number, for case 3 as a representative example. All values were less than 1.7, indicating an excellent quality of fit. More detailed information on the quality of the fit can be obtained by viewing the weighted residuals of each retrieval calculated as $[y(i)-F(i)]/{\sigma }_y(i)$. Figure 9 shows the comparison between measurement and model for the first measurement of case 3 together with the weighted residuals. The comparison highlights the very good agreement, within the experimental uncertainties, of the measurement and model. Similar residual curves have been obtained for all the retrievals of the cases 1 to 6. We have also observed slightly lower values of the ${\widehat{\chi }}^2$ (according to the aforementioned definition) when the OE was used compared to the retrievals with the LM.

Fig. 8

Behavior of ${\widehat{\chi }}^2$ versus the measurement number for the optimal estimation (OE) retrievals of case 3.

Fig. 9

(a) Comparison of time-resolved reflectance from experiment and model. (b) Residuals weighted with the standard deviation of the measured data, pertaining to case 3 and measurement 1.

4.

Discussion

In general, the accuracy of the retrieved optical properties depends on the following factors: (1) the chosen data acquisition technique, (2) the quality of the measured data, (3) the retrieval procedure used in the reconstruction, and (4) the accuracy of the model employed in the retrieval procedure. In this investigation, we have mainly acted on the factors (1) and (3), switching from multidistance to single-distance measurements and from the LM to the OE estimator, in comparison to our previous work.²⁶

The main achievement of the retrievals carried out with the OE method is the extreme reliability and stability of the retrieved absorption coefficient of the second layer ${\mu }_{a2}$ that shows values always consistent and close to the nominal values. This finding is of great importance for many applications that address the absorption in a deep tissue layer. For a better comprehension of this result, we can exploit the information of the mean path length traveled by detected photons inside the two layers of the medium under investigation. Using the DE, we can calculate the TR mean path lengths, $⟨l_1⟩(t)$ and $⟨l_2⟩(t)$, in the first and second layers, respectively. In Fig. 10, these quantities are calculated for measurements 1, 12, and 24 using the nominal values for the geometrical and optical properties of the two-layer phantom. We stress that in these experimental conditions, the diffusion approximation is valid and the mean path lengths calculated with the DE are accurate as shown elsewhere.⁶² Figure 10 illustrates that, except at early times ($t<2000\mathrm{ps}$) and for the measurements with high absorption in the second layer, $⟨l_2⟩(t)$ is always significantly larger than $⟨l_1⟩(t)$. The evidence of this fact seems to suggest that the main reason of the reliable values retrieved for ${\mu }_{a2}$ is related to the larger mean path length inside the second or deep layer. A further justification of this hypothesis can be found in the increase of the error bars of the retrieved ${\mu }_{a2}$ for measurements beyond measurement 17 for which we have lower values of $⟨l_2⟩(t)$ compared to $⟨l_1⟩(t)$. Thus, there is an evident link between the information content of a measurement related to a certain part of the medium and the mean path length spent by the detected photons inside this region of medium. Consequently, we can argue that large values of $⟨l_2⟩(t)$ correspond to a robust reconstruction of ${\mu }_{a2}$ and vice versa. Time-resolved measurements with high dynamic range that can provide data with good SNR at late times facilitate the exploitation of this relationship.

Fig. 10

Time-resolved mean path lengths, $⟨l_1⟩(t)$ and $⟨l_2⟩(t)$, spent in the two layers calculated with the diffusion equation (DE) for measurements 1, 12, and 24 estimated using the nominal values for the geometrical and optical properties of the two-layered medium.

The results presented demonstrate a substantial advantage of the OE method versus the LM method. A single-distance TR measurement is sufficient for the reconstruction of the optical properties of a two-layered medium when the retrieval is carried out with the OE and with an appropriate level of a priori information. In particular, the reliability and stability of the retrieved ${\mu }_{a2}$ of the second layer is excellent. The maximum deviation between retrieved and nominal values of ${\mu }_{a2}$ for cases 1 to 4 remains within about 10%. Previously published results showed a comparable performance obtained with the LM method; however, in those cases, multidistance measurements were employed.^26,27 Indeed, the results of cases 1 and 2, i.e., Figs. 2 and 3, are better but still comparable with those obtained with the LM method and multidistance measurements in Ref. 26 (see Fig. 1). Measurements at multiple distances between the source and detector naturally provide more information to retrieve multiple parameters. However, they are not always practicable. They require several detection channels. Moreover, in in vivo measurements it cannot be assumed that the geometry (e.g., thickness of extra-cerebral tissue) is the same for all source-detector pairs. Furthermore, the optical properties can vary locally, e.g., due to blood vessels. Therefore, single-distance measurements are preferable. The comparison of the results of this study with the previous one suggests that the application of the OE method allows for the successful analysis of single-distance measurements.

The retrievals of cases 3 and 4 obtained with five and six fit parameters suggest some interesting comments on the fit parameters $s_1$ and $t_0$. It is amazing that the retrieved values of $s_1$ and $t_0$ in these cases are close to their nominal values for all the measurements showing a spread well within the error bars of the retrievals. This fact opens the possibility to really treat these parameters as unknowns in a retrieval implemented with in vivo data. Indeed, when dealing with measurements on the adult head and assuming a two-layered tissue model, we have two possible strategies for the retrieval using the OE method: (1) treating $s_1$ and $t_0$ as fit parameters with a reasonable range of variability or (2) including them as fixed forward model parameters with a suitable standard deviation. The results of Figs. 3–5 suggest that the first option seems feasible without deteriorating the correct retrieval of the optical properties. In particular, the retrieved ${\mu }_{a2}$ is rather independent of the way in which $s_1$ and $t_0$ are processed by the inversion procedure. We note that for any inversion procedure implemented with TR data $t_0$ is a key factor for the success of the retrieval. In the history of TR spectroscopy, $t_0$ has not always been considered as fit parameter,⁶³ but as a forward model parameter exactly known¹⁶ (e.g., it can be estimated from a subset of measurements and then be kept fixed for the further analysis). The main reason of this choice was to avoid the crosstalk with other fit parameters of the retrieval, and in particular the crosstalk with ${\mu }_{s1}^{\prime }$ and ${\mu }_{s2}^{\prime }$. The results of Figs. 4 and 5 suggest that the OE retrieval with a priori information included can process $t_0$ as a fit parameter. The information of the retrieved $t_0$ can also be useful to check the stability of this parameter during a set of independent measurements. These figures show that variations of $t_0$ during the set of 24 measurements and crosstalk with the other fit parameters are negligible.

The effect of the systematic errors on the forward model parameters has been investigated with Figs. 6 and 7, where we assumed a systematic error on the time origin $t_0$. Further tests with systematic errors on other forward model parameters such as $s_1$, ${\mu }_{s1}^{\prime }$, and ${\mu }_{s2}^{\prime }$ have been performed (data not shown). The results obtained can be summarized by noting that the distortion introduced by the systematic errors on fixed forward model parameters does not compromise the reconstruction of the absorption in the second layer when the OE method is used. These results further consolidate the test of reliability and stability of the reconstructions for ${\mu }_{a2}$. The results show that the errors on the fixed forward model parameters can affect the retrievals and that the OE allows to account for the effects of these errors in the inversion procedure. This functionality cannot be usually implemented with other procedures such as the LM.

It is interesting to note that the OE retrieval proposed in this work for a two-layered medium may also be used to process in vivo measurements on the adult head. Indeed, given the multicompartment architecture of the head, it is a useful approximation to merge the different compartments into two effective layers, identifying a superficial layer including scalp and skull and a deep layer including cerebrospinal fluid (CSF) and brain. This kind of approximation has been previously discussed, e.g., in Refs. 29 and 64. The effect of the CSF layer on the estimated optical properties, as it is argued in Ref. 64, is likely an underestimation of the reduced scattering coefficient of the second layer, as compared to the true value of brain tissue.

As for the computation time, as also noticed in Ref. 44, the LM is intrinsically faster by roughly a factor of 2 compared to the OE, but the difference between OE and LM also depends on the forward model. In the present contribution, the difference between the LM and OE was lower than a factor of 2 since the total computation time is mainly determined by the calculations of the forward solver for the two-layered medium and by the convolution with the IRF.

Finally, we stress that the presented results, being obtained with a forward solver given by an analytical solution of the DE, in principle may be affected by the approximations of the model. However, according to the results presented in Ref. 44, which were obtained with simulated data generated both with the DE and MC simulations, the approximations of the forward model based on the DE do not influence the overall beneficial effect of a priori information on the retrieval. Thus, this discussion and the conclusions of this work are not affected by the approximations of the forward model, and we do not have any deterioration of the usage of the diffusion approximation compared to the more exact radiative transport equation for the situations considered in this paper.

5.

Conclusions

This work is based on the analysis of single-distance TR measurements on a two-layered tissue simulating phantom. The results presented show that when a priori information on fit parameters and fixed forward model parameters is available and exploited by the OE method a larger number of fit parameters can be retrieved, and/or a better accuracy of the results can be achieved. In a previous work, the OE method was applied to synthetic data⁴⁴, including, in particular, forward simulations based on the DE with added Poisson noise and MC simulations. The present study shows that the optimal estimation method exhibits improved performance, compared to the LM method, also on experimental data. We have demonstrated the evidence of this finding by selecting six types of retrievals that can be considered paradigmatic of most of the retrievals usually carried out in real applications based on TR measurements such as for the reconstructions of the optical properties of muscle and brain tissues. The OE algorithm shows an excellent robustness in the retrieval of the optical properties of a two-layered medium from single-distance TR data. In particular, this characteristic is remarkable for the absorption coefficient of the second layer. For all the retrievals carried out, this parameter was always obtained with an error less than 10%. This is relevant for determining the absorption and oxygenation of the deep tissue compartment in muscle and brain applications. So far, we do not know of other published results with comparable performances.

The results obtained also demonstrate that the use of a priori information both on fit parameters and fixed forward model parameters reduces the crosstalk among the retrieved fit parameters. This is the main reason why the OE approach allows a larger number of fit parameters to be retrieved. For a representative set of retrievals, we have shown that the origin of time $t_0$ and the thickness of the first layer $s_1$ can be retrieved with excellent results in terms of stability and consistency. So far, this is the first time that such results were obtained from time-resolved single-distance measurements.

It is remarkable that in principle, with a single-distance time-resolved measurement and the OE algorithm, it is possible to efficiently reconstruct more than two fit parameters. Thus, the OE algorithm seems to be an appropriate way to exploit the full information content of such measurements.

As a final remark, we note that all the results presented in this work, although quite promising, are a first step only. A real challenge is the application to the analysis of in vivo measurements. The use of the OE method could lead to an improved analysis of the optical properties of the adult head based on TR tissue spectroscopy measurements. In particular, the extreme reliability of the retrieved value of the deep absorption coefficient is a solid ground to improve the robustness of the information retrieved from the cerebral cortex. The exploitation of the OE in the framework of in vivo measurements will be the subject of our future investigations.