## Abstract

Time-resolved near infrared spectroscopy is considered to be a gold standard technique when measuring absolute values of tissue optical properties, as it provides separable and independent information about both tissue absorption and scattering. However, time-resolved instruments require an accurate characterization by measuring the instrument response function in order to decouple the contribution of the instrument itself from the measurement. In this work, a new approach to the methodology of analysing time-resolved data is presented where the influence of instrument response function is eliminated from the data and a self-calibrating analysis is proposed. The proposed methodology requires an instrument to provide at least two wavelengths and allows spectral parameters recovery (optical properties or constituents concentrations and reduced scatter amplitude and power). Phantom and in-vivo data from two different time-resolved systems are used to validate the accuracy of the proposed self-calibrating approach, demonstrating that parameters recovery compared to the conventional curve fitting approach is within 10% and benefits from introducing a spectral constraint to the reconstruction problem. It is shown that a multi-wavelength time-resolved data can be used for parameters recovery directly without prior calibration (instrument response function measurement).

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

## 1. Introduction

The time-resolved (TR) near-infrared spectroscopy (NIRS) is a diagnostic technique [1–3] that uses near-infrared light (600-1000nm) in form of short (femto/pico-seconds in width) and low energy pulses at repetition rates of tens of MHz. Light pulses are delivered to the tissue surface and a TR detection/acquisition electronic chain, often based on photon counting devices measures the time of flight of single photons originating at source pulses and builds a distribution of time of flight of photons (DTOF). The measured DTOF shape depends on both path-lengths (related directly to photon travel time within tissue) and absorption experienced by photons travelling from the source to the detector.

The time-resolved spectroscopy (TRS) has widely become and is considered a gold standard in measurement of tissue optical properties [4,5] as it is able to provide both absorption and reduced scattering coefficient of the tissue being imaged. When used at multiple wavelengths, the spectral information of the absorption provides absolute concentrations of the main tissue constituents such as oxygenated and reduced haemoglobin, water, lipid and in some instances cytochrome c oxidase [6,7]. The spectrally resolved reduced scattering can additionally provide information about the scatter size and density of the medium. Thus, TRS has found many *in-vivo* applications e.g. confirmation of brain death [8], intraoperative brain monitoring during endarterectomy [9], traumatic brain injury monitoring [10], optical mammography [11] or functional brain topography [2]. The possibility to measure the reduced scattering is of high importance on studies combining the TRS and the diffuse correlation spectroscopy [12–14] where the knowledge of the optical properties is required to recover reliable information regarding blood flow/perfusion [13].

The DTOF as measured by a TR instrument is affected by the properties of the system itself, since the measured distribution represents convolution of both tissue and the instrument response to an ideal (infinitely-short) light pulse. The instrument response function (IRF) is often measured independently and used directly within parameters recovery algorithms. The TR NIRS is mainly based on a recovery algorithm which optimises fitting of the theoretical (e.g. semi-infinite medium based) temporal point spread function (TPSF) that is convolved with the IRF to the measured DTOF in order to recover bulk tissue properties [15]. This is the typical approach as the de-convolution of the IRF from the DTOF is undesirable as it is known to amplify the noise within the signal [16].

Multiple channel DTOFs can be also used in a tomographic recovery approach, where the measured distributions are parametrized with photon travel path sensitive parameters such as total intensity, statistical moments [17], time of arrival windows [18] or Mellin-Laplace transform [19,20]. Using this approach, a numerical model of the problem is used to calculate the Jacobian (spatial sensitivity distributions) allowing the mapping of the changes in data measured on the surface to the change of absorption and scattering properties within the discretized tissue model. The utilisation of a gradient based method, such as those involving Jacobians, allow solving of the inverse problem to recover the unknown spatial distributions of the tissue optical properties, and as such must also consider the IRF. As such, in the moments based methods, statistical moments of the DTOF and IRF are subtracted [17], whereas the time windows approach requires convolving of the TPSF with the IRF [18] and the Mellin-Laplace reconstruction uses measurements on a reference optical phantom with known properties to calibrate the tissue data [19].

TR tomography (functional and absolute imaging) and curve fitting methodologies that are not sensitive to the IRF are highly desired, as it will minimise the propagation of error throughout the analysis and parameter recovery. As proposed in [21], bulk tissue optical properties can be recovered without accounting for the IRF using a subtraction time-resolved method where source-detector distance derivative of the statistical moments of a DTOF are sufficient to allow fast and direct calculation of tissue absorption and scatter. A system with a carefully designed source-detector configuration is therefore capable of measuring the derivatives of the moments directly. In systems in which the IRF shape does not change significantly with respect to wavelength, the IRF can be considered as spectrally invariant. Based on this, a new approach is proposed to allow “self-calibrating” DTOFs that are measured at multiple (*N*_{λ} ≥ 2) wavelengths (λ), and utilising Fourier Transformation of TR data to the frequency domain [22]. As the DTOF frequency components can be used directly into well-established frequency domain (FD) based reconstruction [23], this will allow recovery of both absorption and reduced scattering spatial distributions simultaneously. In this work, it is demonstrated that the normalization of the frequency components of TR data by e.g. the first frequency and utilising a spectral derivative approach provides information that are independent on the system characteristics and hence the IRF.

In this work, a new method is presented that is able to self-calibrate the TR data, outlining the utilization of the spectral data within the frequency domain. Comparisons to the curve fitting algorithm show small differences in the recovered parameters between the two methods, albeit removing the need of using the IRF in the parameters recovery process.

## 2. Methodology

This section covers the theory of self-calibrating time-resolved spectroscopy and outlines how self-calibrated data can be used in parameter recovery algorithms. Please refer to [22] for methodology of converting a DTOF into the frequency domain and vice versa. The following notations will be used: *IRF* – measured instrument response function, *DTOF* – measured distribution of time of flight of photons, *TPSF* – theoretical tissue response to the Dirac delta function (temporal point spread function), λ – wavelength, *ω* – angular frequency (*ω* = 2π*f*), *f* – frequency. The IRFs as measured for two spectral TR instruments are used throughout this work [24,25].

#### 2.1 Theory

A key requirement of the proposed method is that the IRFs at multiple wavelengths, for a given system, have the same shape and characteristics; that is the system response function is spectrally independent. As such the time resolved response functions, as transformed into the frequency domain, will have the same frequency components contributing to the signal, up to some frequencies as determined by the noise level. These properties are observed on data from instruments developed at Politecnico di Milano (POLIMI) and Nalecz Institute of Biocybernetics and Biomedical Engineering of the Polish Academy of Sciences (IBIB) [24,25] which are shown in Fig. 1. As evident in Fig. 1(c-f), it is shown that for both systems, the IRF frequency components become noisy and exhibit spectral variation at around 10 GHz.

One condition that should be fulfilled is that the relative positions of IRF maxima is needed to be known a-priori to correct the phase shift accordingly (i.e. align the maxima in time). Otherwise, the phase shift of frequency components as shown in Fig. 1(e-f) will not be consistent. Both of the presented TR instruments utilised in this work, differ by the way which the multi-wavelength measurement is implemented. The POLIMI instrumentation [24] exploits the information and characteristic of the light source, where white laser pulse is first filtered and conventional time-resolved detection is used sequentially on these filtered wavelengths. The IBIB instrument [25] also uses white super-continuum pulse laser but wavelength selection is implemented at the detection side: data at all wavelengths is detected in parallel through a polychromator and multi-channel time-resolved detection system.

To ensure that the IRFs are aligned at maximum values in time, either appropriate compensation on the hardware side can be implemented or by ensuring that the measured DTOFs are shifted in time in post-processing to align the IRFs at maxima. As shown in Fig. 1(c-f) an IRF is wavelength independent and does not undergo phase-wrapping for a wide range of frequencies up to ~10 GHz and ~4 GHz for POLIMI and IBIB instruments respectively. Therefore, within the given frequency range it is assumed that:

*k*(λ) represents hardware dependent amplitude calibration factor and ∗ is the convolution function. This amplitude calibration factor

*k*(λ) can thus be cancelled out by normalizing the DTOF in the frequency domain as follows:

_{N}(i.e. by

*DTOF*(

*ω*,λ

_{N})/

*DTOF*(

*ω*

_{N},λ

_{N})) cancels out the IRFs and leads to:

_{N}is a chosen wavelength, e.g. the shortest wavelength as used in this work.

Equation (4) compares normalized measured data on the left-hand-side with normalized modelled data on the right-hand-side directly. Both sides of the equation are thus normalized in both frequency and wavelength domain and as such provide absolute information about the changes in amplitude and phase with respect to the values at *ω*_{N} and λ_{N}.

#### 2.2 Frequency choice for analysis

The phase shift of the frequency components should not suffer from phase wrapping, i.e. the values are limited between 0 and 2π, otherwise parameter recovery using frequency components of time-resolved data will become non-unique. Considering a source component modulated at some frequency *f*, the detected component of the light signal will not undergo phase wrapping only if the time of flight between the source and the detector is shorter than the modulation period 1/*f*. This time of flight is known as the DTOF, which is always measured, and its frequency components can be analysed. Therefore, the use of the DTOF frequency components up to the wrapping frequency limit is suitable. This frequency wrapping limit will decrease with the time-resolved ‘curve width’ as illustrated in Fig. 2. Moreover, the phase wrapping will be observed at lower frequencies for the DTOF, as compared to the TPSF (since DTOF is normally broader due to the contamination from IRF) and is at its highest for the IRF. Thus, the phase wrapping limit is calculated for the DTOF (left-hand-side of Eq. (4)) and is guaranteed with no phase wrapping for the TPSF (right-hand-side of Eq. (4)) and as such makes Eq. (4) unique in terms of separating the absorption and scattering. TPSFs as show in Fig. 2 have been calculated with NIRFAST [23] for a uniform slab (6 × 6 × 6 cm, 0.23 ± 0.06 mm^{3} tetrahedral element volume) with following optical properties: *μ*_{a} = 0.01 mm^{−1}, *μ*′_{s} = 1 mm^{−1}, *n* = 1.33. Two source-detector separations *r*_{sd}: 15 mm and 30 mm are used to change the TPSF width. The time shift is calculated from the phase shift of the TPSF frequency components.

#### 2.3 Parameters recovery

The normalised spectral time-resolved data in Eq. (4) provides a complex-valued function in **R**^{2} domain (frequency and wavelength). In this work, this function is transformed to two data surfaces (natural logarithm of amplitude and angle of complex values). Thus, the parameters recovery is now a real-valued surface fitting problem discretized at each measured wavelengths and frequencies spanned from the lowest frequency component (greater than 0) up to the phase wrapping frequency as calculated from the measured DTOFs. As such, the recovery problem is now constrained in both wavelength and frequency. The proposed parameters recovery procedure is summarized in Algorithm 1.

The nonlinear fitting problem can be solved for biological tissue using spectrally constrained algorithms [26] to provide parameters such as haemoglobin concentrations, water content and reduced scatter amplitude and power or absorption and reduced scatter at all wavelengths [27].

## 3. Results

In this section, the proposed method is tested against the well-established curve fitting procedure [27] which is conventionally applied separately at individual wavelengths. As for a proof of concept, time-resolved curves (TPSFs) are generated using an analytical solution of the diffusion equation for a semi-infinite model, which are then convolved with measured experimental IRFs from POLIMI as shown in Fig. 1(a-b) and Poisson noise is additionally added considering 10^{7} photon counts per curve. The input data are calculated for a semi-infinite homogenous medium with following properties: oxygenated haemoglobin concentration *C*_{HbO2} = 54.93 μM, reduced haemoglobin concentration *C*_{Hb} = 13.97 μM, water content *W* = 78%, reduced scattering amplitude *S*_{a} = 0.6542 and power *S*_{p} = 0.9260 and the refractive index of *n* = 1.4. The exact partial-flux boundary condition is used with the internal reflectance using the Groenhuis approximation. Haemoglobin extinction coefficients and water absorption spectra are used as integrated in NIRFAST [23] which uses data available at [28]. Data as calculated at 30 mm source-detector distance are shown in Fig. 3. The simulated DTOFs are shifted in time as required to align IRFs maxima.

For the conventional method, the input data (Fig. 3(a)) are analysed in time domain using curve fitting and for extracted frequency domain data (Fig. 3(b-c)) using surface fitting. The curve fitting procedure convolves the IRFs with theoretical TPSFs and calculates the norm of fitting to the input data as limited on the rising and falling edges by 80% and 1% of the maximum value. Curves are fitted at each wavelength separately to get the wavelength dependent absorption *μ*_{a}(λ) and reduced scattering *μ*′_{s}(λ). These spectrally varying optical properties are then used to calculate haemoglobin concentrations (*C*_{HbO2}, *C*_{Hb}), water content (*W*) and scattering amplitude (*S*_{a}) and power (*S*_{p}). Whereas for the proposed self-calibrating method, FD surface fitting is carried out directly for the chromophore concentrations and scattering amplitude and power. No information on IRFs is used in the proposed approach.

Nonlinear curve and surface fitting algorithms are utilised (interior-reflective Newton method) as implemented in MATLAB. In the time domain, the algorithm always starts at optical properties estimated using statistical moments approach [29] and in case of the frequency domain recovery the initial values are: *C*_{HbO2} = *C*_{Hb} = 0.1 mM, *W* = 99% and dimensionless *S*_{a} = *S*_{p} = 1.5. The error function and step tolerances are set to 10^{−6}. Additionally, the algorithm is constrained using the following *min* and *max* values: 0 ≤□*μ*_{a} ≤0.1 mm^{−1}, 0.05 ≤□*μ*′_{s} ≤3 mm^{−1} and for {*C*_{HbO2}; *C*_{Hb}; *W*; *S*_{a}; *S*_{p}} *min* ∈ {0; 0; 0.01; 0.01; 0.01} and *max* ∈ {20; 20; 0.99; 2; 2}. Parameters recovered using both methods are shown in Fig. 4. The Poisson added noise to the data causes variation in the recovered parameters, thus, we show statistics of recovered values based on 30 runs of data generation and recovery procedures.

As shown in Fig. 4 both recovery methods give highly comparable results. As the recovery in time domain is carried out separately in wavelengths, error in the recovered optical properties is observed, whereas this is better constrained in the new method.

#### Phantom data

Fourier transform of a measured time-resolved curve (DTOF) requires the whole curve, including early and late photons. A DTOF undergoing the transform will therefore start and end in noise and due to early photon migration, the diffusion approximation for the generation of the data on the right-hand-side of Eq. (4) is not suitable. As such, for all experimental and *in-vivo* data an open-source Monte Carlo code is utilized [30].

Application of Eq. (4) was tested on a physical phantom as described in [31] where a mechanically switchable phantom allows a change in its global absorption as seen by a pair of source and detector fibres. Data were measured using the POLIMI time-resolved system at a given scattering and absorption of the phantom and a 15 mm source-detector distance. The same algorithms as described above are utilised except that the FD surface fitting is carried out to only recover spectrally varying optical properties. The time domain curve-fitting recovery process still uses the diffusion approximation (since early photons are not considered, i.e. photons on the rising DTOF edge below 80% of the maximum are discarded), whereas the frequency domain recovery uses the whole DTOF and the Monte Carlo method as required by Eq. (4).

As shown in Fig. 5, recovery in time domain using IRFs and in frequency domain (not using the IRFs) gives comparable results with the observed difference being within 10%.

#### In-vivo data

The IBIB time-resolved system is used to measure data on a healthy volunteer during an arm occlusion. The optodes were separated by 30 mm and fixed on a forearm where the cuff placed on the arm was used to change the blood flow. The experiment is composed of three blocks: 0.5 min of rest, 5 min of occlusion (systolic pressure + 50 mmHg) and 1.5 min of rest.

The instrument measured DTOF curves at a sampling rate of 0.3 s. However, since Monte Carlo is used in the recovery loop, only the data at every 30 s is used to minimise the computational burden. To increase the signal to noise ratio the data was post processed using the moving average procedure considering 3 curves (0.9 s) window. The instrument provides measurements at 16 wavelengths from 680 nm to 867.5 nm with the interval of 12.5 nm. Two last wavelengths were discarded as having low signal to noise ratio. The same algorithms settings and recovery approaches as for the phantom experiment are used and the recovery comparisons are shown in Fig. 6.

As shown in Fig. 6 the self-calibrating recovery has successfully been applied to the analyse the *in-vivo* data. Recovered parameters are comparable and follow the same trends. The cuff pressure of systolic pressure + 50 mmHg (usually 170-180 mmHg) has blocked the flow completely resulting in decreased StO2 (tissue oxygenation index) while the *C*_{Hbtot} (total haemoglobin content) has remained constant.

## 4. Discussion and conclusions

It is shown that data from multi-wavelength TR instruments, where the IRF can be considered spectrally invariant, can be utilised as self-calibrating data without the need for measuring and accounting for IRFs. This opens a new way of using the TR instrumentation as multi-frequency high-density diffuse optical tomography devices. Frequency components available for the analysis (up to few GHz) travel at significantly different optical path-lengths (that is, the penetration depth varies), which is a highly desired property to allow tomographic parameter recovery. Moreover, number of available frequency components (typically up to 10 frequencies) is usually greater than the number of data parameters available using other TR based tomographic recovery methods. Additionally, analysing DTOF frequency components up to the phase-wrapping limit supports uniqueness in terms of separating absorption and scattering properties as shown in Fig. 2.

The proposed surface fitting in the frequency domain introduces spectral and frequency constraints regardless of the recovered parameters: optical properties or constituents’ concentrations. Further, normalization in both frequency and amplitude does not require preserving measured DTOFs amplitudes between wavelengths as required by e.g. spectral fitting in time domain [32]. Such instrument calibration would be difficult for the systems as used in the current research. The absorption and scattering spectra and relations between frequencies appropriately constrain the fitting of data as in Fig. 3(b-c). The normalized fitting surfaces are always spanned in frequency and wavelength dimensions regardless of the required recovery parameters. Hence, the parameters recovery should benefit from these constraints as compared to the curve fitting on separate wavelengths as shown in Fig. 4(c). However, for *in-vivo* data as shown in Fig. 6 a variability in the difference of recovered values is observed. More experiments should be carried out to better understand the underlying mechanism responsible for this variation.

Recovery of optical properties introduces number of unknowns equal to two times the number of wavelengths, which in almost all cases will be greater than recovering directly for tissue constituents and scatter amplitude and power. Therefore, spectrally constrained recovery benefits from the fewer degrees of freedom, as is shown, to provide more accurate parameter recovery using FD data [26]. However, the spectral constraint requires a priori knowledge on presence of specific tissue constituents and the extinction spectra of the constituents should be unique within the wavelength frame used.

As shown, Eq. (4) requires exact light transport models and as such, the diffusion approximation cannot be used when analysing measured *in-vivo* data. One option as used in this work is to utilise Monte Carlo, e.g. a GPU-accelerated version [30]. However, one might consider using the direct radiative transport equation (RTE) solutions in semi-infinite space [33] a layered medium [34] or the RTE empirical approximations in time-domain [35] or higher order spherical harmonics approximation [36] in the frequency domain as available e.g. in the FEM-based NIRFAST package [23].

The assumption of the IRF independency of wavelength appears strong. However, lasers can generate different pulse shape (the laser IRF). Additionally, an IRF strongly depends on properties of detecting fibres/fibre bundles (if used) as well as their varying length and numerical aperture (NA) [37]. It may be also dependent on properties of the light coupling optics, which influences the effective NA of the setup. Thus, in case of a time-resolved, multi-wavelength approach in which multiple lasers, independent optical systems and/or fibres are used the method proposed should be used with care. The method can be of use especially when a broadband light source is used for which the IRF is relatively independent on wavelength and in which a single set of fibres/fibre bundles is used for light delivery and detection as in [7,24,25]. However, the usable bandwidth relates directly to the FWHM of a system’s IRF which in turn benefits from a shorter response, and is observed for both tested parameter recovery approaches (proposed one and the curve fitting). It is important to note that both the FWHM of IRF and the time resolution of a system (which is directly related to the detection system), will limit the method. Although there should be very little effect regarding FWHM of the IRF, as long as there is enough temporal sampling within this FWHM, the proposed algorithms are valid.

The relative position of maxima of the IRFs is needed to be known and as such, basic characterisation of the IRFs will be required. However, the time shifts between wavelengths can be fixed in a well-characterised system and the determination of the relative positions of maxima is less challenging as compared to the ‘standard’ IRF collection procedure. Furthermore, the relative position between the IRF and DTOF is no longer an issue as compared to curve fitting.

The coefficient of variation of *IRF*(*ω*,λ) as shown in Fig. 1(c-f) is <3% in amplitude and <1% in phase within the usable frequency range. This translates into <5% variability in recovered haemoglobins concentrations, <10% in water content, <1% in scattering amplitude and <7.5% in scattering power. Comparable variation can also be observed by changing the wavelength used for normalization, λ_{N} in Eq. (4). The same range of variability has also been observed in Fig. 4b where the effects of the added Poisson noise were studied. Therefore, it can be argued that for spectral systems, as used in this work, the effects of variability in IRFs is within the expected noise level in a practical setting.

The parameter recovery error using this algorithm can increase to approximately 15% for very low absorption cases (e.g. *μ*_{a} = 0.0006 mm^{−1} at 730 nm). This is as expected since the low absorption broadens the DTOF, which in consequence limits number of the usable frequencies. Scatter has a minimal effect on the parameters recovery, as changing the reduced scattering coefficient primarily shifts the DTOF left or right on the time axis. Similar trend can also be observed in the curve fitting.

The proposed method can be used almost directly in the frequency domain tomographic approach [23]. Equation (4) transforms raw, IRF-contaminated time-resolved curves into normalized frequency domain data at multiple frequencies and wavelengths. This removes the requirement of using the IRF for data calibration priori to parameter recovery and should increase the fidelity of quantitative imaging using TR NIRS data.

## Funding

Horizon 2020 Framework Programme (688303,675332); Fundació CELLEX (SEV-2015-0522); “la Caixa” Foundation (2014SGR-1555).

## Acknowledgments

This work has received funding from the European Union‘s Horizon 2020 research and innovation programme under grant agreement No 688303 (LUCA project which is an initiative of the Photonics Public Private Partnership), as well as Marie Sklodowska-Curie Innovative Training Networks (ITN-ETN) programme, under grant agreement no. 675332 (BitMap), Fundació CELLEX Barcelona, the “Severo Ochoa” Programme for Centres of Excellence in R&D (SEV-2015-0522), the Obra social “la Caixa” Foundation (LlumMedBcn) and AGAUR-Generalitat (2014SGR-1555)

## Disclosures

The authors declare that there are no conflicts of interest related to this article. LUCA project involves industrial collaboration and, as such, potential conflicts of interest are being monitored by relevant institutional bodies. None has been identified to date.

## References

**1. **A. Pifferi, D. Contini, A. D. Mora, A. Farina, L. Spinelli, and A. Torricelli, “New frontiers in time-domain diffuse optics, a review,” J. Biomed. Opt. **21**(9), 091310 (2016). [CrossRef] [PubMed]

**2. **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**(Pt 1), 28–50 (2014). [CrossRef] [PubMed]

**3. **H. Wabnitz, M. Moeller, A. Liebert, H. Obrig, J. Steinbrink, and R. Macdonald, “Time-resolved near-infrared spectroscopy and imaging of the adult human brain,” Adv. Exp. Med. Biol. **662**, 143–148 (2010). [CrossRef] [PubMed]

**4. **S. Konugolu Venkata Sekar, A. Farina, A. Dalla Mora, C. Lindner, M. Pagliazzi, M. Mora, G. Aranda, H. Dehghani, T. Durduran, P. Taroni, and A. Pifferi, “Broadband (550-1350 nm) diffuse optical characterization of thyroid chromophores,” Sci. Rep. **8**(1), 10015 (2018). [CrossRef] [PubMed]

**5. **A. Farina, A. Torricelli, I. Bargigia, L. Spinelli, R. Cubeddu, F. Foschum, M. Jäger, E. Simon, O. Fugger, A. Kienle, F. Martelli, P. Di Ninni, G. Zaccanti, D. Milej, P. Sawosz, M. Kacprzak, A. Liebert, and A. Pifferi, “In-vivo multilaboratory investigation of the optical properties of the human head,” Biomed. Opt. Express **6**(7), 2609–2623 (2015). [CrossRef] [PubMed]

**6. **A. Rajaram, G. Bale, M. Kewin, L. B. Morrison, I. Tachtsidis, K. St Lawrence, and M. Diop, “Simultaneous monitoring of cerebral perfusion and cytochrome c oxidase by combining broadband near-infrared spectroscopy and diffuse correlation spectroscopy,” Biomed. Opt. Express **9**(6), 2588–2603 (2018). [CrossRef] [PubMed]

**7. **F. Lange, L. Dunne, L. Hale, and I. Tachtsidis, “MAESTROS: A Multiwavelength Time-Domain NIRS System to Monitor Changes in Oxygenation and Oxidation State of Cytochrome-C-Oxidase,” IEEE J. Sel. Top. Quantum Electron. **25**(1), 7100312 (2018). [PubMed]

**8. **W. Weigl, D. Milej, A. Gerega, B. Toczyłowska, P. Sawosz, M. Kacprzak, D. Janusek, S. Wojtkiewicz, R. Maniewski, and A. Liebert, “Confirmation of brain death using optical methods based on tracking of an optical contrast agent: assessment of diagnostic feasibility,” Sci. Rep. **8**(1), 7332 (2018). [CrossRef] [PubMed]

**9. **M. Kacprzak, A. Liebert, W. Staszkiewicz, A. Gabrusiewicz, P. Sawosz, G. Madycki, and R. Maniewski, “Application of a time-resolved optical brain imager for monitoring cerebral oxygenation during carotid surgery,” J. Biomed. Opt. **17**(1), 016002 (2012). [CrossRef] [PubMed]

**10. **W. Weigl, D. Milej, D. Janusek, S. Wojtkiewicz, P. Sawosz, M. Kacprzak, A. Gerega, R. Maniewski, and A. Liebert, “Application of optical methods in the monitoring of traumatic brain injury: A review,” J. Cereb. Blood Flow Metab. **36**(11), 1825–1843 (2016). [CrossRef] [PubMed]

**11. **D. Grosenick, H. Rinneberg, R. Cubeddu, and P. Taroni, “Review of optical breast imaging and spectroscopy,” J. Biomed. Opt. **21**(9), 091311 (2016). [CrossRef] [PubMed]

**12. **T. Durduran and A. G. Yodh, “Diffuse correlation spectroscopy for non-invasive, micro-vascular cerebral blood flow measurement,” Neuroimage **85**(Pt 1), 51–63 (2014). [CrossRef] [PubMed]

**13. **C. Lindner, M. Mora, P. Farzam, M. Squarcia, J. Johansson, U. M. Weigel, I. Halperin, F. A. Hanzu, and T. Durduran, “Diffuse Optical Characterization of the Healthy Human Thyroid Tissue and Two Pathological Case Studies,” PLoS One **11**(1), e0147851 (2016). [CrossRef] [PubMed]

**14. **M. Pagliazzi, S. K. V. Sekar, L. Colombo, E. Martinenghi, J. Minnema, R. Erdmann, D. Contini, A. D. Mora, A. Torricelli, A. Pifferi, and T. Durduran, “Time domain diffuse correlation spectroscopy with a high coherence pulsed source: *in vivo* and phantom results,” Biomed. Opt. Express **8**(11), 5311–5325 (2017). [CrossRef] [PubMed]

**15. **J. A. Guggenheim, I. Bargigia, A. Farina, A. Pifferi, and H. Dehghani, “Time resolved diffuse optical spectroscopy with geometrically accurate models for bulk parameter recovery,” Biomed. Opt. Express **7**(9), 3784–3794 (2016). [CrossRef] [PubMed]

**16. **M. Diop and K. St Lawrence, “Deconvolution method for recovering the photon time-of-flight distribution from time-resolved measurements,” Opt. Lett. **37**(12), 2358–2360 (2012). [CrossRef] [PubMed]

**17. **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] [PubMed]

**18. **Q. Zhu, H. Dehghani, K. M. Tichauer, R. W. Holt, K. Vishwanath, F. Leblond, and B. W. Pogue, “A three-dimensional finite element model and image reconstruction algorithm for time-domain fluorescence imaging in highly scattering media,” Phys. Med. Biol. **56**(23), 7419–7434 (2011). [CrossRef] [PubMed]

**19. **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 transform for experimental detection and depth localization of a single absorbing inclusion,” Biomed. Opt. Express **4**(4), 569–583 (2013). [CrossRef] [PubMed]

**20. **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–4787 (2017). [CrossRef] [PubMed]

**21. **D. Milej, A. Abdalmalak, D. Janusek, M. Diop, A. Liebert, and K. St Lawrence, “Time-resolved subtraction method for measuring optical properties of turbid media,” Appl. Opt. **55**(7), 1507–1513 (2016). [CrossRef] [PubMed]

**22. **S. Wojtkiewicz, T. Durduran, and H. Dehghani, “Time-resolved near infrared light propagation using frequency domain superposition,” Biomed. Opt. Express **9**(1), 41–54 (2018). [CrossRef] [PubMed]

**23. **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. Methods Eng. **25**(6), 711–732 (2009). [CrossRef] [PubMed]

**24. **S. K. V. Sekar, A. D. Mora, I. Bargigia, E. Martinenghi, C. Lindner, P. Farzam, M. Pagliazzi, T. Durduran, P. Taroni, A. Pifferi, and A. Farina, “Broadband (600–1350 nm) Time-Resolved Diffuse Optical Spectrometer for Clinical Use,” IEEE J. Sel. Top. Quantum Electron. **22**(3), 406–414 (2016). [CrossRef]

**25. **A. Gerega, D. Milej, W. Weigl, M. Kacprzak, and A. Liebert, “Multiwavelength time-resolved near-infrared spectroscopy of the adult head: assessment of intracerebral and extracerebral absorption changes,” Biomed. Opt. Express **9**(7), 2974–2993 (2018). [CrossRef] [PubMed]

**26. **S. Srinivasan, B. W. Pogue, S. Jiang, H. Dehghani, and K. D. Paulsen, “Spectrally constrained chromophore and scattering near-infrared tomography provides quantitative and robust reconstruction,” Appl. Opt. **44**(10), 1858–1869 (2005). [CrossRef] [PubMed]

**27. **A. Pifferi, P. Taroni, G. Valentini, and S. Andersson-Engels, “Real-time method for fitting time-resolved reflectance and transmittance measurements with a Monte Carlo model,” Appl. Opt. **37**(13), 2774–2780 (1998). [CrossRef] [PubMed]

**28. **S. Prahl, “Optical Absorption of Hemoglobin”, retrieved https://omlc.org/spectra/hemoglobin/.

**29. **A. Liebert, H. Wabnitz, D. Grosenick, M. Möller, R. Macdonald, and H. Rinneberg, “Evaluation of optical properties of highly scattering media by moments of distributions of times of flight of photons,” Appl. Opt. **42**(28), 5785–5792 (2003). [CrossRef] [PubMed]

**30. **Q. Fang and D. A. Boas, “Monte Carlo Simulation of Photon Migration in 3D Turbid Media Accelerated by Graphics Processing Units,” Opt. Express **17**(22), 20178–20190 (2009). [CrossRef] [PubMed]

**31. **A. Pifferi, A. Torricelli, R. Cubeddu, G. Quarto, R. Re, S. K. V. Sekar, L. Spinelli, A. Farina, F. Martelli, and H. Wabnitz, “Mechanically switchable solid inhomogeneous phantom for performance tests in diffuse imaging and spectroscopy,” in (SPIE, 2015), 10.

**32. **C. D’Andrea, L. Spinelli, A. Bassi, A. Giusto, D. Contini, J. Swartling, A. Torricelli, and R. Cubeddu, “Time-resolved spectrally constrained method for the quantification of chromophore concentrations and scattering parameters in diffusing media,” Opt. Express **14**(5), 1888–1898 (2006). [CrossRef] [PubMed]

**33. **A. Liemert and A. Kienle, “Exact and efficient solution of the radiative transport equation for the semi-infinite medium,” Sci. Rep. **3**(1), 2018 (2013). [CrossRef] [PubMed]

**34. **A. Liemert, D. Reitzle, and A. Kienle, “Analytical solutions of the radiative transport equation for turbid and fluorescent layered media,” Sci. Rep. **7**(1), 3819 (2017). [CrossRef] [PubMed]

**35. **F. Martelli, A. Sassaroli, A. Pifferi, A. Torricelli, L. Spinelli, and G. Zaccanti, “Heuristic Green’s function of the time dependent radiative transfer equation for a semi-infinite medium,” Opt. Express **15**(26), 18168–18175 (2007). [CrossRef] [PubMed]

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

**37. **A. Liebert, H. Wabnitz, D. Grosenick, and R. Macdonald, “Fiber dispersion in time domain measurements compromising the accuracy of determination of optical properties of strongly scattering media,” J. Biomed. Opt. **8**(3), 512–516 (2003). [CrossRef] [PubMed]