Self-calibrating time-resolved near infrared spectroscopy

: 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 multiwavelength time-resolved data can be used for parameters recovery directly without prior calibration (instrument response function measurement).

, 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][13][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. methodology notations will distribution o function (tem -frequency. T work [24,25].  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.

Theory
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: where N ω is a normalization frequency, typically being the first nonzero frequency, as used in this work. The measured DTOF represents the convolved tissue and the instrument response in time and frequency domains such that: where the 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:

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

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 semiinfinite 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  (a) and in freque a)) are analyse Fig. 3(b-c)) us ical TPSFs and falling edges b th separately (λ). These spe entrations (C Hb ereas for the pr the chromoph s is used in the e utilised (inte omain, the algo s approach [2 bO2 = C Hb = 0 nd step tolera ollowing min a S a ; S p } min ∈ using both met in the recovere ns of data ge e highly comp wavelengths, strained in the n in Fig. 3

In-vivo data
The IBIB tim occlusion. Th placed on the blocks: 0.5 m The instru Carlo is used computationa the moving av measurements last waveleng settings and r comparisons a As shown analyse the in The cuff pres . Solid phantom p in (utilising IRF) a n in Fig. 5

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 multifrequency 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.