1.

Introduction

Near-infrared spectroscopy (NIRS) uses light in the red and near-infrared range to monitor concentration changes of oxy- (${\mathrm{HbO}}_2$) and deoxy- (HHb) hemoglobin in the brain.¹ Under the assumption of an intact neurovascular coupling, these can be considered as an indirect measure of brain activity. When multiple sources and detectors are used, arranged in overlapping and multidistance channels, depth-resolved images of these cerebral hemoglobin variations can be reconstructed.^2,3 This technique is usually referred to as diffuse optical tomography (DOT). Eggebrecht et al.² have recently shown that by using a high-density DOT system, DOT can be employed to map distributed brain networks and function. Their maps disclosed a strong correspondence with the same maps obtained with functional magnetic resonance imaging (fMRI). Their results are extremely promising since DOT techniques are characterized by noninvasiveness and portability. They can be easily employed to monitor vulnerable subjects, such as neonates,^4–6 patients with implanted devices,² or critically ill patients who need bedside monitoring⁷ and who could not be previously monitored with fMRI. Both functional NIRS (fNIRS) and fMRI rely on the assumption of intact neurovascular coupling, although this assumption is frequently violated in cases of brain pathology.^8–11 Pathology (e.g., acute brain injury, cerebral ischemia, and neurodegenerative disease) can alter the regional hemodynamic response to localized changes in metabolism and in extreme cases can lead to absence or inversion of the response. Although it is well recognized that an inverted fNIRS hemodynamic response [or an fMRI negative blood-oxygen-level dependent (BOLD) response] can be seen in both healthy volunteers and after brain injury, the exact physiological mechanisms underlying these responses remain unclear. Additional information about colocalized cerebral metabolic state, such as that provided by oxCCO, can help elucidate these pathways. Therefore, there is an unmet clinical need for a direct marker of cerebral metabolism.

Cytochrome c oxidase (CCO), the terminal enzyme in the mitochondrial respiratory chain, is responsible for more than 95% of cellular oxygen metabolism.¹ The redox state of CCO can be measured with broadband NIRS¹² and as such represents an opportunity to directly measure cellular oxidative metabolism. The change in the oxidation state of CCO ($\mathrm{\Delta }[\mathrm{oxCCO}]$) reflects the balance between electron acceptance originating from nicotinamide adenine dinucleotide and donation to oxygen to form water, driving the mitochondrial proton electrochemical gradient and $>90\%$ of adenosine triphosphate production.¹³ The relationship between cerebral oxygen and metabolic substrate delivery and CCO has been extensively investigated in animals,¹⁴ neonatal asphyxia,¹⁵ adult brain injury,¹⁶ and healthy adults,¹⁷ including during functional activation.¹⁸ Importantly, $\mathrm{\Delta }[\mathrm{oxCCO}]$ is concordant with measurements of metabolism, including invasively measured lactate/pyruvate ratio^16,19 and magnetic resonance spectroscopy-derived lactate.¹⁵ However, because the concentration of CCO is an order of magnitude lower than that of hemoglobin, an optimized optical technique to measure $\mathrm{\Delta }[\mathrm{oxCCO}]$ using broadband NIRS^20,21 is typically required.

Present commercial fNIRS systems do not have the capability to measure $\mathrm{\Delta }[\mathrm{oxCCO}]$, and in-house systems developed explicitly to measure $\mathrm{\Delta }[\mathrm{oxCCO}]$ provide only a limited number of channels because of hardware cost and complexity, and therefore cannot perform DOT.¹³ Reconstructing images of $\mathrm{\Delta }[\mathrm{oxCCO}]$ is a natural progression of the technique following the development of DOT for displaying $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ and $\mathrm{\Delta }[\mathrm{HHb}]$. Recovering the spatial distribution of this additional chromophore may well facilitate the investigation of cerebral energy status alongside functional activation, particularly in the context of brain pathology.²² However, the small $\mathrm{\Delta }[\mathrm{oxCCO}]$ signal that occurs in the presence of larger changes in $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ and $\mathrm{\Delta }[\mathrm{HHb}]$ means that imaging it is a significant challenge.

The recovery of images of concentration changes via DOT constitutes a nonlinear, ill-posed inverse problem.^23–25 This process involves three steps. First, a geometric representation of the target object is required. A subject-specific, multilayer, anatomically accurate head model is preferred, but when the magnetic resonance image of the subject is unavailable, a multilayer atlas head model can be employed.²⁶ The atlas head model can be registered to the subject’s cranial landmarks, and wavelength- and tissue-specific optical properties can be assigned to each tissue type [e.g., scalp, skull, cerebrospinal fluid (CSF), gray matter (GM), and white matter (WM)], thus providing an anatomically meaningful solution space. Second, a model of how light propagates through these head tissues, given an array of sources and detectors, is required, so as to relate changes in optical properties to changes in the acquired signals. This forward model can be constructed either with finite-element method (FEM) approaches (e.g., Toast++²⁷ and Nirfast²⁸) or with Monte Carlo approaches (e.g., MMC²⁹ and MCX³⁰). Third, the solution of the forward model must be inverted so as to solve the inverse problem and reconstruct an image of concentration changes. This step is the most computationally expensive, since it is an underdetermined and ill-posed problem. Corlu et al.³¹ demonstrated a multispectral approach, which directly reconstructs changes in concentration (rather than absorption changes at each wavelength). This approach reduces the number of unknowns and better constrains the inverse problem.

Most DOT studies to date have recovered images of $\mathrm{\Delta }[{\mathrm{HbO}}_2]$, $\mathrm{\Delta }[\mathrm{HHb}]$, and of changes in total hemoglobin concentration $\mathrm{\Delta }[\mathrm{HbT}]$ (where $\mathrm{\Delta }[\mathrm{HbT}]=\mathrm{\Delta }[{\mathrm{HbO}}_2]+\mathrm{\Delta }[\mathrm{HHb}])$, since these are the most studied chromophores in clinical and neuroscience settings.^2,4,5 Other chromophores, such as water and lipids, have been successfully reconstructed with DOT,^31,32 mostly for breast imaging applications.

The aim of this work is to provide, for the first time, a proof of concept demonstration of the feasibility of reconstructing images of oxCCO concentration changes, along with $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ and $\mathrm{\Delta }[\mathrm{HHb}]$, in the adult brain. Data were acquired from four healthy volunteers with a broadband system during a visual stimulation paradigm, and images were reconstructed with a multispectral approach. Visual stimulation was employed since it is known to elicit a strong localized activation with a highly repeatable pattern.^33,34 Limitations and assumptions of the present work will be discussed.

2.

Material and Methods

2.1.

Instrumentation

A detailed description of the broadband system used in this study can be found elsewhere.¹⁸ The system has two 50-W halogen light sources, each delivering light to the participant via optical fibers. Only one source can be illuminated at a time. To collect the diffusely reflected light and measure the intensity at each wavelength, the system has eight detector fibers divided into two bundles and two CCD cameras (Pixis 512, Princeton Instruments, Trenton, New Jersey), each connected to a spectrograph. The sampling frequency is 1 Hz.

The optical interface between the system and the participant was based on a three-dimensional (3-D)-printed optode holder. The optode holder incorporated 4 source slots (in the middle row) and 14 detector slots arranged in 2 rows (see Fig. 1).

Fig. 1

Array layout. Sources were positioned in the central row (red dots), while detectors were located in the bottom and top rows (blue dots). Eight channels per source could be recorded, for a total of 32 channels. Note that not all slots could be filled with a source or a detector during a single acquisition, since only one source and eight detectors were available at a time. Each acquisition was repeated four times, with the available source and detector fibers positioned in different slots for each acquisition so as to cover the entire imaging array over four acquisitions in each subject.

2.2.

Participants and Paradigm

Four healthy volunteers (two males and two females; age range 27 to 48 years old) participated in the visual stimulation study after providing written informed consent. Participants were seated in a comfortable chair in a dimly lit room and were asked to focus their attention at the center of the computer monitor, located at a distance of 60 cm in front of the subject. Visual functional activation was achieved with a 4-Hz inverting checkerboard covering the full visual field. The paradigm consisted of 20 s of stimulation followed by 20 s of blank screen repeated over 10 epochs. To monitor the visual cortex, the optode array was fixed horizontally to cover the participant’s left occipital cortex, with the fourth source location positioned over Oz (according to the 10-20 EEG positions) (Fig. 2). The eight nearest detector slots in the array (four above and four below each activated source) were populated with the eight available detector fibers and were used to collect intensity data continuously during the 10 epochs of visual stimulation. The remaining detector slots remained empty during that acquisition. This was repeated for each of the four source slots by translating the fibers over the fixed optode array on the head and repeating the stimulation. A 3-D Patriot^™ Digitizer (Polhemus, Colchester, Vermont) was employed to acquire, for each participant, the coordinates of their cranial landmarks (nasion, inion, preauricular points, and Cz) in addition to the optode locations. The study was approved by University College London's Research Ethics Committee (1133/001).

Fig. 2

Optode holder placed on the left occipital cortex with one source (central row) and eight detectors.

2.5.

Residual Computation

Residuals were computed to investigate the spectral changes due to the oxCCO chromophore. Channel-wise attenuation changes were deduced from the concentration changes in the voxel space (i.e., before their mapping to the volumetric mesh and then to the GM surface mesh, which was performed for visualization purposes only). Attenuation changes were determined for both the concentration changes computed using the three-chromophore model and the changes computed using the two-chromophore model by multiplying the concentration changes voxel-wise data with the multispectral Jacobian.

The difference between the three- and the two-chromophore fits was then investigated. Ideally, if the model explains the data, the residuals should be multivariate Gaussian noise with covariance ${\mathrm{\Gamma }}_e$, i.e., the model accounts for all nonrandom variability in the data. Given that image reconstruction is an ill-posed inverse problem that requires regularization, it is likely that even if the model explains the data, the residuals would not equate to the covariance of the assumed noise model. However, a model that better explains the data will have residuals that tend toward Gaussian noise. If the difference between the residuals of the two models has a nonrandom distribution across wavelengths, with the residuals of the three-chromophore model smaller than those of the two-chromophore models, this would suggest that the two-chromophore model is not accounting for a chromophore that is evident in the data. This analysis of residuals has been previously used in spectroscopy studies recovering $\mathrm{\Delta }[{\mathrm{HbO}}_2]$, $\mathrm{\Delta }[\mathrm{HHb}]$, and $\mathrm{\Delta }[\mathrm{oxCCO}]$.^15,18

Residuals were computed for each channel at the time of the peak of the expected functional activation (i.e., 20 s).

All analyses and reconstruction procedures were performed in MATLAB^® (Mathworks, Natick, Massachusetts). Figure 3 shows the data analysis and image reconstruction procedures.

Fig. 3

Flow diagram summarizing the main steps of the data analysis and image reconstruction procedures.

3.

Results

In Fig. 4, examples of reconstructed $\mathrm{\Delta }[{\mathrm{HbO}}_2]$, $\mathrm{\Delta }[\mathrm{HHb}]$, and $\mathrm{\Delta }[\mathrm{oxCCO}]$ images on the GM surface mesh at five different time points are displayed for one representative participant (participant 1). The selected time points are 2, 10, 20, 30, and 40 s after stimulus presentation. Figure 5 shows the images obtained at 20-s poststimulus onset for the other three participants.

Fig. 4

Examples of reconstructed $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ (first row), $\mathrm{\Delta }[\mathrm{HHb}]$ (second row), and $\mathrm{\Delta }[\mathrm{oxCCO}]$ (third row) images on the GM surface mesh at five different time points for participant 1. The selected time points are 2, 10, 20, 30, and 40 s after stimulus presentation. Images were mapped to the MNI152 GM mesh for visualization purposes.

Fig. 5

Examples of reconstructed $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ (first row), $\mathrm{\Delta }[\mathrm{HHb}]$ (second row), and $\mathrm{\Delta }[\mathrm{oxCCO}]$ (third row) images on the GM surface mesh for participant 2, 3, and 4. The reconstruction for 20-s poststimulus onset is displayed.

Figure 6 shows the grand-average channel-wise recovered $\mathrm{\Delta }[{\mathrm{HbO}}_2]$, $\mathrm{\Delta }[\mathrm{HHb}]$, and $\mathrm{\Delta }[\mathrm{oxCCO}]$ changes obtained using the spectroscopy approach over the whole 740- to 900-nm spectrum. Only some of the channels exhibit a strong activation (defined as an increase in ${\mathrm{HbO}}_2$ with a concomitant decrease in HHb), which is in line with the focal and localized change in chromophores reconstructed in Figs. 4 and 5.

Fig. 6

Grand-average hemodynamic responses in each channel to the functional activation task. Error bars represent the standard deviation across participants. Only the upper error bar for ${\mathrm{HbO}}_2$ and the lower error bar for HHb are reported for visualization purposes. $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ is reported in red, $\mathrm{\Delta }[\mathrm{HHb}]$ in blue, and $\mathrm{\Delta }[\mathrm{oxCCO}]$ in green. These responses were obtained using the data from the whole broadband spectrum (740 to 900 nm) and the UCLn algorithm.

Figure 7 shows the grand-average difference between the residuals of the two- and three-chromophore models in all channels. Residuals were computed as the differences between the measured data and the back-projected attenuation changes after image reconstruction. The spectra in most of the active channels do not have a random distribution. Instead, they show a shape approximating the spectrum of the oxCCO chromophore, with a broad peak at around 810 nm. In most of the channels that do not show activation, the difference between the residuals of the two models is almost nil. There is some variability in the residuals among participants (as there is in the hemodynamic responses, see Fig. 6) likely due to: (1) each participant can have a different timing of brain activation (all residuals are computed at time point 20 s) and (2) the retinotopic mapping of the visual cortex is different in each person, increasing the chance that each channel is probing slightly different brain visual areas in each participant.⁴⁶

Fig. 7

Grand-average residual difference between the two- and three-chromophore model in each channel and for each wavelength at time point 20. Error bars represent the standard deviation across participants.

4.

Discussion

In this study, we explored the feasibility of reconstructing images of $\mathrm{\Delta }[\mathrm{oxCCO}]$, in addition to $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ and $\mathrm{\Delta }[\mathrm{HHb}]$, using data acquired with a broadband NIRS system. This NIRS system allowed us to measure broadband data from overlapping and multidistance channels. With this configuration, we were able to perform DOT on broadband data and reconstruct images of $\mathrm{\Delta }[\mathrm{oxCCO}]$ for the first time. Data were acquired in four healthy volunteers during a visual stimulation task, which is known to elicit a strong activation, with a highly reproducible response.^33,34

We evaluated the reconstructed images both qualitatively and quantitatively. All four participants showed a clear and focused activation for the three chromophores in the left occipital cortex during the period of maximum response to the visual stimulation. These results are in line with the channel-wise spectroscopy data, which show a clear response to the task in a subset of channels. Furthermore, by employing a two-chromophore model, we reconstructed from the same data images of $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ and $\mathrm{\Delta }[\mathrm{HHb}]$ only. We then computed the residuals for both the two- and the three-chromophore models as the difference between the original attenuation changes and the back-projected attenuation changes. The difference between these residuals provides insight on whether a third chromophore is apparent in the data. The residual difference spectra exhibited a shape similar to that of the oxCCO spectrum, a common finding in NIRS studies, thus suggesting that the functionally induced changes in oxCCO concentration are influencing our data and can be imaged via DOT.

It should be noted that while $\mathrm{\Delta }[\mathrm{oxCCO}]$, $\mathrm{\Delta }[{\mathrm{HbO}}_2]$, and $\mathrm{\Delta }[\mathrm{HHb}]$ are normally recovered from NIRS data assuming a linear relationship between the changes in concentration and attenuation, in DOT, a nonlinear, ill-posed, underdetermined inverse problem must be solved. Therefore, the residuals in DOT should be interpreted with caution. Solving the inverse problem involves a balance between our trust in the data and our confidence in the prior knowledge. Residuals can be influenced by the fact that the problem is ill-posed and underdetermined, and by uncertainties in the prior knowledge. Indeed, the back-projected attenuation changes are unlikely to be identical to the measurements. Regularization is required to improve the conditioning of the inverse problem and to cope with noise amplification that can occur in the inversion process due to small singular values in the Jacobian matrix. First-order Tikhonov regularization was utilized here to avoid suppressing the mean value of the images and to prevent the build-up of excessive high-frequency components, which cannot be resolved using diffuse optical methods. By sphering the solution space prior to reconstruction, we ensured that the effect of the regularization was applied equally to each concentration value.

The optimal wavelength range for NIRS studies of functional changes in tissue is still under debate, 740 to 900 nm has been investigated,³⁴ with preference for 780 to 900 nm to recover $\mathrm{\Delta }[\mathrm{oxCCO}]$, $\mathrm{\Delta }[{\mathrm{HbO}}_2]$, and $\mathrm{\Delta }[\mathrm{HHb}]$^15,18 because of evidence of distortion over the peak of the HHb extinction spectrum between 740 and 780 nm. Arifler et al.²¹ reported the best combinations of wavelengths to recover $\mathrm{\Delta }[\mathrm{oxCCO}]$, $\mathrm{\Delta }[{\mathrm{HbO}}_2]$, and $\mathrm{\Delta }[\mathrm{HHb}]$ in spectroscopy when the full broadband spectrum cannot be acquired. With a combination of three wavelengths, the mean recovery error versus the gold standard (i.e., using the full spectrum) was 10%, which was reduced to 2% when using a combination of eight wavelengths. In all combinations, the lowest selected wavelength was around 784 to 785 nm. However, previous DOT studies aimed at recovering only $\mathrm{\Delta }[{\mathrm{HbO}}_2]$, and $\mathrm{\Delta }[\mathrm{HHb}]$ found that the best combination of wavelengths contained at least one wavelength below 780 nm, where measurements are more sensitive to HHb.^37,42 Further studies should be performed to evaluate which wavelength combination is optimum to recover all three chromophores with DOT and to further investigate this issue. In this study, we selected 17 wavelengths from the broadband spectrum, separated at intervals of 10 nm. This choice was a compromise between the computational burden and the accuracy of the reconstructed image^32,47 when using the multispectral approach proposed by Corlu et al.³¹ The computation time for the inversion process for one frame was $\sim 101\mathrm{s}$ with an Intel^® Core^™ processor i7-4790 at 3.60 GHz.

Previous studies have shown the feasibility of recovering DOT images of $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ and $\mathrm{\Delta }[\mathrm{HHb}]$ in real time.^25,48 Real-time imaging could provide valuable clinical information, which can be more naturally interpreted by clinicians. Considering the clinical importance of oxCCO, delivering real-time imaging of $\mathrm{\Delta }[\mathrm{oxCCO}]$, $\mathrm{\Delta }[{\mathrm{HbO}}_2]$, and $\mathrm{\Delta }[\mathrm{HHb}]$ would be a great advance in the optical imaging field. Although from a theoretical perspective this should be readily achievable, there will likely be challenges to overcome when recovering images of all three chromophores in real time, mainly due to the higher number of wavelengths required and the associated increase in computational burden. Further studies should investigate the optimal number of wavelengths that can reliably yield $\mathrm{\Delta }[\mathrm{oxCCO}]$ images, and whether the computational burden of the reconstruction could be reduced to allow images of $\mathrm{\Delta }[\mathrm{oxCCO}]$, $\mathrm{\Delta }[{\mathrm{HbO}}_2]$, and $\mathrm{\Delta }[\mathrm{HHb}]$ to be delivered in real time, perhaps using a dynamic Kalman filter estimator.^49,50

This initial study will hopefully promote further research to optimize image reconstruction methods to deliver oxCCO changes and inform instrument development. Further work is required to explore the validity of the assumptions made in this study (e.g., the prior knowledge on the model covariance matrix, the selected wavelengths, the settings employed to solve the inverse problem, etc.), and the image reconstruction would benefit from further validation with simulations, and with phantom and in vivo studies. Ideally, the procedure should be verified using pathological conditions that induce changes in oxCCO but not in ${\mathrm{HbO}}_2$ or HHb concentration, which would reveal evidence of cross talk during the reconstruction of the three chromophores. Perhaps this could be achieved using an animal model of cardiopulmonary bypass with fluorocarbon perfusion.^51,52

In this study, we have demonstrated for the first time the reconstruction of images of simultaneous changes in concentration of oxCCO, ${\mathrm{HbO}}_2$, and HHb in the healthy adult brain during a visual stimulation task. A clear and focal activation was reconstructed in the left occipital cortex for all three chromophores in all participants. Reconstructing images of $\mathrm{\Delta }[\mathrm{oxCCO}]$ may provide a unique opportunity to investigate cerebral energy status in addition to functional activation, across the cortex, which may have particular implications for studies of brain pathology.