Effect of a thin superficial layer on the estimate of hemodynamic changes in a two-layer medium by time domain NIRS

In order to study hemodynamic changes involved in muscular metabolism by means of time domain fNIRS, we need to discriminate in the measured signal contributions coming from different depths. Muscles are, in fact, typically located under other tissues, e.g. skin and fat. In this paper, we study the possibility to exploit a previously proposed method for analyzing time-resolved fNIRS measurements in a two-layer structure with a thin superficial layer. This method is based on the calculation of the timedependent mean partial pathlengths. We validated it by simulating venous and arterial arm cuff occlusions and then applied it on in vivo measurements. 2015 Optical Society of America OCIS codes: (170.5280) Photon migration; (170.1470) Blood or tissue constituent monitoring; (170.2655) Functional monitoring and imaging; (170.6920) Time-resolved imaging. References and links 1. D. A. Boas, C. E. Elwell, M. Ferrari, G. Taga, “Twenty years of functional near-infrared spectroscopy: introduction for the special issue,” Neuroimage 85(1), 1-5 (2014). 2. M. Ferrari, M. Muthalib, V. Quaresima, “The use of near-infrared spectroscopy in understanding skeletal muscle physiology: recent developments,” Philos Trans A Math Phys Eng Sci. 369(1955), 4577-90 (2011). 3. M. Ferrari, L. Mottola and V Quaresima, “Principles, techniques, and limitations of near infrared spectroscopy,” Can. J. Appl. Physiol. 29(4), 463-487 (2004). 4. F. Scholkmann, S. Kleiser, A. J. Metz, R. Zimmermann, J. M. Pavia, U. Wolf, M. Wolf, “A review on continuous wave functional near-infrared spectroscopy and imaging instrumentation and methodology,” NeuroImage 85(1), 6-27 (2014). 5. O. Pucci, V. and K. St. Lawrence, "Measurement of the optical properties of a two-layer model of the human head using broadband near-infrared spectroscopy," Appl. Opt. 49(32), 6324 (2010). 6. J. T. Elliott, M. Diop, L. B. Morrison, C. D. d'Esterre, T. Y. Lee, K. St. Lawrence, “Quantifying cerebral blood flow in an adult pig ischemia model by a depth-resolved dynamic contrast-enhanced optical method,” 94, 303 (2014). 7. A. Torricelli, D. Contini, A. Pifferi, M. Caffini, R. Re, L. Zucchelli, L. Spinelli, “Time domain functional NIRS imaging for human brain mapping,” NeuroImage 85(1), 28-50(2014). 8. J. Steinbrink, H. Wabnitz, H. Obrig, A. Villringer, and H. Rinneberg, “Determining changes in NIR absorption using a layered model of the human head,” Phys. Med. Biol. 46(3), 879–896 (2001). 9. S. Del Bianco, F. Martelli, and G. Zaccanti, “Penetration depth of light re-emitted by a diffusive medium: theoretical and experimental investigation,” Phys. Med. Biol. 47(23), 4131–4144 (2002). 10. D. Milej, D. Janusek, A. Gerega, S. Wojtkiewicz, P. Sawosz, J. Treszczanowicz, W. Weigl and A. Liebert, “Optimization of the method for assessment of brain perfusion in humans using contrast-enhanced reflectometry: multidistance time resolved measurements,” J. Biomed. Opt., 20(10), 106013 (2015). 11. A. Liebert, H. Wabnitz, J. Steinbrink, H. Obrig, M. Moller, R. Macdonald, A. Villringer and H. Rinneberg, “Time-resolved multidistance nearinfrared 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). 12. L. Hervé, A. Puszka, A. Planat-Chrétien and J.M. Dinten, “Time-domain diffuse optical tomography processing by using the Mellin–Laplace transform,” Appl. Opt. 51(25), 5978–5988 (2012). 13. L. Zucchelli, D. Contini, R. Re, A. Torrcelli and L. Spinelli, “Method for the discrimination of superficial and deep absorption variations by time domain fNIRS,” l4(12), BOE (2013). 14. M. Diop and K. St. Lawrence, "Improving the depth sensitivity of time-resolved measurements by extracting the distribution of times-of-flight," BOE 4(3), 447 (2013). 15. J. R. Mourant, T. Fuselier, T., J. Boyer, T. M. Johnson, I. J. Bigio, “Predictions and measurements of scattering and absorption over broad wavelength ranges in tissue phantoms,” Applied Optics 36(4), 949-957 (1997). 16. M. Niwayama, L. Lin, J. Shao, N. Kudo, and K. Yamamoto, “Quantitative measurement of muscle hemoglobin oxygenation using near-infrared spectroscopy with correction for the influence of a subcutaneous fat layer,” Rev. Sci. Instrum. 71, 4571-4575 (2000). 17. S. Koga, T. J. Barstow, D. Okushima, H.B. Rossiter, N. Kondo, E. Ohmae and D. C. Poole, “Validation of a high-power time resolved near-infrared spectroscopy sytem for measurement of superficial and deep muscle deoxygenation during exercise,” Journal of Applied Physiology, (2015). 18. S. L. Jacques, “Optical properties of biological tissues: a review,” Phys. Med. Biol. 58 (2013). 19. S. Ferrante, D. Contini, L. Spinelli, A. Pedrocchi, A. Torricelli, F. Molteni and G. Ferrigno, “Monitoring muscle metabolic indexes by time-domain near-infrared spectroscopy during knee flexextension induced by functional electrical stimulation,” J. Biomed. Opt. 14(4), 044011-1 (2009). 20. F. Martelli, S. Del Bianco, A. Ismaelli and G. Zaccanti, Light Propagation through Biological Tissue and Other Diffusive Media: Theory, Solutions and Software (SPIE Press, 2010). Chapter 6 pp 109-129. 21. S. H. Hyttel-Sorensen, T. W. Hessel, A. la Cour and G. Greisen, "A comparison between two NIRS oximeters (INVOS, OxyPrem) using measurement on the arm of adults and head of infants after caesarean section," BOE 5(10), 3671-3683 (2014). 22. F. Martelli, S. Del Bianco, L. Spinelli, S. Cavalieri, P. Di Ninni, T. Binzoni, A. Jelzow, R. Macdonald, H. Wabnitz, “Optimal estimation reconstruction of the optical properties of a two-layered tissue phantom from time-resolved single distance measurements,” J. Biomed. Opt. 20(11), 115001 (2015). 23. S. Prahl,, “Optical absorption of hemoglobin”, (2013), http://omlc.ogi.edu/spectra/hemoglobin/index.html 24. M. Muthalib, M. Jubeau, G. Y. Millet, N. A. Maffiuletti, M. Ferrari and K. Nosaka, “Biceps brachii muscle oxygenation in electrical muscle stimulation,” Clin Physiol Funct Imaging 30(5), 360–368 (2010). 25. M. Ferrari, M. Muthalib and V. Quaresima, “The use of near-infrared spectroscopy in understanding skeletal muscle physiology: recent developments,” Phil. Trans. R. Soc. A 369(1955), 4577–4590 (2011). 26. G. A. Tew, A. D. Ruddock and J.M. Saxton, “Skin blood flow differentially affects near-infrared spectroscopyderived measures of muscle oxygen saturation and blood volume at rest and during dynamic leg exercise,” Eur J Appl Physiol 110(5), 1083-1089 (2010). 27. R. C. Mesquita, S. S. Schenkel, D. L. Minkoff, X. Lu, C. G. Favilla, P. M. Vora, D. R. Busch, M. Chandra, J. H. Greenberg, J. A. Detre and A. G. Yodh, “Influence of probe pressure on the diffuse correlation spectroscopy blood flow signal: extra-cerebral contributions,” BOE 4(7), 978-994 (2013). 28. A. Pifferi, D. Contini, L. Spinelli, A. Torricelli, R. Cubeddu, F. Martelli, G. Zaccanti, A. Dalla Mora, A. Tosi, F. Zappa, “The spread matrix: a method to predict the effect of a non time-invariant measurement system,” in Biomedical Optics, OSA Technical Digest (Optical Society of America, 2010), paper BSuD22R. 29. M. S. Patterson, B. Chance and B. C. Wilson, “Time resolved reflectance and trasmittance for the non-invasive measurement of tissue optical properties,” Applied Optics 28(12) (1989). 30. Re, D. Contini, M. Caffini, R. Cubeddu, L. Spinelli and A. Torricelli, “A compact time-resolved system for near infrared spectroscopy based on wavelength space multiplexing,” Rev. Sci. Instrum 81, 113101 (2010).


Introduction
Near infrared spectroscopy (NIRS), is a non-invasive optical technique for monitoring hemodynamic parameters (e.g.oxy-and deoxy-hemoglobin concentration) of biological tissues [1][2][3].In most applications the target tissues are the brain cortex and the muscle.Both targets cannot be directly investigated by NIRS since they are naturally covered by other tissue types: the scalp, the skull, and the cerebral spinal fluid layer shelter the brain cortex, while the skin and the fat layer cover the muscle.To a first approximation, from a physical point of view, these target tissues can be modelled as a layered medium with a superficial layer and a deep layer.Depending on the thickness and on the optical properties of the superficial layer, the estimate of hemodynamic parameters in the deep layer by NIRS measurements can be greatly affected.In the typical configuration for NIRS measurements, optical fibers for light injection and collection are placed on the surface of the superficial medium at a relative distance of about 2 cm or greater.As a result, the NIRS signal receives contributions from photons that have travelled in both layers.There is therefore the need of extracting from the NIRS signal the information regarding the deep layer, discriminating the contribution coming from the superficial one.A multi-distance approach has been proposed to tackle this issue in the continuous wave (CW) [4,5], in the frequency domain (FD) or in the time domain (TD) [6].However, with TD NIRS it is possible to discriminate signals coming from different depths without the need of implementing a multi-distance approach [7][8][9][10].In fact, in TD NIRS there is a direct and intuitive link between the time spent by photons inside the tissue and the probed depth.Therefore, it is possible to qualitatively discriminate superficial and deep changes by selecting early and late time-windows in the recorded distribution of photon time-of-flights (DTOF).However, the accurate quantification of the hemodynamic changes in the different compartments of tissues is still an open issue, and different solutions were proposed [11][12][13][14].In particular, a method, based on the refined computation of the pathlengths traveled by photons in a two-layer medium, was recently proposed by our group [13].This method was applied for modelling the adult head, where the thickness of the upper layer was varied in the range 0.5 -2 cm in 0.5 cm steps.No studies were performed for other cases (e.g.neonatal head) or other tissues (e.g.muscle) in other body compartments (e.g.arm, leg and abdomen).In these cases, we have to take into account that the superficial layer is much thinner than the one we encounter on average during measurements on the adult head.Our aim is to test if the proposed method can be applied with success also when the superficial layer is in the range 0.2-1 cm.
In the first part of the work we simulate a realistic situation where the concentrations of oxy-[O2Hb], deoxy-[HHb], and total-hemoglobin [tHb], and tissue oxygen saturation [SO2] vary as during an arterial or a venous occlusion of the arm muscle.We implemented different conditions in terms of geometry (layer thickness) and instrument set-up (instrument response function) and investigated different ways to analyze data (homogenous or bi-layer model).In the second part, we applied the results on data acquired during in vivo measurements.

Simulations
We modeled the arm as a two-layer structure composed of a superficial (UP) and a deep (DW) layer, respectively (see Fig. 1).The UP layer represents the part of the arm which is not involved in muscle activity (i.e.skin, lipid and surface capillary bed), while the DW layer refers to the muscular tissue.The thickness of the UP layer (SUP) was varied in the range 2-10 mm, in 2 mm steps, while the thickness of the DW layer (SDW) was kept fixed to 100 mm.The refractive index was set at n = 1.40 for both the UP and DW layers.The spectral dependence of the reduced scattering coefficient was modeled applying an approximation to the Mie theory [15] as , with 0 λ = 600 nm, a = 12 cm -1 and b = 0.5 for the UP layer, a = 12 cm -1 and b = 1 for the DW layer.The absorption coefficient was obtained from the concentrations Ci of the tissue constituents by means of the Beer's law, ( ) ( ) , where ( ) i ε λ is the specific absorption of the i-th component.We choose 690 nm and 821 nm as operating wavelengths, being two wavelengths in the range employed by most TD NIRS instruments [7].As main tissue components, we considered HHb, O2Hb, lipid and water.Starting from some studies [16][17][18] and some our previous experience [19] we set the baseline concentrations as shown in Table 1.The values for total-hemoglobin ( ) and tissue oxygen saturation ( )

SO
O Hb tHb = are also shown, since they are the relevant hemodynamic parameters for clinical applications.In Table 1 the corresponding values for the optical properties in the two layers are also reported.
Given the value of the optical parameters at the two wavelengths, the time-resolved reflectance (TRR) curves were simulated , with an interfiber distance ρ = 20 mm and for an  ideal time-resolved set-up with a δ-like instrument response function (IRF).The TRR curves were calculated by computing the solution of the diffusion equation in a two-layer geometry [20] by means of a homemade software.The total number of collected photons was set at 10 8 photons, and Poisson noise was added to the simulated curves to mimic real measurements.We did not implement repetitions since it was already demonstrated by Zucchelli et al. [13] that the error introduced by repetitions is just a negligible stochastic one.In the previously mentioned paper was also already discussed the influence of different counts rates for the TRR curves.
Starting from the baseline values reported in Table 1 different scenarios were simulated to mimic typical hemodynamic changes in the arm.The first scenario simulated is an arterial occlusion (AO): when a pneumatic cuff is inflated around the arm at a pressure higher than 180 mmHg, both the arteries and the veins are closed and the typical response is an increase of HHb and a decrease of O2Hb, with a corresponding decrease of SO2, while tHb stays fixed (apart from a slight increase during the initial cuff inflation phase, given that arteries are closed only after veins occlusion and a limited amount of blood can enter the measured volume before the complete AO starts).If the cuff is then released, a typical reactive hyperemia peak is observed, eventually followed by a recovery to the baseline values.
The second simulated scenario is a venous occlusion (VO) obtained by inflating a pneumatic cuff at a pressure below 100 mmHg that is sufficiently mild to keep the arteries open.In this case both HHb and O2Hb increase (as well as tHb), while SO2 shows only a slight decrease.For both scenarios the amount of the variations were tuned according to literature data [21] For representation purposes we have decided to simulate a time course of hemodynamic changes in which the baseline values were kept fixed for 20 s, the occlusion lasted 60 s and the recovery phase 80 s.The sampling time was 10 s per point, so as to create 2 points in the baseline, 6 points in the occlusion phase and 8 points in the recovery phase.With these settings, we are aware of the fact that the rate of change of the hemodynamic parameters might be distorted with respect to the natural case (e.g. in AO the maximum change in the hemodynamic parameters may require much more than 1 min).Nonetheless, we were not interested in monitoring in details the time course, rather in estimating the maximal errors when large variations occur.In Table 2 and Table 3 the simulated variations of the hemodynamic parameters during AO and VO are reported for both UP and DW layers, respectively.
In order to thoroughly test the performances of the analysis method (described in the next section) we have simulated for the AO (VO) scenarios, the sub-cases AOUP (VOUP) where the changes in the hemodynamic parameters occur in the UP layer only (with DW layer kept constant at the baseline values) and AODW (VODW), where the changes occur in the DW layer only (with UP layer kept constant at the baseline values).While these sub-cases can be far from real physiological situations, they are useful to test the ability of the method in discriminating changes that occur in the UP or DW layer.

Data analysis
The method, we have used to estimate the absorption changes in the UP and DW layers, was described in details in Zucchelli et al. [13].With a custom-made software we calculated the time dependent mean partial pathlength (TMPP), traveled in the j th layer by photons detected at time t, as: where is the unperturbed TRR curve convoluted with the IRF, and 0 a j µ is the absorption coefficient of the baseline period for the j th -layer.This calculation was done for both layers and wavelengths.
The simulated TRR curves were then divided into 10 consecutive temporal gates (time windows) with the same width (200 ps).As starting and ending points, we choose the time where the counts were one hundredth of the peak counts on the leading and trailing edge of the TRR curve.
Once the TMPPs are calculated, we can estimate for each wavelength the absorption changes in each layer ( ) aj µ λ ∆ with respect to the baseline, by inverting the following linear system: where ( ) The estimate for ( ) 0 a j µ λ can be obtained by fitting of the baseline TRR curve ( ) 0 R λ  with a proper model for photon migration (e.g.two-layer, with a priori knowledge of SUP).The performances of this fitting procedure were discussed elsewhere [22].Conversely in this work, a priori knowledge of the values of ( ) 0 a j µ λ is assumed (e.g. in Eq. 1 and 3) because we are interested in testing the method against changes in the hemodynamic parameters only.From the absorption coefficient ( ) aj µ λ , HHb and O2Hb can be obtained by means of the Beer's law [23], assuming a known background absorption due to lipid and water.

Results
Starting from the calculation of the TMMP and following the analysis, as explained in the previous sections, we analyzed the different experiments we have simulated.In the following paragraphs the results are shown.

Arterial occlusion (AO)
In this section we show the results for the simulated AO, obtained assuming a priori knowledge of ( ) 0 a j µ λ .In Figure 2 the simulated and reconstructed time courses of the hemodynamic parameters for UP and DW layers are shown for the different values of SUP.In both the cases, we can first notice that, as expected, the calculated baseline values well fit with the nominal ones for all the SUP.We can observe just some negligible fluctuations around the nominal value, due to the Poisson noise.For what concern the whole time course, we can observe that, qualitatively, the reconstructed values closely match the simulated ones.As expected, for the UP layer, a discrepancy is observed when SUP < 4 mm (i.e. when the corresponding TMMPs are very short), while for the DW layer the larger is SUP (i.e. the smaller are the corresponding TMMPs), the larger are the deviations.In both layers, the deviations are larger when the changes with respect to the baseline are larger as well.This is in accordance with the fact that our model works in a better way for small variations [13].
In order to better quantify the differences between simulated and reconstructed values, the percentage relative error ε(t) for all the hemodynamic parameters, with respect to the nominal simulated value x0(t), was calculated as ε(t) = [x(t) -x0(t)]/x0(t) for both the UP and DW layers.
In order to be able to evaluate the dependence of ε from the change in the hemodynamic parameters and from SUP, Figure 3 reports radar graphs of the relative error as a function of the time (1-160 s) for UP and DW layer.The different colors represent the different SUP.The biggest error (−7 % for HHb) in the UP layer is found when the UP layer thickness is small (2 mm).In this case, the UP layer is so thin that the path traveled by photons inside it is much smaller than the total pathlength, causing an error in the estimation of the optical properties and then of the hemodynamic parameters.Looking at thicker UP layers (6, 8, 10 mm) the error absolute value decreases under the 1 % and can be then considered as negligible.Conversely in the DW layer: when the UP layer is 10 mm, the error is bigger (about 12 % for SO2) because the path traveled by photons in the DW layer is much smaller than the one traveled in the UP layer.Instead, when the UP layer is 2 mm the error absolute value is much smaller (< 2 %).

Arterial occlusion with one layer kept constant
In the AOUP sub-case the hemodynamic changes are applied to the UP layer only, while the baseline values are kept constant in the DW layer.The absolute value of the relative error ε(t) is always smaller than 2% in both layers for all SUP values (see Fig. 4).This means that in case the hemodynamic parameters in the DW layer are constant, the variations in the UP layer can be determined with a good accuracy and that they have no significant influence in the estimation of the DW layer parameters.Conversely, in the AODW sub-case the hemodynamic changes are applied to the DW layer only, while the baseline values are kept constant in the UP layer.The relative error in the two layers are not always negligible (see Fig. 5).In particular, we can observe that the hemodynamic parameters in the DW layer can be calculated with an error, which is bigger if the layer is thicker (max ε(t) around 19% for SO2).Those variations affect also the capability of calculating the UP layer hemodynamic parameters, above all if the UP layer is thinner (max ε(t) around 16% for HHb).Anyway, the errors found in this particular case are higher than the ones found when both the layers are varying.

Venous occlusion
We calculated the hemodynamics parameters for VO in the same way as described in Sec 3.1 for AO.In Figure 6 we show their time courses for respectively the UP and DW layer at the different SUP .All the considerations done for AO concerning the baseline period and the increasing in the error with the increasing of the absolute values of the variations, can be repeated for VO as well.During the occlusion we can notice an underestimation of the nominal values, as happens during the AO.In particular, looking at the radar graphs in figure 7, we can observe that for the UP (DW) layer the maximum ɛ(t) found was about 8% (−15%) for HHb, when SUP is small (thick).

Venous occlusion with one layer kept constant
As done for AO, we considered, also for VO, just the sub-cases where the only UP (VOUP) or DW (VODW) layer is varying.In figures 8 and 9, the relative errors are respectively shown.

Discussion
The results, shown in Secs.3.1 and 3.3, clearly identify the potentiality of the analysis method when applied to tissues, such as muscle, where SUP is thin (<= 10 mm).In fact, if we contextualize these results in real measurements, like monitoring the arm muscle during a voluntary movement or an electrical muscle stimulation, the maximal SO2 variation we expect is around 20% [24].In this case, the relative error is below 3% for both the UP and DW layers and for all SUP values.For all the applications where we have a comparable amount of SO2 variations or lower (e.g.knee flex-extension during electrical stimulations on the calf muscle [19]) the error, which this method introduces, is negligible.Higher variations can be reached, for example, in the biceps brachii muscle during sustained isometric contraction at the maximal voluntary contraction [25].In this case, we can reach the higher value we found for the ε(t) which is around 11% for the DW layer when SUP is thicker, for SO2 parameter.One example of real experiments where it is interesting to study physiological events occurring only in the UP layer is the study of the tissue heating when light is injected in [26], or the study of the influence of the probe pressure on the measurements [27].Also in this case (Secs.3.2 and 3.4) the error is negligible.(below 3%).
On the other hand, when we are really performing measurement and not simulations, we have to take into account a series of non-idealities and unknown parameters that are discussed in the following.

On the effect of the instrument response function
The reference dataset was simulated considering a δ-like temporal response.In this way, the errors we found were referred only to the data analysis method.However, when we perform real measurements we have to consider the instrument response function (IRF).Generally speaking, the broadening of the IRF affects negatively the depth sensitivity and the contrast [28].Zucchelli et al. [13] have already shown that the worsening in the results is due not only to the broadening of IRF (i.e. its full width at high maximum -FWHM), but also to the value of the IRF decay time constant τ.For this reason, we considered two different IRF with different τ and FWHM, and we convoluted them with the theoretical model before fitting the experimental DTOF.The first IRF (IRF1) we considered has τ = 200 ps and FWHM = 700 ps, while the second IRF (IRF2) has τ = 100 ps and FWHM= 160 ps.IRF1 has been obtained with a system employing photomultiplier tubes and 3 mm custom-made glass optical fiber, while IRF2 has been obtained employing a hybrid detector and 1 mm plastic optical fibers.We calculated the percentage relative error ( ) with respect to the δ-like impulse as . Since the introduction of the instrument contribution causes a decreasing in the contrast, an increasing of the errors in the hemodynamic variation calculation is introduced as well.With the increasing of the hemodynamic variation with respect to the baseline values, we have an increasing in the error due to the non-ideality of the pulse shape.If we focus on the SO2 values obtained with the IRF1, which is often the goal of in-vivo measurements on muscle, we obtain a negligible error for all the percentage variations of the UP layer only for the 10 mm SUP case.If the thickness is smaller, we have an increase in ( ) IRF t ε which is negligible only if the hemodynamic variation is low (< 10 %).When SUP is thin (2 mm) the relative error is much higher and can reach values around 90%.For the DW layer, the ( ) IRF t ε found when the variations are high are around 20 % (40 %) when SUP is 2 mm (4 mm).
If we consider the IRF2, we have almost the same results when SUP is thin, but the relative error ( ) IRF t ε decreases with the increasing of SUP.In particular the relative error occurring when SUP = 10 mm at the maximum of the variation is only 6%.This data underline the huge increase in the error due to the non-idealities introduced by the IRF and in particular that the IRF1 cause errors larger than the IRF2 principally because of the larger τ.

On the effect of source detector distance
While in CW measurements the source-detector distance ρ affects the depth sensitivity, in TD measurements it is principally related with the lateral resolution [7], since it is the time that affects depth sensitivity: then TMPPs should not be affected by a change in the source detector distance.This assumption is true when the arm is model as a homogeneous tissue.We wanted to understand if it can be considered true also when the arm is modeled as a bilayer structure and the SUP is small.Therefore, starting from the reference dataset, we changed the sourcedetector distance from 20 to 40 mm.If we look at the UP (DW) layer, the smallest (biggest) SUP does not show significant differences for the two different source-detector distances.On the other hand, in the UP (DW) layer, the biggest (smallest) thickness shows higher value for TMMPs when the source-detector distance is 40 mm.The increase in the TMPPs values can be higher than 70 %.We can then study how these differences affect the retrieve of the hemodynamic parameters.Total hemoglobin is never affected by the source-detector distance.HHb, O2Hb and SO2, show differences only for the UP layer, when the SUP is thin.In this case, the relative errors ( ) , when the amplitude variation is at maximum 9%, 11% and 14% respectively.

On the uncertainty of unknown parameters and fitting procedures
If we perform real measurements, typically we do not have the a priori information about the value of the absorption coefficient during the baseline period: ( ) 0 a µ λ .We have then to perform a fit procedure for the baseline period, in order to extrapolate this value.The finding of the best procedure lies outside the aim of this study, and it was recently addressed in the context of optimal estimation by using a bilayer model [22].Here we just observe that fitting the baseline with a homogeneous model causes not negligible error when the SUP is higher than 2 mm (data not shown).
Another point of interest can be to understand if, once we have calculated ( ) 0 a µ λ , since the SUP was kept very thin in our simulations (2-10 mm thickness), it is possible, in terms of acceptable error, to analyze data with an homogenous model, which can be of course easier to implement because we do not need the SUP as a priori information.We therefore fit the TRR curves, simulated as a bilayer medium with the different UP layer thicknesses, with a semiinfinite homogeneous model [29] as if the arm was composed only of muscular tissue (DW layer).We found that for all the SUP, the error is much smaller if we fit the TRR curves with the new method with respect to the homogeneous model.The error carried out fitting with the homogeneous model is acceptable only when a 2 mm SUP is considered, but anyway is higher than the one we found with the method proposed here.For bigger SUP layer thicknesses, the error is unacceptable not only in terms of maximum variations, but also in terms of baseline value determination.We can then affirm that is not possible to implement a homogeneous fit if the muscle is a bilayer structure even if the UP layer is thin.
The knowledge of the SUP values is then important.During muscle measurement it can be determined with techniques such as ultrasounds or with the employment of much easier to use and cheaper instrument such as skin folder calipers.The problem can be how much confident we need to be with this measurement.We performed some trials on the reference dataset, where we analyzed data referring to SUP= 2mm with 10 mm of thickness and vice-versa.We calculate the percentage relative error x t is the value obtained fitting with SUP equal to 2 mm (case 1) or 10 mm (case 2).If we look at the errors calculated for the SO2 parameter we obtain additional errors that can reach the 33 % (case one) and the 224 % (case 2) in the point of the maximum variation for the DW layer.This discussion underline the necessity to have a really accurate measurement of the SUP value, because small errors (< 2 mm) in its determination cause huge errors in the estimation of the hemodynamic parameters.
Finally, it is important also to take into consideration that some parameters we set during simulations are unknown or not ideal, such as the photons number per second that can be lower and the concentration of lipid and water in both the layers.In particular, the consideration about the number of photons and its relation with errors in the estimation of the hemodynamic time courses because of the worsening of the signal to noise ratio was already discussed by Zucchelli et al. [13].

In vivo measurements
Finally, we performed in vivo measurements in order to test the applicability of our method on muscle during a real experiment.We performed an arterial (180 mmHg) and a venous (100 mmHg) occlusion on the left arm of two different adult male volunteers.The TD instrument employed is described in Re et al. [30 ].The measurement protocol consisted of 1 min of baseline, 2 min of occlusion, and 3 min of recovery.The measurements were performed in a reflectance geometry, with an interfiber distance of 20 mm, with an acquisition rate of 400 ms, (200 ms for each wavelength).The optical fibers were fixed to the arm by means of a black rubber pad connected to brand fasteners (ONE-WRAP ® , Velcro Italia Srl, Italy).SUP was measured with a skin folder caliper and set at 2 mm for the subject undergoing the arterial occlusion and at 5 mm for the one undergoing the venous occlusion.Data were analyzed as described in Secs.2.2 and 4.3 with the two-layer model.
In figure 10, the time courses for the absolute values of O2Hb, HHb, SO2 and tHb, for the DW layer of the left arm, during the arterial and a venous occlusion, respectively, are presented.During the arterial occlusion, the deoxygenated (oxygenated) blood through the veins (arteries) cannot go outside (inside) the occluded part, so that we have an increase of the HHb and a decrease of the O2Hb during the entire occlusion task, as shown in figure 10 left panel, for the DW layer.When the cuff is deflated, we observe the typical reactive hyperemia peak and then a slow return to the baseline values.We notice a delay in the response of HHb compared to O2Hb, as expected, since the veins are released after the arteries, this effect being enhanced by the slow release of the cuff.No significant changes were measured in the right arm.In figure 10 right panel, we can notice the typical time course of a venous occlusion: the oxygenated blood can flux inside the occluded part through the arteries (increase in O2Hb) but the deoxygenated one cannot go outside them (increase in HHb).

Conclusion
In this paper, we have shown with success the application of an analysis method for TD NIRS on muscle tissues.This method, based on the calculation of the TMMPs exploiting a bilayer model for the photon propagation in turbid media, allows to decouple the path traveled by photons in the most superficial part of the tissue from the path spent in the deeper part, where the muscular metabolism is located.In this way, it is possible to calculate the concentrations of the hemodynamics parameters in both layers.This method was at first tested on a simulated dataset and then applied with success on in vivo measurements.The parameters' estimation of the DW layer is always well achieve in particular when the UP layer is small (2mm).Conversely, the estimation in the UP layer can be less accurate if the UP layer is small.
are the intensities in the time gate g in the reference TRR curve (average of the values over the baseline period) and in a generic TRR curve, respectively.The absorption changes from the baseline values to obtain the absolute estimate for the absorption coefficient in the two layers

Fig. 2 .Fig. 3 .
Fig. 2. AO simulated and reconstructed hemodynamic parameters in the UP (a) and DW (b) layer.The different colors represent the SUP and the initial simulated value.

Fig. 4 .
Fig. 4. Relative error with respect to the baseline of the hemodynamic parameters during AOUP for the UP and DW layer calculated for different SUP (different colors).

Figure 5 :
Figure 5: Relative error with respect to the baseline of the hemodynamic parameters during AODW for the UP and DW layer calculated for different SUP (different colors).

Figure 6 :
Figure 6: VO Simulated and reconstructed hemodynamic parameters in the UP (a) and DW (b) layer.The different colors represent the SUP and the initial simulated value.

Figure 7 :
Figure 7: Relative error with respect to the baseline of the hemodynamic parameters during VO for the UP and DW layer calculated for different SUP (different colors).

Figure 8 :
Figure 8: Relative error with respect to the baseline of the hemodynamic parameters during VOUP for the UP and DW layer calculated for different SUP (different colors).

Figure 9 :
Figure 9: Relative error with respect to the baseline of the hemodynamic parameters during VODW for the UP and DW layer calculated for different SUP (different colors).

Figure 10 .
Figure 10.In-vivo AO (left panel) and VO (right panel) hemodynamics variations during an arterial cuff occlusion of the left arm.Results refer to DW layer.