**Advances in Hyperspectral and Multispectral Optical Spectroscopy and Imaging of Tissue**

Editors

**Vladislav Toronov Mamadou Diop Angelo Sassaroli Ilias Tachtsidis**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editors* Vladislav Toronov Toronto Metropolitan University Canada Ilias Tachtsidis University College London UK

Mamadou Diop Western University Canada

Angelo Sassaroli Tufts University USA

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Applied Sciences* (ISSN 2076-3417) (available at: https://www.mdpi.com/journal/applsci/special issues/tissue optical Imaging).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-4477-9 (Hbk) ISBN 978-3-0365-4478-6 (PDF)**

© 2022 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**


### *Editorial* **Advances in Hyperspectral and Multispectral Optical Spectroscopy and Imaging of Tissue**

**Vladislav Toronov 1,2**


#### **1. Introduction**

Optical imaging and characterization of tissue has become a huge applied field due to the advantages of the optical analysis methods, which include non-invasiveness, portability, high sensitivity, and high spectral specificity. This research field continues to grow and spread in many different directions due to the development of new light sources and detectors, such as, for example, the supercontinuum and tunable lasers, portable highly sensitive spectrometers, multiwavelength photoacoustic imagers, due to novel methods of data analysis, such as the machine-learning methods of spectral analysis, and due to novel applications, such as the imaging of embryogenesis or monitoring of the cerebral oxygen metabolism. Polarization analysis of the anisotropy of optical, mechanical, and electrical properties of materials is another active research direction in hyperspectral imaging from THz to IR spectral ranges.

The purpose of this Special Issue is to provide an overview of recent advances in the methods of tissue imaging and characterization, which benefit from using large numbers of optical wavelengths.

#### **2. Review of Issue Contents**

The human brain is an enchanting object of studies by near-infrared spectroscopy and imaging. In this Special Issue, Guerouah et al. [1] has contributed a profound study of the responses of adult human brain to breath-holding challenges by the hyperspectral near-infrared spectroscopy (hNIRS). This work examined brain signals across the entire near-infrared therapeutic spectral window and showed the critical role of the short-distance channel for the robust measurements of the hemodynamic and metabolic changes in the brain.

Hyperspectral spectroscopy and imaging requires broadband light sources. In the last 20 years, supercontinuum laser sources covering wide spectral bandwidths at high spectral power and fast switching have been developed. Lange et al. [2] contributed a timely and comprehensive review of the features and biomedical and clinical applications of supercontinuum laser sources.

Near-infrared spectroscopy is a quantitative tool that allows for the inspection and characterization of biological tissues and other turbid media such as wood, food and pharmaceutical products, soil, marble, etc. The quantitative accuracy of the medium characterization can be improved using hNIRS. In this Special Issue, Blaney et al. [3] reported a development of a calibration-free hNIRS system that can measure the absolute and broadband absorption and scattering spectra of turbid media.

The ability of the optical imaging techniques to provide real-time quantitative assessment of the biological tissue urges research on their use for tissue monitoring during surgeries. Slooter et al. [4] studied the utility of the measuring multiple tissue parameters

**Citation:** Toronov, V. Advances in Hyperspectral and Multispectral Optical Spectroscopy and Imaging of Tissue. *Appl. Sci.* **2022**, *12*, 3543. https://doi.org/10.3390/app12073543

Received: 21 March 2022 Accepted: 29 March 2022 Published: 31 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

simultaneously by four optical techniques operating at different wavelengths of light: optical coherence tomography (1300 nm), sidestream darkfield microscopy (530 nm), laser speckle contrast imaging (785 nm), and fluorescence angiography (~800 nm) of the gastric conduit during esophagectomy.

Neurosurgical procedures on the open human brain require localization of the functional cortical areas in real time. A Monte Carlo simulation study by Caredda et al. [5] showed feasibility to accurately quantify the oxy- and deoxy-hemoglobin and cytochromec-oxidase responses to neuronal activation and to obtain the spatial maps of these responses using a setup consisting of a while light source and a hyperspectral or a standard RGB camera.

#### **3. Conclusions**

This Special Issue is of interest for the developers and potential users of clinical brain and tissue optical monitors, and for the researchers studying brain physiology and functional brain activity.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


### *Review* **The Use of Supercontinuum Laser Sources in Biomedical Diffuse Optics: Unlocking the Power of Multispectral Imaging**

**Frédéric Lange 1,\*, Luca Giannoni 1,2,3 and Ilias Tachtsidis <sup>1</sup>**


**Abstract:** Optical techniques based on diffuse optics have been around for decades now and are making their way into the day-to-day medical applications. Even though the physics foundations of these techniques have been known for many years, practical implementation of these technique were hindered by technological limitations, mainly from the light sources and/or detection electronics. In the past 20 years, the developments of supercontinuum laser (SCL) enabled to unlock some of these limitations, enabling the development of system and methodologies relevant for medical use, notably in terms of spectral monitoring. In this review, we focus on the use of SCL in biomedical diffuse optics, from instrumentation and methods developments to their use for medical applications. A total of 95 publications were identified, from 1993 to 2021. We discuss the advantages of the SCL to cover a large spectral bandwidth with a high spectral power and fast switching against the disadvantages of cost, bulkiness, and long warm up times. Finally, we summarize the utility of using such light sources in the development and application of diffuse optics in biomedical sciences and clinical applications.

**Keywords:** supercontinuum laser; NIRS; tissue optics; diffuse optics

#### **1. Introduction**

Non-invasive optical monitoring of the human physiology has been a great asset in the medical diagnosis over the last decades [1], and this field is still expanding as technological advances drives the development of new instruments for use in various medical applications [2]. One of the main techniques of optical monitoring of biological tissues is near-infrared spectroscopy (NIRS). It is based on the use of light between 650 and 900 nm, which is within the so-called "transparency window". Within this light window, the optical absorption of biological tissue is at its lowest and light scattering is sufficiently high so the propagation of light in tissue is facilitated. Therefore, NIRS provides a way to gather information from deep into the tissue [3]. As one of the main absorbers in this window is the hemoglobin, in its oxygenated and deoxygenated forms (oxy- and deoxyhemoglobin ([HbO2] and [HHb], respectively), NIRS has been broadly used to quantify the changes in oxygenation of various tissues like muscle and brain [4]. This ability to quantify oxygenation non-invasively, and in real time has helped NIRS to become a standard technique in the labs. However, a few pitfalls are still holding back its wide adoption at a clinical level. Indeed, most of the widespread techniques used for NIRS are based on the socalled continuous wave (CW) technique where the tissue is illuminated with a non-coherent light and the variation in light intensity is collected at the detector [5]. In this configuration, it is harder to estimate absolute quantities, as the underlying optical properties of the tissue (mainly absorption and scattering coefficient) cannot be measured, and the scattering properties are assumed to be fixed (i.e., in time and for every subject). Therefore, most of the CW-NIRS systems only quantify the changes in hemoglobin concentration. This can be a drawback for clinical use, as an absolute comparison between two patients

**Citation:** Lange, F.; Giannoni, L.; Tachtsidis, I. The Use of Supercontinuum Laser Sources in Biomedical Diffuse Optics: Unlocking the Power of Multispectral Imaging. *Appl. Sci.* **2021**, *11*, 4616. https:// doi.org/10.3390/app11104616

Academic Editor: Byoung-Kwan Cho

Received: 11 April 2021 Accepted: 7 May 2021 Published: 18 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

cannot be done unless you involve a well-controlled physiological stimulus. The other drawback is that the measurement is often performed with only two wavelengths, the minimum necessary to retrieve two absorbing compounds [HbO2] and [HHb]. However, the contribution of other chromophores is either neglected (like collagen), or assumed and fixed, like water. This results in an increase in the physiological noise, which reduces the accuracy of the measurement.

These drawbacks have been tackled over the recent years, and some new techniques and instrumentation enabled researchers to work around these drawbacks. Firstly, one can use more advanced CW-NIRS techniques in order to extract absolute information about the tissue. We can cite the broadband NIRS that also has the capability to extract absolute information of tissue optical properties [6]. The added benefits of this technique are that, as it uses a large number of wavelength (typically more than a hundred), it can perform a true spectroscopic analysis, hence also gathering information about other chromophores like cytochromes, water, and fat [6]. This extra information enables having a complete optical view of the tissue under investigations. Moreover, on top of providing a more precise view of the oxygenation, the use of many wavelengths enables targeting novel optical contrasts, like the oxidative state of the cytochrome-c-oxidase (oxCCO) which is a biomarker of metabolism [7].

However, even though these improvements are promising, performing a true spectroscopic measurement on tissue using NIRS is still challenging. One of the main issues is the penetration depth of the light and the contamination of the tissues of interest by the upper layers (i.e., skin) [3]. Indeed, the tissues are often modelled as homogeneous for simplicity, even though this hypothesis is not true. Therefore, new methodologies are required to enhance the discrimination between different tissue types, thus improving the accuracy of the measurement.

It is well known that the time-domain NIRS (TD-NIRS) is the most advanced and accurate type of measurement, capable of retrieving absolute information and to discriminate between shallow and deep tissue [8]. Briefly, in TD-NIRS an ultra-short light impulse (typically in the picosecond range) illuminates the tissue. Tissue has the effect of attenuating and broadening the re-emitted pulse (the resulting pulse-width is typically few nanoseconds wide). The measured signal is typically called a temporal point spread function (TPSF) or distribution of time of flights (DTOFs) of photons. Indeed, what is measured here is the histogram of the arrival time of photons. Thus, the amount of information acquired by TD-NIRS systems outstrips the classical CW-NIRS, which only acquires the intensity. From these DTOF, various data processing methods can be employed to extract information about the tissue. It is out of the scope of this paper to detail these methods, but the interested reader can refer to the following reviews for a more in-depth description of TD-NIRS methodologies [9,10], and for a special focus on broadband TD-NIRS [11], which will be largely discussed in the present work.

The extraction of more accurate information with TD-NIRS comes at the cost of a more complex instrumentation, requiring pulsed-laser source (pulse in order of MHz) and fast electronic and sensitive detectors in order to be able to unlock single photon counting. Thus, even though TD-NIRS has been around for many decades now and has been developed in parallel to CW-NIRS, these more specific requirements have hindered its development. Nevertheless, much progress has been made over the last two decades, and TD-NIRS is closing the gap with CW-NIRS in terms of user's adoption.

One of the key developments that has promoted the development of TD-NIRS is the supercontinuum laser (SCL) sources, which provide a white laser source, with a high power (several watts) and a very short pulse width (typical less than 10 ps) at several tens of MHz. Indeed, these characteristics are essential in order to develop reliable TD-NIRS systems, that started to be developed in the nineties. During these early days of development, no commercial pulsed laser diodes were available, and one had to rely on homemade SC generation or Ti-Sapphire lasers in order to be able to produce a coherence pulsed light with a spectral power density sufficient to be detected by the photon counting electronic. Therefore, SCL sources have been a true driver of TD-NIRS development in the early days, and have more recently unlocked new horizons like spectroscopy, which have broadened the fields of applications.

Finally, besides traditional sparse measurement of tissue properties with diffuse optics, Diffuse Optical Tomography (DOT) has also been rapidly expending in the past decades. This imaging technique is based on the acquisition of a large number of channels, coupled to a more refined data analysis, based on anatomically realistic models of tissue. This allows to reconstruct 3D spatial distribution of optical/physiological parameters of the tissues, rather than providing tissue information on a limited number of points. An extension of this technique is called fluorescence DOT, that uses extrinsic near-infrared fluorescent dyes which targets specific tissues to harvest extra information and enable molecular imaging. The interested reader can refer to reviews [12,13] for more information about these techniques. We will see that SCL sources have also a role to play in order to improve the performances of the DOT instruments.

In this work, we will review the use of SCL in the field of biomedical diffuse optics. Thus, we will begin by a short overview of the history and basic principle of supercontinuum generation. Then, we will review works that have been using SCL in biomedical diffuse optics, from the instrument development to the application. Finally, we will summarize and highlight the strength and weakness of SCL for these applications and show what areas can benefits from its use.

#### **2. History and Basic Principles of Supercontinuum Lasers**

Supercontinuum (SC) lasers are broadband, coherent light sources that are capable of combining the wide spectral range of conventional white lamps (including UV, visible and IR light) with the optical performances of single-mode lasers, in particular their high output intensity, the capability of fast, pulsed emission, and the fine collimation of the beam down to diffraction limit. Such characteristics make this type of photonic technology highly versatile for a number of applications, in particular in biomedical optics and biophotonics, thanks also to recent advancement towards making SC lasers more compact and affordable [14].

The physical process of optical SC generation was discovered and studied in the 60 s and 70 s with milestones publications by Alfano and Shapiro [15,16]. It is based on the propagation of extreme non-linear effects from a narrow-bandwidth pulsed laser beam in order to generate a broadening of the overall spectrum of the emitted light. These nonlinear effects typically include self-phase modulation, four-wave mixing, Raman scattering, and solitons dynamics, among others [14,17]. The propagation of these effects is performed by dispersing the pulsed laser beam through a non-linear material. These can range from bulk materials, such as simple water cells, silicas, or gases (as in the earliest attempts of SC generation), to more advanced crystal-structured optical fibres [17,18].

Early SC lasers using bulk materials for non-linear propagation were limited in the efficiency of coupling the input laser light with these propagating material (coupling efficiencies stood at below 10%) and thus they were able to only achieve ranges of about 200–300 nm, with pulses duration in the order of tens of nanoseconds. Further developments in SC generation efficiency and coupling using optic fibres lead to enlarging the range up to 1–2 μm and increasing the pulse duration to the range of picoseconds (for commercial, compact lasers), or even femtoseconds (for research-grade, benchtop lasers) [18]. A critical turning point for modern SC lasers to achieve such high performances was the invention and application of photonic crystal fibres (PCF), a particular type of microstructured optic fibre (Figure 1).

**Figure 1.** Example of commercial SC laser sources and their typical spectra. (**a**) Fianium source, (**b**) Leukos source, (**c**) FYLA source, (**d**) YSL source. (**e**) Schematic of the classical triangular cladding single-core photonic crystal fiber in which light is guided in a solid core embedded in a triangular lattice of air holes. The fiber structure is determined by the hole-size, d, and the hole-pitch, Λ. Like standard fibers, the PCF is coated with a high index polymer for protection and to strip off cladding-modes. Extracted from the application note on supercontinuum generation in photonic crystal fibers by NKT. (**f**) Example of full spectrum generated by a SC laser, showing the width of the broadband emission and the central peak wavelength that is a residual of the single-mode laser used for SC generation. The spectrum is taken from the SuperK COMPACT SC laser manufactured by NKT Photonics (https://www.nktphotonics.com/lasers-fibers/product/superkcompact-supercontinuum-lasers/, accessed date: 1 November 2021).

PCFs, invented in 1996 [19], are currently the most used non-linear propagating materials for SC generation [17,18]. As the name suggests, this class of optic fibres is based on recreating the optical properties of crystalline structures, in particular their photonic bandgap effect, by combining various high and low refractive indexes in their inner structure. Such features can be achieved in different manners, which typically define the given category of PCFs [20]: (1) hollow-core PCFs present a central hole in their structures, filled with air; (2) holey PCFs are composed of a matrix of several air-filled holes in their cross-sections that is interspersed with the core material; (3) solid-core PCFs involves a solid core surrounded by a matrix of air-filled holes.

PCFs offer high coupling efficiencies and low power losses, thus allowing ultrabroadband, high-brightness SC generation at cost-effective access [17,21]. SC lasers based on PCFs used in biomedical optics applications can typically generate broad spectrum spanning from 350 to 2400 nm, generally showing a central sharp peak which is a residual of the original shape of the pulsed beam. An example of this is shown in Figure 1. Most commercial SC lasers show performances similar to pulsed laser diodes, with integrated output power up to tens of watts and the capability of generating pulses in the order of picoseconds and repetition rates of tens of MHz [22].

Finally, we can report that in the recent year, more companies are able to provide ready-to-use SC solution. A search from the website https://www.rp-photonics.com (accessed date: 1 November 2021), displays at least 4 companies manufacturing SCL: NKT photonics, who recently acquired the former UK company Fianium (www.nktphotonics. com, accessed date: 1 November 2021), LEUKOS (www.leukos-systems.com, accessed date: 1 November 2021), FYLA (www.fyla.com, accessed date: 1 November 2021) and YSL photonics (www.yslphotonics.com, accessed date: 1 November 2021). Figure 1 displays pictures of different SCL offered by these companies.

#### **3. The Use of Supercontinuum Lasers in Biomedical Diffuse Optics**

The focus of the present review was to evaluate the use of SCL in biomedical diffuse optics, from the systems and methodological developments to their medical applications.

Therefore, papers were identified using PubMed, Scopus, and Web of Science, searching for a combination of keywords including (supercontinuum laser|white laser) and (diffuse optics|NIRS). Moreover, a manual search from articles' references was performed and the personal reference library of the authors of the current papers has also been searched. Papers were rejected if no SCL was used in the method section of the reported work, if the application was not medically related, if they were reviews or if the method was not purely optical (i.e., no photoacoustic papers). A total of 95 publications covering a variety of focus were identified. We have identified 4 main categories of publications: system development (31 studies), optical properties estimation (16 studies), methodology (43 studies) and application (5 studies). Table 1 provides a summary of the main characteristics of the SCL used in the studies reported in this review, Figure 2 shows pictures of the instrumentation typically used, and Figure 3 presents typical output in terms of datatypes of such systems. Moreover, a summary of the studies included in this review is presented in Table 2 and includes the system used, the year, the category of the study and the target of the study. Finally, Figure 4 summarizes the main topics and characteristic of the systems used in the reviewed publications and shows the number of publications over the last 30 years that report the use of SCL in biomedical diffuse optics.

#### *3.1. Novel Instrument Developement*

The first paper reporting the use of a supercontinuum laser in biomedical diffuse optics was reported in the early nineties [29–31] by a group based at the Lund Laser Centre (LLC, Lund, Sweden). The femtosecond white light supercontinuum was generated by a Terawatt laser running at 10 Hz. The spectral decomposition was performed on the detection side, by coupling an imaging spectrometer and a streak camera, enabling to acquire parallelly the DTOF at all wavelengths. To improve the signal-to-noise ratio, the measurement had to be repeated for several laser shots, but the overall acquisition time was of just a few minutes, which is acceptable when no dynamic contrast is targeted. Using such a system, the very first in vivo absorption and scattering spectra of the female breast were acquired from 650 to 850 nm. This system was then updated in the early 2000s with a supercontinuum generation system based on a crystal fibre, which allowed to increase the optical power of the SCL generation [32].

*Appl. Sci.* **2021**, *11*, 4616

**Table 1.** List of the SC laser-based systems used in the publications reviewed, with their principal characteristics. SC: SuperContinuum, FWHM: Full Widht at Half Maximum, CCD:Charge-Coupled Device, ICCD: Intensified Charge-Coupled Device, PMT: PhotoMultiplier Tube, SPAD: Single-Photon Avalanche Diode, SiPM: Silicon PhotoMultiplier, HPM: Hybrid Photodetector Module, sCMOS: scientific Complementary Metal–Oxide–Semiconductor, NIR: Near InfraRed.


*Appl. Sci.* **2021**, *11*, 4616

**Table 1.**

*Cont.*




**Table1.***Cont.*

10

**Figure 2.** Pictures and schematic of various instruments based on SCL as a source. (**a**) TD-NIRS instrument aimed at broadband characterization of tissues. Figure reproduced from [23]. (**b**) TD-NIRS instrument aimed at monitoring both oxygenation and metabolism functional changes. Figure reproduced from [24]. (**c**) Preclinical instrument aimed at monitoring the exposed brain cortex of rodents. (reproduced with permission from [25], Copyright IEEE, 2021). (**d**) DOT system based on a SPAD camera aimed at functional brain monitoring in neonates. Figure reproduced from [26,27].

**Figure 3.** Example of typical data extracted form system based on SCL as sources. (**a**) absorption and reduce scattering coefficient of several tissues acquired the TD-NIRS system developed in ref [23]. Figure reproduced from [23]. (**b**) 2D maps of the changes in optical properties of a phantom using MONSTIR II. (reproduced with permission from [28], Copyright AIP Publishing, 2021). (**c**) Changes in [HbO2] (red), [HHb] (blue) and [oxCCO](green), in the occipital lobe, acquired with the MAESTROS system [24], during a typical brain activation. Data presented at the fNIRS conference in Tokyo in 2018. (**d**) Raw image and reconstructed variation in [HbO2], [HHb], and [oxCCO] acquired on the exposed cortex during anoxia. (reproduced with permission from [25], Copyright IEEE, 2021).

**Figure 4.** Summary of the literature review results. (**a**) Graph showing the number of publications using SCL in biomedical diffuse optics since the first study in 1993. (**b**) Summary of the main focus of the papers reviewed.

Parallel to these initial developments, a team based at Politecnico di Milano (POLIMI) produced a system able to characterise tissue between 610 and 1050 nm [33]. The system was however not based on a supercontinuum laser generation but on two tuneable lasers for the source, and the detection scheme was based on time-correlated single-photon counting (TCSPC), enabling light detection down to the single photon, which drastically improves the sensitivity of the system. This increased sensitivity is of great advantage to reduce the acquisition time and improve the light collection in the region where water absorption is high (900–1050). However, handling two lasers drastically increases the complexity and bulkiness of the system. Nevertheless, with that system, that team was able to record the optical properties of absorption and scattering of various tissues invivo [34–37]. The availability of commercial SCL sources gave rise to a new version of that system, where the 2 lasers were replaced by a SCL source coupled with a prism for the spectral selection [38]. This new source allowed to make the system more compact and transportable making it more suitable for the clinic. Moreover, the detector was also replaced (previously photomultiplier (PMT)) by a Single-Photon Avalanche Diode (SPAD) which enable to increase the spectral coverage, allowing a measurement between 600 and 1000 nm. This system could be used to characterise the optical and physiological properties of the tissue of the female breast. Another version of that system was developed, combining the SCL with a Ti-Sapphire laser to expand the spectral coverage to 1100 nm [39]. The SCL was used to cover the 600–900 nm and the Ti-Sapphire was used between 900 and 1100 nm. Indeed, even though the SCL could cover the entire considered wavelength band, the spectral resolution provided by the couple SCL/dispersion prism was too low to provide accurate results, as the source bandwidth increased with wavelength from 3 nm at 600 nm to 10 nm at 900 nm, and to approximately 15 nm at 1100 nm with that solution. This instrument has indeed been used to demonstrate that the broadening of the spectral bandwidth effect could impair the results, as the spectral region within the bandpass that exhibits the lowest absorption will dominate, leading to significant underestimations of the absorption and distortions in the spectral shape [40]. The spectral quality of the acquired data being extremely important, having the narrowest bandwidth is a great advantage. To do so, using a SCL is thus a true advantage as it provides a high power across the entire spectra. However, the method used to filter the white light is very important. As it has been seen above, several techniques can be used, either filtering on the detection side with a spectrometer, or on the source side with a dispersive prism or using an acousto-optic tuneable filter (AOTF), as we will see in systems described below. The POLIMI team has compared the wavelength selection from a SCL with an AOTF and a prism. They have

shown that one had to be careful when using AOTF, as they can produce side lobes in the selected bandpass, which can drastically impair the retrieval of the optical properties [41].

Nevertheless, the extension of the spectral coverage offered by this newly developed system was of particular interest, as several components, like lipid and collagen, have a distinct signature above 1000 nm. The LLC group also reported a similar system, based on a commercial SCL coupled with 2 AOTF and a detection scheme based on TCSPC with a coverage between 650 and 1400 nm [42]. The POLIMI group has then extended this work by focusing on the spectral window between 900 and 1700 nm, by using an InGaAs/InP SPAD as a detector [43–45]. Finally, they have produced a system able to cover the entire bandwidth from 600 to 1350 nm using a single SCL coupled with a prism for the wavelength selection, and 2 TCSPC cards couples with 2 detectors covering different bandwidth [23]. They managed to fix the issue of the spectral broadening induced by the prism by refining the optical chain. Additionally, this system was mounted on a 19-inch medical cart, making it easily used in a clinical environment. This system can be seen in Figure 2 and in an example of data in Figure 3. Moreover, by adjusting the detector used, the same instrument can be used to cover an even larger bandwidth, between 500 to 1700 nm [46].

With all these developments, POLIMI has been a major contributor to multispectral TD-NIRS and were able to characterize various tissue over a large bandwidth. They were able to characterize bones [47–49], thyroid tissues [50], adipose tissues [51], collagen [46,52], and elastin [53]. This groundwork was crucial in order to investigate the tissue composition of breast for example, which can help for diagnosis of cancerous tissue [54,55]. Finally, we can report that the latest reported system has also been used to monitor the thermal effect of Radio Frequency (RF) on tissues [56,57].

The previously reported systems were only composed of one channel, which can be a drawback when an image is required, like for performing mammography. Therefore, parallel to these developments, another system aimed at optical mammography was developed, based on a supercontinuum light generated in PCF and a 32-channel parallel detection [58]. This allowed to acquired image of the breast between 606 and 885 nm, and to extract hemoglobin, tissue saturation, lipid and water content, and scattering parameters of the tissue.

All the systems reported above were designed to collect spectral information about tissue, with a first application in breast cancer detection. Therefore, the acquisition time was not limited, and functional contrast was not explored. However, the potential of monitoring real time change in optical properties over a large bandwidth is quite appealing for various application, like brain monitoring [59]. Therefore, other systems were developed in order to achieve that goal. A system was developed at POLIMI based on a supercontinuum source generated in a PCF by a self-locking mode-locked Ti:Sapphire oscillator and a detection scheme based on an imaging spectrometer coupled to a 16-channel multianode PMT connected to a TCSPC board [60]. The fact that all the 16 wavelengths were captured at once enabled a fast acquisition and permitted to track dynamic variations in absorption and scattering spectra following hemodynamic changes [61]. The limitation of the detections scheme constraint the bandwidth of this system between 520 and 850 nm.

A similar approach was taken by a group based in Poland [62,63]. The system was based on a commercial SCL for the source, coupled to a 16-channel multianode PMT connected to a TCSPC board for the detection part. This system was used to follow a contrast agent (Indocyanine green (ICG)), which allows to estimate the cerebral blood perfusion. The same system can be used without ICG to monitor the intrinsic optical contrast [64]. Similarly, the author and colleagues used a TD-NIRS system based on a SCL filtered at 808 nm and on a detection scheme based on a single channel TCSPC to follow ICG bolus [65]. The novelty of this study was the fact that the TD-NIRS was coupled with a diffuse correlation spectroscopy system (DCS) [66]. The coupling of these two techniques allowed to retrieve absolute cerebral blood flow measurement in patients in the intensive care unit, showing that these systems could be used in a particularly difficult clinical environment. Moreover, the possibility to also measure blood perfusion is of great interest and might help in the adoption of the technique in the clinic.

Another potential detection scheme in order to collect parallel data is to use an Intensified CCD (ICCD) camera. Selb and collaborators used a commercial SCL combined with a ICCD camera to perform functional brain imaging with 2 wavelengths at a frequency of 3 Hz [67]. This system was an upgrade of a previously developed instrument based on pulsed Ti:Sapphire laser [68]. This source was tuneable but the switching time between wavelength was too slow to enable dual-wavelength functional imaging. Therefore, the authors report that the use of SCL unlocked this possibility with their detection scheme.

Similarly, Lange and collaborators described a TD-NIRS system developed to measure the human brain tissue physiology [69,70]. This system was based on commercial SCL for the emission and on an ICCD camera coupled with an imaging spectrometer in order to perform the spectral decomposition. This system allowed to retrieve the typical hemodynamic response, and preliminary results on the detection of oxCCO were also performed. The ability to retrieve oxCCO is exiting for brain applications, as it has been shown that this contrast can be more specific than the oxygenation, and that it can be an early biomarker of neonatal encephalopathy for example [7]. This possibility has been explored by the UCL team by developing a 4-channel 16 wavelengths TD-NIRS system in order to monitor both the oxygenation and oxCCO signal [24,71]. This system was based on time multiplexing, with a TCSPC detection scheme, but the fast-switching capabilities of the SCL coupled with AOTF allowed the acquisition of all the 16 wavelengths in less than 2 s. This system has also been used to evaluate the long-term reproducibility of the evaluation of cerebral tissue saturation in healthy adults with good reproducibility [72]. This system can be seen in Figure 2, and in an example of data in Figure 3.

Another approach that has been used using SCL sources is the development of noncontact systems. This can present several advantages to avoid issues attributed to sensortissue contact like skin compression and to obtain a reduction in measurement preparation time. A team at the Physikalisch-Technische Bundesanstalt (PTB) has developed a noncontact TD-NIRS system based on a commercial SCL and fast-gated SPAD to explore the human brain physiology [73–76]. At the time, the authors reported that the use of SCL permitted a shorter pulse width (<100 ps), an achievable output power considerably larger than picosecond diode lasers, and the possibility to quickly switch between two different wavelengths (i.e., on the μs to ms time scale). Moreover, the use of the SPAD in that system enables to reject the burst of early photons, which can improve the spatial resolution of the system at very short detector distances [77]. This process had been previously demonstrated by the POLIMI team that had develop a single channel TD-NIRS system based on SCL and a fast gated SPAD detector [78]. The POLIMI team has also worked on the modelling of the parameters of such detection schemes, in order to optimise their design [79].

The non-contact scheme is also very useful for pre-clinical applications. Indeed, several systems, based on SCL for the source, aiming at collecting physiological information of small animals, have been proposed. In the early 2000s, a system based on a home-made SC generation associated with an imaging spectrometer coupled with a streak camera was developed and used to monitor the cerebral hemodynamic of songbirds [80,81]. More recently, the UCL team have developed a system aiming at investigating both the cerebral hemodynamic and metabolic response (via oxCCO) of the exposed cortex of rodents [25]. In this work, a commercial SCL was used in combination with a dispersive prism, for the emission, and a sCMOS camera, as a detector, to record 11 images (11 wavelengths) of the exposed cortex in less than 0.3 s/wavelength. Here the recording was not time resolved, but the power of the laser and the ability to quickly switch between the wavelength had dictated the choice of this source. This system can be seen in Figure 2 and an example of data in Figure 3. Another approach can be to use spatial frequency domain imaging (SFDI), a wide-field, noncontact, model-based technique that can quantitatively assess absorption and scattering, on a pixel-by-pixel basis by calculating the modulation transfer function

of structured light projected onto the tissue. This technique can be hyperspectral if a proper source is used. However, the system developed thus far was based on conventional broadband sources coupled with tunable filters, which limited their acquisition speed and spectral resolution. Torabzadeh and colleagues [82] applied this technique using the principles of the single-pixel compressive sensing. They used a commercial SCL coupled with a dispersive prism to achieve a final bandwidth 1000 spectral bands between 580 and 950 nm, which was larger than the previously developed systems. With this system, the authors were able to characterize the optical properties of a beef sample ex-vivo in 2D, and extract physiological information (i.e., hemoglobin content and fat and water content). The authors noted that the wavelength selection technique used was still slow (i.e., 150 s of acquisition time for a 4 cm × 6 cm), but that switching to an AOTF solution could greatly improve this speed. Using the same principle of the single-pixel compressive sensing, Farina and colleagues reported a TD-DOT instrument based on time-resolved single pixel camera [83]. Their approach was based on a SCL source filtered at 620 nm coupled with a digital micromirror device (DMD) to generate the spatial pattern. The time domain acquisition was based on TCSPC. Pian and colleagues used a similar approach but expended it to a broadband illumination and detection [84]. The authors reported that their system, aimed at fluorescence molecular tomography, benefitted for the commercial SCL source to perform hyperspectral excitation of the sample which offers a spectral information boost within the same data acquisition time compared to single wavelength excitation situations.

Another work aiming at improving the spatial resolution of TD-DOT was to investigate the effect of only considering early photons, which undergo fewer scattering events [85]. In this work, this approach has been tested using a commercial SCL together with a bandpass filter at 670 nm. It worth noting that the commercial version of this laser was modified by the manufacturer so that it had a shortened fibre compared to standard models, which reduced the pulse width to approximately 300 ps. This is important as the benefit of the gating of the early photons is highly dependent on the instrument response function (IRF) of the system. The authors reported that with this technique, and this instrument, a 64 to 84% image resolution could be achieved, which could be a great benefit for preclinical imaging with diffuse optics.

Novel algorithms can also be used in combination to TD-DOT instruments to recovers 3D images of the optical properties and physiological parameters. A team in Grenoble has explored the use of the Mellin-Laplace transform to reconstruct 3D images of the chromophores inside phantoms and for the non-invasive assessment of flap viability in rats' models. There instrumentation was based on commercial SCL, coupled with either an AOTF or an interferential filter, for the source and a TCPSC scheme for the detection [86,87].

TD-DOT has also been explored in humans, in order to reconstruct 2D or 3D maps of the optical properties and hemodynamic parameters. To this end, the UCL team had design MONSTIR II to perform tomographic brain imaging [28]. This system is based on a commercial SCL with an AOTF filtering and has 32 channels based on TCSPC. It was successfully used to scan healthy and ill infants in a clinical environment [88,89], an example of data produced by this system is provided in Figure 3. More recently, a new TD-DOT system based on a commercial SCL filtered by an AOTF, and a 32 × 32-pixel SPAD camera, has been developed by a team at the Biomedical Optics Research Laboratory at the University of Zurich [26,27,90]. This system operates at 2 wavelengths and at 11 sources points, which can be scanned in about 3 s. It is a big step forward compared to previously developed system, especially considering its size, as it is handheld. A picture of this instrument is reported in Figure 2.

Finally, all the instrument mentioned above were mostly based on the TD technique, or traditional white light illumination. We can note that we found one frequency domain (FD) system, that is based a commercial SCL to collect our system collects FD data at 170 wavelengths from 680 to 850 nm and allows to accurately measure absorption and scattering information about the tissue [91].

#### *3.2. Methodology*

We have seen that the SCL has been used to develop several systems, mainly aiming at characterising the optical properties of tissues over a large bandwidth or to non-invasively follow the oxygenation and/or metabolic changes both in humans and in preclinical models. Both of these capacities are extremely appealing for clinical applications. However, the transition to the clinical application of system based on SCL and more generally TD-NIRS/DOT is still happening at a slow pace. One of the reasons of the slow adoption of optical techniques is the lack of standards in the community. Indeed, standardisation is key in order to test all the systems on the same ground, and reenforce the trustworthiness of the technique. In the recent years, standardisation have been a main goal of the NIRS community, and several protocols have been developed in order to compare and test the capabilities of NIRS systems. The main produced protocols were the MEDPHOT [92], BIP [93], and nEUROpt [94]. It is out of the scope of this paper to describe these protocols in detail, but the cornerstones of all of these are to produce standard phantoms, well characterized spectrally in terms of absorption and scattering coefficients, in order to be able to quantitatively compare different systems or methods. In this matter, SCL based systems were extensively used to characterize the material used to produce these phantoms and the actual final phantoms used as standard. These optical phantoms can either be solid (easier to store and stable over a long period of time) or liquid (easier to precisely tune the optical properties of the solution and to adjust it dynamically). For the liquid phantoms, the core principle is to use a fatty solution to control the scattering coefficient of the solution, and ink to control their absorption properties. Several recipes have been tested, like using Agar and Triton [95], but the most common recipe to date is to use a combination of water, intralipid and Indian ink. This choice has arisen after a careful characterization of both intralipid and Indian ink, principally with instrument based on SCL [96–100]. We can note that for the precise optical characterisation of liquid, techniques based on integrating sphere are often used. However, this methodology cannot be applied for large solid phantoms. In this matter, TD-NIRS has proved a good methodology to accurately estimates the optical properties. For example, in 2010, Bouchard and colleagues used a TD-NIRS setup based on a SCL source filtered at 610 nm in order to estimate the optical properties of phantom with an absolute error estimates of 0.01 cm−<sup>1</sup> (11.3%) and 0.67 cm−<sup>1</sup> (6.8%) for the absorption coefficient and reduced scattering coefficient, respectively [101].

On top of designing phantom able to mimic the background optical properties of tissues, some methods have emerged to mimic a perturbation of the absorption coefficient within a homogeneous medium. To that end, one can either use a liquid phantom for the homogeneous medium, and a small black PVC target for the perturbation. A protocol has been designed and it is possible to precisely set the absorption perturbation value based on the size of the black PVC inclusion [94]. This protocol has been validated with Monte Carlo simulations and in-situ using a system based on SCL [102]. A similar phantom based on a completely solid phantom, with a movable solid rod, has also been designed and characterised [103]. Lastly, some phantoms have also been characterized in order to be a tool to evaluate quantitatively the responsivity of NIRS systems [104].

Finally, we can also mention new approaches based on digital phantoms. Basically, rather than physically simulating the process of light propagation in tissue, the digital phantom provides the detector with light signals mimicking the signals obtained in-vivo. Such systems based on an integrating sphere and several light sources, including a commercial SCL, have been proposed to spectrally characterise and calibration optical systems [105]. More recently, a digital phantom has been developed to mimic TD-NIRS signal acquisition, and the authors used various sources and detectors in order to test this approach, including a SCL [106]. One other important component to optically characterize is the optode holders used to secure the optical fibres to the patient. In a recent study, Amendola and colleagues used a system previously developed to evaluate the optical properties of fruit [107] (based on a commercial SCL and a filter wheel to select 14 spectral band in the range 550 to 940 nm) to spectrally characterised poly lactic acid (PLA), a material used for

3D printing [108]. As 3D printing is more and more common in order to produce probe holders, it is important to optically characterise the material used in order to avoid any unwanted effects on the measurement provoked by light interaction with the material. In that study, the authors concluded that the different material tested were not optically the same, and that characterisation was important before their use.

The development of all these phantoms has provided a good framework for the developers to compare their systems and methods in order to provide robust clinical information. Regarding the instrumentation, a lot of effort are currently being made in order to improve the detection scheme of TD-NIRS. To that end, SCL have been a source of choice to test new detectors like the promising fast silicon photomultiplier (SiPM) which drastically improve light harvesting compared to the old PMTs [109–111]. SCL based systems have also been used to characterise novel bio-compatible fibres that can be implanted in the patient with potential use in the monitoring of the evolution of the physiology after a surgical intervention or in photodynamic therapy (PDT) [112,113].

Regarding the methodology, several new data processing techniques have been explored and tested in conjunction with SCL. Indeed, system based on SCL can collect multidimensional data that can be coupled together in order to improve the accuracy of the data, reduce the noise or explore new contrast. One of the first work regarding this matter dates from 2006, where D'Andrea and colleagues, showed that one could increase the robustness of the determination of the chromophore's concentrations in tissues (hemoglobin, water, lipid) by using a spectral constraint during the fitting procedure [114]. This showed the strength of using the spectral dimension coupled with the DTOF, which is accessible when using SCL. This method has also been used to monitor absorption changes in a layered diffusive medium [115], which is relevant in the monitoring of brain activation for example. Similar approaches were explored in the same group in order to accurately monitor the spectral changes in absorption coefficients of turbid media [116].

One can also use the spatial dimension in order to enhance the amount of information recorded, by recording data from multiple source-detector distances. This approach has been recently tested and it has been shown that coupling the spatial dimension with the DTOF could increase the accuracy of the estimation of the optical properties [117–119]. Another interesting use of the spatial dimension is the self-calibrating method for TD-NIRS, which uses the spatial dimension in order to avoid the measurement of the IRF, that is needed to remove the effect of the instrument from the response of the tissue [120]. The characterisation of the IRF necessitates an extra step that can be a burden in a clinical environment. Therefore, this method could facilitate the use of TD-NIRS in the clinic. Finally, we can also report that an approach combining the spatial and spectral dimension proved to also improve the accuracy of the results of CW-NIRS. This has been tested recently, especially with blood phantoms and compared to the results of a TD-NIRS instrument based on a SCL source [121].

The large amount of data collected by TD-NIRS also enables different data processing scheme to extract the relevant information. Therefore, there is a need to properly compare all the available methods in order to evaluate their strength and weaknesses. This work has recently been undertaken in order to compare the moment and time windowing techniques [122,123]. Having a deep understanding of all the data processing method is also crucial so the proper technique can be used depending on the application.

Several multi-centre initiatives have recently promoted this standardisation effort to compare instruments and algorithms. To that end, partners from the BitMap network (http://www.bitmap-itn.eu/, accessed date: 1 February 2021) have started comparing instruments from different institutions all over Europe, by using the same sets of solid phantoms from the 3 main protocols [124,125]. Moreover, several centres in Europe have used their instruments (several based on SCL) on 9 subjects, to study the inter-subject variability of the optical properties of the human head measured by several instruments [126]. They could show that the inter-subject variability was significant, with a big effect of the technique used (either CW or TD) because of their different depth sensitivity, which in

turn affects the averaged optical properties retrieved. Therefore, it is evident that all the efforts to precisely characterized instruments and methods will have a huge impact on the precision of the optical measurement of tissue parameters.

Finally, the use of SCL can help to explore new avenues in diffuse optics. The main one is to go above the "therapeutical window" in the second (1100 nm to 1350 nm) and the third (1600 to 1870 nm) near-infrared optical window. On top of the previously reported works that expand the bandwidth of TD-NIRS systems up to 1700 nm, Sordillo and colleagues reported the use of a commercial SCL source, filtered with a bandpass and long pass filters to target the second and third NIR windows, coupled with an IR-CCD camera, to acquire transmission images of 3 targets embedded in chicken tissue of various thickness [127]. They were able to distinguish the targets up to a depth of 10 mm, which was not possible when using the same setup but with a traditional lamp source. These results are promising in order to implement transmission measurement, especially for preclinical imaging.

As we have seen previously, the use of these long NIR wavelength is useful to resolve more chromophores like collagen. It is also promising to resolve glucose, which could be a true benefit for several clinical applications like monitoring diabetes [128], or monitoring brain injuries [129] or neonatal encephalopathy [130]. Once again, in this wavelength range, the power delivered by SCL has proved efficient to monitor the glucose contrast in the NIR [131,132].

Finally, a tentative to non-invasively probe the lungs with a TD-NIRS system based on SCL source was reported [133]. The measurements were performed on three subjects, and spectra from 600 to 1100 nm could be acquired. Here, the depth sensitivity of a TD-NIRS is essential and this application shows that instruments based on that technology could be used to non-invasively probes tissues previously thought out of reach.


*Appl. Sci.* **2021**, *11*, 4616








*Appl. Sci.* **2021**, *11*, 4616

#### **4. Discussion**

We have summarized the use of SCL sources in biomedical diffuse optics. Historically, these sources were allowed to develop TD-NIRS system as they could deliver a pulsed coherent light at several MHz with narrow pulse width in order to reduce the IRF [137], and high power, of several mW/Wavelength, which were characteristically difficult to achieve in the early 1990s. Moreover, the development of commercially available compact SCL in the 2000s also helped the development of the use of these sources, and we have reported that more companies (see Figure 1) are now able to provide such sources.

The main advantage of these sources is that they make it easier to explore more wavelength, enabling a multispectral or even broadband or hyperspectral measurement. Several measurement strategies can be implemented in order to do so. The first one is to use the white light as a source and achieve the spectral unmixing at the detection stage. The first available option to do so was to use an imaging spectrometer coupled with a streak camera, which achieve very high temporal resolution (in terms of arrival time of photons). However, this solution was bulky, and the sensitivity of these camera was low compared to the newest detection technologies. The same method can be used by either coupling the spectrometer with an ICCD camera, or a 16-channel TCSPC. This approach of performing the spectral unmixing on the detection can achieve faster measurement since no spectral switching mechanism is employed. However, the resolution of these approaches, either spectrally or temporally, can be limited. For example, typical ICCD cameras have a gate width of 200 ps, which limits the temporal resolution of the DTOF. On the other end, the 16-channel TCSPC system used by Sudakou and colleagues [64] allows spectral bands of 12.5 nm, limiting the spectral resolution. The other option to perform the spectral unmixing is to do a sequential acquisition, by filtering the white light on the source side. To do so, two main approach are used, (1) using an AOTF, with the ability of fast switching (few microseconds) but with a potential degradation of the spectral response (i.e., side lobs), and (2), using a dispersing prism which can be mounted on a rotating axis, which can achieve better spectral performances but at the cost of a slower switching time. Therefore, the choice of the instrumentation depends on the end application, and trade-off needs to be made between spectral quality and acquisition speed. Whatever the requirement, we have seen that SCL can be used and coupled with various type of instrumentation in order to accommodate to the requirements of the applications.

From an application point of view, the ability to monitor a large number of several wavelengths has several advantages. First of all, it enables to provides a true spectroscopic analysis, both in terms of absorption and scattering. This is extremely appealing for cancer application, where this type of "optical biopsy", enables to distinguish between healthy and cancerous tissue. This type of application has been extensively used in breast cancer [138,139], and more recently, thyroid applications [50,140]. In this type of application, where the spectral quality is the most important criteria, the second instrumental approach, based on spectral unmixing with a dispersing prism, was used. Secondly, the first instrumental approach, based on spectral unmixing at the detection stage, or the second instrumental approach using ATOF, give the possibility to acquire several wavelengths, either parallelly, or with a fast switching, in a time scale compatible with functional imaging. With these approaches, several authors have explicitly noted that the adoption of SCL enabled them to design or update existing systems, in order to retrieve dynamic physiological information [67,75].

The large number of wavelengths enabled by SCL also give the possibility to go beyond the traditional contrasts, like the hemodynamic one, in order to get information about metabolism for example, by extracting information about more chromophores, like oxCCO [7,24]. We can also note that beyond the intrinsic optical signal, a contrast agent, like ICG, can also be used in order to get information on tissue perfusion. This has notably been extensively used to monitor cerebral perfusion, and the SCL has also proved to be a suitable light source for these applications, even used in extremely complicated clinical settings like ICUs [65].

From a historical point of view, we can note that the development of commercial fibre also largely helped with the adoption of that technology. Indeed, out of the 35 different systems used in the reviewed paper, only 4 (less than 12%) are home-made. They correspond to the system developed before 2010. Past that date, all systems reviewed used a commercial source. We can also see the increase number of papers since the adoption of the commercial sources, with a number of papers published per year roughly doubling from that date. One of the main advantages reported by the authors when they adopted commercial SCL was their relative compactness, compared to the initial development of SC generation. This can be seen in Figure 2, with systems able to be mounted on trolley compatible with the clinical environment. However, the use of these lasers in actual clinical studies remains marginal. In our recent review [141], which has explored the use of TD-NIRS for clinical applications, less than 4% of the 52 papers reviewed were using a system based on a SCL source. Indeed, despite their advantages, a few drawbacks still hindered the adoption of SCL as a standard source. The first one is the cost which, despite a significant decrease in the recent years, remains high (i.e., more than 20,000 pounds for a SCL with several watts of power, enabling a final power of at least 1 mW per wavelength). The second main disadvantage is the warmup time, and the overall stability of these sources. Indeed, the stability of the output can be a concern, but the use of a reference arm can solve this issue [23]. This of course comes at the cost of a more complicated and bulkier instrument. Regarding the warmup time, typically, 1 h is required so the instrument reach thermal equilibrium, which can be too long for certain clinical applications. Some recent TD-NIRS systems based on pulsed laser diode have reported a warmup time of only 15 min, which is way more advantageous in a clinical setting [142]. Finally, even though the novel commercial SCL are more compact than in the early days, their size is still significantly bigger than traditional pulsed laser diode, which limit the reduction of the overall instrument size. Therefore, this can also slow down the adoption of these systems, as space are often very limited in busy clinical environment like ICUs.

These drawbacks are less of an issue for the preclinical settings as the environment is way less constrained. Therefore, the great potential of applications based on SCL outweighs its current main drawbacks, and several instruments based on SCL have been developed for preclinical imaging. The interested reader can refer to the recent review by Bérubé-Lauzière and colleagues, which explore more particularly the use of DOT in preclinical settings [143]. It worth noting that these systems are non-contact, and that this particular detection approach has also been explored for human imaging. Here, having an extremely narrow pulse, like the one provided by SCL, is of great interest as one can take advantage of the null-source approach [77]. Even if these approaches are less common, they can have potential great applications for perioperative measurements, or when the patient is difficult to access, and represent a promising avenue.

The above-mentioned drawbacks are also less problematic in a typical lab setting. SCL are a great tool for teams aiming at instrumental and methodological developments. In the reviewed papers, 33% were aimed at the initial system developments, and 18% and 44% were focused on optical tissue estimations or methodological developments, respectively. Only 5% of the reported studies were focused on pure applications, like tissue oxygenation estimation or cancer monitoring. Indeed, the availability of these SCL sources in teams focused on methodological and instrumental developments have pushed their used to test new possibilities and enhance the reliability of the system by promoting standardization. It is worth mentioning that having a SCL is not mandatory in a lot of settings, and that these sources were used as they were available in most of the groups working on instrumental and methodological developments. Indeed, in the case of a CW system, we have seen that SCL systems were used to provide a high number of wavelengths. In that case, other broadband sources are available, like Tungsten Halogen Light Sources (THLS) or like the Laser-Driven Light Sources (LDLS). The main disadvantage of the THLS is the low power density/wavelength, compared to coherent sources, which can limit their possibilities in biomedical diffuse optics. However, we can still report that THLS were successfully used

to developed simple CW instruments [144,145]. On the contrary, the novel LDLS are of interest as they can produce a stable, high power broadband light. This technology is still relatively new and, to our knowledge, has not been used in biomedical diffuse optics yet. However, its attractive characteristics could be of interest for developing CW biomedical diffuse optics systems. The interested reader can refer to review [146] for more details. Finally, we can note that designing a multi-wavelengths CW system by adding up several single wavelengths light sources, like LEDs, is possible, but it is limited to a low number of wavelengths as the bulkiness and cost of the systems would rapidly increase with the number of wavelengths. In the case of TD systems, the choice of the light source is also largely dictated by the number of wavelengths needed, as can be seen in Table 3, that summarizes the relative advantages and disadvantages of single wavelength laser diode sources versus SCL in two cases: (1) low and (2) high number of wavelengths for TD systems. The use of single wavelength laser diode sources is indeed well adapted when a low number of wavelengths is required, as these sources are compact and relatively not expensive regarding the SCL ones. However, if one wants to design a truly spectroscopic system, with a high number of wavelengths required, these advantages became less evident as, similarly to the CW case, the bulkiness and cost of the system increase as function as the number of wavelengths. This is even more pronounced for the TD cases as the electronics are more complex. Finally, we can mention that if one wants to focus on methodological development, when specific wavelengths are still not defined, the flexibility of SCL is then evident as this choice does not have to be made from the start of the designing process, and not all the wavelengths are available commercially with laser diodes. Therefore, even though other light sources with characteristics compatible with diffuse optics are available, the high power, broadband spectrum, and fast switching capabilities of SCL makes them a source of choice, sufficiently versatile to be able to explore various avenues without source limitations.

**Table 3.** Main relative advantages and drawbacks of single wavelength laser diode sources and SCL as function as the number of wavelengths used for TD systems. (+ advantages /− disadvantages). FW: Few Wavelength, HNW: High number of Wavelengths.


This spectral flexibility offered by SCL has indeed been a true advantage regarding the standardization efforts, giving the ability to optically characterize over a large bandwidth material and tissue. Phantom materials, like Intralipid or Indian ink, have been able to be precisely characterized using SCL, enabling to produce reproducible recipes for liquid phantoms, or to produces reliable solid phantoms. This now gives the possibility to accurately and quantitatively test and compare newly developed instruments and algorithms. We can note the effort undertaken by the BitMap project that have started comparing instruments, from different institutions all over Europe, using the same sets of solid phantoms [124,125].

SCL can also be used as companions to develop new electronical components. Indeed, one of the big instrumental development in the recent year have been the refinement of the detection scheme. Especially, a lot of improvement have been seen in novel solid state detectors [147] which enable to boost light harvesting detectors electronics. The development of gated detector [148] also enables techniques like null source distance to be practically implemented, which can improve the resolution and sensitivity of the

technique [77]. In these works, the availability of SCL sources enabled to perform a precise characterisation of these detectors in situ, both in terms of temporal resolution, thanks to the short pulse width of the SCL, and in terms of spectral responsivity over a large bandwidth. The combination of the SCL capabilities with these new detectors can notably boost the capabilities of TD-NIR instrument, overcoming some current limitations like, for example, the penetration depth [22,149].

The inherent ability of SCL sources to construct systems gathering multidimensional information, i.e., spatial/spectral/temporal, also makes them a perfect tool to test novel algorithms able to exploit this large amount of data. The recent works by Yang and colleagues [117–119] show the advantages of using multi-dimensional datasets, by enhancing the accuracy of the optical properties estimation of tissues by using simultaneously the temporal and spatial dimensions. Kovacsova and colleagues [121] also showed this utility, here combining the spectral and the spatial dimensions. It worth noting that this work is based on broadband CW-NIRS, but that the results of this algorithm has been tested against a TD-NIRS instrument based on SCL as a source. Combining all the dimensions available is certainly a really exiting avenue for the field of diffuse optics, to make the system even more robust and precise.

Finally, we can mention that in the recent years, progress has also been made in individual pulsed laser source, which now achieve several 10th of mW and a pulse width shorter than 100 ps. A good example of the uses of these sources can be found in [150], where an 8 wavelengths TD-NIRS system, compatible for breast imaging for example, has been developed. This paves the way on the use of these sources to enable miniaturisation of TD-NIRS instruments, making them truly compatible with the clinic. These sources could then replace the SCL in clinical settings. However, depending on the application and the chromophores targeted, the wavelength selection is a process that requires some investigations in order to balance the system complexity, and the possible crosstalk effects between the chosen wavelengths. Therefore, systems based on SCL able to easily explore various combinations of wavelengths are once more extremely useful. Lastly, the flexibility of these systems also allows to explore new horizons and applications area, like we have seen with exploring contrast above the classical therapeutic windows (i.e., above 1000 nm) [46], or to explore new organs, like lungs [133]. Therefore, even though the main drawbacks of SCL in terms of absolute cost, size, and stability might prevent their use on a wide scale in the clinical world, their main advantages in terms of power, wavelength choice, and relatively low cost for high number of wavelengths is likely to still push their use in the labs to drive forward the clinical innovation in biomedical diffuse optics.

#### **5. Conclusions**

We have seen that the use of SCL drove the developments in biomedical diffuse optics, and particularly in TD-NIRS. First of all, it is a tool of choice to estimate the optical properties of tissues and phantom material over a large bandwidth. This knowledge is crucial in order to develop a standardized method, and to provide a-priori information used in the calculation of in vivo chromophores concentration for example. Moreover, the high spectral power combined with the ability to easily and quickly select wavelength make SCL suitable not only for in vivo tissue characterization, and application like breast cancer monitoring, but also to follow fast dynamic physiological processes like cerebral hemodynamic. The large number of wavelength available in the systems based on SCL (several systems reported parallel or quasi-parallel acquisition of up to 16 wavelengths) not only enables to refine the precision of the measurement, but also provides extra information about different processes like metabolism. These sources are also relevant to design instrument aiming at non-contact scanning, which make them also a good candidate for preclinical applications. In conclusion, SCL are a valuable tool at every step of the developmental and translational work, from fundamental characterization to preclinical and clinical use, and we have no doubt that it will still be an important brick to keep driving the innovation in biomedical diffuse optics.

**Author Contributions:** F.L., L.G. and I.T. conceived and designed the study. F.L. conducted the literature search and drafted the first version of the paper. L.G. created Tables 1 and 2. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Medical Research Council (MR/S003134/1).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


## *Article* **Measurement of Adult Human Brain Responses to Breath-Holding by Multi-Distance Hyperspectral Near-Infrared Spectroscopy**

**Zahida Guerouah 1, Steve Lin 2,3 and Vladislav Toronov 1,2,\***

	- Toronto, ON M5G 2C4, Canada

**Abstract:** A major limitation of near-infrared spectroscopy (NIRS) is its high sensitivity to the scalp and low sensitivity to the brain of adult humans. In the present work we used multi-distance hyperspectral NIRS (hNIRS) to investigate the optimal source-detector distances, wavelength ranges, and analysis techniques to separate cerebral responses to 30 s breath-holds (BHs) from the responses in the superficial tissue layer in healthy adult humans. We observed significant responses to BHs in the scalp hemodynamics. Cerebral responses to BHs were detected in the cytochrome C oxidase redox (rCCO) at 4 cm without using data from the short-distance channel. Using the data from the 1 cm channel in the two-layer regression algorithm showed that cerebral hemodynamic and rCCO responses also occurred at 3 cm. We found that the waveband 700–900 nm was optimal for the detection of cerebral responses to BHs in adults.

**Keywords:** near-infrared spectroscopy; brain; BOLD signal; breath-holding; cytochrome C oxidase

#### **1. Introduction**

Near-infrared spectroscopy (NIRS) was proposed for human brain measurements in 1970s [1]. It has also been considered for clinical monitoring of cerebral status during various medical conditions in adults, such as cardiac surgeries, cardiac arrest, and traumatic brain injury [2]. However, in such conditions significant circulatory and metabolic changes occur in the entire body, including the scalp, where NIRS sensitivity is maximal [3,4]. In general, high sensitivity to the scalp and low sensitivity to the brain of adult humans remains a main problem of NIRS in spite of numerous attempts to resolve it using different continuous-wave [5], time-domain [6], and frequency-domain approaches [7]. In particular, in several papers the combination of long and short source-detector channels was proposed and investigated [5,8–12]. Other aspects of cerebral NIRS requiring investigation include optimization of the spectrum of wavelengths [13–15], and of the range of source-detector distances [14–16] for specific categories of subjects and patients—children, adults, and seniors. In the present work we approached the above aspects of cerebral NIRS using multi-distance hyperspectral NIRS (hNIRS).

Apart from the ability to measure the optical properties of tissue at all NIR wavelengths simultaneously, hNIRS is the most robust technique to spectrally separate hemoglobin and redox cytochrome C oxidase (rCCO) changes [4,17,18]. An increase in rCCO signal corresponds to an increase in oxidized CCO and an equal decrease in reduced CCO. While rCCO is a direct marker of cellular oxygen metabolism, the robust measurement of rCCO is much more challenging than measurement of the blood parameters due to

**Citation:** Guerouah, Z.; Lin, S.; Toronov, V. Measurement of Adult Human Brain Responses to Breath-Holding by Multi-Distance Hyperspectral Near-Infrared Spectroscopy. *Appl. Sci.* **2022**, *12*, 371. https://doi.org/10.3390/ app12010371

Academic Editor: Hidenao Fukuyama

Received: 26 September 2021 Accepted: 24 December 2021 Published: 31 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the much lower concentration of CCO than of hemoglobin. Therefore, the development and validation of the methodology for rCCO measurements remains a hot topic of NIRS research. Reference [18] provide comprehensive reviews of the physiological significance, and of the history and various aspects of measurements of cerebral rCCO. A review of multiple possible factors that can change rCCO in the brain (changes in oxygen tension, ATP, NADH, etc.) is available in [19].

In this work we used hNIRS to study the possibility to measure specific cerebral autoregulation changes concomitant with systemic changes at source-detector distances of 1 cm, 3 cm, and 4 cm in healthy adults during breath-holding (BH) respiratory challenges. Since BH was proposed as a clinical paradigm to assess cerebral status in various clinical conditions such as stroke, concussion, etc. [20–22], the dynamics of cerebral responses to BH have been studied both by NIRS [23,24] and functional magnetic resonance imaging (fMRI) [25,26], which allows for a direct comparison of our results with other studies.

#### **2. Materials and Methods**

The study was conducted according to the guidelines of the Declaration of Helsinki, and approved by the Research Ethics Board of Ryerson University (REB: 2008-003-1, 4 May 2015).

#### *2.1. Measurement Setup*

An hNIRS custom sensor was placed on the left forehead over the Fp1 position according to the International 10–20 system [27] for the entire measurement procedure (see Figure 1). The distance between two sources was 8 mm (Figure 1b). According to the Monte-Carlo simulation of the light propagation in the adult human head [3] and to the fMRI studies [25,26], with such a geometry the light from both sources collected at 3 cm and 4 cm interrogated the cortical region which was expected to show the same dynamics during BH.

**Figure 1.** hNIRS measurement setup: (**a**) approximate interrogated volumes in the scalp, skull, and brain; (**b**) optical sensor layout; and (**c**) subject position.

The spectra were collected at the sampling rate of 2 Hz by three fiber optic spectrometers— AvaSpec (Avantes, CO, USA) and QE 65000 and USB 4000 (Ocean Optics, Dunedin, FL, USA) at 4, 3, and 1 cm, respectively, to separate the extracerebral and cerebral measurements. The spectrometers were tested on the same phantom and on a human forearm muscle to ensure that at the same distance they measured the same absorbance within the 700–1000 nm wavelength range. AvaSoft-Full software (Avantes, CO, USA) was used to collect data from the AvaSpec spectrometer (4 cm), and The Spectra Suite (Ocean Optics, FL, USA) software was used to collect the broadband continuous-wave hNIRS data from both Ocean Optics spectrometers with dark-signal correction. Data acquired at 1 cm represented the extra-cerebral layer (mostly scalp). Data acquired at 3 cm and 4 cm channels represented a combined extracerebral and cerebral tissue volume. The common spectral range of these spectrometers was from 650 to 1024 nm. The spectrometers at 4 cm and 3 cm had a high signal-to-noise ratio (over 1000:1 single acquisition) and the slit width of 0.5 mm to provide

the sensitivity required to measure light at large distance from the source. Three custommade 2-m-long optical fiber bundles (each made of seven 0.5 NA, 400 μm core-diameter multimode polymer-clad fibers with broad UV/VIS/NIR spectral range of 400 to 2200 nm Thorlabs, NJ, USA) connected spectrometers with the patient's head. Two other optical fiber bundles were used to connect the probe with a halogen lamp light source (Fiber-Lite Dc 950H Fiber Optic Illuminator, Dolan-Jenner, MA, USA). The source light was injected into the tissue at two symmetric scalp locations (Figure 1b) in order to increase the total light power without exceeding the maximum permissible exposure.

#### *2.2. Breath-Holding Paradigm*

Data were obtained from 12 healthy adult participants (seven males, five females, 25–55 years old). All participants gave informed consent before participation and the experiment was performed according to Ryerson University Research Ethics protocol. Participants were audibly cued to perform a 30 s BH at the end of expiration. BH was repeated three times (60–90 s, 180–210 s, and 300–330 s from the beginning of the recording) with 90 s rest intervals. Long time intervals between BHs were used to avoid resonance induction of systemic Mayer waves in arterial blood pressure [28,29]. All data sets including the baseline periods before and after BHs were acquired during 10 min.

#### *2.3. Data Processing*

All our custom signal processing methods were implemented in MATLAB (Mathworks, Natick, MA, USA, Version R2020b). At the pre-processing stage the data from all spectrometers were resampled in both time and spectral domains to the time step of 0.5 s and wavelength step of 1 nm between 650 and 1024 nm, filtered in the time domain using a band-pass filter with the window of 0.01–0.3 Hz, and smoothed in the spectral domain with the median filter of the 20 nm width. To analyze the chromophores, we used two custom data processing approaches. Our first algorithm was based on the analytical solution to the diffusion equation [30,31] for the semi-infinite homogeneous medium. This allows for measuring the bulk absolute concentrations of tissue HbO2 and HHb, and changes in rCCO concentration without accounting for the layered tissue structure. The baseline concentrations of HbO2 and HHb were calculated by performing a non-linear least square fitting [30,31] of the measured absorbance spectrum at each moment of time by the analytical solution to the diffusion equation, in which the optical absorption coefficient *μa*(*λ*) was modeled as

$$
\mu\_{\rm \mu}(\lambda) = [\rm Hb]\varepsilon(\lambda)\_{\rm Hb} + [\rm HbO\_2]\varepsilon(\lambda)\_{\rm HbO\_2} + \eta(\lambda) \tag{1}
$$

According to Beer–Lambert law, the reduced scattering coefficient *μ<sup>s</sup>* (*λ*) as a function of wavelength *λ* was modeled using the power law as described in [30]. In Equation (1) and below the square brackets denote concentrations measured in micromoles (μM), *η*(*λ*) is the absorption by water assuming 80% water concentration. The temporal changes in the hemoglobin concentrations HbO2, Hb, and rCOO were resolved using a two-step data-fitting algorithm [30,31] (also based on the same analytical solution to the diffusion equation) by relating changes in HbO2, Hb, and rCOO to the changes in the optical absorbance as

$$
\Delta\mu\_a(\lambda, t) = \Delta[\text{Hb}](t)\varepsilon(\lambda)\_{\text{Hb}} + \Delta[\text{HbO}\_2](t)\varepsilon(\lambda)\_{\text{HbO}\_2} + \Delta[\text{COO}](t)\varepsilon(\lambda)\_{\text{CO}},\tag{2}
$$

where *ε*(*λ*)*<sup>x</sup>* were the spectra of the extinction coefficients of HbO2, Hb, and rCCO (see Figure 2) [30,31].

**Figure 2.** Absorption spectra of HHb, HbO2, H2O, and rCCO.

The data-fitting was performed in two steps: first Δ[HbO2] and Δ[HHb] were calculated assuming Δ[CCO] = 0, and after that Δ[CCO] was calculated and retained non-zero only if the addition of Δ[CCO] resulted in an improvement in the fit quality of at least 20% in terms of the norm of residuals. Such a stepwise regression worked as an adaptive filter reducing the noise in the time-course of rCCO without suppressing changes recognized from significant spectral distortions corresponding to *ε*(*λ*)cco. Additional details on recovering the absolute values and changes of chromophore concentrations from the hNIRS data can be found in [32,33].

As the fraction of oxygenated hemoglobin relative to the total hemoglobin in the blood, the cerebral tissue saturation of oxygen was calculated as:

$$\text{tfSO}\_2 = \frac{\text{[HbO}\_2\text{]}}{\text{[HbO}\_2\text{]} + \text{[Hb]}} \text{(\text{\textquotedblleft}o\text{)}.} \tag{3}$$

In order to assess differences between the cerebral and scalp responses to BH we used another common NIRS model—the modified Lambert–Beer law for a two-layer medium [32,33]. This model assumes that at any moment of time, *t*, the absorbance measured at 1 cm and at the wavelength λ can be related to the changes of the absorption coefficient of the scalp only Δ*μa*,*s*(*λ*, *t*):

$$
\Delta OD\_{1\text{cm}}(\lambda, t) = L\_{\text{s}\ 1\text{cm}} \Delta \mu\_{\text{a}, \text{s}} \ (\lambda, t), \tag{4}
$$

where *Ls* 1cm is the optical pathlength in the scalp at 1 cm, and the absorbance at 3 cm and 4 cm could be expressed as

$$
\Delta OD\_{\text{Sym}}(\lambda, t) = L\_{\text{c3cm}} \Delta \mu\_{a,s}(\lambda, t) + L\_{\text{c2cm}} \Delta \mu\_{a,c}(\lambda, t) \tag{5}
$$

$$
\Delta OD\_{\text{4cm}}(\lambda, t) = L\_{\text{s 4cm}} \Delta \mu\_{a,s}(\lambda, t) + L\_{\text{c 4cm}} \Delta \mu\_{a,c}(\lambda, t) \tag{6}
$$

where Δ*μa*,*c*(*λ*, *t*) is the change in the cerebral absorption coefficient, *Ls* 3cm and *Ls* 4cm are the scalp partial pathlength at 3 cm and 4 cm, respectively, and *Lc* 3cm and *Lc* 4cm are the cerebral partial pathlength at 3 cm and 4 cm, respectively. Since the pathlengths in Equations (4)–(6) were specific for every individual and unknown, Δ*μa*,*<sup>s</sup>* could not be algebraically excluded to find Δ*μa*,*c*. However, since Equations (4) and (5) were linear with respect to the common time-dependent Δ*μa*,*s*(*λ*, *t*) and Δ*μa*,*c*(*λ*, *t*), for every wavelength, *λ*, one could exclude the time-course common with Δ*μa*,*s*(*λ*, *t*) from Δ*OD*3cm(*λ*, *t*) and Δ*OD*4cm(*λ*, *t*). Indeed, one could rewrite Equations (5) and (6) in the form

$$
\Delta O D\_{\text{3cm}}(\lambda, t) = L\_{\text{c.3cm}} \delta \mu\_{\text{a.c}}(\lambda, t) + (L\_{\text{c.3cm}} + L\_{\text{s.3cm}}) \left(L\_{\text{(s.1cm)}}\right)^{-1} \Delta O D\_{\text{1cm}}(\lambda, t) \tag{7}
$$

$$
\Delta OD\_{\text{4cm}}(\lambda, t) = L\_{\text{c 4cm}} \delta \mu\_{\text{c}, \text{c}}(\lambda, t) + \left(L\_{\text{c 4cm}} + L\_{\text{s 4cm}}\right) \left(L\_{\text{(s 1cm)}}\right)^{-1} \Delta OD\_{\text{1cm}}(\lambda, t) \tag{8}
$$

where

$$
\delta \mu\_{a,\varepsilon}(\lambda, t) = \Delta \mu\_{a,\varepsilon}(\lambda, t) - \Delta \mu\_{a,\varepsilon}(\lambda, t) \tag{9}
$$

represents the difference between the cerebral and scalp responses to BH. Note that Equations (7) and (8) are linear with respect to Δ*OD*1cm(*λ*, *t*) and *δμa*,*c*(*λ*, *t*), and also we expect that the latter functions of time and wavelength will be uncorrelated in the time domain. Therefore the time-domain linear regression coefficients between Δ*OD*1cm(*λ*, *t*) and Δ*OD*3cm(*λ*, *t*) and between Δ*OD*1cm(*λ*, *t*) and Δ*OD*4cm(*λ*, *t*) should be close to *β*1,3 cm ≈ (*Lc* <sup>3</sup>*cm* + *Ls* <sup>3</sup>*cm*) - *L*(*<sup>s</sup>* <sup>1</sup>*cm*) −<sup>1</sup> and *β*1,4 cm ≈ (*Lc* <sup>4</sup>*cm* + *Ls* <sup>4</sup>*cm*) - *L*(*<sup>s</sup>* <sup>1</sup>*cm*) −<sup>1</sup> , respectively. Then, by subtracting *β*1,3 cm(*λ*)Δ*OD*1cm(*λ*, *t*) and *β*1,4 cm(*λ*)Δ*OD*1cm(*λ*, *t*) from Equations (7) and (8) we obtain:

$$
\delta OD\_{3cm}(\lambda, t) \approx L\_{\text{c } 3cm} \delta \mu\_{a\varepsilon}(\lambda, t), \text{ and } \delta OD\_{4cm}(\lambda, t) \approx L\_{\text{c } 4cm} \delta \mu\_{a\varepsilon}(\lambda, t) \tag{10}
$$

Note that at all wavelengths the true cerebral response to BH could have a component close to the scalp response. Therefore, *δOD*3,4 cm(*λ*, *t*) represent not the pure time-course of the cerebral absorption coefficient, but rather the difference between the cerebral and scalp absorption time-courses.

In accordance with Equation (2) one could model *δOD*3,4 cm(*λ*, *t*) as a linear combination of the contributions from *δ*[Hb](*t*), *δ*[HbO2](*t*), and *δ*[rCCO](*t*):

$$\delta O D\_{3,4\text{ cm}}(\lambda, t) \sim \delta [\text{Hb}](t) \varepsilon(\lambda)\_{\text{Hb}} + \delta [\text{HbO}\_2](t) \varepsilon(\lambda)\_{\text{HbO}\_2} + \delta [\text{rCOO}](t) \varepsilon(\lambda)\_{\text{rCOo}}.\tag{11}$$

By applying the stepwise linear regression spectral unmixing to the linear model (Equation (9)) (where the regressors are the extinction coefficients *ε*(*λ*)*x*) one could obtain the specific cerebral changes *δ*[Hb](*t*), *δ*[HbO2](*t*), and *δ*[rCCO](*t*), whose time-courses represent the differences between the cerebral and scalp responses in HHb, HbO2, and rCCO, respectively. In this spectral unmixing step the initial regression model did not include *δ*[rCCO](*t*). It was added to the model at the second step of the stepwise regression only if the *p*-value corresponding to non-zero *δ*[rCCO](*t*) was greater than 0.05.

In addition, we also used Equation (10) to calculate the specific cerebral response in the total hemoglobin *δ*[*t*Hb]. However, *δ*[*t*Hb] was calculated not by using the spectral unmixing, but by directly using Equation (10) around 800 nm since at this wavelength *ε*Hb = *ε*HbO2 . Since near 800 nm *ε*(*λ*)*r*cco is greater than *ε*Hb, *δ*[*t*Hb] could include a crosstalk with *δ*[rCCO](*t*).

Since the values of cerebral pathlengths *Lc* 3cm and *Lc* 4cm were unknown, for the quantitative comparison in μM with the responses obtained using the homogeneous tissue model (explained above) we normalized *δ*[Hb](*t*), *δ*[HbO2](*t*) , and *δ*[rCCO](*t*) to the pathlength corresponding to the homogeneous tissue equal to the product of the differential pathlength factor (known be close to 6 for the mid-age human forehead [34]) and the sourcedetector distance.

#### *2.4. Statistical Analysis*

In the individual subject data the peak magnitudes of responses in Δ[Hb], Δ[HbO2], and Δ[CCO], and in *δ*[Hb], *δ*[HbO2], and *δ*[CCO] were measured as the difference between the time average over ±5 s around each peak time (see Table 1) and the average over 20 s just before BHs. These samples were checked for normality (positively) and further tested

by the one-sample *t*-test (*ttest* MATLAB function) and *p*-values reported in Tables 1 and 2. In addition, the differences in response magnitudes among three BH episodes were analyzed using *ranova* MATLAB function (repeated measures analysis of variance, RM ANOVA).

**Table 1.** Cross-subject average peak magnitudes, *p*-values, and peak times obtained from the homogeneous semi-infinite diffusion analysis. The peak magnitudes were measured from the 20 s averaged values before BHs. The peak times measured from the beginning of BHs were found in the averaged traces shown in Figure 4 with the errors estimated as the mean absolute deviation from the individual data.


**Table 2.** Cross-subject average peak magnitudes and *p*-values obtained from the two-layer regression analysis. The peak magnitudes were measured from the 20 s averaged values before BHs.


#### **3. Results**

Figure 3 shows the assessment of de-noised signals in all three channels and at all wavelengths in terms of the cross-subject averaged temporal correlation (first and second rows), temporal logarithmic signal means (third row, used for the data quality assessment), and temporal standard deviations (fourth row).

**Figure 3.** Hyperspectral signal analysis: cross-subject averaged correlation, temporal means, and temporal standard deviations for signals from 1 cm, 3 cm, and 4 cm channels. Blue color shows the cross-subject standard deviations.

Since the noise was filtered out, the latter represent the average magnitude of changes at different wavelengths. The correlation coefficient maps (first row) show that signals at all wavelengths from 700 nm to over 900 nm and in all channels were positively correlated. The latter means that there were no pairs of signals at any wavelengths from 700 nm to 900 nm and in any channels where temporal changes occurred simultaneously in opposite directions. Also from the correlation maps and same-wavelength correlation plots (second row) one could see that all signals between 750 nm and 900 nm were highly correlated (correlation coefficient > 0.6). This high correlation occurred because, as shown in [3], at both 3 and 4 cm the optical pathlength in the extracerebral tissue was several times longer than in the brain, and therefore the scalp contribution was high at both 4 and 3 cm channels. Correlations were low for the wavelengths longer than 900 nm in the 4 cm channel. While the mean signal value in the 4 cm channel at these wavelengths was high (third row), the amplitude of changes (fourth row) was small. This lack of changes should be due to fact that at longer wavelengths the partial optical pathlengths in the tissues where the changes occurred (scalp and brain) was shorter than at wavelengths below 900 nm, and at the 4 cm source-detector separation the light of wavelengths greater than 900 nm mostly interrogated the volume which belongs to the skull where no changes occurred. For the above reasons the further analysis was performed using the waveband from 750 nm to 900 nm.

Figure 4 shows the cross-subject averaged changes in [HbO2], [HHb], [rCCO], and [rSO2] measured at 1 cm, 3 cm, and 4 cm using the homogeneous tissue model and 750–900 nm waveband. One could see that the hemoglobin signals showed responses to BHs in all channels, while rCCO exhibited clear responses only at 4 cm. Table 1 presents the average amplitude, p-values, and temporal values characteristic for the chromophore responses to BHs. In addition to Figures 4 and 5 provides a direct comparison of the hemoglobin time-courses at 1 cm and 4 cm and shows a return to the baseline 4 min after BHs.

As shown in Figures 4 and 5, and also from the peak times for HbO2 in Table 1 the time-courses of HbO2 responses appeared similar in all channels. The common features of HHb responses to BHs in 1 cm and 3 cm channels were the initial dips followed by quick increases peaking up 4–10 s after the end of BH. The depth of initial dips was much greater in the 1 cm channel than in 3 cm and 4 cm channels (Figure 5b). However, the magnitudes of these dips measured from the HHb levels just before BH were statistically insignificant (*p* > 0.05). The HHb peak times were longest for the 1 cm channel (about 9 s after the end of BH, see Table 1). Much earlier peaks and faster falls in the HHb responses at 4 cm compared to those features at 1 cm (Figure 5b, Table 1) might indicate a prevalence of the cerebral component in the 4 cm HHb responses. In particular, the faster falls of HHb at 4 cm than at 1 cm could be due to the cerebral autoregulation restoring the oxygen concentration in the brain faster than in extracerebral tissues. In Figure 4 one can also see that at 3 cm the HHb responses were intermediate between those measured at 1 cm and 4 cm.

The rCCO signals from the 1 cm and 3 cm channels showed no clear responses to BHs. In contrast, at 4 cm the rCCO responses were very clear (Figures 4 and 5c), which also might indicate that these responses belonged solely to the cerebral tissue. Each rCCO fall after the end of BH was quickly followed by a rebound which also could be a result of the faster restoration of the oxygen supply in the brain than in other tissues. Note that rCCO peak changes had about 10 times smaller magnitudes than HbO2 and HHb changes (Table 1), which was in agreement with the fact that the concentration of CCO in the brain is approximately 10 times smaller than the concentration of hemoglobin.

**Figure 4.** Cross-subject averaged changes in [HbO2], [HHb], and [rCCO], and regional SO2 measured at 1 cm, 3 cm, and 4 cm using the homogeneous tissue model and 750–900nm waveband. Blue color shows the cross-subject standard deviations.

**Figure 5.** Cross-subject averaged temporal traces of [HbO2] (**a**), [HHb] (**b**), and [rCCO] (**c**) showing a comparison of 1 cm and 4 cm time-courses and a return to the baseline 4 min after BHs.

In Figure 4 the SO2 time-course at 4 cm was also very different from those at 1 cm and 3 cm, which were very similar to each other by showing SO2 increases characteristic for the scalp tissue. However, the quick SO2 falls at 4 cm at the end of each BH should not be interpreted as the failure of cerebral autoregulation but rather, these falls resulted from the smaller magnitude of the HbO2 increase than of the HHb increase due to the difference in the partial volume of the skull interrogated by light at shorter and longer wavelengths. From Table 1 one can see that the magnitudes of HbO2 and HHb responses to BHs decreased with the distance. The reason for this was the application of the homogeneous model to the measurements of the inhomogeneous tissue, which included bone, where almost no blood was present. With the increase of the distance the partial volume of the bone increased leading to the underestimation of the magnitudes of HbO2 and HHb responses to BHs at 3 cm and 4 cm.

The cumulative *p*-values in Table 1 show that all included responses were statistically significant at *p* < 0.05 when the results of all three BHs were combined. HbO2 and HHb responses to each individual BH were significant at 1 cm and at 3 cm. Unlike HbO2 and HHb, at 4 cm, rCCO responses to each individual BH were statistically significant. Although the average peak values for different BHs in Table 1 were different, RM ANOVA analysis did not show statistical significance of these differences.

As explained above, in order to reveal differences between responses to BHs measured in different channels we performed time-spectral linear regression analysis. For the regression analysis the range of wavelengths was limited to 750–900 nm where the highest correlation between different channels was detected (see Figure 3). Figure 6 shows the cross-subject average correlation and standard deviations spectra of signals from the 3 cm and 4 cm channels after regressing out changes measured at the 1 cm channel. One could see that the correlation coefficient for all wavelengths in both channels was greater than 0.6 between 750 nm and 900 nm. The correlation coefficient for the same wavelengths in both channels (corresponding to the diagonal of the correlation map in Figure 6) was particularly high for all wavelengths. This high correlation indicated that after regressing out the superficial changes, both 3 cm and 4 cm channels exhibited very similar time-courses at all wavelengths. Also in Figure 6 one can see that the spectra of average amplitudes of changes (measured by SDs) were similar in shape to those before regression (shown in Figure 3). They showed peaks at the HHb and rCCO absorption maxima around 760 nm and 830 nm, respectively, (compare with Figure 2), and decreased at wavelengths longer than 870 nm, corresponding to maximum HbO2 absorption. Therefore, after the regression the HbO2 changes should have smaller magnitudes than the HHb and rCCO responses.

**Figure 6.** Cross-subject average correlation and standard deviation spectra of signals from the 3 cm and 4 cm channels after regressing out changes measured at the 1 cm channel.

In Figure 7 the green and blue curves show the cross-subject averaged changes in HbO2, HHb, tHB, and rCCO calculated by applying the spectral-domain regression to the residuals after regressing out the in time-domain the changes measured at the 1 cm channel from the signals measured at the 3 cm and 4 cm channels, respectively. The tHB changes were computed by averaging changes in the 5 nm waveband around 800 nm where the HHb- and HbO2-specific absorptions were equal. Note that the changes in Figure 7 show the differences in the time-courses of changes in the deep and superficial tissues. (Note that signals measured at different channels never changed in opposite directions, which follows from the correlation analysis shown in Figure 3 and explained above.) In particular, the rising slopes seen on the curves in Figure 7 correspond to the time intervals when changes measured at 3 or 4 cm increased faster than changes measured at 1 cm. Conversely, the curves in Figure 6 show negative slopes during the times when changes measured at long-distance channels decreased faster than changes measured at the 1 cm channel. The fact that the HbO2 responses to BHs were not very clear means that the time-courses of [HbO2] changes in the superficial and cerebral tissues were similar. (The "noise" seen in Figure 7 was mostly due to the respiratory sinus arrhythmia between BHs). In Figure 7 one can also see that the peak times of [HHb], [tHB], and [rCCO] responses during all three BH episodes were between 15 and 25 s after the BH onsets, which was similar to the peak times of rCCO responses shown in Figure 4, but very different from the peak times of HbO2 and HHb responses shown in Figure 4, which never occurred earlier than 26 s from the beginning of a BH episode (see also Table 1). Some details of the HHb, tHB, and rCCO responses in Figure 7 were quite different. In particular, during BH episodes, the HHb curves showed initial dips during the first 10–15 s turning into the fast increases after that. During the same initial periods of BHs the tHB curves showed steady positive changes, while rCCO curves showed almost no changes during the first 10 s. These differences indicated that the responses to BHs in Figure 7 indeed corresponded to different chromophores and did not result from a crosstalk between different NIR wavelengths. The fact that in Figure 7 the curves corresponding to the 4 cm channel showed similar changes to the curves corresponding to the 3 cm channel but with much smaller magnitudes was due to the overall smaller magnitudes of changes at 4 cm as shown in Figure 3 and explained above.

In addition, the red curves in Figure 7 correspond to regressing out the 3 cm changes from the 4 cm changes. The fact that these curves show almost no changes means that after regressing out the changes due to the contribution from the superficial tissues, the time-courses of changes at all wavelengths in both 3 cm and 4 cm channels were very similar, which was also in agreement with the high correlation between signals from these channels at all wavelengths as shown in Figure 6.

Table 2 shows that HHb-, tHb-, and rCCO-specific cerebral responses to each individual BH were significant at both 3 cm and at 4 cm. rCCO responses showed highest significance (*p* < 0.03) among all chromophores. HbO2-specific cerebral responses to individual BHs were not significant at 3 cm or at 4 cm. RM ANOVA analysis did not reveal statistical significance of differences in responses to different BH trials.

**Figure 7.** Cross-subject averaged changes in [HbO2], [HHb], [tHB], and [rCCO] calculated from the time-domain regression residuals using stepwise linear regression in spectral domain.

#### **4. Discussion**

Our results showed that BHs induced significant hemodynamic responses in all three channels including the 1 cm channel, which mainly interrogated the scalp and partially the skull. Hemodynamic responses measured by the 3 cm channel were quite similar to those measured at 1 cm. In Figures 4 and 5 only HHb and rCCO signals measured at 4 cm showed different responses to those of the extracerebral tissue. The clearest difference from the scalp response was observed at 4 cm in the rCCO response (Figures 4 and 5). No rCCO responses were observed at 1 cm and 3 cm using the homogeneous tissue model. Clear HHb-, tHB-, and rCCO-specific cerebral responses seen in Figure 7 and Table 2 indicate that the dynamics of changes measured at 3 cm and 4 cm channels were different from that measured at 1 cm, and this difference can be attributed to the differences between the brain and scalp responses to BHs. Note that hemodynamic responses at 3 cm and 4 cm shown in Figure 7 and Table 2 had significantly smaller magnitudes than those shown in Figures 4 and 5 and in Table 1. However, rCCO responses in both Figures 4 and 7 had their magnitudes close to 0.1 μM. This result supports the conclusion that we succeeded in decoupling rCCO changes from changes in other chromophores, and that the measured rCCO changes occurred in the brain. We would like to underline that in both homogeneous and two-layer algorithms we used the most conservative treatments of rCCO changes such as thresholding by the quantitative improvement in the goodness of fit, which suppressed spurious rCCO changes.

The time-courses of cerebral responses to BHs in healthy adult humans were studied in detail using blood-oxygen-level-dependent (BOLD) fMRI in [25,26]. The quantitative relationship between the BOLD changes measured by fMRI and hemodynamic signals measured by NIRS was investigated in [35]. Figure 8 (adopted from [25]) shows averaged BOLD signal time-courses with a 15 s BH after the full expiration on the areas supplied by the middle cerebral artery (MCA, red), anterior cerebral artery (ACA, blue) and posterior cerebral artery (PCA, black). Since our hNIRS probe was positioned on the left side of the forehead near the midline, it could interrogate cortical regions supplied by both ACA and MCA. The main difference between ACA and MCA regional responses was the early rise of the BOLD signal response in the ACA region and a 10 s delay in the MCA region. In our results the initial dynamics (during the first 15 s of BH) of cerebral HbO2 and HHb responses rather corresponded to the MCA type of the BOLD response [25,33]. A delayed rCCO response after the regression (Figure 7) corresponded to the delayed ACA response. Such delayed cerebral responses were most different from the scalp hemodynamic responses to BHs measured at 1 cm (Figure 4), which were also very significant.

**Figure 8.** (Adapted from [23]) Averaged BOLD signal time-course with breath holding of 15 s after expiration on areas supplied by the MCA (red), ACA (blue), and PCA (black).

Another significant difference between the long- and short-distance channels was the longer peak times of HHb responses in the short-distance channel (longer than BH duration, Table 1). The peak times of responses to BHs after expiration measured by hNIRS at 4 cm were close to those measured by BOLD fMRI in [26] (see Figure 2 in [26]).

Our results show that without using a short-distance channel to remove biasing by the scalp, the cerebral responses to BHs were clearest in the rCCO time-course at 4 cm. Some features different from the scalp responses could also be seen in the HHb time-course at 4 cm. At 3 cm, some features of the cerebral response could be poorly recognized in the HHb time-course. The HbO2 time-course did not show clear cerebral responses to BH at 3 cm or at 4 cm, in particular due to the high inter-subject variability. However, measurement of rCCO changes required hyperspectral or many-wavelength multispectral technology.

We found the waveband 750–900 nm be optimal for the detection of cerebral responses to BHs in adults. The wavelength that exhibited the largest fractional change in the detected optical signal with respect to the baseline value both at 3 and 4 cm was between 800 and 850 nm, which was in good agreement with the results of [13] (830 nm).

#### **5. Conclusions**

Cerebral changes could be efficiently separated from the extra-cerebral biasing using the time-domain linear regression of hNIRS signals measured at 3 cm and 4 cm sourcedetector separations by the signals from the short-distance channel. Without using the short-distance channel, none of the signals measured at 3 cm clearly reflected cerebral changes. At 4 cm source-detector separation, cerebral changes could be detected without using short-distance channels in the [HHb] and [rCCO] time-courses. The clearest cerebral responses were detected in the [rCCO] time-course at 4 cm. The optimal waveband for the rCCO measurements was between 750 nm and 900 nm. The [HbO2] time-course at all source-detector distances was most dominated by the extracerebral biasing and therefore it was least suitable for cerebral signal detection. Further clinical studies are required to investigate the ability of the NIRS breath-holding paradigm to detect cerebral circulation disorders.

**Author Contributions:** Conceptualization, S.L. and V.T.; methodology, V.T.; software, V.T.; validation, S.L., V.T. and Z.G.; formal analysis, Z.G.; investigation, Z.G.; resources, V.T.; data curation, V.T.; writing—original draft preparation, V.T.; writing—review and editing, S.L.; visualization, Z.G.; supervision, V.T.; project administration, V.T.; funding acquisition, V.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** The study was conducted according to the guidelines of the Declaration of Helsinki, and approved by the Research Ethics Board of Ryerson University (REB: 2008-003-1, 4 May 2015).

**Informed Consent Statement:** Informed consent was obtained from all subjects involved in the study.

**Acknowledgments:** This research was supported by the Dean's Research Fund, Faculty of Science, Ryerson University.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


### *Article* **Dual-Slope Diffuse Reflectance Instrument for Calibration-Free Broadband Spectroscopy**

**Giles Blaney \*, Ryan Donaldson †, Samee Mushtak †, Han Nguyen †, Lydia Vignale †, Cristianne Fernandez, Thao Pham, Angelo Sassaroli and Sergio Fantini**

> Department of Biomedical Engineering, Tufts University, Medford, MA 02155, USA; Ryan.Donaldson@tufts.edu (R.D.); Samee.Mushtak@tufts.edu (S.M.); Han.Nguyen@tufts.edu (H.N.); Lydia.Vignale@tufts.edu (L.V.); Cristianne.Fernandez@tufts.edu (C.F.); Thao.Pham@tufts.edu (T.P.); Angelo.Sassaroli@tufts.edu (A.S.); sergio.fantini@tufts.edu (S.F.)

**\*** Correspondence: Giles.Blaney@tufts.edu

† These authors contributed equally to this work.

**Abstract:** This work presents the design and validation of an instrument for dual-slope broadband diffuse reflectance spectroscopy. This instrument affords calibration-free, continuous-wave measurements of broadband absorbance of optically diffusive media, which may be translated into absolute absorption spectra by adding frequency-domain measurements of scattering at two wavelengths. An experiment on a strongly scattering liquid phantom (milk, water, dyes) confirms the instrument's ability to correctly identify spectral features and measure absolute absorption. This is done by sequentially adding three dyes, each featuring a distinct spectral absorption, to the milk/water phantom. After each dye addition, the absorption spectrum is measured, and it is found to reproduce the spectral features of the added dye. Additionally, the absorption spectrum is compared to the absorption values measured with a commercial frequency-domain instrument at two wavelengths. The measured absorption of the milk/water phantom quantitatively agrees with the known water absorption spectrum (*R*<sup>2</sup> = 0.98), and the measured absorption of the milk/water/dyes phantom quantitatively agrees with the absorption measured with the frequency-domain instrument in six of eight cases. Additionally, the measured absorption spectrum correctly recovers the concentration of one dye, black India ink, for which we could accurately determine the extinction spectrum (i.e., the specific absorption per unit concentration). The instrumental methods presented in this work can find applications in quantitative spectroscopy of optically diffusive media, and particularly in near-infrared spectroscopy of biological tissue.

**Keywords:** broadband diffuse reflectance spectroscopy; frequency-domain near-infrared spectroscopy; dual-slope; absorption spectra

#### **1. Introduction**

Diffuse optics is concerned with the propagation of light in highly-scattering or diffusive media. One notable application is Near-Infrared Spectroscopy (NIRS) of biological tissue, which is typically performed in the wavelength range from about 600 nm to about 1000 nm [1–3]. Diffuse optics finds a variety of applications in several fields of study. In food science [4], it may be applied for inspection [5] or evaluation [6]. In pharmaceutical manufacturing, diffuse optics may join other process analytical technologies, for example to analyze and characterize particles [7] or powders [8]. A few other examples include archaeological soil analysis [9], dendrology (study of wood) [10], and art authentication analysis [11]. In the study of biological tissue, diffuse optics finds applications in basic research, medical diagnostics, and physiological monitoring. Examples include clinical brain monitoring [12], the study of brain activation [13], breast imaging [14], and muscle measurements in sports science [15]. This is by no means an exhaustive list. The work

**Citation:** Blaney, G.; Donaldson, R.; Mushtak, S.; Nguyen, H.; Vignale, L.; Fernandez, C.; Sassaroli, A.; Pham, T.; Fantini, S. Dual-Slope Diffuse Reflectance Instrument for Calibration-Free Broadband Spectroscopy. *Appl. Sci.* **2021**, *11*, 1757. https://doi.org/10.3390/ app11041757

Academic Editor: Johannes Kiefer

Received: 29 December 2020 Accepted: 12 February 2021 Published: 16 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

that we present in this article aims to improve the accuracy and robustness of quantitative broadband spectroscopy of optically turbid media, and it is relevant to a variety of applications of diffuse optics, even though we mostly focus on its role in the area of biomedical optics.

Quantitatively, light propagation in diffusive media is characterized by a lower probability of absorption (related to the absorption coefficient; *μa*) compared to the probability of effectively isotropic scattering (related to the reduced scattering coefficient; *μ <sup>s</sup>*). These two optical properties are the chief quantities of interest in the field of diffuse optics and diffuse biomedical optics [2,3,16,17], and they feature a wavelength dependence that is often of crucial importance. In the case of most biological tissues, the dominating scattering condition of diffuse optics is realized in the NIR optical window. Measurements of the wavelength dependent absorption of tissue yields information about the concentration of chromophores with known extinction spectra, while measurements of scattering spectra yields structural information related to the size and density of scattering centers [2].

These properties are often obtained in the reflectance geometry, where light is delivered onto the sample surface and detected at some distance away from the source on the same surface. The primary difficulty in these measurements is the decoupling of absorption and scattering contributions to the measured optical signal. One method to achieve this is by using a source light with a temporally varying intensity. These time-resolved methods fall into two categories: time-domain methods where the light emission has an impulse profile, and Frequency-Domain (FD) methods where the light emission has a sinusoidal profile [1,18]. Focusing on FD-NIRS, the measurement of both absorption and scattering can be achieved via a calibrated measurement of the spatial dependence (as a function of the distance from the source on the tissue surface) of the amplitude and phase of the detected modulated intensity (which may be represented as a complex reflectance) [19–21].

This spatially dependent measurement of FD-NIRS data to achieve absolute optical properties is often summarized as measuring the slope (versus source-detector distance; *ρ*) of a modified intensity amplitude and phase (i.e., modified complex reflectance). The chief complication of this technique is calibration, since each source-detector pair may have differing instrumental contributions (because of the individual source emission and detector sensitivity properties) and differing coupling factors between source/detector and sample. A method which compensates for the differing instrumental contributions is the so-called Multi-Distance (MD) scanning, where either the source or the detector is moved across the surface of the optical medium [22]. In doing so, the instrumental factors remain the same for each distance (same source and same detector used for each distance); thus, the measurement of slope is independent of these factors. However, differing coupling factors may still be present if during the scan the source (or detector) moves away or toward the tissue surface (in a non-contact case), or experiences a variable contact pressure (in a contact case). Aside from MD scanning, a potentially more effective and more easily implemented calibration free method may be employed. This method is the Self-Calibrating (SC) method, where a symmetric optical probe that features two sources and two detectors is used such that instrumental and coupling factors cancel in the calculation of slopes [23]. This method is also effective at suppressing motion artifacts [24] and allows calibrationfree saturation measurements with Continuous-Wave (CW) methods after assuming the wavelength dependence of scattering [25]. The original purpose of the SC method was to achieve measurements of absolute optical properties with FD-NIRS. However, recently the optode geometry inspired by the SC method has been extended to the use of only the amplitude or phase of the complex reflectance in FD-NIRS [26,27]. This extension was named Dual-Slope (DS) since it relies on the average of two slopes measured with this special optode configuration.

One of the limitations of time-resolved NIRS methods, such as FD-NIRS, is the added instrumental complexity compared to CW. There has been work in broadband time-resolved spectroscopy [28]; however, such instruments require complex instruments (super-continuum lasers, single photon counting, and instrument response calibration).

Because of this, the norm in time-resolved NIRS methods is the use of two or a few wavelengths as opposed to a continuous broadband spectrum [1]. In contrast, CW methods may use White-Light or broadband light sources, such as halogen lamps, and spectrometer detectors, since there is no need for measurements of fast temporal characteristics. These spectroscopic methods in CW are so-called Diffuse Reflectance Spectroscopy (DRS) due to their collection geometry and use of a spectrometer as a detector (instead of avalanche photodiodes or photomultiplier tubes typical in time-resolved methods). Therefore, CW broadband DRS (CW-bDRS) has the advantage of collecting data over many wavelengths but the disadvantage that absorption and scattering are coupled and inseparable due to the use of CW illumination. A solution that has become more and more common is to combine a time-resolved NIRS instrument at few wavelengths with a CW-bDRS system [29–34]. Such a technique allows extrapolation of scattering from few to many wavelengths, by assuming a power law decay of the reduced scattering coefficient with wavelength [2], as well as decoupling of absorption from scattering in the CW-bDRS data.

At this point, a small note on nomenclature is valuable. The acronym MD FD-NIRS refers to a frequency-domain method capable of measuring absolute optical properties (absorption and scattering) at discrete wavelengths from measurements at multiple distances between source and detector. The acronym CW-bDRS refers to a continuous-wave technique capable of measuring wavelength-resolved data (which depends on absorption and scattering) over a range of many wavelengths. Further, to achieve a calibration free full spectral technique, DS CW-bDRS is introduced. It is noted that the distinction between SC (self-calibrating) and DS (dual-slope) is the use of complex reflectance in SC FD-NIRS versus continuous-wave reflectance in DS CW-bDRS, for example. The main focus of this work will be to develop the DS CW-bDRS technique in combination with MD FD-NIRS to decouple absorption and scattering contributions to the absorbance spectra.

In this paper, a new DS CW-bDRS/MD FD-NIRS instrument is described and validated by measuring absolute absorption spectra of highly scattering media, or phantoms. The novelty of this work is the design of a CW instrument for broadband spectroscopy that implements DS methods for robust and calibration free measurements of absorbance spectra. The addition of MD FD-NIRS measurements at two wavelengths allows for the translation of absorbance spectra into quantitative absorption spectra. The organization of the paper is as follows. The methods, in Section 2, are split into techniques, in Section 2.1, where FD-NIRS and CW-bDRS are described, and is followed by experiments, Section 2.2, which lays out a phantom experiment, and then by an analysis, Section 2.3, which explains methods to calculate optical properties. Next, the results, Section 3, report measured phantom absorption spectra. Finally, the discussion, Section 4, elaborates on the validation of the DS CW-bDRS/MD FD-NIRS instrument.

#### **2. Materials and Methods**

#### *2.1. Techniques*

2.1.1. Frequency-Domain Near-Infrared Spectroscopy

Frequency-Domain Near-Infrared Spectroscopy (FD-NIRS) was implemented with the purpose of measuring absolute optical properties of highly scattering media. Those being the absorption coefficient (*μa*) and the reduced scattering coefficient (*μ <sup>s</sup>*). A commercial FD-NIRS instrument was used (Imagent V2, ISS, Champaign, IL, USA) operating with a modulation frequency of 140.625 MHz and optical wavelengths of 690 and 830 nm.

FD-NIRS was implemented in a Multi-Distance scan (MD FD-NIRS). To achieve this, a single detector fiber bundle (-3 mm) was held at a fixed location and two co-localized source fibers (-600 μm; two wavelengths) were scanned via a linear stage (Figure 1a). In doing so, the complex reflectance (amplitude and phase) was measured at eleven distances (from 15 to 25 mm spaced by 1 mm; Figure 1b). This MD scan allowed for measurements of complex reflectance slopes (amplitude and phase) versus source-detector distance without the need for calibration (assuming unchanging coupling with the sample since the same fibers are used for each distance and the fiber/sample contact remains about the same

during the linear scan). These measurements were used to calculate the absolute optical properties of the diffuse medium as described in Section 2.3.

**Figure 1.** Frequency-domain near-infrared spectroscopy methods to achieve measurements of absolute optical properties. (**a**) Render of multi-distance scan done on diffuse optical phantoms. (**b**) Schematic of eleven different source positions realized during the multi-distance scan.

#### 2.1.2. Broadband Diffuse Reflectance Spectroscopy

Continuous-Wave broadband Diffuse Reflectance Spectroscopy (CW-bDRS) was implemented with the purpose of measuring absolute absorption spectra of optically diffusive media (Figure 2a). This was realized using a Dual-Slope (DS) optode geometry which allowed for calibration free measurements of the CW reflectance slope (versus sourcedetector distance). The DS optode configuration (Figure 2c) used was the same as described previously [26]; it contains two source positions and two detector positions. The linearly symmetric arrangement resulted in a calibrated measurement of the slope (versus source-detector distance) of CW reflectance from 25 to 35 mm.

The DS CW-bDRS system was custom built to achieve the multiplexing needed for DS measurements. This measurement requires the reflectance signal be acquired from all combinations of sources (named 1 and 2) and detectors (named A and B; Figure 2c). To do so, both the sources and detectors must be multiplexed (Figure 2b). Sources were multiplexed by using two shuttered light sources (AvaLight-HAL-S-Mini, Avantes, Louisville, CO, USA), each connected to one source position. Shutter state was controlled via a Transistor-Transistor Logic (TTL) signal from a micro-controller (Uno R3, Arduino, Ivrea, Italy) which was connected to the control computer via Universal Serial Bus (USB). Light was delivered from the sources to the probe using -2 mm fiber bundles. The light sources output an approximate black-body spectrum with a temperature of about 2650 K (peak at about 1000 nm). On the detector side, diffuse CW light was collected using a -600 μm fiber at each of the two detector positions. Multiplexing of the detector positions was done using a 1 × 2 fiber switch (LBMB, Photonwares, Woburn, MA, USA) where the common end was connected to a spectrometer (AvaSpec-HERO, Avantes) via a -600 μm fiber. The spectrometer was configured with a 500 μm slit and collected 1024 wavelengths between about 498 nm and about 1064 nm.

The bottleneck in the multiplexing sequence is the switching time of the fiber switch. Because of this, the cycle of source-detector position combinations was optimized to minimize fiber switch actuation. Naming source-detector combinations as source number (1 or 2), then detector letter (A or B; Figure 2c) as an example measurement sequence could be:

. . .1A]→[1A2A2B1B]→[1B2B2A1A]→[1A. . .

where "" is a source switch, "" is a detector switch, and the square brackets show one full DS set acquisition. The reflectance from each source-detector combination is linearly interpolated to the average absolute time for their corresponding DS set. The four reflectance measurements that contribute to a DS measurement are synchronous, thus preventing potential temporal artifacts when studying dynamic samples or biological tissue. The specific acquisition parameters (including sampling rate and wavelength range considered) are stated in Section 2.2, and the analysis of the DS data (slope of CW reflectance versus source-detector distance) to yield absolute absorption spectra is described in Section 2.3.

**Figure 2.** Broadband diffuse reflectance spectroscopy methods to achieve measurements of absolute absorption spectra. (**a**) Render of diffuse reflectance spectroscopy probe on a diffuse optical phantom. (**b**) Schematic of dual-slope diffuse reflectance spectroscopy device. Acronyms: Universal Serial Bus (USB) and Transistor-Transistor Logic (TTL). (**c**) Schematic of the source (1 and 2) and detector (A and B) positions on the dual-slope diffuse reflectance spectroscopy probe.

#### *2.2. Experiment*

The purpose of the experiment was to validate the Dual-Slope broadband Diffuse Reflectance Spectroscopy (DS CW-bDRS) instrument. This was done using a liquid optical phantom which was measured at different dyes concentrations to confirm the DS CWbDRS's ability to distinguish spectral features and measure absolute absorption. A total of 5 L of the phantom was made and placed in a cylindrical tank (Figures 1a and 2a). The base of the phantom was made of 2% (reduced fat) Milk and Water (MW; 43% milk and 57% water volume fraction) such that the scattering properties were similar to those of biological tissue in the near-infrared. The scattering was assumed to not change as dyes were added since the addition of the dyes is expected to have a negligible effect on scattering. Each time the phantom was measured, the measurement was done by both Multi-Distance Frequency-Domain Near-Infrared Spectroscopy (MD FD-NIRS) and DS CW-bDRS.

Three dyes were added in the following order: black India-Ink (II; Higgins, Leeds, MA, USA), NIR746A (N7; QCR Solutions, Palm City, FL, USA), and NIR869A (N8; QCR Solutions). II is expected to have a relatively flat or decreasing absorption spectrum with wavelength in the Near-Infrared (NIR) range [35,36]. By contrast, N7 and N8 are expected to have absorption peaks at around 746 nm and 869 nm, respectively (Figure 3). However, the actual peak wavelength may shift depending on the chemical properties of the solvent [37,38].

**Figure 3.** Expected spectral features of near-infrared dyes in water, provided by the manufacturer (QCR Solutions, Palm City, FL, USA). [37] (**a**) Normalized Absorbance (A) versus wavelength (*λ*) of NIR746A dye (N7). (**b**) Normalized A of NIR869A dye (N8).

The measurements and dye additions proceeded as follows (measurements meaning both MD FD-NIRS and DS CW-bDRS). First, the MW base phantom was mixed, then measured. Then, II was mixed in and the phantom measured again. Next, N7 was mixed and again the phantom measured. Finally, N8 was also mixed into the phantom and the measurement was done for a final time. The amount of each addition of dye was such that all absorption spectra were expected to stay within typical absorption values of biological tissue (approximately 0.005 to 0.020 mm<sup>−</sup>1) [2].

The temporal sampling parameters of the DS CW-bDRS system were set to yield a relatively fast sampling rate while still measuring a relevant wavelength range. Along these lines, DS CW-bDRS data were analyzed between 600 nm and 900 nm and the sampling rate was set to 2.5 Hz. For DS CW-bDRS 5 min of data were averaged for analysis. The MD FD-NIRS measurement took 30 s at each step for a total measurement time of 5.5 min (eleven steps; Section 2.1.1), and had a sampling rate of 2.8 Hz. These parameters were chosen such that the rate and time of DS CW-bDRS and MD FD-NIRS were approximately equal.

This phantom experiment resulted in MD FD-NIRS measurements at two wavelengths and DS CW-bDRS measurements between 600 nm and 900 nm for four phantoms (MW, MW+II, MW+II+N7, and MW+II+N7+N8). Section 2.3 describes the analysis of these data to result in absolute optical properties at two wavelengths and absolute absorption spectra.

#### *2.3. Analysis*

In Section 2.1, two techniques were described: Multi-Distance Frequency-Domain Near-Infrared Spectroscopy (MD FD-NIRS; Section 2.1.1) and Dual-Slope Continuous-Wave broadband Diffuse Reflectance Spectroscopy (DS CW-bDRS; Section 2.1.2). MD yields the spatial dependence of reflectance across multiple distances, while DS yields two symmetric spatial dependencies of reflectance. Additionally, FD-NIRS yields complex reflectance (at a few wavelengths), while CW-bDRS yields CW reflectance (at many wavelengths). However, regardless of these differences, the same basic analysis was used to find absolute optical properties (with the process repeated for each separate wavelength). In the following, we first describe the analysis of complex reflectance across multiple distances (MD FD-NIRS). Then, we explain how the method is modified to handle DS data and what additional assumptions are needed to use CW data (DS CW-bDRS).

Considering a semi-infinite medium with extrapolated-boundary conditions, the complex reflectance (*R*), defined as the net optical flux exiting the tissue per unit source power, may be written as [2,39]:

$$\widetilde{\mathcal{R}} = \frac{1}{4\pi t} \left( \overline{\mathcal{C}\_1} e^{-\overline{\mu\_{eff}}r\_1} + \overline{\mathcal{C}\_2} e^{-\overline{\mu\_{eff}}r\_2} \right), \tag{1}$$

$$
\overline{\mathcal{C}\_1} = z\_0 \left( \frac{1 + \overline{\mu\_{eff}} r\_1}{r\_1^3} \right) \qquad \overline{\mathcal{C}\_2} = -z\_0' \left( \frac{1 + \overline{\mu\_{eff}} r\_2}{r\_2^3} \right) , \tag{2}
$$

$$r\_1 = \sqrt{\rho^2 + z\_0^2} \qquad r\_2 = \sqrt{\rho^2 + z\_0^{'}^2} \tag{3}$$

where *z*<sup>0</sup> is the depth of the real isotropic source (*z*<sup>0</sup> = 1/*μ <sup>s</sup>*), *z* <sup>0</sup> is the height of the imaginary isotropic source (*z* <sup>0</sup> = −*z*<sup>0</sup> + 2*zb*), *ρ* is the source-detector distance on the surface of the semi-infinite medium, and *<sup>μ</sup>eff* is the complex effective attenuation coefficient. The extrapolated boundary height (*zb* = 2*A*/(3(*μ <sup>s</sup>* + *μa*))) and the complex effective attenuation coefficient (*μeff* <sup>=</sup> 3(*μ <sup>s</sup>* + *μa*)(*μ<sup>a</sup>* − *ωnini*/*c*) ) are both expressed in terms of the absorption coefficient (*μa*), reduced scattering coefficient (*μ <sup>s</sup>*). The remaining parameters are the angular modulation frequency (*ω*), the reflection parameter (*A*), the index of refraction inside the medium (*nin*; assumed to be 1.4, typical for biological tissue [2]), and the speed of light in vacuum (*c*). Additionally, *ω* = 2*π fmod*, where *fmod* is the modulation frequency. Finally, *A* is a function of *nin* and the index of refraction outside the medium (*nout*; assumed to be 1.0) [39]. To yield a linear relationship, Equation (1) can be rewritten as:

$$\ln\left[\frac{4\pi\bar{\mathcal{R}}}{\overline{\mathcal{C}\_1} + \overline{\mathcal{C}\_2}e^{\overline{\mu\_{eff}}(r\_1 - r\_2)}}\right] = -\overline{\mu\_{eff}}r\_{1\prime} \tag{4}$$

which shows a linear dependence of a modified complex reflectance (*y* <sup>=</sup> *ln*[4*πR*/ (*C* <sup>1</sup> + *C* 2*e <sup>μ</sup>eff*(*r*1−*r*2) )]) on *<sup>r</sup>*<sup>1</sup> (i.e., *<sup>y</sup>* <sup>=</sup> <sup>−</sup>*μeffr*1).

The linear complex slope (−*μeff*) in Equation (4) is used to calculate the desired optical properties (*μ<sup>a</sup>* and *μ <sup>s</sup>*). However, *<sup>y</sup>* is dependent on *<sup>μ</sup><sup>a</sup>* and *<sup>μ</sup> <sup>s</sup>*, and *r*<sup>1</sup> is dependent on *μ s*. Thus, an iterative method is used to recover the *<sup>μ</sup>eff* :

$$\widetilde{y}\left[\mu\_{a,k'}\mu\_{s,k}^{'}\right] = -\overline{\mu\_{eff,k+1}}r\_1\left[\mu\_{s,k}^{'}\right] + \widetilde{y}\_{k+1'} \tag{5}$$

where *k* is the iteration number, and *b* is the complex intercept (dependent on instrumental factors not considered in Equation (4)). For the iterative method, an initial guess of *μ<sup>a</sup>* and *μ <sup>s</sup>* is needed for *k* = 0. This guess was determined based on the linear slopes of complex reflectance amplitude and phase as previously described [19,21]. At each iteration, the current values of *<sup>y</sup>* and *<sup>r</sup>*<sup>1</sup> were used to fit a complex linear slope (−*μeff*):

$$
\begin{bmatrix} -\widehat{\mu\_{\ell f}} \\ \widetilde{b} \end{bmatrix}\_{k+1} = \begin{bmatrix} r\_{1,1} & 1 \\ \vdots & \vdots \\ r\_{1,n} & 1 \end{bmatrix}\_{k}^{+} \begin{bmatrix} \widetilde{y}\_{1} \\ \vdots \\ \widetilde{y}\_{n} \end{bmatrix}\_{k} \tag{6}
$$

where the + superscript is the Moore-Penrose inverse, *b* is the complex intercept (ignored), and there are *n* measurements (distances). In the case of MD FD-NIRS, *n* = 11 since there were eleven measurements of complex reflectance at eleven different source-detector distances. To convert *<sup>μ</sup>eff* to *<sup>μ</sup><sup>a</sup>* and *<sup>μ</sup> <sup>s</sup>* the following expressions were used:

$$
\mu\_t' = \frac{-2c\Re\left[\overline{\mu\_{eff}}\right] \Im\left[\overline{\mu\_{eff}}\right]}{3n\omega},
\tag{7}
$$

$$\mu\_a = \frac{\Re\left[\overline{\mu\_{eff}}\right]^2 - \Im\left[\overline{\mu\_{eff}}\right]^2}{\Im\mu\_t^{'}},\tag{8}$$

where *μ <sup>t</sup>* is the total reduced attenuation coefficient, and *μ <sup>s</sup>* = *μ <sup>t</sup>* − *μa*. At each iteration, the new fitted <sup>−</sup>*μeff* yielded the *<sup>μ</sup><sup>a</sup>* and *<sup>μ</sup> <sup>s</sup>* used to calculate *<sup>y</sup>* and *<sup>r</sup>*<sup>1</sup> for the next iter-

ation (Equation (5)). The procedure terminated when the condition <sup>|</sup>*μeff <sup>k</sup>* <sup>−</sup> *<sup>μ</sup>eff <sup>k</sup>*−<sup>1</sup> | < 10−<sup>4</sup> L/mm was met.

The above expressions show how MD FD-NIRS data were used to find absolute *μ<sup>a</sup>* and *μ <sup>s</sup>* at two wavelengths. To extend this to DS CW-bDRS, first, Equation (6) was modified to use DS data. DS yields two symmetric measurements of *R* versus *ρ* (where *R* is now real not complex since CW data is used). Thus, instead of Equation (6), the following was used for DS CW-bDRS:

$$
\mu\_{eff\_{k+1}} = \frac{-1}{2} \left( \frac{y\_{2,1} - y\_{1,1}}{r\_{1,2,1} - r\_{1,1,1}} + \frac{y\_{2,2} - y\_{1,2}}{r\_{1,2,2} - r\_{1,1,2}} \right)\_k \tag{9}
$$

where the first subscript of *y* and second subscript of *r*<sup>1</sup> corresponds to source-detector distance, and the second and third subscript, respectively, corresponds to the symmetric slope in the DS set. This expression is simply the average of the slopes of *y* versus *r*<sup>1</sup> between the two DS symmetric sets. In addition, note that the tildes were removed since *y* is now real (from real *R*). This results in a real effective attenuation coefficient (*μeff* = 3*μa*(*μ <sup>s</sup>* + *μa*) ) instead of a complex one. Equations (1)–(5) may all be used as is, with removing the tilde hats to now represent real data.

Since the DS CW-bDRS method yields a real *μeff* , *μ <sup>s</sup>* must be known to find *μ<sup>a</sup>* using the following expression instead of Equation (8):

$$
\mu\_a = \sqrt{\frac{\mu\_s'^2}{4} + \frac{\mu\_{eff}^2}{3}} - \frac{\mu\_s'}{2}. \tag{10}
$$

To address the need to know *μ <sup>s</sup>*, the data from the MD FD-NIRS measurements were used. FD-NIRS yielded a measurement of *μ <sup>s</sup>* at two wavelengths (*μ <sup>s</sup>*,*FD*(*λ*1) and *μ <sup>s</sup>*,*FD*(*λ*2); where *λ* is wavelength). These *μ <sup>s</sup>* measurements are then extrapolated using the following expression to find *μ <sup>s</sup>* at all the wavelengths used by bDRS:

$$
\mu\_s'(\lambda) = \mu\_{s,FD}'(\lambda\_1) \left(\frac{\lambda}{\lambda\_1}\right)^{-\ln\left[\frac{\mu\_{s,FD}'(\lambda\_2) + \mu\_{s,FD}'(\lambda\_1)}{\ln\left[\lambda\_1 + \lambda\_2\right]}\right]}\tag{11}
$$

Therefore, when analyzing DS CW-bDRS data, *μ<sup>a</sup>* is found for each wavelength using the assumed *μ <sup>s</sup>* expressed above for each wavelength.

In summary, the same basic method, based on iteratively fitting to Equation (4), is used for MD FD-NIRS and DS CW-bDRS to find absolute optical properties. DS differs from MD in the calculation of slope, and the *μ <sup>s</sup>* from MD FD-NIRS is extrapolated to find an assumed *μ <sup>s</sup>* for each DS CW-bDRS wavelength.

#### **3. Results**

The phantom experiment described in Section 2.2 consisted of measuring four phantoms, each with both Multi-Distance Frequency-Domain Near-Infrared Spectroscopy (MD FD-NIRS; Figure 1a) and Dual-Slope Continuous Wave broadband Diffuse Reflectance Spectroscopy (DS CW-bDRS; Figure 2a). The phantoms were Milk and Water (MW), that with India Ink added (MW+II), then that again with NIR746A added (MW+II+N7), and then that yet again with NIR869A added (MW+II+N7+N8).

Figure 4 shows the results of this experiment. The absolute absorption (*μa*) spectrum is shown across wavelengths (*λ*) for all four phantoms and two measurement methods (Figure 4a). Additionally, Figure 4a shows a dashed line representing the modeled (as 99.1% water [40,41] and 0.9% lipid [42] volume fractions) absorption of the MW phantom. The coefficient of determination (*R*2) was calculated for the MW data and model, yielding a value of 0.98, implying that 98% of the variance in the MW data is explained by the modeled MW absorption. In Figure 4b–d, the change in absorption (Δ*μa*) as a result of

adding the three different dyes (II, N7, and N8) is shown. The error regions in all plots are dominated by the systematic uncertainty in the distances on the DS-CW-bDRS probe (estimated at ±0.1 mm for chained dimensions) and within the MD-FD-NIRS measurement (estimated at ±0.5 mm for initial position and ±0.01 mm for scan pitch). Note the flat absorption contribution from II and the spectral features from N7 and N8. Particularly, a peak between 700 nm and 800 nm for N7 and a peak between 800 nm and 900 nm for N8. How these results serve to validate the instrument will be discussed in Section 4. Briefly, notice the agreement between the MD FD-NIRS and DS CW-bDRS measurements and the agreement between the expected MW spectrum and the DS CW-bDRS measured MW spectrum.

**Figure 4.** Results from phantom experiment. (**a**) Absolute absorption (*μa*) spectra as a function of wavelength (*λ*). Showing results from Dual-Slope Continuous-Wave broadband Diffuse Reflectance Spectroscopy (DS CW-bDRS) and Multi-Distance Frequency-Domain Near-Infrared Spectroscopy (MD FD-NIRS) measurements. Spectra shown for the following phantoms: Milk and Water (MW), MW plus India Ink (II), MW plus II plus NIR746A (N7), and, finally, MW plus II plus N7 plus NIR869A (N8). DS CW-bDRS points show individual wavelength measurements and lines show smoothed (moving average) spectra for visualization. Dashed line shows the expected spectrum for MW modeled as water and lipid. (**b**) Change in absorption (Δ*μa*) from adding N8 (i.e., *<sup>μ</sup>MW*+*I I*+*N*7+*N*<sup>8</sup> *<sup>a</sup>* <sup>−</sup> *<sup>μ</sup>MW*+*I I*+*N*<sup>7</sup> *<sup>a</sup>* ), (**c**) <sup>Δ</sup>*μ<sup>a</sup>* from adding N7 (i.e., *<sup>μ</sup>MW*+*N*<sup>7</sup> *<sup>a</sup>* <sup>−</sup> *<sup>μ</sup>MW*+*I I <sup>a</sup>* ), and (**d**) <sup>Δ</sup>*μ<sup>a</sup>* from adding II (i.e., *<sup>μ</sup>MW*+*I I <sup>a</sup>* <sup>−</sup> *<sup>μ</sup>MW <sup>a</sup>* ).

#### **4. Discussion**

The experiment involved a liquid phantom and sought to validate the Dual-Slope Continuous Wave broadband Diffuse Reflectance Spectroscopy (DS CW-bDRS) instrument (Sections 2.2 and 3). This was done by measuring a highly scattering phantom as dyes were progressively added (Figure 4). The expectation was that the DS CW-bDRS (in combination with Multi-Distance Frequency-Domain Near-Infrared Spectroscopy; MD FD-NIRS) instrument would be able to measure the known spectral features of the dyes and the absolute absorption of the phantom.

First, before any dyes were added (just Milk and Water; MW), the phantom was measured to yield an absorption spectrum dominated by water (Figure 4a). This MW spectra had the expected spectral features of water (hump at about 750 nm and sharp increase starting at about 800 nm) [41]. Further, the expected MW spectrum agrees (within error) with the experimental data for the spectral range between 690 nm and 830 nm, where the reduced scattering coefficient is interpolated (not extrapolated) from the MD FD-NIRS measurement. Below 690 nm, the agreement is lost likely due to the very low water absorption, such that, even low absorption, contributions from fat, proteins, or other milk constituents in the MW medium may become detectable. Above 830 nm, the agreement also degrades possibly because of incorrect reduced scattering values (from extrapolation) or contributions from other absorbers in the milk other than water. Overall, the quantitative agreement between the expected and measured MW spectra are quite good with a coefficient of determination (*R*2) of 0.98, indicating an accurate measurement of absolute absorption by DS CW-bDRS.

Further, the absorption measured by MD FD-NIRS agrees within error with DS CWbDRS at all points of comparison except for NIR869A (N8). For N8, the difference between the measurements is about 10%, and the combined errors for the two measurements amount to about 8% (dominated by uncertainty in the distances; Section 3). Admittedly, the reduced scattering found by MD FD-NIRS is used to find the absorption for DS CWbDRS; thus, the measurements are related (Section 2.3). However, the observed agreement demonstrates a reliable absolute measurement of the slope of diffuse reflectance with the DS CW-bDRS instrument, given the accurate measurements with the scanning MD method used by MD FD-NIRS in controlled laboratory conditions. In the case of N8, where the two instruments do not agree, the most likely explanation is instrumental error (distance uncertainty and boundary conditions). The distance uncertainty accounts for most of the difference between the measurements, and the remaining unaccounted disagreement may be explained by uncertainties in boundary conditions (true depth of fibers in phantom, existence of meniscus, etc.). These distance uncertainties and boundary effects may have changed for each measurement, since after the addition of each dye the instruments needed to be re-setup and placed on the phantom. This may explain why instrumental errors impacted differently the measurements of different dyes. Additionally, there is evidence that the NIR dyes are not stable over time (discussed below), and the two measurements are not simultaneous (5–10 min between measurements). Any of these considerations may have led to the lack of agreement for N8. However, given the existence of these considerations and the close agreement in six out of eight cases, we find that the results serves to validate the absolute absorbance measured by DS CW-bDRS.

Focusing on the individual dye additions allows for validation of the measurement of spectral features. The first dye added was India Ink (II), which is expected to have a flat or decreasing spectral dependence with wavelength [35]. The change in the absorption spectrum observed confirmed this wavelength dependence and, again, agreed with the measurements at two wavelengths with MD FD-NIRS (Figure 4d). Further evaluation was carried out to estimate the recovered concentration of II given the measured change in absorption. The true concentration (in volume fraction) of II in the phantom was 5.7 × <sup>10</sup><sup>−</sup>6, given the phantom volume and the amount of ink added. The same II was measured in transmission in non-diffuse solutions of known concentrations to yield spectra of the total attenuation coefficient (*μt*; assuming only unscattered light is detected). To yield the

spectrum of absorption (*μa*; needed to recover concentration in the diffuse experiment) for II, the single scattering albedo (*a*) must be assumed (*μ<sup>a</sup>* = *μt*(1 − *a*)). This must be considered for II since it is not in solution in water but instead a suspension of carbon particles. Wavelength independent values of *a* = 0 to *a* = 0.15 were assumed which yielded recovered concentrations of 5.2 × <sup>10</sup>−<sup>6</sup> to 6.1 × <sup>10</sup>−6, respectively. This range of low albedo values of pure II indicates that its scattering coefficient is much smaller than the absorption coefficient, which is consistent with the literature [35,36]. The range resulted in measurement errors of −8.8% to 7.0% (given true value of 5.7 × <sup>10</sup>−6). Therefore, it is expected that DS CW-bDRS is capable of recovering accurate chromophore concentrations given that the extinction spectra is known (not the case here since the albedo for II was not measured). For NIR746A (N7) and N8, the concentration recovery was not carried out due to temporal and spectral instability of the dyes (discussed below). Future experiments will be undertaken to validate the accurate chromophore concentration recovery of DS CW-bDRS using soluble and stable dyes.

Moving on from II, the next two dyes were expected to feature a more interesting spectrum (Figure 3). NIR746A (N7) was added after II and the change of absorption presented in Figure 4c. The expected peak between 700 nm and 800 nm was present. However, upon closer examination, the exact peak location is shifted about 12 nm higher for the DS CW-bDRS measurement compared to the expected spectrum (Figure 5a). When the next and final dye, NIR869A (N8), was added, yet again there was a peak present in approximately the expected location (Figure 4b). But, as with N7, there was an approximate 12-nm shift to longer wavelengths in the DS CW-bDRS measurement. This is consistent with a bathochromic shift, which is possible for these dyes given the manufacturer information and previous studies [37,38]. The hypothesis is that, when in milk, the dyes exhibit a bathochromic shift due to the different bulk polar nature of milk versus water. To test this possibility, two further measurements were done on N8. N8 was measured in transmission (in a semi-micro cuvette), both in water and in the MW mixture (Figure 5b). The results show that the transmission spectral absorbance for N8 in milk matches the DS CW-bDRS peak location, whereas the transmission spectral absorbance of N8 in water matches the manufacturer provided spectra. Thus, this transmission experiment supports the hypothesis that these dyes exhibit a bathochromic shift of about 12 nm when in the MW mixture. However, it is unknown if the amplitude of the absorption peak is also effected (hyperchromic or hypochromic shift) for these dyes. The dyes were also found to be temporally unstable, as extended time in solution caused N7 to lose its near-infrared absorption peak, and N8's peak shifted roughly 200 nm to the blue. These spectral and temporal instabilities stopped the analysis of the concentrations of these dyes in the diffuse phantom. But, despite this, DS CW-bDRS was still capable of distinguishing the dye's spectral features.

Summarizing the discussion above, these results on the liquid phantom serve to validate the DS CW-bDRS instrument's ability to measure spectral features and absolute values of absorption in a diffuse medium. In six of eight cases (all except 690 nm and 830 nm for N8), the absolute absorption measured with DS CW-bDRS agreed within error with MD FD-NIRS. DS CW-bDRS also accurately measured the expected absorption spectrum of MW (*R*<sup>2</sup> = 0.98). Additionally, the DS CW-bDRS instrument correctly measured the flat spectra of II and the peaks of N7 and N8. The concentration of II was estimated, suggesting the ability to recover chromophore concentrations given the extinction. Finally, it was confirmed that the recovered peak locations of the N7 and N8 dyes are what is expected for these dyes in milk. Future work will be done to validate the overall methods ability to accurately recover chromophore concentrations.

**Figure 5.** Comparison of wavelength (*λ*) normalized Absorbance (A) peak locations for two dyes. Four types of spectra shown: Diffuse Reflectance (Diff Refl, i.e., in milk), provided by the Manufacturer (Manu; (QCR Solutions, Palm City FL, USA) [37], i.e., in water), Diffuse Transmission (Diff Tran, i.e., in milk), and Transmission (Tran, i.e., in water). (**a**) NIR746A (N7) dye. (**b**) NIR869A (N8) dye.

#### **5. Conclusions**

In this article, we have presented a new Dual-Slope Continuous Wave broadband Diffuse Reflectance Spectroscopy (DS CW-bDRS) instrument. Experiments on highly scattering liquid phantoms demonstrated the instrument's capability of measuring absolute absorbance spectra without any need for instrumental calibration. By combining DS CW-bDRS and Frequency-Domain Near-Infrared Spectroscopy (FD-NIRS; to account for scattering contributions to the absorbance spectrum), we were able to measure absolute absorption spectra of the liquid phantoms that contained a combination of three dyes. These experiments demonstrated the technique's ability to perform absolute spectral absorption measurements and retrieve the correct spectral features of various dyes. Future work will focus on further validation of chromophore concentration measurements. The importance of this work lies in the development of the DS CW-bDRS instrument. This DS CW-bDRS instrument is novel in that it combines DS and bDRS to achieve calibration-free measurements and provides a valuable tool for absolute spectral measurements of highly scattering media, including biological tissues.

**Author Contributions:** Conceptualization, G.B. and S.F.; methodology, G.B., R.D., S.M., H.N., L.V., C.F., A.S., T.P., and S.F.; software, G.B., R.D., S.M., H.N., and L.V.; validation, G.B., R.D., S.M., H.N., L.V., and C.F.; formal analysis, G.B.; investigation, G.B., R.D., S.M., H.N., L.V., and C.F.; resources, A.S. and S.F.; data curation, G.B., A.S., and S.F.; writing—original draft preparation, G.B., A.S., and S.F.; writing—review and editing, G.B., C.F., A.S., T.P., and S.F.; visualization, G.B.; supervision, G.B., C.F., A.S., and S.F.; project administration, G.B. and S.F.; funding acquisition, S.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was funded by the National Institutes of Health (NIH) grant R01 NS095334.

**Data Availability Statement:** Data and supporting codes are available upon request.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


## **Comparison of Optical Imaging Techniques to Quantitatively Assess the Perfusion of the Gastric Conduit during Oesophagectomy**

**Maxime D. Slooter 1, Sanne M. A. Jansen 2, Paul R. Bloemen 2, Richard M. van den Elzen 2, Leah S. Wilk 2, Ton G. van Leeuwen 2, Mark I. van Berge Henegouwen 1,\*, Daniel M. de Bruin <sup>2</sup> and Suzanne S. Gisbertz <sup>1</sup>**


Received: 5 July 2020; Accepted: 5 August 2020; Published: 10 August 2020

**Abstract:** In this study, four optical techniques—Optical Coherence Tomography, Sidestream Darkfield Microscopy, Laser Speckle Contrast Imaging, and Fluorescence Angiography (FA)—were compared on performing an intraoperative quantitative perfusion assessment of the gastric conduit during oesophagectomy. We hypothesised that the quantitative parameters show decreased perfusion towards the fundus in the gastric conduit and in patients with anastomotic leakage. In a prospective study in patients undergoing oesophagectomy with gastric conduit reconstruction, measurements were taken with all four optical techniques at four locations from the base towards the fundus in the gastric conduit (Loc1, Loc2, Loc3, Loc4). The primary outcome included 14 quantitative parameters and the anastomotic leakage rate. Imaging was performed in 22 patients during oesophagectomy. Ten out of 14 quantitative parameters significantly indicated a reduced perfusion towards the fundus of the gastric conduit. Anastomotic leakage occurred in 4/22 patients (18.4%). At Loc4, the FA quantitative values for "T1/2" and "mean slope" differed between patients with and without anastomotic leakage (*p* = 0.025 and *p* = 0.041, respectively). A quantitative perfusion assessment during oesophagectomy is feasible using optical imaging techniques, of which FA is the most promising for future research.

**Keywords:** optical imaging; fluorescence imaging; fluorescence angiography; indocyanine green (ICG); optical coherence tomography (OCT); laser speckle contrast imaging (LSCI); esophagectomy; gastric conduit; Sidestream Darkfield Microscopy (SDF); multispectral

#### **1. Introduction**

Surgical morbidity related to inadequate perfusion is a major challenge after oesophagectomy with gastric conduit reconstruction for oesophageal cancer. This type of surgical morbidity includes anastomotic leakage and graft necrosis, and remains high even after recent developments in minimally invasive surgery, enhanced recovery programs, and the centralisation of cancer surgery, occurring in up to 20% of the patients [1–4].

To achieve digestive continuity after oesophagectomy, the replacement of the native oesophagus by a constructed gastric conduit is a potential solution. The perfusion of the gastric conduit is mainly based on the right gastroepiploic artery, which usually terminates before the future anastomotic site, causing the future anastomotic site to be at risk for inadequate perfusion. Surgeons assess the perfusion of the gastric conduit by observing the colour of the serosal surface, the bleeding of resection margins, or arterial palpation. To aid surgeons' decision-making, different optical imaging techniques have been introduced to assess perfusion intraoperatively [5,6].

Optical imaging, which is established by the interaction between light and tissue, is potentially able to image changes in perfusion real-time and at high resolution. These characteristics of optical imaging have the potential to analyse a perfusion intraoperatively and in a quantitative manner [7,8]. Some optical imaging techniques are already Food and Drug Administration (FDA)-approved and Conformité Européenne (CE)-marked systems, but others lack clinical evaluation [5,9]. Not all techniques have an adequate specificity for patient outcomes and are so far unable to objectify perfusion in quantitative values for clinical usage [5]. The comparison of techniques in the literature is hampered by the large range of available imaging systems and the heterogeneity of quantitative parameters and perfusion endpoints [5]. Therefore, there is scant evidence on which optical imaging technique and related quantitative parameter is the most predictive for inadequate perfusion and patient outcomes.

In this explorative study (IDEAL phase 2a), we aim to compare four optical techniques using different wavelengths of light on performing a quantitative perfusion assessment of the gastric conduit in patients undergoing oesophagectomy with gastric conduit reconstruction for oesophageal cancer. The four optical techniques are Optical Coherence Tomography (OCT), Sidestream Darkfield Microscopy (SDF), Laser Speckle Contrast Imaging (LSCI), and Fluorescence Angiography (FA). We hypothesised that the quantitative parameters of the optical techniques show decreased perfusion (A) towards the fundus in the gastric conduit and (B) in patients with anastomotic leakage. By comparing these four techniques, we aimed for a recommendation of an optical imaging technique for a quantitative perfusion assessment of the gastric conduit for future studies and clinical intraoperative use.

#### **2. Materials and Methods**

#### *2.1. Study*

This is a single-centre, prospective, observational study in patients undergoing elective oesophagectomy with gastric conduit reconstruction for oesophageal cancer. This study adhered to the STARD (Standards for Reporting of Diagnostic Accuracy Studies) [10] and STROBE (STrengthening the Reporting of OBservational studies in Epidemiology) [11] guidelines.

This study was approved by the Institutional Review Board of the Amsterdam UMC, location AMC, (NL52377.018.15) and submitted to the clinicaltrials.gov database (NCT02902549) [12]. Written informed consent was obtained from all patients.

#### *2.2. Procedure*

All the procedures were performed totally minimally invasively with laparoscopy and thoracoscopy. A McKeown (cervical anastomosis) or Ivor Lewis (intrathoracic anastomosis) procedure was performed by two specialised upper-gastrointestinal surgeons (MIvBH, SSG), as previously described [13]. The stomach was reconstructed into a 3–4 cm wide gastric conduit using a powered linear stapler. Directly after gastric conduit reconstruction, perfusion measurements were performed using four optical techniques in the following order: OCT, SDF, LSCI, and FA. Perfusion measurements were performed at four locations (Loc1–4): 3 cm proximal of the level of (Loc1), at the level of (Loc2), 3 cm distal to (Loc3) at the level of the watershed between the right and left gastroepiploic artery, and at the level of the gastric fundus (Loc4) (Figure 1). A sterile gauze was placed to point out the watershed area (the end of the right gastroepiploic artery). Subsequently, a cervical or intrathoracic anastomosis to the proximal oesophagus was constructed. The anastomotic site of the gastric conduit was determined by a conventional white light assessment, independent of imaging outcomes, and was located between Loc3 and Loc4.

**Figure 1.** Qualitative images by all four imaging techniques: upper left Optical Coherence Tomography (OCT), lower left Sidestream Darkfield Imaging (SDF), upper right Laser Speckle Contrast Imaging (LSCI), and lower right Fluorescence Angiography (FA). Perfusion measurements were performed at four locations (Loc1–4): 3 cm proximal of the level of (Loc1), at the level of (Loc2), 3 cm distal to (Loc3) the level of the watershed between the right and left gastroepiploic artery, and at the level of the gastric fundus (Loc4).

#### *2.3. Imaging Modalities and Software*

OCT was performed by a commercial 50 kHz IVS 2000 swept source system (Santec, 5823 Ohkusa-Nenjozaka, Komaki-City, Aichi, Japan), an FDA-approved and CE-marked system, using a centre wavelength of ∼1300 nm. Measurements were performed as previously described [14]. Data analysis was performed using custom-made scripts written in MATLAB (Mathworks, Natick, MA, USA) [14].

SDF is a hand-held, in-contact, green light (530 nm) imaging system that visualises microcirculation. Image contrast is based on the absorption of green light by oxyhaemoglobin in red blood cells (RBCs) [15]. Therefore, RBCs will appear as dark spots in the image. SDF measurements were performed by the Microscan (Microvision Medical, Amsterdam, The Netherlands), an FDA-approved system to monitor the microcirculation of critically ill patients. To stabilise the image, a custom-made stabiliser based on the design by Balestra et al. [16] was used to adhere the tissue to the imaging system. The SDF obtained images were analysed using the AVA3.2 software (Microvision Medical BV, Amsterdam, The Netherlands).

LSCI is a near infrared imaging system with a laser operating at 785 nm and measures speckle fluctuations induced by the motion and change in number of RBCs. The LSCI was used as previously described [17]. Imaging was performed by the MoorFLPI-1™ (Moor instruments, Devon, UK), an FDA-approved and CE-marked system, positioned at 40 cm above the gastric conduit, as measured by a laser distance meter (Leica Geosystems D110, Leica Camera, Wetzlar, Germany).

FA is a near infrared imaging modality that images the spatial distribution of the contrast agent indocyanine green (ICG), detecting at an emission wavelength of approximately 800 nm. FA was performed by the Artemis system® (Quest Medical Imaging, Middenmeer, The Netherlands), an FDA-approved and CE-marked system, after the intravenous administration of 2.5 mg of ICG. The Artemis system was placed under a 90-degree angle and also at a distance of 40 cm above the gastric conduit, as measured with the laser distance meter. A metric ruler was placed in the field of view for image calibration. The anaesthesiologist administered 2.5 mg of ICG intravenously. The perfusion was imaged in continuous video recordings of 2–3 min. In-house software was developed to plot the ICG in- and outflow in a fluorescence–time curve. A ruler in the field of view was used for calibration, and for every location a circular region of interest (ROI) with a diameter of 1 cm was used to measure the ICG signal.

#### *2.4. Outcomes*

The primary objective of this study was the comparison of the ability of the four optical imaging systems—OCT, SDF, LSCI, and FA—to obtain the quantitative perfusion parameters during surgery, and to explore the relation between the four locations (Loc1–4) and the occurrence of anastomotic leakage. The primary outcome included 14 quantitative parameters outlined below in "quantitative parameters" at four locations and anastomotic leakage. Quantitative assessment was deemed successful per technique when one or more measurements could be obtained per patient. Anastomotic leakage was defined by the Esophageal Complications Consensus Group (ECCG) classification and diagnosed by a CT-scan with contrast, upper gastrointestinal endoscopy, drainage of saliva, gastric content via wound/drains, and/or by reoperation [18].

The secondary outcomes were imaging characteristics and hemodynamic parameters. The imaging characteristics included the field of view, imaging depth, ease of use, practical limitations (overall and for quantitative analysis), invasiveness, time consumption, and quality. The quality was deemed high when no technical issues occurred to obtain qualitative images. The relations between the quantitative and hemodynamic parameters were explored at Loc1. Hemodynamic parameters were recorded at the start of imaging and included the mean arterial pressure (MAP; mmHg), heart rate (beats/min), systolic and diastolic blood pressure (mmHg), temperature nasopharynx (◦C), cardiac output (L/min), stroke volume (mL), stroke volume variation (%), and cardiac index (L/min/m2).

#### *2.5. Quantitative Parameters*

Fourteen quantitative parameters were measured (OCT *n* = 1, SDF *n* = 6, LSCI *n* = 2, and FA *n* = 5). The quantitative parameter for OCT was [14]:

1. The percentage of speckle contrast pixels (speckle %) indicative of the flow obtained in the M-mode of the scans;

Quantitative parameters derived for SDF included [15]:


The quantitative parameters calculated for LSCI were [17]:


FA was quantified by the assessment of the time-dependent change in the fluorescence signal at the ROIs in the form of the software-derived fluorescence–time curves [8,19]. The quantitative parameters for FA included:


#### *2.6. Statistics*

The quantitative perfusion parameters were shown as the median and interquartile (IQR) or total range. All the available data were used for the analysis. The comparison of quantitative perfusion parameters between Loc1 and Loc4 were analysed using the Wilcoxon Signed Ranks Test. The results at Loc3 and Loc4 were compared between patients with or without anastomotic leakage using the Mann–Whitney U test. A *p*-value of below 0.05 was considered significant. Correlations between the quantitative values and hemodynamic parameters were assessed by calculating the Spearman's rank correlation coefficient ρ. A *p*-value ≤ 0.05 was considered statistically significant. The imaging characteristics were summarised using simple descriptive statistics or in quantitative terms. The added surgical time was shown in the median and IQR. All the analyses were performed using IBM® SPSS® Statistics (version 26.0.0, IBM Corporation, NY, USA).

#### **3. Results**

Between October 2015 and June 2016, in total 26 patients that underwent minimally invasive oesophagectomy with gastric conduit reconstruction were included after written informed consent was obtained. Four patients did not receive imaging due to extended operation time and were not taken into the analysis. The baseline characteristics of 22 patients are shown in Table 1 [17]. The gastric conduit was connected to the proximal oesophagus by a cervical (*n* = 2) or intrathoracic (*n* = 20) anastomosis.


**Table 1.** Baseline characteristics.

\* Only diabetes mellitus type 2. COPD = chronic obstructive pulmonary disease.

#### *3.1. Imaging Characteristics*

Qualitative imaging was successful in all the patients (*n* = 22) with all four techniques. In Figure 1, qualitative images by all four imaging techniques are shown for one patient. The imaging characteristics are shown in Table 2. OCT and SDF had a narrow field of view, of which SDF assessed the microvasculature for a region of 1 mm × 1 mm. LSCI and FA were wide field techniques and showed similarity in the qualitative images (Figure 1). All the practical limitations are summarised in Table 2, of which artefacts, shadowing, and focus were limitations for quantitative analysis. The motion artefacts for OCT and SDF were secondary to the respiratory rate, heart rate, peristalsis, or probe handling, and other artefacts were due to air bubbles in the surgical drape. The use of OCT, SDF, LSCI, and FA added 6 (5–7), 5 (4–7), 3 (2–4), and 4 (2–4) minutes to the surgical time, respectively.



FA Fluorescence Angiography, ICG Indocyanine Green, LSCI Laser Speckle Contrast Imaging, OCT Optical Coherence Tomography, OR Operating Room, ROI Region-of-interest, SDF Sidestream Darkfield Microscopy. \* Not when used laparoscopically

#### *3.2. Quantitative Parameters*

Of 22 patients, quantitative assessment was successful in 15, 22, 20, and 20 patients for OCT, SDF, LSCI, and FA, respectively. At the four locations (Loc1/Loc2/Loc3/Loc4), the number of obtained measurements for OCT were 12/10/8/10, for SDF 21/21/21/21, for LSCI 20/20/20/20, and for FA 17/20/20/18, respectively. The measurements for the SDF parameters were available in 21 patients at Loc4, except for MFI (*n* = 14).

Comparing Loc1 with Loc4, a significant difference was noticed for 10 out of 14 parameters (Figure 2): the SDF parameters PVD, MFI, PPV, and velocity; the LSCI parameters Flux and Δ%Flux; and the FA parameters influx time point, Fmax, T1/2, and mean slope. Of those, PVD, MFI, PPV, velocity, Flux, Fmax, and mean slope decreased towards the fundus, while the Δ%Flux, influx time point, and T1/<sup>2</sup> increased. No significant differences were noticed between the locations for the OCT parameter speckle %, the SDF parameters TVD and DBS, and the FA parameter Tmax.

Anastomotic leakage occurred in 4/22 patients (18.2%). Of the measurements available at Loc3, anastomotic leakage occurred in 1/8 (12.5%), 3/21 (14.3%), 4/20 (20.0%), and 2/20 (10.0%) for OCT, SDF, LSCI, and FA, respectively. No significant differences were observed for the quantitative parameters between the patients with or without anastomotic leakage at Loc3 (data not shown). Of the measurements available at Loc4, anastomotic leakage occurred in 1/10 (10.0%), 3/21 (14.3%), 4/20 (20.0%), and 2/18 (11.1%) for OCT, SDF, LSCI, and FA, respectively. A comparison of quantitative parameters at Loc4 between patients with and without anastomotic leakage is shown in Table 3. Of all the quantitative parameters, the FA values for the T1/<sup>2</sup> and mean slope significantly differed between patients with and without anastomotic leakage. The T1/<sup>2</sup> was 13.0 (total range 6.3–21.3) seconds in patients without anastomotic leakage, and increased to 29.3 (total range 27.8–30.8) seconds in the case of anastomotic leakage (*p* = 0.025). For T1/2, the total range of values showed no overlap between patients with or without anastomotic leakage. The mean slope decreased from 1.3 (total range 0.1–11.0) to 0.2 (total range 0.1–0.2) in patients without and with anastomotic leakage (*p* = 0.041), respectively. All the values for the SDF parameters and for Flux were non-significantly higher (Table 3).

**Figure 2.** Fourteen quantitative parameters at four locations of the gastric conduit. Data are shown as medians and interquartile ranges. *p*-value is shown for the difference between Loc1 and Loc4. OCT Optical Coherence Tomography: Speckle in %. SDF Sidestream Darkfield Imaging: TVD total vessel density in mm/mm2, PVD perfused vessel density in mm/mm2, MFI microvascular flow index, DBS de Backer score, PPV proportion (%) of perfused vessels. LSCI Laser Speckle Contrast Imaging: Flux in LSPU (laser speckle perfusion units), Δ%Flux percentage difference in Flux relative to Location 1 (right vertical axis). \* Not calculated, *p*-value calculated on absolute data (see Flux). FA Fluorescence Angiography: Influx time point in sec, Fmax in arbitrary units, Tmax in sec, T1/<sup>2</sup> in sec, mean slope in arbitrary units per sec.


**Table 3.** Comparison of the perfusion-related parameters at location 4 between patients with and without anastomotic leakage.

Data are shown as median and total range. OCT Optical Coherence Tomography, SDF Sidestream Darkfield Imaging, TVD total vessel density in mm/mm2, PVD perfused vessel density in mm/mm2, MFI microvascular flow index, DBS de Backer score, PPV proportion (%) of perfused vessels, LSCI Laser Speckle Contrast Imaging, Flux in LSPU (laser speckle perfusion units), Δ%Flux percentage difference in Flux relative to Location 1, FA Fluorescence Angiography, Fmax in arbitrary units, mean slope in arbitrary units per sec.

At Loc1, the heart rate was correlated with the OCT parameter speckle % (ρ = 0.669, *p* = 0.017). Systolic blood pressure was inversely correlated with LSCI parameter Flux (ρ = −0.451, *p* = 0.046), and stroke volume with SDF parameter MFI (ρ = −0.528, *p* = 0.017). No other significant correlations were observed for the quantitative and hemodynamic parameters.

#### **4. Discussion**

Four optical techniques were tested intraoperatively for a quantitative perfusion assessment of the gastric conduit in oesophageal cancer patients undergoing resection. The perfusion restriction towards the fundus was quantitatively assessed and was significantly objectified in 10 out of 14 parameters. At the fundus of the gastric conduit (Loc4), the FA parameters T1/<sup>2</sup> and mean slope significantly differed between the patients with or without anastomotic leakage.

The gastric conduit is mainly vascularised by the right gastroepiploic artery, which usually terminates before the future anastomotic site, making the future anastomotic site prone to vascularisation problems. Therefore, the gastric conduit is an adequate model for perfusion, as perfusion restriction is expected towards the gastric fundus. Quantitative data from this study indeed showed perfusion restriction towards the fundus. Ten quantitative parameters significantly differed between Loc1 and Loc4, of which seven decreased towards the fundus. Comparing those to the literature, the Flux and Fmax are described to indicate perfusion restriction towards the fundus [19,20]. Three parameters (Δ%Flux, influx time point (τ), and T1/2) increased towards the fundus, indicating decreased perfusion as well. Influx time point (τ) and T1/<sup>2</sup> measured the time-dependent change in fluorescence. In the literature, the time-dependent change in fluorescence intensity is a potential method for FA quantification, and increased τ and T1/<sup>2</sup> might indicate inadequate perfusion [19,21–23]. An increase in these values might be due to either diffusion through tissue instead of an intact arterial route or to resistance in the vascular network in the case of venous outflow obstruction.

Quantitative parameters were compared between patients with and without anastomotic leakage at Loc4. In patients with anastomotic leakage, a significantly higher T1/<sup>2</sup> and significantly lower mean slope were found. T1/<sup>2</sup> and mean slope provide information about the ICG influx velocity, amount of vessel density, and potentially the outflow. T1/<sup>2</sup> and mean slope are potentially hands-on parameters for intraoperative perfusion evaluation and anastomotic leakage prediction [19]. At Loc4, in contrast to an expected decrease, the SDF parameters and Flux appeared to be higher for patients with anastomotic leakage compared with patients without leakage. This increase might be explained by venous outflow obstruction. SDF measures the velocity of RBCs using a space–time diagram. The Flux measured by LSCI is dependent on the amount of blood particles and movement of these particles. In the case of venous outflow obstruction at the fundus of the gastric tube, veins swell to 5–10 times their original size [24,25]. The SDF parameters and Flux are possibly influenced by the number of RBCs or blood particles in swollen veins and therefore could be higher at Loc4 in the case of anastomotic leakage secondary to venous outflow obstruction. No differences were observed at Loc4 for the OCT parameter speckle %.

This explorative study allowed the comparison of four optical techniques in terms of their quantitative assessment and imaging characteristics. Both OCT and SDF have a limited field of view focusing on microperfusion, while LSCI and FA are wide-field (considering the heterogeneity of the capillary network). Wide-field techniques allow the immediate overview of entire gastric conduit, and potentially show vascular confines visually, which is valuable for clinical determination of the anastomotic site. Besides, the time consumption for imaging seemed lower for wide-field techniques. OCT, LSCI, and FA are all non-contact techniques, thereby avoiding stabilising problems. Only OCT is depth-resolving (depicting individual overlying vessels). In terms of quantification, quantitative SDF, LCSI, and FA had acceptable success rates. The quantitative parameter for OCT did not significantly differ between Loc1 and Loc4, nor in the case of anastomotic leakage. This observation might be explained by the method of analysis that assessed only one transection of 2.5 mm depth for quantification, indicating that limited field of view techniques might be prone to sampling errors. SDF and LSCI show perfusion restriction towards the fundus, but both show non-significant higher values at Loc4 in the

case of anastomotic leakage. As mentioned above, the SDF and LSCI trends might be comparable due to focus on the same perfusion units and might be higher in the case of anastomotic leakage due to venous outflow obstruction. FA visualises the spatial distribution of ICG that is bound to plasma proteins, and thus assesses the patency of vessels and shows vascular confines. FA parameters were able to indicate perfusion restriction towards the fundus as well as in patients with anastomotic leakage. In the case of venous congestion, FA parameters might show resistance in the vasculature, which leads to prolonged time values and a reduced rate in the increase in fluorescence intensity. Furthermore, FA time values might show what degree of arterial diffusion (instead of direct arterial flow) or venous outflow obstruction might lead to anastomotic leakage. Although FA use is dependent on ICG injection, ICG is a safe dye and can be applied at low dosages. Moreover, FA imaging systems are available for laparoscopy, which is clinically beneficial. Therefore, and on the basis of this study, FA quantitative parameters are most promising for future research.

T1/<sup>2</sup> and mean slope will be promising quantitative parameters for intraoperative clinical decisionmaking if a threshold can be derived for patient outcomes, including anastomotic leakage. According to the threshold, change in management can include the administration of medication, to shorten the gastric conduit at the maximum extent, or no anastomotic reconstruction and oesophagostomy. The latter option is a very radical solution with a decreased quality of life and need for surgery in the future to make a new reconstruction. Most surgeons would probably choose to make an extra effort to shorten the gastric conduit and make an anastomosis at increased risk instead of choosing no reconstruction as an option. Furthermore, the threshold can be used to identify high-risk patients for anastomotic leakage.

The quantitative parameters for FA may be biased by the duration of imaging, as this might not have been long enough to reach the peak value (Fmax). This study was further limited by a relatively small number of patients with a low absolute number of anastomotic leakages, and particularly measurements at Loc4 were infrequent. However, this pilot study can guide future studies, by for example serving as a base for sample size calculations. Furthermore, the exact anastomotic site was not recorded, but the anastomosis was constructed between Loc3 and Loc4. At Loc3, no differences were observed for any quantitative parameter, however at Loc4 differences were found. The site of anastomotic reconstruction may have led to bias in the observations. In addition, the system used for SDF measurements involves an older model and older software, and possibly new models (e.g., Cytocam) allow higher resolution, and new software (AVA 5, MicroVision Medical, Amsterdam, The Netherlands) might perform better. However, the SDF technique is still employed, and sampling error and increased values in the case of anastomotic leakage can still be expected. In addition, the parameters were not assessed in relation to patient risk factors. Co-morbidities such as diabetes mellitus and cardiovascular disease affect the (micro) circulation [26,27], and might therefore influence the read-out of the parameters. Furthermore, diabetes mellitus and chronic obstructive pulmonary disease are known risk factors for anastomotic leakage [28]. Owing to the small sample size, correction for co-morbidities was not performed. Nevertheless, it is difficult to compare different optical imaging systems in the literature, and this study allows comparison for the same group of patients. Furthermore, this study evaluated quantitative parameters for both perfusion and patient outcomes.

In conclusion, for future studies, FA and software-derived fluorescence–time curves are most promising to quantitatively assess decreased perfusion (A) towards the fundus in the gastric conduit and (B) in patients with anastomotic leakage. Future FA ideally has a real-time perfusion map with analysis per pixel showing perfusion thresholds. In future studies, FA data should be combined with risk factors for patient outcomes to provide additional information on FA thresholds.

**Funding:** This research received no external funding.

**Author Contributions:** Conceptualization, M.D.S., S.M.A.J., P.R.B., T.G.v.L., M.I.v.B.H., D.M.d.B., and S.S.G.; Data curation, S.M.A.J.; Formal analysis, M.D.S., S.M.A.J., P.R.B., R.M.v.d.E., and L.S.W.; Investigation, M.I.v.B.H. and S.S.G.; Methodology, S.M.A.J., L.S.W., T.G.v.L., M.I.v.B.H., D.M.d.B. and S.S.G.; Software, P.R.B. and R.M.v.d.E.; Supervision, M.I.v.B.H., D.M.d.B., and S.S.G.; Writing—original draft, M.D.S.; Writing—review and editing, S.M.A.J., P.R.B., R.M.v.d.E., L.S.W., T.G.v.L., M.I.v.B.H., D.M.d.B., and S.S.G. All authors have read and agreed to the published version of the manuscript.

**Conflicts of Interest:** M.I.v.B.H.: unrestricted grant and materials Stryker European Operations B.V.; unrestricted grant Olympus; consultant for Medtronic, Johnson & Johnson and Mylan. The other authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## **Optimal Spectral Combination of a Hyperspectral Camera for Intraoperative Hemodynamic and Metabolic Brain Mapping**

**Charly Caredda 1,\*, Laurent Mahieu-Williame 1, Raphaël Sablong 1, Michaël Sdika 1, Jacques Guyotat <sup>2</sup> and Bruno Montcel 1,\***


Received: 24 June 2020; Accepted: 21 July 2020 ; Published: 27 July 2020

**Abstract:** Intraoperative optical imaging is a localization technique for the functional areas of the human brain cortex during neurosurgical procedures. These areas are assessed by monitoring the oxygenated (HbO2) and deoxygenated hemoglobin (Hb) concentration changes occurring in the brain. Sometimes, the functional status of the brain is assessed using metabolic biomarkers: the oxidative state of cytochrome-c-oxidase (oxCCO). A setup composed of a white light source and a hyperspectral or a standard RGB camera could be used to identify the functional areas. The choice of the best spectral configuration is still based on an empirical approach. We propose in this study a method to define the optimal spectral combinations of a commercial hyperspectral camera for the computation of hemodynamic and metabolic brain maps. The method is based on a Monte Carlo framework that simulates the acquisition of the intrinsic optical signal following a neuronal activation. The results indicate that the optimal spectral combination of a hyperspectral camera aims to accurately quantify the HbO2 (0.5% error), Hb (4.4% error), and oxCCO (15% error) responses in the brain following neuronal activation. We also show that RGB imaging is a low cost and accurate solution to compute Hb maps (4% error), but not accurate to compute HbO2 (48% error) or oxCCO (1036% error) maps.

**Keywords:** hemodynamic brain mapping; metabolic brain mapping; Monte Carlo simulations; intraoperative imaging; optical imaging; hyperspectral imaging; RGB imaging

#### **1. Introduction**

Non-invasive functional brain mapping is an imaging technique used to localize the functional areas of the patient brain. This technique is used during brain tumor resection surgery to indicate to the neurosurgeon the cortical tissues that should not be removed without cognitive impairment. Functional magnetic resonance imaging (fMRI) [1] is the preoperative gold standard for the identification of the patient brain functional areas. However, after patient craniotomy, a brain shift invalidates the relevance of neuro-navigation to intraoperatively localize the functional areas of the patient brain [2]. In order to prevent any localization error, intraoperative MRI has been suggested, but it complicates the surgery gesture, which makes it rarely used. For these reasons, electrical brain stimulation [3] is preferred during neurosurgery. However, this technique suffers from limitations because the measurements could trigger epilepsy seizures. Since optical imaging combined with a quantitative modeling of brain hemodynamic biomarkers could evaluate in real time the functional areas during neurosurgery [4–6], this technique could serve as a tool of choice to complement the electrical brain stimulation.

Hyperspectral imaging allows the in vivo monitoring of the hemodynamic and metabolic status of an exposed cortex. Hyperspectral imaging provides spatially and spectrally resolved images using numerous and contiguous spectral bands [7]. In comparison, a standard color camera (or RGB camera) acquires three colors (red, green, and blue) using broad and overlapping spectral detectors. Both techniques have the ability to measure the oxygenation changes in the tissue using the modified Beer–Lambert law [5,8–12]. In functional brain mapping studies, the concentration changes of oxy- (Δ*CHbO*<sup>2</sup> ) and deoxy-hemoglobin (Δ*CHb*)can be analyzed to identify the activated cortical areas [5,13–21]. The acquisition of the intrinsic signal in the near-infrared range offers the potential to monitor the brain metabolism with the quantification of the concentration changes of the oxidative state of cytochrome-c-oxidase (Δ*CoxCCO*) [22–25]. Hyperspectral and color cameras combined with a white light illumination are simple and powerful tools for the computation of intraoperative functional brain maps. The objective is to guide the neurosurgeon during brain surgery to prevent any functional impairments after surgical procedures (tumor resection).

In the literature, all wavelength bands acquired by hyperspectral imaging setups are used to measure the hemodynamic and metabolic changes in the brain [10,24–26]. However, there are some studies in which the choice of the selected spectral bands is discussed. Bale et al. [23] showed that tens to hundreds of spectral bands acquired with a broadband near-infrared spectroscopy setup (780 nm to 900 nm) can be used to measure the oxCCO concentration changes. Arifler et al. [27] showed that eight wavelength combinations between 780 nm and 900 nm give rise to the least possible estimation errors for the deconvolution of Δ*CHbO*<sup>2</sup> , Δ*CHb*, and Δ*CoxCCO* when compared to a gold standard (121 wavelengths included between 780 nm and 900 nm). Giannoni et al. [28] proposed a Monte Carlo framework to investigate the performances of broadband spectroscopy to quantify the brain hemodynamic and metabolic responses. The results of this study indicated that eight wavelength between 780 nm and 900 nm should be selected to provide minimal differences in quantification compared to a gold standard of 121 wavelengths (780 nm to 900 nm). Sudakou et al. [29] proposed a method based on the error propagation analysis and Monte Carlo simulations (three layer model: scalp, skull, and brain) allowing the estimation of the cytochrome-c-oxidase uncertainty in data measured with a multispectral time-resolved near-infrared spectroscopy device. The results of this study indicated that 16 wavelengths between 688 and 875 nm could be used to minimize the standard deviation of the cytochrome-c-oxidase concentration changes in the brain layer. Wavelength optimization problems have also been studied for other optical imaging techniques such as near-infrared optical tomography. Chen et al. [30] identified seven laser diodes among 38 commercially available diodes in the range of 633–980 nm to estimate four chromophores (HbO2, Hb, water, and lipids) and the scattering prefactor in breast tissue. Chen et al. used the condition number and the residual norm to identify the optimal matrix used for the resolution of a linear system and thus to estimate four chromophores and the scattering prefactor in breast tissue. The optimal wavelengths were identified for a large residual norm and small condition number. The residual norm and the condition number can be interpreted as parameters representing the uniqueness and stability of the solution, respectively.

The commercial hyperspectral cameras have limited choices in the available spectral bands and do not have a spectral resolution as high as the broadband spectroscopy devices used by Bale et al. and Arifler et al. Therefore, the optimal spectral bands identified in these studies may not be available with the commercial cameras. Moreover, the more spectral bands are used, the more time is needed to compute functional brain maps. Since time is the key factor in intraoperative imaging, the smallest number of spectral bands must be acquired while ensuring minimal quantification errors.

First, we propose in this study a method to define the optimal spectral combinations of a commercial hyperspectral camera for intraoperative hemodynamic and metabolic brain mapping. This method could be used with any hyperspectral or standard RGB camera to evaluate its ability to compute accurate hemodynamic (Δ*CHbO*<sup>2</sup> and Δ*CHb*) and/or metabolic (Δ*CoxCCO*) brain maps following neuronal activation. The method is based on the Monte Carlo simulations of the acquisition of the intrinsic signal acquired by a camera. All spectral combinations of the hyperspectral camera are tested to evaluate the optimal spectral configuration that minimize the quantification errors in

Δ*CHbO*<sup>2</sup> , Δ*CHb*, and Δ*CoxCCO*. In this work, we also show that a spectral correction [10] of the reflection spectra acquired by a mosaic hyperspectral sensor is mandatory to minimize the chromophores' quantification errors. Finally, we compare standard RGB imaging and hyperspectral imaging for hemodynamic and metabolic brain mapping. We demonstrate that RGB imaging is a low cost, but not an accurate solution to identify the functional areas in a patient brain based on the analysis of the cortical hemodynamics. Hyperspectral imaging is the ideal solution for an accurate computation of hemodynamic and metabolic brain maps.

#### **2. Materials and Methods**

We simulated the acquisition of the intrinsic reflection spectra of a patient exposed cortex using hyperspectral imaging to determine the optimal spectral configuration for the computation of hemodynamic and/or metabolic brain maps. We also simulated the acquisition of the intrinsic reflection spectra based on standard RGB imaging to evaluate the potential benefit of using hyperspectral imaging.

#### *2.1. Simulated Setup*

In this study, the imaging system represented in Figure 1 was simulated. This setup was used in our previous study [5] for the identification of functional areas based on RGB imaging. A continuous wave white light source illuminated the patient exposed cortex. The light source was a halogen bulb (OSRAM Classic 116W 230V) whose spectra was measured with a USB 2000 spectrometer [31]. The intrinsic reflection spectra were acquired by a hyperspectral camera (XIMEA MQ022HG-IM-SM5X5-NIR) and standard RGB camera (BASLER acA2000-165uc). Both RGB and hyperspectral cameras acquire spectrally and spatially resolved images, but these two systems differ in their spectral resolution. The RGB camera acquires three broad and overlapping spectra over 300 nm to 700 nm using a Bayer filter [32] mounted on a CMOS sensor. The hyperspectral camera acquires 25 narrow and contiguous spectral bands in the red in the near-infrared range. A 5 × 5 mosaic filter (25 Fabry–Perot filters) mounted on a CMOS sensor was used for the acquisition of the 25 spectral bands. Each Fabry–Perot filter has two transmission peaks. The first or the second peak can be selected by mounting an interference filter on the camera lens. In our study, a long pass filter mounted on the hyperspectral camera lens selected the wavelength included between 675 nm and 975 nm.

**Figure 1.** Schematic of the simulated imaging system.

The intensity *R* measured by the spectral channel *k* of the camera (*k* ∈ [1; 3] for RGB imaging and *k* ∈ [1; 25] for hyperspectral imaging) is expressed by the following equation:

$$R\_k = \int\_{\lambda\_{\min}}^{\lambda\_{\max}} D\_k(\lambda) . \Phi(\lambda) . d\lambda. \tag{1}$$

The integral runs from *λmin* = 400 nm to *λmax* = 700 nm for RGB imaging and from *λmin* = 675 nm to *λmax* = 975 nm for hyperspectral imaging. These wavelengths delimit the ranges of the spectral sensitivity profiles provided by the camera manufacturers. *Dk* is the spectral sensitivity of spectral channel *k* of the camera, provided by the camera manufacturer. Φ is the diffuse reflection spectra of the illuminated tissue.

#### *2.2. Hyperspectral Correction*

In hyperspectral imaging based on mosaic filters, the use of the interference filter does not completely eliminate the crosstalk between the camera spectral sensitivities. To approach the ideal spectral sensitivity *Dideal* (provided by the camera manufacturer), Pichette et al. [10] proposed to compute a linear combination of the spectral sensitivities:

$$D\_m^{Corr}(\lambda) = \sum\_{k=1}^{25} X\_{m,k}.D\_k(\lambda) \tag{2}$$

with *DCorr <sup>m</sup>* the corrected spectral sensibility of the spectral channel *m* (*m* ∈ [1; 25]) and *X* the 25 × 25 correction matrix. This correction matrix *X* was computed with a linear least squares fitting:

$$\min\_{\mathbf{X}} \|D^{\text{id}\alpha \text{al}} - \mathbf{X}.D\|^2\_{\text{\textquotedblleft}}\tag{3}$$

The normalized spectral sensitivities ( *D*(*λ*).*dλ* = 1) are represented in Figure 2.

**Figure 2.** Spectral sensitivities of the hyperspectral camera XIMEA MQ022HG-IM-SM5X5-NIR. The uncorrected spectral sensitivities are plotted in solid lines on the left side of the figure and the corrected ones on the right side. The ideal spectral sensitivities are plotted in green dashed lines.

As the normalized spectral sensitivity ( *D*(*λ*).*dλ* = 1) represents a probability density function, the Kullback–Leibler divergence can be used to quantify the performance of the spectral correction [33]:

$$KL\_k(D\_k \| D\_k^{ideal}) = \int\_{\lambda\_{\rm min}}^{\lambda\_{\rm max}} D\_k(\lambda) . \log \left( \frac{D\_k(\lambda)}{D\_k^{ideal}(\lambda)} \right) d\lambda. \tag{4}$$

*KLk*(*Dk Dideal <sup>k</sup>* ) represents the Kullback–Leibler divergence (*KL* divergence) computed between the spectral sensitivity of the channel *k* (corrected or uncorrected) and the ideal spectral sensitivity of the channel *k*. The *KL* divergence represents a measure of dissimilarity between *Dk* and *Dideal <sup>k</sup>* . A value equal to 0 indicates the equality of the two probability density functions. The *KL* divergence values computed for the corrected and uncorrected spectral sensitivities are represented in Figure 3.

**Figure 3.** Measure of dissimilarities between the corrected (in orange) or uncorrected (in blue) spectral sensibilities and the ideal spectral sensitivities.

The *KL* divergence computed between the uncorrected spectral sensitivities and the ideal spectral sensitivities is on average four times higher than those computed between the corrected spectral sensitivities and the ideal spectral sensitivities.

#### *2.3. Sensor Signal-to-Noise Ratio*

The signal-to-noise ratios (*SNR*) of the spectral channel *k* of the RGB and hyperspectral cameras were measured by acquiring the white light reflected at the surface of a calibration target (plain white sheet of paper) during two minutes with an integration time set to 33 ms:

$$SNR(k) = \frac{\mu\_k}{\sigma\_k}.\tag{5}$$

*μ<sup>k</sup>* denotes the spectral channel *k* mean value obtained by computing the mean value of the temporal intensity averaged over the sensor area. *σ<sup>k</sup>* is the spectral channel *k* standard deviation value obtained by computing the standard deviation value of the temporal intensity averaged over the sensor area. The illumination source power, quantum efficiency and integration time of the cameras directly impact the SNR values. The SNR values of the RGB and hyperspectral cameras we used in this study are represented in Figure 4. The SNR of the red channel of the RGB camera is 1.1 times greater than the one of the green channel and 1.4 times greater than the one of the blue channel. The mean value of the SNR of the hyperspectral is 3.9 times lower than the SNR value of the red channel of the RGB camera. Note that these values are dependent on our experimental configuration and directly impact the amount of noise and the quantification error in the simulated concentration changes.

**Figure 4.** SNR values of each spectral channel of the RGB and hyperspectral cameras.

#### *2.4. Simulation of the Patient Cortical Activation*

Following the patient physiological activity, a hemodynamic change occurs in the activated cortical areas with an increase of the oxygenated hemoglobin concentration (*CHbO*<sup>2</sup> ) and a decrease of the deoxygenated hemoglobin concentration changes (*CHb*). These hemodynamic changes are not the only ones that occur during a patient physiological activity. The concentration of the oxidative state of cytochrome-c-oxidase (*CoxCCO*) also varies [22]. Cytochrome-c-oxidase (CCO) is an enzyme in the mitochondria that is involved in the aerobic metabolism of glucose. The total concentration of CCO does not change over a short time period (in the order of hours). However, it is possible to assess the differences between the oxidized and reduced states of CCO to obtain an indicator of the change in the CCO redox state.

The patient cortical activation was simulated by Monte Carlo simulations [34]. A volume of 60 mm × 60 mm × 60 mm of grey matter was modeled under a homogeneous white light illumination. Each voxel of the modeled tissue included the information of optical parameters (absorption *μa*, scattering *μs*, anisotropy *g* coefficients, and refractive index *n*). A white light illumination was simulated by scanning the optical parameters along the entire illumination spectrum (from 400 nm to 1000 nm in steps of 1 nm). The size of the modeled tissues was chosen in accordance with the photon sensitivity profile [5,35] computed for the detector situated at the center of the top face of the volumes. To avoid any photon loss and inexact results due to the boundary conditions (the simulation of the travel of a packet of photon stops when this packet of photon leaves the volume), the size of the models was set to 60 × <sup>60</sup> × 60 voxels with a resolution of 1 mm3. The 108 packets of photons were homogeneously illuminating the modeled surface. The optical mean path length and the diffuse reflection spectra were measured at the detector position with an integration time set to 33 ms, and the diffuse reflection spectra was normalized to a unitary source.

The optical parameters were taken from the literature and correspond to a nominal physiological condition [36–41]. The anisotropy, reduced scattering coefficients, and the refractive index used in the Monte Carlo simulations are summarized in Table 1. The scattering coefficient *μ<sup>s</sup>* was calculated from the reduced scattering coefficient (*μ*, *<sup>s</sup>*) and the anisotropy coefficient (*g*): *μ<sup>s</sup>* = *μ*, *<sup>s</sup>*/(1 − *g*). In Table 1, *λ* denotes the wavelength dependence of the reduced scattering coefficient (in nm).

**Table 1.** Anisotropy, reduced scattering coefficients, and refractive index used in the Monte Carlo simulations.


The absorption coefficients were calculated using the chemical composition of the chromophores of the simulated tissues, which are summarized in Table 2.


**Table 2.** Chemical composition of the chromophores of the simulated volumes of grey matter.

The chemical differences between the activated and non-activated grey matter are represented by a 5 <sup>μ</sup>mol·L−<sup>1</sup> increase of *CHbO*<sup>2</sup> , a 3.75 <sup>μ</sup>mol·L−<sup>1</sup> decrease of *CHb*, and a 0.5 <sup>μ</sup>mol·L−<sup>1</sup> increase of *CoxCCO*. These concentration changes are consistent with those defined in the literature. The reflection spectra and the optical mean path length were measured at the center of the top face of the volume.

#### *2.5. Modified Beer–Lambert Law*

The modified Beer–Lambert law is used to convert the acquired intrinsic reflection spectra into concentration changes. Optical functional brain maps are computed by assessing the oxyand deoxy-genated hemoglobin concentration changes (Δ*CHbO*<sup>2</sup> and Δ*CHb*) [9,10,43]. Using the near-infrared range, the oxidative state of cytochrome-c-oxidase concentration changes (Δ*CoxCCO*) can also be quantified [23]. The modified Beer–Lambert law can be expressed as a linear system [9]:

$$
\begin{bmatrix}
\Delta A\_1 \\ \vdots \\ \Delta A\_N
\end{bmatrix} = \begin{bmatrix}
E\_{1,HbO\_2} & E\_{1,Hb} & E\_{1,ox\,\text{CCO}} \\ \vdots & \vdots & \vdots \\ E\_{N,HbO\_2} & E\_{N,Hb} & E\_{N,ox\,\text{CCO}}
\end{bmatrix} \times \begin{bmatrix}
\Delta C\_{HbO\_2} \\ \Delta C\_{Hb} \\ \Delta C\_{w\,\text{CCO}}
\end{bmatrix} \tag{6}
$$

with

$$E\_{k,n} = \int \epsilon\_n(\lambda) . D\_k(\lambda) . S(\lambda) . L(\lambda) . d\lambda. \tag{7}$$

Δ*Ak* is the absorbance changes measured by the spectral channel *k* of the camera (*k* ∈ [1; *N*] with *N* ∈ [3; 25] for hyperspectral imaging and *k* ∈ [1; 3] for RGB imaging):

$$
\Delta A\_k = \log\_{10} \left( \frac{R\_k^{GM\_1}}{R\_k^{GM\_2}} \right). \tag{8}
$$

*RGM*<sup>1</sup> *<sup>k</sup>* is the intensity of the non-activated grey matter (see Section 2.4) acquired by the spectral channel *k* of the camera. *RGM*<sup>2</sup> *<sup>k</sup>* is the intensity of the activated grey matter (see Section 2.4) acquired by the spectral channel *k* of the camera. The intensities were calculated with Equation (1) using the reflection spectra simulated by Monte Carlo simulations. Δ*CHb* is the deoxygenated hemoglobin molar concentration changes (in mol·L<sup>−</sup>1). <sup>Δ</sup>*CHbO*<sup>2</sup> is the oxygenated hemoglobin molar concentration changes (in mol·L<sup>−</sup>1) and <sup>Δ</sup>*CoxCCO* the oxidative state of cytochrome-c-oxidase concentration changes (in mol·L−1). *<sup>n</sup>* is the extinction coefficient of the chromophore *<sup>n</sup>* (in L·mol−1·cm−1). The spectral sensitivity of the channel *k* is represented by *Dk*(*λ*), and *S*(*λ*) is the normalized intensity spectrum of the light source. *L*(*λ*) is the wavelength dependent mean optical path length of the photons traveling in tissue estimated by Monte Carlo simulations (in cm); see Section 2.4. The concentration changes were obtained by matrix inversion once the matrix *E* was calculated.

#### *2.6. Determination of the Optimal Spectral Configuration of the Hyperspectral Camera*

The quantification of Δ*CHb*, Δ*CHbO*<sup>2</sup> , and Δ*CoxCCO* can be achieved using at least three wavelengths [44]. With our hyperspectral camera, twenty-five spectral bands were acquired. Therefore, we can then ask ourselves what is the optimal spectral configuration for the assessment of Δ*CHb*, Δ*CHbO*<sup>2</sup> , and Δ*CoxCCO*? If *N* spectral bands were used (*N* ∈ [3; 25]), *P* combinations may be investigated:

$$\forall N \in [3; 25] \quad P = \frac{25!}{(25-N)! \times N!} \tag{9}$$

#### 2.6.1. Noise in Simulations

In order to evaluate the robustness of the chromophores' quantification, zero mean Gaussian noises were added to the simulated quantities. The Monte Carlo noises (noise on diffuse reflection spectra *φ* and mean path length *L*) were estimated for each simulated wavelength by launching the MCX software 10<sup>3</sup> times using random seeds. The seed is a number used to generate a random number with a pseudorandom number generator. For each wavelength, the standard deviations of the Monte Carlo noises were measured on these 10<sup>3</sup> simulations. The Monte Carlo noise depends on the number of packets of photons used in the simulation. In addition, the noise of the hyperspectral and RGB cameras that were experimentally measured (see Section 2.3) were added to the simulated intensities.

#### 2.6.2. Simulation Based Method

For all the possible *P* combinations for a group of *N* spectral bands (see Equation (9)), the quantification errors in Δ*CHb*, Δ*CHbO*<sup>2</sup> , or Δ*CoxCCO* were computed:

$$E\_{\Lambda\text{C}\_n}(N, p) = \left| \frac{\Delta\mathcal{C}\_n^{\text{Expected}} - \Delta\mathcal{C}\_n^{\text{Estimated}}(N, p)}{\Delta\mathcal{C}\_n^{\text{Expected}}} \right| \times 100. \tag{10}$$

Note that this quantification error is represented in absolute values. Δ*CExpected <sup>n</sup>* is the concentration changes of the chromophore *n* that were expected to be measured in the tissue; see Section 2.4. Δ*CEstimated <sup>n</sup>* (*N*, *p*) is the concentration changes of the chromophore *n* estimated from the simulated data using the combination *p* (*p* ∈ [1; *P*]) for a group of *N* spectral bands; see Equation (9). The computation of Δ*CEstimated <sup>n</sup>* (*N*, *p*) and *E*Δ*Cn* (*N*, *p*) was repeated 10<sup>3</sup> times to get their mean and the standard deviation values. Different noises were added to the simulated reflection spectra *φ*, mean path length *L*, and camera intensities *R* for each iteration. For each group of *N* spectral bands, the optimal spectral combination among the *P* possibilities was determined using Equation (11):

$$\min\_{p} \left( \sqrt{m \left( E\_{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\cdots}}}}}}}}}}}}}}}}}}}}}} \right}} \right)} \right}} \right)} \right} \right} \right} \right} \right} \right} \right} \right} \right} \right} \frac \frac{\left( \left( E\_{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\rm{\cdots}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}} \right}}} \right}}} \right}}} $$

*m* and *σ* designate the mean and standard deviation functions, respectively. This metric was chosen to identify the spectral combination of *N* spectral bands that minimize the average quantification errors while ensuring the reproduction of measurements. Equation (11) leads to the determination of the combination *p<sup>N</sup> <sup>n</sup>* , which is the best combination of *N* spectral bands for the deconvolution of the chromophore *n*. Finally, the optimal spectral configurations of the hyperspectral camera for the quantification of Δ*CHbO*<sup>2</sup> , Δ*CHb*, and Δ*CoxCCO* were obtained using Equation (12):

$$\min\_{p\_n^N} \left( \sqrt{m \left( E\_{\Lambda \mathbb{C}\_n} (N, p\_n^N) \right)^2 + \sigma \left( E\_{\Lambda \mathbb{C}\_n} (N, p\_n^N) \right)^2} \right). \tag{12}$$

Note that the spectral configurations determined by Equation (12) are not necessarily the same for the deconvolution of Δ*CHbO*<sup>2</sup> , Δ*CHb*, and Δ*CoxCCO*. In other words, this means that three distinct configurations may be used.

#### *2.7. Application to a Functional Brain Mapping Study*

In functional brain mapping studies, the measured hemodynamic concentration changes' time courses are compared to the patient theoretical hemodynamic responses. These theoretical responses are obtained by convolving the hemodynamic impulse response function to the window function representing the patient physiological stimulus. The functional areas may be identified by testing the Pearson correlation coefficients (computed between the concentration changes' time courses and the theoretical hemodynamic responses) and the concentration changes averaged over the patient activity period with several *T*-tests [5,25]. Another technique derived by BOLD fMRI (blood oxygenation level dependent, functional magnetic resonance imaging) is the NIRS SPM (near-infrared spectroscopy statistical parametric mapping) analyses [13,14] that use a general linear model and random field theory to identify the functional areas in a pixel-wise manner.

We built a simulation of the grey matter hemodynamic and metabolic responses following a patient physiological stimulus based on the Monte Carlo framework described in Section 2.4. The intrinsic signal was acquired during 60 s at 0.5 frames per second with the integration time set to 33 ms. The following paradigm was simulated: 18 s of rest, 20 s of neuronal stimulation, followed by 22 s of rest. We incorporated the temporal hemodynamic response, as well as the temporal oxCCO response following neuronal activation. The hemodynamic responses were obtained by convolving the hemodynamic impulse response function to the window function representing the patient physiological stimulus. The oxCCO response function was introduced by Wobst et al. [45] to describe the metabolic response in the primary and adjacent visual cortex for a stimulus of 6 to 24 s. The theoretical time courses can be roughly represented by the function *h*(*t*) = *a*.*t* 2.*e*−*<sup>t</sup>* convolved with the window function representing the patient physiological stimulus. In this function, *t* represents the time (in s) and *a* is the characteristic amplitude of Δ*CoxCCO*. The theoretical temporal variation of the Δ*CHbO*<sup>2</sup> , Δ*CHb*, and Δ*CoxCCO* values are represented in Figure 5. Note that these values varied homogeneously in the modeled tissue. During these absorption perturbations, the other optical coefficients were kept constant.

**Figure 5.** Theoretical Δ*CHbO*<sup>2</sup> (in red), Δ*CHb* (in blue), and Δ*CoxCCO* (in green) responses following the patient physiological activity. Δ*CHbO*<sup>2</sup> and Δ*CHb* are expressed by the left vertical axis. Δ*CoxCCO* is expressed by the right vertical axis. The blue rectangle represent the "patient" physiological stimulus.

With these simulations, we can evaluate the potential benefit from using hyperspectral imaging for the computation of hemodynamic brain maps instead of using a standard RGB camera (see our previous studies [5,11,12]). This approach aims to simulate the quantitative values that are computed in clinical functional brain mapping applications. This makes it possible to test the influence of the optimal spectral band configuration on the quantities measured in clinical applications. For this purpose, the Pearson correlation coefficient is computed. The correlation coefficient was calculated between the

theoretical hemodynamic time courses (see Figure 5) and the hemodynamic time courses simulated with the optimal spectral configuration of the hyperspectral camera (see Section 2.6). The correlation coefficient was also computed with the data simulated with RGB imaging. We can also investigate the ability of hyperspectral imaging to monitor the brain metabolism. For this purpose, the Pearson correlation coefficient values were computed between the theoretical oxCCO time course (see Figure 5) and the Δ*CoxCCO* time courses simulated with the optimal spectral configuration of the hyperspectral camera (see Section 2.6).

We also compared our results with the measurements computed with the gold standard of 121 wavelengths included between 780 nm and 900 nm [23] and with the 8 wavelengths identified by Arifler et al. [27]. For this purpose, we simulated the hemodynamic and/or metabolic monitoring using two or three chromophore systems that used these wavelengths. Note that we did not incorporate noises in these simulations of references.

#### **3. Results**

#### *3.1. Determination of the Optimal Spectral Configuration of the Hyperspectral Camera*

In Figure 6, the quantification errors in Δ*CHbO*<sup>2</sup> , Δ*CHb* and Δ*CoxCCO* (see Equation (10)) obtained for the best simulation-based deconvolution systems of *N* spectral bands (*N* ∈ [3; 25]) are represented. The solid lines stand for the mean quantification error. The colored areas represent the dispersion range of the quantification errors. The vertical dashed lines indicate the optimal number of spectral bands for the deconvolution of each chromophore. Both corrected (in green) and uncorrected spectral configuration (in red) are represented. Note that these data were obtained using noisy measurements; zeros mean Gaussian noises were added to the camera intensities and to the Monte Carlo quantities (mean path length and reflection spectra); see Section 2.6.1.

**Figure 6.** Quantification errors (see Equation (10)) in Δ*CHbO*<sup>2</sup> , Δ*CHb*, and Δ*CoxCCO* obtained for the best simulation-based deconvolution systems of *N* spectral bands (*N* ∈ [3; 25]). The solid lines stand for the mean quantification error. The colored areas represent the dispersion range of the quantification errors. The vertical dashed lines indicate the optimal number of spectral bands for the deconvolution of each chromophore. Both corrected (in green) and uncorrected spectral configuration (in red) are represented.

For the three chromophores, the quantification errors are important when a small number of spectral bands are used, but the quantification errors decrease with a larger number of bands. The quantification errors in Δ*CHbO*<sup>2</sup> are on average 30% lower when the spectral correction is applied than when it is not applied. The quantification errors in Δ*CHb* are on average 20% lower when the spectral correction is applied than when it is not applied. The quantification errors in Δ*CoxCCO* are

equivalent with and without spectral correction. We can notice that the quantification errors reach a plateau when 15 corrected and uncorrected spectral bands are used. These spectral bands are indicated by violet points in Figure 7 and will be named optimal reduced spectral bands in the rest of the paper.

**Figure 7.** Optimal spectral configuration of the hyperspectral camera XIMEA MQ022HG-IM-SM5X5- NIR obtained with the simulation-based method. The peaks of the selected spectral bands are represented as vertical lines in the molar extinction coefficients of the deconvolved chromophores. The reduced optimal spectral bands are indicated by violet points in the molar extinction coefficients of the deconvolved chromophores.

For each chromophore, different optimal spectral combinations are found. For the quantification of Δ*CHbO*<sup>2</sup> , the best spectral configuration of our hyperspectral camera is composed of 25 spectral bands if the spectral correction is applied and of 16 spectral bands if the spectral correction is not applied. These spectral bands are indicated by black vertical lines in the HbO2 molar extinction spectra in Figure 7. For these configurations, the quantification error (with Monte Carlo and intensity noise addition) is equal to 28.7% ± 21.4% (mean ± standard deviation of the 103 measurements) when the spectral correction is applied and is equal to 38.5% ± 26.9% when the spectral correction is not applied.

For the quantification of Δ*CHb*, the best spectral configuration of our hyperspectral camera is composed of 24 spectral bands if the spectral correction is applied and of 21 spectral bands if the spectral correction is not applied. These spectral bands are indicated by black vertical lines in the Hb molar extinction spectra in Figure 7. For these configurations, the quantification error (with Monte Carlo and intensity noise addition) is equal to 39% ± 29.7% (mean ± standard deviation of the 103 measurements) when the spectral correction is applied and is equal to 49.7% ± 36.6% when the spectral correction is not applied.

For the quantification of Δ*CoxCCO*, the best spectral configuration of our hyperspectral camera is composed of 21 spectral bands if the spectral correction is applied and of 25 spectral bands if the spectral correction is not applied. These spectral bands are indicated by black vertical lines in the oxCCO molar extinction spectra in Figure 7. For these configurations, the quantification error (with Monte Carlo and intensity noise addition) is equal to 251.1% ±195.4% (mean ± standard deviation of the 10<sup>3</sup> measurements) when the spectral correction is applied and is equal to 237.6% ± 183.1% when the spectral correction is not applied.

#### Quantification Performances of Deconvolution Systems

When using the optimal or the reduced optimal spectral configurations of the hyperspectral camera (see Figure 7), three linear systems have to be defined to independently quantify Δ*CHbO*<sup>2</sup> , Δ*CHb*, and Δ*CoxCCO*. In the rest of the paper, these systems will be named multiple deconvolution

systems. When inspecting the spectral band used in the optimal reduced multiple systems, twenty-one different spectral bands are used with and without spectral correction. Thus, a single system can be designed using these spectral bands to quantify Δ*CHbO*<sup>2</sup> , Δ*CHb*, and Δ*CoxCCO* in a single measurement. In the same way, when inspecting the spectral variables used in the optimal multiple systems, all spectral bands of the hyperspectral camera are used. The accuracy of the quantification of the chromophore *n* can be measured using the mean and the standard deviation of the quantification error (see Equation (10)). According to ISO 5725-1standards, the term "accuracy" stands for the trueness of the measurement (mean quantification errors close to 0%) and the repeatability of the measurement (standard deviation of quantification errors close to 0%). The mean and standard deviation values were calculated on 10<sup>3</sup> noisy measurements. These metrics were calculated using the multiple and single deconvolution systems of the hyperspectral camera; see Table 3.


**Table 3.** Quantification performances of the multiple and single deconvolution systems of the hyperspectral camera.

For each chromophore, the lowest values of the mean and standard deviation of the quantification errors in Δ*CHbO*<sup>2</sup> and Δ*CHb* are obtained with the optimal multiple deconvolution system using corrected spectral bands. When the spectral bands are not corrected, the mean and standard deviation values of the quantification errors in Δ*CHbO*<sup>2</sup> and Δ*CHb* are higher than the values obtained with corrected spectral bands. However, we can notice that the mean and standard deviation values of the quantification errors in Δ*CoxCCO* computed without spectral correction are lower than the values obtained with corrected spectral bands. The quantification errors computed with the optimal multiple system are similar to those obtained with the single system composed of all the spectral bands of the hyperspectral camera. The quantification errors computed with the optimal reduced multiple system are similar to those obtained with the single optimal reduced system.

#### *3.2. Impact of the Signal-to-Noise Ratio in the Measurements*

As mentioned in Section 2.3, the SNR values of the imaging system directly impact the amount of noise and the accuracy of the simulated quantities. To illustrate the effect of SNR on the measurements, the mean and standard deviation of the quantification errors in Δ*CHbO*<sup>2</sup> , Δ*CHb* and Δ*CoxCCO* (see Equation (10)) obtained for the optimal spectral configuration of the hyperspectral camera (see Figure 7) and the RGB camera are represented as a function of SNR; see Figures 8 and 9. For each chromophore system, the mean and standard deviation values of the quantification errors are high for low SNR values and decrease with increasing SNR values.

**Figure 8.** Mean quantification errors in Δ*CHbO*<sup>2</sup> (b), Δ*CHb* (blue curves), and Δ*CoxCCO* (green curves) (see Equation (10)) obtained for the optimal spectral configuration of the hyperspectral camera (see Figure 7) and the RGB camera as a function of SNR. The solid lines were computed using the RGB camera; the curves with circles and with triangles were computed using the hyperspectral camera with and without spectral correction, respectively. Red, green, and blue vertical lines indicate the SNR values of the R, G, and B channels of the RGB camera and black vertical lines the SNR values of the spectral channels of the hyperspectral camera.

**Figure 9.** Standard deviation of the quantification errors in Δ*CHbO*<sup>2</sup> (red curves), Δ*CHb* (blue curves), and Δ*CoxCCO* (green curves) (see Equation (10)) obtained for the optimal spectral configuration of the hyperspectral camera (see Figure 7) and the RGB camera as a function of SNR. The solid lines were computed using the RGB camera; the curves with circles and with triangles were computed using the hyperspectral camera with and without spectral correction, respectively. Red, green, and blue vertical lines indicate the SNR values of the R, G, and B channels of the RGB camera and black vertical lines the SNR values of the spectral channels of the hyperspectral camera.

When a two-chromophore system is considered with the RGB camera (SNR = 10), the quantification errors in Δ*CHbO*<sup>2</sup> and Δ*CHb* are equal to 236% ± 187% and 226% ± 179%, respectively. When the SNR value is equal to 1000, the quantification errors in Δ*CHbO*<sup>2</sup> and Δ*CHb* are equal to 52% ± 3% and 6% ± 3%, respectively. For high SNR values, it is interesting to note that the

quantification of Δ*CHbO*<sup>2</sup> is not accurate (low *σ*(*EHbO*<sup>2</sup> ) values and high *m*(*EHbO*<sup>2</sup> ) values) and the quantification of Δ*CHb* is accurate. Moreover, the mean quantification errors in Δ*CHbO*<sup>2</sup> reach a plateau for an SNR value equal to 110. For the hyperspectral camera with an SNR value equal to 10, the quantification errors in Δ*CHbO*<sup>2</sup> are equal to 232% ± 173% when the spectral correction is applied and 259% ± 193% when it is not applied. The quantification errors in Δ*CHb* are equal to 235% ± 176% when the spectral correction is applied and 267% ± 201% when it is not applied. When the SNR value is equal to 1000, the quantification errors in Δ*CHbO*<sup>2</sup> are equal to 8% ± 6% when the spectral correction is applied and 28% ± 6% when it is not applied. The quantification errors in Δ*CHb* are equal to 8% ± 6% when the spectral correction is applied and 52% ± 5% when it is not applied. For high SNR values, the quantification of Δ*CHbO*<sup>2</sup> and Δ*CHb* is accurate when the spectral correction is applied and not accurate when the spectral correction is not applied (low *σ*(*En*) values and high *m*(*En*) values). When the spectral correction is not applied, the mean quantification errors in Δ*CHbO*<sup>2</sup> and Δ*CHb* reach a plateau for an SNR value equal to 200 and 110, respectively.

When a three-chromophore system is considered with the RGB camera (SNR = 10), the quantification errors in Δ*CHbO*<sup>2</sup> , Δ*CHb*, and Δ*CoxCCO* are equal to 712% ± 531%, 822% ± 617%, and 9427% ± 7011%, respectively. When the SNR value is equal to 1000, the quantification errors in Δ*CHbO*<sup>2</sup> , Δ*CHb*, and Δ*CoxCCO* are equal to 164% ± 10%, 138% ± 12%, and 1431% ± 138%, respectively. The mean quantification errors in Δ*CHbO*<sup>2</sup> , Δ*CHb*, and Δ*CoxCCO* reach a plateau for an SNR value equal to 110, 130, and 130, respectively. For high SNR values, the quantification of Δ*CHbO*<sup>2</sup> , Δ*CHb*, and Δ*CoxCCO* is not accurate (low *σ*(*En*) values and high *m*(*En*) values). For the hyperspectral camera with an SNR value equal to 10, the quantification errors in Δ*CHbO*<sup>2</sup> are equal to 345% ± 267% when the spectral correction is applied and 328% ± 248% when it is not applied. The quantification errors in Δ*CHb* are equal to 486% ± 379% when the spectral correction is applied and 496% ± 382% when it is not applied. The quantification errors in Δ*CoxCCO* are equal to 2908% ± 2190% when the spectral correction is applied and 2951% ± 2201% when it is not applied. When the SNR value is equal to 1000, the quantification errors in Δ*CHbO*<sup>2</sup> are equal to 10% ± 8% when the spectral correction is applied and 33% ± 8% when it is not applied. The quantification errors in Δ*CHb* are equal to 14% ± 11% when the spectral correction is applied and 41% ± 11% when it is not applied. The quantification errors in Δ*CoxCCO* are equal to 90% ± 68% when the spectral correction is applied and 91 % ±70% when it is not applied. For high SNR values, the quantification of Δ*CHbO*<sup>2</sup> and Δ*CHb* is accurate when the spectral correction is applied and not accurate (low *σ*(*En*) values and high *m*(*En*) values) when the spectral correction is not applied. The quantification of Δ*CoxCCO* is not accurate with and without spectral correction (low *σ*(*EoxCCO*) values and high *m*(*EoxCCO*) values). We can notice that the measurement uncertainty of Δ*CoxCCO* is similar for all SNR values with and without spectral correction. When the spectral correction is not applied, the mean quantification errors in Δ*CHbO*<sup>2</sup> and Δ*CHb* reach a plateau for an SNR value equal to 200.

#### *3.3. Hemodynamic Monitoring*

The hemodynamic monitoring (Δ*CHbO*<sup>2</sup> and Δ*CHb* measurements) following a simulated patient neuronal stimulation (see Section 2.7) is represented In Figures 10 and 11. In Figure 10, the Δ*CHbO*<sup>2</sup> and Δ*CHb* values were computed with the optimal spectral combination of the hyperspectral camera; see Figure 7. In Figure 11, the Δ*CHbO*<sup>2</sup> and Δ*CHb* values were computed with the RGB camera. In these figures, the modified Beer–Lambert law was computed 103 times, using different Monte Carlo and intensity noise occurrences for each time iteration. When using the 121 wavelengths included between 780 nm and 900 nm [23] to quantify the concentration changes during the stimulation period (from *t* = 18 s to *t* = 38 s), there is a 12% overestimation of the Δ*CHbO*<sup>2</sup> values and a 34% underestimation of the Δ*CHb* values compared to the theoretical measurements. When using the eight wavelengths identified by Arifler et al. [27] to quantify the concentration changes during the stimulation period, there is a 13% overestimation of the Δ*CHbO*<sup>2</sup> values and a 31% underestimation of the Δ*CHb* values compared to the theoretical measurements.

**Figure 10.** Hemodynamic monitoring following neuronal activation computed with the spectral configuration of the hyperspectral camera obtained with the simulation-based method (see Figure 7). The dashed lines represent the theoretical hemodynamic responses following the neuronal stimulation. The dispersion ranges of the measurements are represented by colored areas. The concentration changes averaged over the 10<sup>3</sup> noisy measurements are represented in solid lines. The concentration changes' time courses computed with the 121 wavelengths included between 780 nm and 900 nm [23] are represented with yellow points. The ones computed with Arifler et al.'s wavelengths [27] are represented with black triangles. The blue rectangle represents the "patient" physiological stimulus.

In Figure 10, the hemodynamic monitoring following a simulated neuronal stimulation was computed using hyperspectral imaging. The optimal corrected and uncorrected spectral configurations of the hyperspectral camera were used to quantify the Δ*CHbO*<sup>2</sup> and Δ*CHb* values. The quantification dispersion ranges of the corrected and uncorrected spectral configurations have approximately the same range of magnitude. The standard deviation averaged over all time measurements for Δ*CHbO*<sup>2</sup> and <sup>Δ</sup>*CHb* are equal to 0.9 <sup>μ</sup>mol·L−<sup>1</sup> and 0.7 <sup>μ</sup>mol·L−1, respectively. When the spectral bands are corrected, the <sup>Δ</sup>*CHbO*<sup>2</sup> and <sup>Δ</sup>*CHb* values averaged over the 103 noisy measurements have a good match with the theoretical hemodynamic responses. During the stimulation period, there is on average a 14% overestimation of the Δ*CHbO*<sup>2</sup> values and a 18% underestimation of the Δ*CHb* values compared to the theoretical measurements. When the spectral bands are not corrected, there is on average a 17% underestimation of the Δ*CHbO*<sup>2</sup> values and a 58% underestimation of the Δ*CHb* values compared to the theoretical measurements. The correlation coefficients computed between the theoretical HbO2 response and the Δ*CHbO*<sup>2</sup> time courses are higher when the spectral bands are corrected (*r* = 0.951±0.013) than when they are not (*r* = 0.904±0.026). In the same way, the correlation coefficients computed between the theoretical Hb response and the Δ*CHb* time courses are higher when the spectral bands are corrected (*r* = 0.932 ± 0.018) than when they are not (*r* = 0.783 ± 0.059). The correlation coefficient values are summarized in Table 4.

In Figure 11, the hemodynamic monitoring following a simulated neuronal stimulation was computed using RGB imaging. The standard deviation averaged over all time measurements is equal to 0.35 <sup>μ</sup>mol·L−<sup>1</sup> and 0.20 <sup>μ</sup>mol·L−<sup>1</sup> for <sup>Δ</sup>*CHbO*<sup>2</sup> and <sup>Δ</sup>*CHb*, respectively. The <sup>Δ</sup>*CHbO*<sup>2</sup> values values averaged over the 10<sup>3</sup> noisy measurements do not have a good match with the theoretical HbO2 response; however, the match is rather good between the Δ*CHb* values and the theoretical Hb response. During the stimulation period, there is on average a 48% underestimation of the Δ*CHbO*<sup>2</sup> values and a 4% underestimation of the Δ*CHb* values compared to the theoretical measurements. The correlation coefficients computed between the theoretical HbO2 response and the Δ*CHbO*<sup>2</sup> time courses are equal to

*r* = 0.949 ± 0.013. The correlation coefficients computed between the theoretical Hb response and the Δ*CHb* time courses are equal to *r* = 0.989 ± 0.002. The correlation coefficient values are summarized in Table 4.

**Figure 11.** Hemodynamic monitoring following neuronal activation computed with RGB imaging. The dashed lines represent the theoretical hemodynamic responses following the neuronal stimulation. The dispersion ranges of the measurements are represented by colored areas. The concentration changes averaged over the 10<sup>3</sup> noisy measurements are represented in solid lines. The concentration changes' time courses computed with the 121 wavelengths included between 780 nm and 900 nm [23] are represented with yellow points. The ones computed with Arifler et al.'s wavelengths [27] are represented with black triangles. The blue rectangle represents the "patient" physiological stimulus.

**Table 4.** Pearson correlation coefficients computed between the theoretical hemodynamic responses and the simulated concentration changes' time courses; see Figures 10 and 11. *μ*(*rn*) and *σ*(*rn*) designate the mean and standard deviation, respectively, of the Pearson correlation coefficients computed with the 103 noisy time courses of the chromophore *n*.


#### *3.4. Hemodynamic and Metabolic Monitoring*

The hemodynamic and metabolic monitoring (Δ*CHbO*<sup>2</sup> , Δ*CHb*, and Δ*CoxCCO* measurements) following a simulated patient neuronal stimulation (see Section 2.7) are represented in Figures 12 and 13. In Figure 12, the Δ*CHbO*<sup>2</sup> , Δ*CHb*, and Δ*CoxCCO* values were computed with the optimal spectral combination of the hyperspectral camera; see Figure 7. In Figure 13, the Δ*CHbO*<sup>2</sup> , Δ*CHb*, and Δ*CoxCCO* values were computed with the RGB camera. In these figures, the modified Beer–Lambert law was computed 103 times, using different Monte Carlo and intensity noise occurrences for each time iteration. When using the 121 wavelengths included between 780 nm and 900 nm [23] to quantify the concentration changes during the stimulation period (from *t* = 18 s to *t* = 38 s), there is a 1.2% underestimation of the Δ*CHbO*<sup>2</sup> values, a 0.4% overestimation of the Δ*CHb* values, and a 1.3% overestimation of the Δ*CoxCCO* values compared to the theoretical measurements. When using the eight wavelengths identified by Arifler et al. [27] to quantify the concentration changes during the stimulation period, there is a 1.7% underestimation of the Δ*CHbO*<sup>2</sup> values, a 1.91% overestimation of the Δ*CHb* values, and a 6% overestimation of the Δ*CoxCCO* values compared to the theoretical measurements.

**Figure 12.** Hemodynamic and metabolic monitoring following neuronal activation computed with the optimal spectral configuration obtained with the simulation-based method (see Figure 7). The dashed lines represent the theoretical hemodynamic and metabolic responses following the neuronal stimulation. The dispersion ranges of the measurements are represented by colored areas. The concentration changes averaged over the 10<sup>3</sup> noisy measurements are represented in solid lines. The concentration changes' time courses computed with the 121 wavelengths included between 780 nm and 900 nm [23] are represented with yellow points. The ones computed with Arifler et al.'s wavelengths [27] are represented with black triangles. The blue rectangle represents the "patient" physiological stimulus.

In Figure 12, the hemodynamic and metabolic monitoring following a simulated neuronal stimulation was computed using hyperspectral imaging. The optimal corrected and uncorrected spectral configurations of the hyperspectral camera were used to quantify the Δ*CHbO*<sup>2</sup> , Δ*CHb*, and Δ*CoxCCO* values. The quantification dispersion ranges of the corrected and uncorrected spectral configurations have approximately the same range of magnitude. The standard deviation averaged over all time measurements for <sup>Δ</sup>*CHbO*<sup>2</sup> , <sup>Δ</sup>*CHb*, and <sup>Δ</sup>*CoxCCO* are equal to 1.23 <sup>μ</sup>mol·L<sup>−</sup>1, 1.27 <sup>μ</sup>mol·L<sup>−</sup>1, and 1.07 <sup>μ</sup>mol·L<sup>−</sup>1, respectively. When the spectral bands are corrected, the <sup>Δ</sup>*CHbO*<sup>2</sup> and <sup>Δ</sup>*CHb* values averaged over the 10<sup>3</sup> noisy measurements have a good match with the theoretical hemodynamic responses. The Δ*CoxCCO* values averaged over the 103 noisy measurements have a good match with the theoretical metabolic response with and without spectral correction. When the spectral bands are corrected, there is on average a 0.5% underestimation of the Δ*CHbO*<sup>2</sup> values, a 4.4% overestimation of the Δ*CHb* values, and a 15% overestimation of the Δ*CoxCCO* values compared to the theoretical measurements. When the spectral bands are not corrected, there is on average a 25% underestimation of the Δ*CHbO*<sup>2</sup> values, a 42% underestimation of the Δ*CHb* values, and a 19% underestimation of the Δ*CoxCCO* values compared to the theoretical measurements. The correlation coefficients computed between the theoretical HbO2 response and the Δ*CHbO*<sup>2</sup> time courses are higher when the spectral bands are corrected (*r* = 0.843 ± 0.043) than when they are not (*r* = 0.749 ± 0.068). The correlation coefficients computed between the theoretical Hb response and the Δ*CHb* time courses are higher when the spectral bands are corrected (*r* = 0.780 ± 0.063) than when they are not (*r* = 0.577 ± 0.112). The correlation coefficients computed between the theoretical oxCCO response and the Δ*CoxCCO* time

courses are higher when the spectral bands are corrected (*r* = 0.256 ± 0.177) than when they are not (*r* = 0.177 ± 0.175). The correlation coefficient values are summarized in Table 5.

**Figure 13.** Hemodynamic and metabolic monitoring following neuronal activation computed with RGB imaging. The dashed lines represent the theoretical hemodynamic and metabolic responses following the neuronal stimulation. The dispersion ranges of the measurements are represented by colored areas. The concentration changes averaged over the 10<sup>3</sup> noisy measurements are represented in solid lines. The concentration changes' time courses computed with the 121 wavelengths included between 780 nm and 900 nm [23] are represented with yellow points. The ones computed with Arifler et al.'s wavelengths [27] are represented with black triangles. The blue rectangle represents the "patient" physiological stimulus.

In Figure 13, the hemodynamic and metabolic monitoring following a simulated neuronal stimulation was computed using RGB imaging. The standard deviation averaged over all time measurements is equal to 1.24 <sup>μ</sup>mol·L−1, 1.04 <sup>μ</sup>mol·L−1, and 1.54 <sup>μ</sup>mol·L−<sup>1</sup> for <sup>Δ</sup>*CHbO*<sup>2</sup> , <sup>Δ</sup>*CHb*, and Δ*CoxCCO*, respectively. The computed Δ*CHbO*<sup>2</sup> , Δ*CHb*, and Δ*CoxCCO* values do not have a good match with the theoretical hemodynamic and metabolic responses. During the stimulation period, there is on average a 169% underestimation of the Δ*CHbO*<sup>2</sup> values, a 147% underestimation of the Δ*CHb* values, and a 1036% overestimation of the Δ*CoxCCO* values compared to the theoretical measurements. We can notice that there is an important crosstalk between the chromophores. The Δ*CHbO*<sup>2</sup> values were mainly interpreted as Hb variations, the Δ*CHb* values as oxCCO variations, and the Δ*CoxCCO* as HbO2 variations.

**Table 5.** Pearson correlation coefficients computed between the theoretical hemodynamic and metabolic responses and the simulated concentration changes' time courses; see Figures 12 and 13. *μ*(*rn*) and *σ*(*rn*) designate the mean and standard deviation, respectively, of the Pearson correlation coefficients computed with the 103 noisy time courses of the chromophore *n*.


#### **4. Discussion**

Matcher et al. [46] showed that the performance of spectroscopic analysis can be improved by increasing the number of wavelengths of illumination. The results of our study are consistent with Matcher et al.'s results; see Figure 6. The simulation-based method presented in this paper aimed to identify the optimal spectral combination of our hyperspectral camera for the brain hemodynamic and metabolic monitoring. When the spectral bands are corrected, the optimal spectral configuration is composed of 21 and 22 spectral bands; see Figure 7. This configuration could be however reduced from 10 to 12 spectral bands while keeping fairly constant performances. We can observe a plateau in the quantification errors of Δ*CHbO*<sup>2</sup> when 10 or more spectral bands are used. We also can observe a plateau in the quantification errors in Δ*CHb* and Δ*CoxCCO* when 12 or more spectral bands are used; see Figure 6. This reduction of the number of the spectral bands could be interesting for the real-time computation of hemodynamic and metabolic brain maps in the operative room using commercial hyperspectral cameras and for the conception of dedicated cameras used in functional brain mapping studies.

In this study, the optimal spectral configurations were obtained by searching the three chromophores deconvolution systems that minimize the quantification errors in Δ*CHbO*<sup>2</sup> , Δ*CHb*, and Δ*CoxCCO*. Therefore, three different spectral configurations were obtained for the deconvolution of the three chromophores. These spectral configurations, as well as the reference configuration identified by Bale et al. [23] and Arifler et al [27] are represented in Figure 14.

**Figure 14.** Peaks of spectral bands that minimize the quantification errors in Δ*CHbO*<sup>2</sup> , Δ*CHb* and Δ*CoxCCO*. The spectral configuration identified by Bale et al. [23] and Arifler et al. [27] are also indicated.

The comparison between the spectral configuration identified by Arifler and al. and ours is not trivial since our hyperspectral camera does not acquire exactly the same wavelengths; see Figure 2. There are some similarities between the wavelengths used in our system and in the system of Arifler et al. Indeed, between 780 nm and 900 nm, some wavelengths correspond. Our system also used others wavelengths inferior to 780 nm (703, 754, and 767 nm) and superior to 900 nm (916, 929, 937, 947, 953, and 958 nm). The quantification of Δ*CHbO*<sup>2</sup> , Δ*CHb*, and Δ*CoxCCO* computed with the optimal spectral bands of our hyperspectral camera are consistent with those computed with 121 wavelengths included between 780 nm and 900 nm [23] and with those computed with the eight wavelengths identified by Arifler et al. [27]. For a two-chromophore deconvolution system, the quantification error in Δ*CHbO*<sup>2</sup> measured with our system is equivalent to those measured with 121 wavelengths. However, the choice of our spectral bands aims to reduce the quantification error in Δ*CHb*; see Figure 10. For a three-chromophore deconvolution system, the quantification errors measured with our system are a little higher than those measured with 121 wavelengths or with Arifler et al.'s wavelengths; see

Figure 12. This difference may be explained because the illumination and the acquisition of the 121 wavelengths of reference, as well as the eight wavelengths identified by Arifler et al. were simulated as ideal sources and detectors. This means that for these simulations, the term *D* × *S* in Equation (7) is equal to one. To efficiently compare our spectral configuration with these spectral configurations of reference, the spectral sensitivities of the detectors and the illumination sources have to be considered. However, the acquisition condition is different since we acquired the intrinsic signal of an exposed cortex, whereas the wavelengths identified in Bale et al.'s and Arfiler et al.'s studies were identified for functional near-infrared spectroscopy devices.

We showed that a spectral correction (see Section 2.2) is required when a hyperspectral mosaic sensor is used [10]. The spectral correction aims to reduce the quantification errors in the measurements and allows a better follows-up of the temporal hemodynamic and metabolic variations; see Table 5.

We also compared hyperspectral imaging to RGB imaging for the computation of hemodynamic and metabolic brain maps. Hyperspectral imaging is the suitable solution to compute hemodynamic maps and metabolic maps thanks to its ability to acquire the intrinsic signal in the near-infrared range [23]. A very important crosstalk between HbO2, Hb, and oxCCO can be observed when the RGB camera is used for the computation of hemodynamic and metabolic brain maps; see Figure 13. Therefore, RGB imaging is not a suitable solution to compute metabolic brain maps. An RGB camera could be however a low-cost solution to compute hemodynamic maps; see Figure 11. This solution is not accurate to quantify Δ*CHbO*<sup>2</sup> , but is very accurate for the quantification of Δ*CHb*. This original result is interesting because most surgical microscopes used in the operating room are equipped with standard RGB cameras. It is known that the BOLD signal used in fMRI studies is predominately due to the paramagnetic properties of deoxygenated hemoglobin [47]. This result indicates that the Δ*CHb* quantified with RGB imaging can be used in a robust way for intraoperative functional mapping based on SPM analyses.

We incorporated the camera noises in our simulation to find the most robust and reliable spectral configuration. In Figures 8 and 9, we show that the SNR of the imaging system has a drastic impact on quantification performances. An accurate quantification could be obtained with a high SNR value. Therefore, the light source and the camera specifications and settings have to be carefully chosen in order to guarantee an optimal SNR. This simulation framework could be a great tool for industry or researchers working on intraoperative functional brain mapping solutions to help them in the choice of commercial camera. However, our simulation framework needs to be improved. For the moment, a homogeneous volume of grey matter was considered. A realistic mapping of the exposed cortex could be simulated as suggested by Giannoni et al. [28]. This will be considered in future studies. Moreover, the hemodynamic and metabolic changes following neuronal activation were homogeneously simulated in the volume of grey matter. These events are obviously not consistent with those appearing in a real cortical tissue. This modeling also has to be taken into account to improve our method of identification of the optimal spectral bands of a hyperspectral camera for brain hemodynamic and metabolic monitoring.

#### **5. Conclusions**

We present in this paper a method for the identification of the optimal spectral bands of a commercial camera for the intraoperative monitoring of the brain hemodynamic and metabolic responses following neuronal activation. The method described in the report is based on Monte Carlo simulations of the light propagation in a volume of grey matter and incorporates a realistic modeling of the camera acquisition with the addition of Gaussian noises (experimentally measured with the cameras). We identified that an optimal spectral combination of our hyperspectral camera composed of 21 to 22 spectral bands can be used to compute accurate hemodynamic and metabolic brain maps. This configuration could be however reduced from 10 to 12 spectral bands while keeping fairly constant performances, which is consistent with the spectral configurations proposed in the literature. We also showed that RGB imaging is not a suitable technique to compute metabolic brain maps, but is very accurate to compute hemodynamic maps with the quantification of deoxygenated hemoglobin concentration changes. Our Monte Carlo framework needs to be improved, namely with the consideration of the perfusion of grey matter by blood capillaries.

**Author Contributions:** Conceptualization: C.C., L.M.-W., and B.M. Methodology: C.C., L.M.-W., R.S., M.S., and B.M. Software: C.C., L.M.-W., and M.S. Writing, original draft: C.C. Writing, review and editing: C.C., L.M.-W., R.S., M.S., J.G., and B.M. Supervision: R.S. and B.M. Funding acquisition: R.S. and B.M. Investigation: J.G. Resources: J.G. Project administration: B.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** These works were funded by LABEX PRIMES (ANR-11-LABX-0063) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007), operated by the French National Research Agency (ANR); Cancéropôle Lyon Auvergne Rhône Alpes (CLARA) within the program "OncoStarter"; Infrastructures d'Avenir en Biologie Santé (ANR-11-INBS-000), within the program "Investissements d'Avenir" operated by the French National Research Agency (ANR) and France Life Imaging.

**Acknowledgments:** We want to acknowledge the PILoT facility for the support provided for the image acquisition.

**Conflicts of Interest:** No conflicts of interest, financial or otherwise, are declared by the authors.

#### **Abbreviations**

The following abbreviations are used in this manuscript:

HbO2 Oxygenated hemoglobin Hb Deoxygenated hemoglobin oxCCO Oxidative state of cytochrome-c-oxidase BOLD Blood oxygenation level dependent fMRI functional magnetic resonance imaging NIRS Near infrared spectroscopy SPM Statistical parametric mapping SNR Signal-to-noise ratio

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Applied Sciences* Editorial Office E-mail: applsci@mdpi.com www.mdpi.com/journal/applsci

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 www.mdpi.com