## Abstract

We report a method for deriving the absolute value of absorption coefficients at depth in bilayered media. The method was simplified from that of time-resolved diffuse optical tomography (TR-DOT) into one dimension to validate and set up the main parameters with the help of simulations, and to test it in an easy preclinical model. The method was applied to buried flaps as used in reconstructive surgery, and absolute chromophore concentrations in the flap and in the upper (skin and fat) layer were derived. The encouraging results obtained lay a foundation for developing more complex multidimensional models.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Near infrared spectroscopy (NIRS) approaches have been developed for years to characterize the optical properties of diffusing biological tissue. Time-resolved (TR) NIRS approaches focus on the use of a pulsed laser for excitation and single-photon detectors with photon counting for detection. The photon distribution resulting from time of flight highlights the scattering and absorption contribution of the scanned medium, as well as the influence of depth: the longer the time of flight, the longer the photon has traveled in the medium, and thus the deeper it may have penetrated. These two properties - reconstruction of the scattering and absorption contributions, and sensitivity at depth leading to separation of the contributions of various layers- can be studied even in the case of a null source detector distance approach [1], allowing probe geometries with a small distance between the source and detector. Based on these properties, we propose a method for the noninvasive bedside monitoring of buried flaps in the context of reconstructive surgery. A major complication of such a surgery is thrombosis of the pedicle, which leads to flap necrosis. Since the flap is buried under a superficial layer of skin and fat, clinical examination with color observation and flap refill analysis is impossible. Moreover, the first layer may evolve as well - because of surgery - leading to a possible misinterpretation of the observations. Thus, there is a need for bedside noninvasive monitoring of hemodynamic parameters at depth. To achieve such a buried flap assessment, we developed a time-resolved system [2] that we tested on a pig model [3]. We implemented two main analysis approaches. One approach is based on a windowing strategy that has been described in the literature [4]; it provides oxy- and deoxyhemoglobin- concentration variations and allows the detection and identification of all arterial and venous occlusion events. We developed a new correction approach to determine the contribution of the superficial layer to the measured reflectance signals. The results are detailed in [5]. The second approach (object of this paper) is based on an underconstrainted reconstruction of a 3D map to obtain the absolute concentration of oxyhemoglobin and deoxyhemoglobin at depth. Indeed, despite promising results in different application domains (brain [6], breast [7], prostate [8]), facilitated by the use of free open-source software [9–13]], the accurate quantification of these absolute concentrations at depth is still an open issue in TR-NIRS. We aim to contribute to new results with the present proposal.

Using a time-resolved multidistance/multiwavelength data set to reconstruct a 3D map of the absorption parameter of the medium makes it possible to determine the local concentrations of targeted chromophores. Different approaches that focus on different data types, discretization schemes and analysis methods to optimize computational effort and memory issues have been proposed in the literature, as briefly discussed in [14]. Some years ago, we proposed a 3D reconstruction approach based on the use of Mellin-Laplace transforms, with a finite volume discretization scheme and an iterative optimization process built on the perturbation approach [15]. Some improvements have been proposed [16] by considering the biological sample as a set of distinct homogeneous areas; the resulting method uses a dictionary with labeled areas to induce spatial a priori constraints and reduce the solution space dimensions. This approach leads to better quantification of chromophore concentrations at depth [17]. In the particular case of buried flaps, one can assume a homogeneous evolution in both layers (the buried flap and superficial layer). Moreover, the observed signal variations are due to hemoglobin variations in the tissue; thus, only the absorption evolves, and the scattering is assumed to be constant during the observation time. Therefore, we optimized our 3D reconstruction method in the present paper with a focus on a bilayered model: because of surgery, the upper layer may itself be subjected to infusion hazards; these surface variations are likely to mask variations induced by a more severe deep occlusion problem. The challenge is to reconstruct both the absorption properties of the upper surface layer and those of the lower layer (the flaps itself), assuming the scattering is known. The theory of our method is resumed in section 2 of this paper. We apply the whole reconstruction method to two buried flaps with occlusion in pigs in section 3, and conclude by discussing the results in section 4.

## 2. Model development

#### 2.1 Main points of the algorithm

For several years, our group has worked to develop a time-resolved diffuse optical tomography (TR-DOT) approach based on the absence of the need for the instrument response function (IRF). The IRF is replaced by a reference measurement, acquired in a medium whose optical properties are known or assumed to be known, as will be developed in section 2.2. This choice was made in our group considering the difficulty of measuring the IRF with some devices, particularly built-in optical fiber probes or probes containing many fibers. We recall the main points of the approach to highlight the role of the reference. The interested reader will find more details in [15].

**Main principles**

- • We have access to a limited number of distributions of photon times of flight (DTOFs), measured by a limited number of spatially distributed detectors and originating from a limited number of spatially distributed sources. If we note $\overrightarrow {r_S}$ and $\overrightarrow {r_D}$ as the source and detector positions in space, the signal issued from $\overrightarrow {r_S}$ and measured at $\overrightarrow {r_D}$ is denoted as $u_S(\overrightarrow {r_{D}},t)$.
- • By using a reference measure in a medium whose absorption and scattering coefficients are assumed as known (hereafter denoted as the medium "ref") and performing the same measurements in that medium as in the current one, the cost function to minimize is:$$\chi^2 =\lVert \Sigma_{ \left\{\overrightarrow{r_{S}},\:\overrightarrow{r_{D}},\: t\right\} } Y_{SD}(t) \rVert^2$$with$$Y_{SD}(t)=u_S(\overrightarrow{r_{D}},t)\ast G^{ref}(\overrightarrow{r_{D}},\overrightarrow{r_{S}},t)-u_S^{ref}(\overrightarrow{r_{D}},t)\ast\tilde{G}(\overrightarrow{r_{S}},\overrightarrow{r_{D}},t)$$where the $\widetilde {}$ superscript is applied to the values predicted by the model. $G(\overrightarrow {r}_{1},\overrightarrow {r}_{2},t)$ is the Green function solution of the diffusion equation at position $\overrightarrow {r_{2}}$, time t for a Dirac source at position $\overrightarrow {r_{1}}$, and time 0. The symbol $\ast$ represents the temporal convolution. To resolve the model, we need to estimate the optimal Green function $G(\overrightarrow {r_{S}},\overrightarrow {r_{D}}, t)$ via the absorption distribution $\mu _{a}(\overrightarrow {r})$ at any position $\overrightarrow {r}$ in the medium of interest.

To solve Eq. (2), we use a reconstruction approach based on the Mellin-Laplace transform of the DTOF with a finite volume discretization scheme. To keep the perturbations small, we solve the problem through an iterative approach over steps denoted as $k$. This leads to Eq. (3):

with $\mathcal {W}_{SD}^{^{k-1}}$ indicating the Jacobian matrix in the Mellin-Laplace space, $\mathcal {Y}_{SD}^{k}$ indicating the measurements vector for source S and detector D in the Mellin-Laplace space, and $X^{k}$ indicating the vector of unknowns, of size the number of volumetric elements.**Labeled area dictionary approach**

Conditioning matrices [19] and/or regrouping voxels by interest areas [16,17] and layering can be introduced to compensate for the ill-posedness of the system. In the case of buried flaps, as explained in the introduction, it is possible to implement a bilayered model to monitor both the upper layer covering the flap and the flap itself. To this end, we introduced a matrix, called D, that is the number of voxels by the number of unknowns after the bilayered reshaping in size, which is 2. Matrix D consists of two columns of 1 and 0, mapping each voxel to its corresponding layer. The unknown vector becomes $2\times 1$ in size with $X^{k}=D\times X'^{k}$. The resolution of Eq. (3) becomes straightforward: all the signal measures (for every source and detector position SD) can be concatenated in a single vector $\mathcal {Y}$. For each iteration step $k$, and until convergence, the direct model is solved by finite volume modeling, and the global matricial system turns into:

with Jacobian $\mathcal {W}^{^{k-1}} \times D$, of a small size (number of measures $\times$ number of MLT) $\times$ 2, that can be inversed analytically quickly at each iteration step k.#### 2.2 Absorption initialization and reference medium

Although the method has been widely detailed in some of our group’s papers [15,16,20], the influence of the reference and of the absorption initialization step has been less well described. We wish to address these two points in the particular context of the two-layered model applied here for buried flaps. By means of numerical simulations, we will see that the choice of these parameters is fundamental for assessing the convergence and accuracy of the quantified absorption of both the surface layer and the deep layer.

In this section, we focus on absorption reconstruction, so only one wavelength is used in the signals. Multiwavelength measurements and chromophore reconstruction will be mentioned in section 2.3.

### 2.2.1 Simulation parameters

We simulated signals backreflected from surface source points onto surface detector points through a homogeneous bilayered medium. The original probe geometry tested for the simulations is that of our experimental probe, consisting of 5 sources and 4 detectors, with 3 different source/detector distances available: 6, 8 and 13 mm (see Fig. 1).

The signals are generated by the convolution of a typical instrument response function and the Green function of the simulated media (obtained by finite-volume approach from the diffusion equation in a regularly 3D voxelised media). The geometry of the medium voxels in the reconstruction is chosen to be identical to that of voxels in the simulation.

The algorithm parameters were optimized using preliminary simulations and reconstructions. In particular, as an outcome of these preliminary simulations, the Mellin-Laplace parameters are set to n=15 orders and a precision of p=4 $ns^{-1}$. Preliminary reconstructions to find the depth that can be reconstructed as a function of medium absorption showed that noise is the only factor limiting the reconstruction if the whole medium absorption initialization is that of the top layer, and if the reference medium has the same absorption and scattering coefficients as the top layer. In this case, the algorithm perfectly reconstructs both the upper and bottom layers, regardless of the bottom-layer absorption (variations in the bottom-layer absorption between -75% and +250% that of the top layer were tested).

This supports the robustness of the method to probe a lower layer consisting of varying media. In the following, we will use $10^8$ photons per signal and simulate a 10-mm-thick layer on top of a 16-mm-thick layer. The longitudinal dimensions are 28 by 30 mm. The acquisition geometry corresponding to 1 source point and 2 detection points 13 mm from the source was selected for all simulations in this paper. For each reconstruction point, 20 similar noise variant-only reconstructions were made, and the values were averaged and the standard deviation tagged. This choice of 20 tests was validated by comparison to a set of 100 tests and presented the same average and standard deviation values. The noise model is the sum of the Poisson noise applied to the number of photons detected and a detector dark count rate. Each reconstruction is made with a different noise reference signal.

### 2.2.2 Absorption initialization impact

To test the effect of the absorption initialization parameter, a bilayered medium is simulated, with a constant reduced scattering (set at 10 $cm^{-1}$), a top layer varying slightly in absorption around physiological values in humans (from 0.1 to 0.3 $cm^{-1}$), and a bottom layer varying widely in absorption from smaller to larger than that of the top layer (from 0.05 to 0.5 $cm^{-1}$), to mimic flap occlusion (see Fig. 2). The reconstruction algorithm is set on a bilayered model with homogeneous values and a known thickness for each layer. Two variables are then reconstructed: the absorption of the upper and bottom layers.

The reduced scattering parameter in the reconstruction model is set to its simulation value (i.e., no error is introduced in the scattering estimation). The reference medium is taken as a uniform medium identical to that of the top layer (its variation effects will be described below in section 2.2.3). The absorption initialization is tested with values varying from -30% to +30% that of the upper layer.

The results are as follows:

- • Whatever the initialization values, the top layer is reconstructed correctly with a relative error $\frac {\mu _{a}\textrm {reconstructed}-\mu _{a}\textrm {simulated}}{\mu _{a}\textrm {simulated}}$ less than 0.1%.
- • The bottom-layer reconstruction presents a larger error, but this error is always lower than 16% (obtained for high absorption in the bottom layer), including the error bar at 1$\sigma$.
- • As mentioned in 2.2.1, if the initialization absorption is exactly that of the top layer, the error on the reconstruction is negligible for both layers and only due to photonic and detector noise, whatever the bottom layer absorption.

These results confirm the preference to initialize absorption with values not too far from those of the media to reconstruct absorption, so that the assumption of small variations is true at the first iterative step of the algorithm. However, subsequent variation is tolerable without major error in the reconstruction, and globally, initialization close enough to the top-layer absorption will give satisfactory results.

### 2.2.3 Impact of the reference medium

Numerical simulations using a reference medium with known absorption and scattering coefficients (such as a phantom medium) indicated that the reference has to be identical in scattering to the medium of interest. Indeed, in our iterative reconstruction scheme, the Jacobian matrix is computed according to Green functions with a fixed $\mu _s'$ and evolutive $\mu _a$: thus, an incorrect $\mu _s'$ will introduce major crosstalk between scattering and absorption. For example, simulations showed that a 10% difference in $\mu _s'$ between the medium of interest and the reference medium leads to significant errors in the reconstructed absorption parameters (>100%).

As a solution, in the context of the monitoring of buried flaps addressed in this paper, we decided to take for reference the medium itself. Thus, the scattering coefficient of the reference equals that of the reconstructed volume. The reference can be taken at the beginning of the experiment (at rest, just after flap surgery), or be chosen differently for every reconstructed point.

The reference absorption will then have to be estimated, along with its reduced scattering. Sensitivity to these two parameters for the reference measure is explored here by numerical simulations and reconstructions.

In these simulations, the top layer is homogeneous and constant (with absorption set at 0.2 $cm^{-1}$ and reduced scattering set at 10 $cm^{-1}$), the bottom-layer absorption varies among the simulations, and its scattering is the same as that of the first layer. The optical coefficients of the reference medium are those of the first layer. To mimic incorrect estimation of those coefficients, the reference scattering and absorption coefficients used in the model are assumed to be 10% higher or lower than their actual values. The initialization absorption coefficient is taken to be the same as that of the top layer. This is summarized in Fig. 3.

For the 10% reduced scattering error in the estimation (Fig. 3(a)), the results are as follows:

- • The top-layer absorption is reconstructed with an error less than 0.01%, regardless of the bottom-layer absorption.
- • The bottom-layer absorption is reconstructed well, with a relative error less than 6% (achieved for high absorption values), including error bars at 1$\sigma$.

The absorption reconstructions are good for this test scheme, which is explained by the fact that scattering is incorrectly estimated for both the reference medium and the medium of interest, so that an error in its value has only a limited impact.

For the 10% absorption error in the estimation (Fig. 3(b)), the results are as follows:

- • The top-layer absorption reconstructed is identical to that of the estimated reference $\pm 0.1\%$.
- • The bottom-layer reconstructions are not good, with relative errors up to 40%, including error bars at 1$\sigma$.

These simulations physically correspond to incorrect calibration of the top layer, with all of the calibration error being deferred to the bottom layer. Top-layer overestimation implies bottom-layer underestimation, and vice versa. The reference acts as a calibration tool whose relevance depends on the accuracy in absorption of the initial estimates.

#### 2.3 Estimating the absorption of a reference based on the medium of interest

To estimate absorption at any wavelength used for measurement, we use the comparison between the total number of photons from measurements at several wavelengths. We deduce the partial oxygen saturation ($StO_2$) and the hemoglobin fraction ($f$) according to a Lookup Table (LUT) approach, as described in [21, 22]. Then, given the $StO_2$ and $f$ parameters, it is possible to determine the absorption coefficient at any wavelength. Indeed, the $StO_2$ is linked to the ratio of absorption coefficients at various wavelengths, while the hemoglobin fraction $f$ is linked to a global absorption multiplication factor. Several measurements (from at least three wavelengths) need to be available for that estimation method. The interest in that method compared to classical fitting methods, such as presented in [23–25], lies in the search for a working point coherent for all wavelengths and not in an independent estimation of absorption at every wavelength.

**LUT approach**

We established a Lookup Table approach according to the following steps:

- • The signals coming from a uniform medium at various hemoglobin fractions and $StO_2$ values are simulated at three wavelengths (750 nm, 800 nm and 850 nm, which are the instrument wavelengths we use for in vivo and phantom measurements) for the following reduced scattering and absorption coefficients:
- – Absorption: Absorption is calculated through the hemoglobin fraction $f$ and $StO_2$ based on Eq. (5), considering hemoglobin the main contributor to tissue absorption at our wavelengths:$$\mu_{a}(\lambda_{k})=H{}_{rate}\times(StO_{2}\times\varepsilon_{HbO_{2}}(\lambda_{k})+(1-StO_{2})\times\varepsilon_{Hb}(\lambda_{k}))\times f\,,$$where $H_{rate}$ is the hematocrit of the blood (estimated through typical values), and $\varepsilon _{HbO_{2}}(\lambda _{k})$ and $\varepsilon _{Hb}(\lambda _{k})$ are the oxy- and deoxyhemoglobin molar extinction coefficients at wavelength $\lambda _{k}$. They are derived from tabulated molar extinction coefficients for hemoglobin in water from OMLC.
The calculations are made with the typical blood values of 150 grams of hemoglobin per liter and 50% hematocrit. Melanin and water are not taken into account in this model.

- • The total number of photons of each simulated signal is calculated and the ratios ($N_{800}/N_{750}$, $N_{75}/N_{850}$, $N_{850}/N_{800}$) of the outcoming photons are 3D mapped as a function of the $StO_{2}$ and hemoglobin fraction.
- • Instrument-derived attenuation, which is different at the three wavelengths and detector dependent, is calibrated by measures of phantoms with our instrument and compensated.

**Fitting step and absorption estimation**

Once these LUTs were established, we considered the medium as homogeneous at the start of the experiment (resting state). We then fitted the ratio of photons at the three wavelengths to a hemoglobin fraction and $StO_2$. These values were then converted into absorption coefficients through Eq. (5).

**Remarks**

The definition of the LUT requires the knowledge or estimation of the reduced scattering coefficient at each wavelength. An error in this scattering coefficient will lead to an error in the estimated absorption coefficient. Here, we used sr-DRS to obtain this coefficient, but in the case of a known IRF, the TR signal itself can be used for this objective. A specific fitting approach can be developed to consider $\mu _s'$, $StO_2$, and $f$ as parameters to estimate given the three wavelengths.

Note that if water can be assumed to be constant during the whole experiment (which is not always the case in buried flaps), we could introduce a baseline for taking it into account in the method. This would need an evaluation of the amount of water in the tissue, for example, by a table of values published in the literature. The water concentration could also be reconstructed apart from the hemoglobin concentration by adding a measurement at a higher wavelength more sensitive to its variation.

## 3. Reconstructions in a porcine bilayered flap model

To validate this original method for absorption coefficient evaluation and associated multilayered reconstruction in a time-evolving medium, we performed experiments in porcine models of flaps buried 1 cm deep with venous or arterial occlusion. As part of our multilayered reconstruction project [26], two flaps buried 1 cm deep were created on 15 pigs [3]. Venous or arterial occlusion was then induced in each buried flap. During the entire operation, the flaps were measured with our time-resolved optical instrument to quantitatively reconstruct the absorption at any time. A disposable Licox probe, as an invasive nonquantitative reference for blood oxygenation estimation, was placed inside the flap to provide a standard for comparison with our measures. The clinical setup is depicted in Fig. 4.

The measurement prototype instrument used a supercontinuum fiber laser (Fianium WL-SC-480-10) generating picosecond laser pulses at 80 MHz. A first wavelength range selection was performed by two longpass filters, after which specific wavelengths were selected by means of interferential filters mounted on a motorized wheel and centered at 750, 800 and 850 nm with a bandwidth of 10 nm. Light exiting from the interferential filter was then attenuated by means of a motorized variable optical attenuator (at optical densities from 0 to 4). Before a 1x6 optical fiber switch sending light to one of the five probe sources, half of the signal was collected by a dual-window fiber coupler to check the stability of the laser. The mean source power was below 1 mW at the excitation fiber output, with a numerical aperture of 0.22, meeting security standards. For detection, we used 4 hybrid photomultipliers (Becker-Hickl GmbH, HPM 100-50) coupled to time-correlated single-photon-counting electronics. A synchronous signal was collected at the highpass filter to trigger the TCSPC boards. A longpass filter was added in front of the sensors to optimize the useful signal while reducing noise interference below 690 nm. A probe was designed for the experiments (Fig. 1) as an optical stethoscope built with materials suitable for clinical use. The shape of the probe allows a good grip and maintenance. The distances between the sources and detectors (SD) were 6, 8 and 13 mm. One source (source 5 in Fig. 1) and two detectors (detectors 1 and 2 in Fig. 1) were used for the measurements.

We chose here to report the bilayered absorption reconstruction and derived chromophore concentrations of one venous and one arterial flap as an illustration of the method detailed above. In that respect, for each flap, the reference measure was the first measure after the surgery (in a medium considered uniform, before any time evolution of the flap should occur). The reference absorption coefficients were estimated by the method described in section 2.3. The reduced scattering coefficients were estimated through sr-DRS, which is a standard validated method for estimating absolute superficial optical coefficients [22]. We thus considered the value measured through this method as a precise estimation of the reduced top-layer scattering and extended that scattering value to the whole sample as a first approximation of the reference medium (considered uniform).

#### 3.1 General parameters

Reduced scattering coefficients were measured by sr-DRS in several pigs, for an average of $15\,cm^{-1}$ at 750nm, $14\,cm^{-1}$ at 800nm and $13\,cm^{-1}$ at 850nm. The hematocrit in pig blood was fixed at 45%, an average for pigs given data in the literature [27]. The flap and top-layer widths were measured by ultrasound immediately after flap surgery. The two pigs chosen for the experiments weighed 88 kg and 92 kg. Flaps were collected from the thigh and buried in the right flank. Total arterial occlusion was induced in the first pig, and total venous occlusion was induced in the second pig. The general protocol consisted of an experiment lasing 1 h 20 minutes divided into a resting period of approximately 30 minutes after flap surgery, a clamping period lasting approximately 20 minutes, a declamping period and a second resting period of approximately 30 minutes. We measured a complete reconstruction point every 1.5 minutes, and all the reconstruction points are shown in the figures. The initialization parameters were those estimated for the reference. The reconstructions used the layer widths given by the ultrasound measurements. The Mellin-Laplace parameters, as set for the simulations, were chosen as 15 orders and a precision of 4 $ns^{-1}$. The medium was regularly 3D-voxelized, the spacing between the summits of two voxels in Cartesian coordinates was 2 mm in the X and Y directions and 1 mm in the Z direction, with Z being the vertical direction, and X and Y being the surface directions. The medium dimensions were set to 30$\times$28$\times$26 mm.

To retrieve the tissue blood oxygenation information, absorption was reconstructed at three wavelengths and hemoglobin chromophores were deduced for each flap. The reconstruction was based on Beer’s law of mixtures, which states that the absorption spectrum is a weighted average of its main chromophores.

#### 3.2 Arterial occlusion

### 3.2.1 Absorption reconstruction

The flap width, measured by ultrasound after surgery, was 1.15 cm. The flap was buried under a top layer with an average width of 0.95 cm that was assumed to be mainly composed of dermal, fat and muscle tissue. The bilayered model thus appears here obviously as an approximation, the absorption evolution of each layer being representative of the evolution of an equivalent homogeneous tissue. This model was chosen here since the chromophore evolution in the whole flap is expected to be very homogeneous at the time scale of a measurement, although the physiological composition of the flap may not be. The first measurement after surgery was taken as reference. Applying the method described in section 2.3, fitting using the number of photons at each wavelength gives an $StO_2$ of 82% and a hemoglobin fraction of 2.2%. This corresponds to absorption coefficients of 0.071 cm-1, 0.085 cm-1 and 0.105 cm-1 at 750, 800 and 850 nm, respectively (according to Eq. (5) with the parameters specified in section 3.1). Reconstructions for all measurement points after surgery are given in Fig. 5, with a moving average of each three data points applied. The instant of arterial clamping, at 10 h 33 minutes, is marked in Fig. 5 by a vertical dashed line. This clamping is followed at 10 h55 minutes by declamping, which is also marked by a vertical dashed line. The Licox data (in black) show a progressive decrease and increase after the points of clamping and declamping, respectively; the units of these data are arbitrary. For the top layer, the reconstructed absorption at the three wavelengths continuously increases by approximately 9% between the beginning and end of the measures, with no significant manifestation of the clamping and declamping in the bottom layer. The bottom-layer absorption shows a decreasing slope after surgery before the clamping. After the clamping, to the absorption decreases more consistently at 800 and 850 nm and then increases at 750 nm. The three absorption values return to their regular slope trends after declamping.

### 3.2.2 Chromophore retrieval

Figure 6 shows the absolute oxy- and deoxyhemoglobin chromophore values retrieved from the absorption values for the top layer (a) and for the bottom layer (b), i.e., the flap, which underwent clamping. A moving average over 3 reconstruction points is applied.

Before commenting on the trend in the concentration of the two chromophores, we wish to comment on their absolute values. The concentration of each chromophore is calculated in $mol.L^{-1}$ in the tissue as a whole. Considering the average hematological value of 116 g of total hemoglobin per liter of blood of a pig [28], given $64.5\times 10^{3}$ $g.mol^{-1}$ as the molar weight of hemoglobin and an average blood content of 7% in bulk tissue, we calculate $12.6\times 10^{-5\,\,}mol.L^{-1}$ as the hemoglobin concentration in the tissue. Here, the total hemoglobin concentration in the top and bottom layers before clamping was approximately $11\times 10^{-5\,\,}mol.L^{-1}$. The ratio of oxyhemoglobin to deoxyhemoglobin in the top and bottom layers at rest was close to 4.5, corresponding to an $StO_2$ of 82% (green dashed curve), as was set by the reference $StO_2$ calculation method. This means that 82% of the hemoglobin contained in the blood in the tissue is in its oxygenated state, while 18% is in its deoxygenated state.

Regarding the variation in the chromophores during the clamping period, in terms of absorption, in the bottom layer, both chromophores show a decreasing baseline after surgery. After clamping, oxyhemoglobin absorption decreases, while deoxyhemoglobin absorption increases, and after declamping, both return to their respective preclamping values. These variations are coherent pathophysiologically. Indeed, during arterial clamping, oxygenated blood flow is disrupted while deoxygenated blood flow is maintained. The blood volume also decreases slightly. The $HbO_2$ concentration should fall sharply, while the Hb concentration should remain initially stable. During clamping, the tissue consumes oxygen, so the $HbO_2$ concentration continues to decrease progressively as oxyhemoglobin is transformed into Hb [29,30]. In the top layer, a slightly increasing baseline is present. In the bottom layer, the drop in $StO_2$ during arterial clamping (down to 58%) corresponds to values reported for total arterial clamping in rats [31].

### 3.2.3 Influence of top layer geometry

The decrease in $HbO_2$ and increase in $Hb$ may be minimized by our reconstruction algorithm because of the underestimation of the top-layer width. In this way, some invariant top-layer medium is included in what is reconstructed as the bottom layer, minimizing its variations during the clamping period. Therefore, we calculated the same reconstruction but with the top layer assumed to be 1 mm thicker. To complete the analysis, we also performed the calculations with a top layer 1 mm thinner. Observing the impact of a small error in the top layer estimation is also interesting since reconstructions without US measures may be desired in the future. Figure 7 shows the absorption reconstructions during arterial clamping for a top-layer width assumed to be 1 mm thinner and 1 mm thicker. The thicker the assumed top-layer width was, the more pronounced the absorption variations of the bottom layer. The top-layer absorption coefficients were always well reconstructed.

#### 3.3 Venous occlusion

### 3.3.1 Absorption reconstruction

The flap width, measured by ultrasound after surgery, was approximately 1.2 cm. It was buried under a top layer of 1.2 cm. As for arterial occlusion, the first measure after surgery was taken as a reference. Fitting the ratio of the measured photon number here gives an $StO_2$ of 72% and a hemoglobin fraction of 2.1%, corresponding to absorption coefficients of 0.078 cm-1, 0.083 cm-1 and 0.098 cm-1 at 750, 800 and 850 nm, respectively.

Reconstructions for all the measurement points after surgery are given in Fig. 8, with a moving average of each three data points being applied. The instant of venous clamping, at 15 h 24 minutes, is marked in Fig. 8 by a vertical dashed line. Venous clamping was followed by declamping at 15 h 49 minutes, also marked by a vertical dashed line in the figure. The Licox data (in black) show a 2-minute delay after the instant of clamping. In the top layer, the reconstructed absorption at the three wavelengths remains stable between the beginning and end of the measurement period. In the bottom layer, absorption slightly decreases before clamping and increases consistently immediately after clamping at the three wavelengths, reaching 3 times its initial value at 750 nm. Just after declamping, absorption at the three wavelengths decreases immediately and tends to stabilize at roughly the same time as the Licox data but at an absorption level higher than that before clamping.

### 3.3.2 Chromophore retrieval

Figure 9 shows the absolute oxy- and deoxyhemoglobin chromophore values retrieved for the absorption values for the top layer (a) and the bottom layer (b), i.e., the flap, which underwent clamping. A moving average over 3 reconstruction points is applied. The absolute concentration of the two chromophores in the top and bottom layers during the resting period was close to its value at rest in the other pig (Fig. 6) but with a slightly higher $Hb$. In the top layer (Fig. 9(a)), these values remained stable throughout the experiment, as expected. In the bottom layer (Fig. 9(b)), clamping induced a strong rise in the total hemoglobin concentration (dashed brown curve). That rise corresponded to an immediate moderate increase in oxyhemoglobin (red curve) after clamping, followed by a greater increase in deoxyhemoglobin (blue curve) some minutes later. Despite the noise, it can be observed that the oxyhemoglobin level slowly returned to its preclamping value, even before declamping, as the deoxyhemoglobin level increased. After declamping, the deoxyhemoglobin level decreased again to its preclamping value. These trends are coherent pathophysiologically, since clamping is expected to induce tissue congestion, resulting in a strong increase in oxy- and total hemoglobin. This increase is followed by a decrease in the level of deoxyhemoglobin as tissues consume oxygen [29,30]. Additionally, the $StO_2$ evolves from 80 to 90% before clamping down to 40% during clamping.

#### 3.4 In-vivo reconstruction outcome

Evolution of the absorption coefficient showed a behavior coherent with clinical expectations and absolute reconstructed values in accordance with expectations in biological tissues. The absorption values were nearly stable in the top and bottom layers at rest, decreased during arterial clamping at two out of the three wavelengths, and increased during venous clamping at the three wavelengths. For both flaps, chromophores in the top and bottom layers before clamping showed nearly unvarying behavior, and the absolute concentrations were consistent with expectations. Their behavior during arterial and venous clamping was coherent pathophysiologically, although noisy. In both cases, clamping could be easily detected and identified by analyzing the evolution of the chromophore concentration.

## 4. Discussion

Given the importance of using initialization parameters as close as possible to those of the medium of interest (see the simulations in section 2.2), in a continuously time-evolving medium, an intuitive choice for initialization would be the reconstructed absorption value for each previous point in the medium. In this way, a strong variation in absorption during clamping or even at baseline could be smoothed in terms of initialization errors. That initialization choice was tested on the pig flaps, and the reconstruction results were very similar to those obtained with the first point initialization. This is because the top-layer absorption variation is quite small in experimental flaps, so even by initializing all the points with the first measure, the initialization absorption is close to those of the top layer at any point. These considerations validate the use of a single initialization value for all points in one experiment.

The same kind of considerations could be applied to the reference measures, which could consist of the previous measure for each point. Indeed, an outcome from the calculations carried out in section 2.2 was that the closer the reference medium is from the top layer, the smaller the error gap is if badly estimated in absorption. Then, taking as a reference each previously reconstructed measure could be interesting, while remaining aware that this would propagate the reconstruction error from one point to the other. Absorption reconstruction for both flaps was tested with such a reference scheme (the initialization value in these tests was taken as identical to the reference value), and again nearly no difference emerged in the reconstructed points. This is because although the reconstruction error between the reference and medium of interest was minimized for each point, the error for the medium of interest propagated from one point to the other instead of lying only in the difference between the reference and the reconstructed medium.

## 5. Conclusion

We developed a reconstruction method based on TR-DOT that retrieves the absorption coefficients of a medium composed of homogeneous area. The aim of this work was to explore and validate the 3D reconstruction algorithm we developed in a simplified context of buried flaps. In the context of total clamping (arterial occlusion and venous occlusion) as performed in this paper, we considered a bilayered model with homogeneous layers. We proceeded to the absolute absorption estimation by simplifying the 3D-DOT method into a 1D TR-DOT method through a strong a priori spatial constraint. A reference medium, similar in reduced scattering to the medium and with absorption estimated as precisely as possible, is necessary to our method. Thus, in biological tissue experiments, this reference has to be taken in the tissue. A method for calculating absorption coefficients for that reference measure was developed.

The whole method was tested in vivo in pigs with buried flaps. Absolute absorption concentration was reconstructed on arterial and venous flap clamping. Absolute oxy- and deoxyhemoglobin chromophore concentrations were then reconstructed as well. These absolute values and the evolution of these values, though noisy, were as expected pathophysiologically.

Once simplified for two layers, our method can be compared to spectroscopic approaches, such as time-windowing TR-NIRS or moments TR-NIRS. Although these methods focus mainly on absorption variations in a time-evolving medium [32–35], some of them report absolute values by fitting the initial absorption value with the solution of the diffusion equation for semi-infinite homogeneous uni- or bilayered media [24,25,36]. A quantitative comparison of TR-NIRS methods with our TR-DOT method would be interesting.

The flap animal model we developed is a powerful model for validating our algorithm. However, the algorithm is designed to be used in more complex situations (in the brain, for example, with 4 areas and several sources and detectors), that cannot be addressed by curve-fitting approaches. Moreover, in a real clinical situation, the occlusion may be partial, leading to the need for a 3D map of perfusion. In the specific context of this paper (buried flaps with total occlusion), a simpler approach such as the windowing approach may be sufficient to provide diagnostic information in this specific application [5].

Finally, we wanted to emphasize that by providing absolute values of chromophore concentrations, we aim to enable easy plug-and-play monitoring without the need for the instrument to stay on the patient throughout the monitoring period. This would also enable, for example, a unique instrument for use by a whole hospital. Clinical translation of the method and the instrument is now under study.

## Funding

Agence Nationale de la Recherche (ANR-155-CE19-0010); Union des Blessés de la Face et de la Tête.

## Disclosures

The authors declare no conflicts of interest.

## References

**1. **A. Torricelli, A. Pifferi, L. Spinelli, R. Cubeddu, F. Martelli, S. Del Bianco, and G. Zaccanti, “Time-resolved reflectance at null source-detector separation: improving contrast and resolution in diffuse optical imaging,” Phys. Rev. Lett. **95**(7), 078101 (2005). [CrossRef]

**2. **A. Planat-Chrétien, A. Dot, M. Perriollat, M. Berger, R. Lartizien, J.-L. Coll, and G. Bettega, “Non-invasive assessment of deep buried flap viability with time-resolved optical monitoring: results on pigs,” in Biophotonics Congress: Biomedical Optics Congress 2018 (Microscopy/Translational/Brain/OTS), (Optical Society of America, Hollywood, Florida, 2018), OSA Technical Digest, p. JW3A.27.

**3. **R. Lartizien, M. Bouyer, L. Vinit, and G. Bettega, “Porcine abdominal cutaneous flap model: A simple and reliable method,” J. Stomatol. Oral Maxillofac. Surg. **119**(4), 301–303 (2018). [CrossRef]

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

**5. **R. Lartizien, A. Planat-Chrétien, M. Berger, M. Henry, J.-L. Coll, A. Dot, and G. Bettega, “Noninvasive monitoring of deep tissue oxygenation in buried flaps by time-resolved near-infrared spectroscopy in pigs,” Plast Reconstrc. Surg. **146**(5), 565e–577e (2020). [CrossRef]

**6. **R. J. Cooper, E. Magee, N. Everdell, S. Magazov, M. Varela, D. Airantzis, A. P. Gibson, and J. C. Hebden, “MONSTIR II: A 32-channel, multispectral, time-resolved optical tomography system for neonatal brain imaging,” Rev. Sci. Instrum. **85**(5), 053105 (2014). [CrossRef]

**7. **J. M. Cochran, D. R. Busch, L. Lin, D. L. Minkoff, M. Schweiger, S. Arridge, and A. G. Yodh, “Hybrid time-domain and continuous-wave diffuse optical tomography instrument with concurrent, clinical magnetic resonance imaging for breast cancer imaging,” J. Biomed. Opt. **24**(05), 1–11 (2019). [CrossRef]

**8. **J. He, R. Weersink, I. Veilleux, K. Mayo, A. Zhang, D. Piao, A. Alam, J. Trachtenberg, and B. C. Wilson, “Development of transrectal diffuse optical tomography combined with 3D-transrectal ultrasound imaging to monitor the photocoagulation front during interstitial photothermal therapy of primary focal prostate cancer,” Proc. SPIE **8578**, 85781J (2013). [CrossRef]

**9. **H. Dehghani, M. E. Eames, P. K. Yalavarthy, S. C. Davis, S. Srinivasan, C. M. Carpenter, B. W. Pogue, and K. D. Paulsen, “Near infrared optical tomography using NIRFAST: Algorithm for numerical model and image reconstruction,” Commun. Numer. Meth. Engng. **25**(6), 711–732 (2009). [CrossRef]

**10. **M. Jermyn, H. R. Ghadyani, M. A. Mastanduno, W. D. Turner, S. C. Davis, H. Dehghani, and B. W. Pogue, “Fast segmentation and high-quality three-dimensional volume mesh creation from medical images for diffuse optical tomography,” J. Biomed. Opt. **18**(8), 086007 (2013). [CrossRef]

**11. **M. Schweiger and S. R. Arridge, “The Toast++ software suite for forward and inverse modeling in optical tomography,” J. Biomed. Opt. **19**(4), 040801 (2014). [CrossRef]

**12. **Q. Fang, “Mesh-based Monte Carlo method using fast ray-tracing in Plücker coordinates,” Biomed. Opt. Express **1**(1), 165–175 (2010). [CrossRef]

**13. **F. Hecht, “New development in freefem++,” J. Numer. Math. 20 **20**(3-4), 251–266 (2012). [CrossRef]

**14. **D. Orive-Miguel, L. Hervé, J. Mars, L. Condat, and P. Jallon, “Time-resolved diffuse optical tomography: a novel method to compute datatypes allows better absorption quantification,” in Clinical and Preclinical Optical Diagnostics II, vol. EB101 of *SPIE Proceedings* (Optical Society of America, Munich, 2019), p. 11074_8.

**15. **L. Hervé, A. Puszka, A. Planat-Chrétien, and J.-M. Dinten, “Time-domain diffuse optical tomography processing by using the MLT,” Appl. Opt. **51**(25), 5978–5988 (2012). [CrossRef]

**16. **J. Zouaoui, L. Di Sieno, L. Hervé, A. Pifferi, A. Farina, A. D. Mora, J. Derouard, and J.-M. Dinten, “Quantification in time-domain diffuse optical tomography using Mellin-Laplace transforms,” Biomed. Opt. Express **7**(10), 4346 (2016). [CrossRef]

**17. **J. Zouaoui, L. Di Sieno, L. Hervé, A. Pifferi, A. Farina, A. D. Mora, J. Derouard, and J.-M. Dinten, “Chromophore decomposition in multispectral time-resolved diffuse optical tomography,” Biomed. Opt. Express **8**(10), 4772 (2017). [CrossRef]

**18. **S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. **15**(2), R41–R93 (1999). [CrossRef]

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

**20. **A. Puszka, L. Hervé, A. Planat-Chrétien, A. Koenig, J. Derouard, and J.-M. Dinten, “Time-domain reflectance diffuse optical tomography with Mellin Laplace transorm for experimental detection and depth localization of a single absorbing inclusion,” Biomed. Opt. Express **4**(4), 569–583 (2013). [CrossRef]

**21. **V. Sorgato, M. Berger, C. Emain, C. Vever-Bizet, J.-M. Dinten, G. Bourg-Heckly, and A. Planat-Chrétien, “ACA-Pro: calibration protocol for quantitative diffuse reflectance spectroscopy. Validation on contact and noncontact probe- and CCD-based systems,” J. Biomed. Opt. **21**(6), 065003 (2016). [CrossRef]

**22. **V. Sorgato, M. Berger, C. Emain, C. Vever-Bizet, J.-M. Dinten, G. Bourg-Heckly, and A. Planat-Chrétien, “Validation of optical properties quantification with a dual-step technique for biological tissue analysis,” J. Biomed. Opt. **23**(09), 1–14 (2018). [CrossRef]

**23. **L. Zucchelli, D. Contini, R. Re, A. Torricelli, and L. Spinelli, “Method for the discrimination of superficial and deep absorption variations by time domain fNIRS,” Biomed. Opt. Express **4**(12), 2893–2910 (2013). [CrossRef]

**24. **H. Auger, L. Bherer, É. Boucher, R. Hoge, F. Lesage, and M. Dehaes, “Quantification of extra-cerebral and cerebral hemoglobin concentrations during physical exercise using time-domain near infrared spectroscopy,” Biomed. Opt. Express **7**(10), 3826–3842 (2016). [CrossRef]

**25. **V. Quaresima, M. Ferrari, A. Torricelli, L. Spinelli, A. Pifferi, and R. Cubeddu, “Bilateral prefrontal cortex oxygenation responses to a verbal fluency task: a multichannel time-resolved near-infrared topography study,” J. Biomed. Opt. **10**(1), 011012 (2005). [CrossRef]

**26. **A. Dot, A. Planat-Chrétien, M. Perriollat, M. Berger, R. Lartizien, M. Henry, G. Bettega, and J.-L. Coll, “Blood oxygenation in buried flaps: a bi-layer reconstruction,” Proc. SPIE **11074**, 1107427 (2019). [CrossRef]

**27. **A. V. Cardoso and A. O. Camargos, “Geometrical aspects during formation of compact aggregates of red blood cells,” Mater. Res. **5**(3), 263–268 (2002). [CrossRef]

**28. **J. Ježek, J. Starič, and M. Nemec, “The influence of age, farm, and physiological status on pig hematological profiles,” J. Swine Health Prod. **26**, 72–78 (2018).

**29. **M. S. Irwin, M. S. Thorniley, C. J. Doré, and C. J. Green, “Near infra-red spectroscopy: a non-invasive monitor of perfusion and oxygenation within the microcirculation of limbs and flaps,” Br. J. Plast. Surg. **48**(1), 14–22 (1995). [CrossRef]

**30. **M. S. Thorniley, J. D. S. Sinclair, N. J. Barnett, C. B. Shurey, and C. J. Green, “The use of near-infrared spectroscopy for assessing flap viability during reconstructive surgery,” Br. J. Plast. Surg. **51**(3), 218–226 (1998). [CrossRef]

**31. **Y. Kagaya, N. Ohura, M. Kurita, A. Takushima, and K. Harii, “Examination of tissue oxygen saturation (StO2) changes associated with vascular pedicle occlusion in a rat Island flap model using near-Infrared spectroscopy,” Microsurgery **35**(5), 393–398 (2015). [CrossRef]

**32. **J. J. Selb, J. J. Stott, M. A. Franceschini, A. Gregory Sorensen MD, and D. A. Boas, “Improved sensitivity to cerebral hemodynamics during brain activation with a time-gated optical system: analytical model and experimental validation,” J. Biomed. Opt. **10**(1), 011013 (2005). [CrossRef]

**33. **D. Contini and L. Spinelli, “Novel method for depth-resolved brain functional imaging by time-domain NIRS,” Eur. Conf. on Biomed. Opt. **6629**, 6629_7 (2007). [CrossRef]

**34. **A. Liebert, H. Wabnitz, J. Steinbrink, H. Obrig, M. Möller, R. Macdonald, A. Villringer, and H. Rinneberg, “Time-resolved multidistance near-infrared spectroscopy of the adult head: intracerebral and extracerebral absorption changes from moments of distribution of times of flight of photons,” Appl. Opt. **43**(15), 3037–3047 (2004). [CrossRef]

**35. **A. Liebert, H. Wabnitz, and C. Elster, “Determination of absorption changes from moments of distributions of times of flight of photons: optimization of measurement conditions for a two-layered tissue model,” J. Biomed. Opt. **17**(5), 057005 (2012). [CrossRef]

**36. **G. Giacalone, M. Zanoletti, D. Contini, R. Re, L. Spinelli, L. Roveri, and A. Torricelli, “Cerebral time domain-NIRS: reproducibility analysis, optical properties, hemoglobin species and tissue oxygen saturation in a cohort of adult subjects,” Biomed. Opt. Express **8**(11), 4987–5000 (2017). [CrossRef]