**Functional Near Infrared Spectroscopy and Diffuse Optical Tomography in Neuroscience**

Matteo Caffini1, Davide Contini1, Rebecca Re1, Lucia M. Zucchelli1, Rinaldo Cubeddu1, Alessandro Torricelli1 and Lorenzo Spinelli2 *1Dipartimento di Fisica - Politecnico di Milano, Milano 2CNR - Istituto di Fotonica e Nanotecnologie, Milano Italy* 

#### **1. Introduction**

50 Advances in Brain Imaging

Tsujii, T., Masuda, S., Akiyama, T. & Watanabe, S. (2010a). The role of inferior frontal cortex

Tsujii, T., Okada, M. & Watanabe, S. (2010b). Effects of aging on hemispheric asymmetry in

Tsujii, T., Yamamoto, E., Ohira, T., Takahashi, T. & Watanabe, S. (2010d). Antihistamine

Tsujii, T., Sakatani, K., Masuda, S., Akiyama, T., & Watanabe, S. (2011a). Evaluating the roles

Tsujii, T., Sakatani, K., Nakashima, E., Igarashi, T., & Katayama, Y. (2011b). Characterization

Zhu, Z., Zhang, J.X., Wang, S., Xiao, Z., Huang, J. & Chen, H.C. (2009). Involvement of left

May; Vol.67, No.1 (May 2010), pp.80-85, ISSN: 0168-0102

pp.2005-2008, ISSN: 0028-3932

pp.178-183, ISSN: 0166-4328

No.2 (August 2009), pp.756-763.

8119

1932.

in belief-bias reasoning: an rTMS study. *Neuropsychologia*, Vol.48, No.7 (June 2010),

inferior frontal cortex activity during belief-bias syllogistic reasoning: a nearinfrared spectroscopy study. *Behavioural Brain Research*, Vol.210, No.2 (July 2010),

effects on prefrontal cortex activity during working memory process in preschool children: a near-infrared spectroscopy (NIRS) study. *Neuroscience Research*. 2010

of the inferior frontal gyrus and superior parietal lobule in deductive reasoning: An rTMS study. *Neuroimage*, Vol.58, No.2 (September 2011), pp.640-646, ISSN: 1053-

of the acute effects of alcohol on asymmetry of inferior frontal cortex activity during a Go/No-Go task using functional near-infrared spectroscopy. *Psychopharmacology*, Vol.217, No.4 (October 2011), pp.595-603, ISSN: 0033-3158 Xue, G., Aron, A.R. & Poldrack, R.A. (2008). Common neural substrates for inhibition of

spoken and manual responses. *Cerebral Cortex*, Vol.18, No.8 (August 2008), 1923-

inferior frontal gyrus in sentence-level semantic integration. *Neuroimage*, Vol.47,

Neuroimaging is the branch of medicine whose purpose is to provide visual information about the structure and the anatomy of the brain. The main techniques in clinics are: *Computed Tomography* (CT), *Magnetic Resonance Imaging* (MRI), *Diffuse Optical Tomography* (DOT), *Positron Emission Tomography* (PET) and *Single Photon Emission Computed Tomography* (SPECT).

CT scanning uses X-rays crossing the sample to image sections of the specimen in study. Specimen can be a living being, a part of if (e.g. the abdomen, a knee, the head, ...) or whatever non-living object. X-rays travel ballistically inside most of the materials (living tissue included), so measuring absorption of X-rays we can guess the composition of the sample we are measuring. Changing the direction of injection of the X-rays and merging absorption data coming from multiple directions makes a planar reconstruction of the examined section of the sample possible. CT is invasive in the sense that irradiates the patient with ionizing radiation. A little dose is given to the patient during a single CT session, anyway. CT scanning of the head is typically used to detect skull fractures, brain injuries, aneurysms, strokes, brain tumors and arteriovenous malformations in the brain.

MRI is based on nuclear magnetic resonance principles. It uses a strong static magnetic field to align the nuclear magnetization of hydrogen atoms of water in the body and then radio frequency fields are generated to alter this magnetization alignment. Several coils mounted on the scanner are then able to detect the magnetic field produced by the altered hydrogen atom magnetization and to relate the recovery time of these short-lasting magnetic fields to the environment in which resonant hydrogen atoms lie. MRI, using non-ionizing radiation, is generally considered non-invasive for the patient and provides greater contrast between different soft tissues than CT does. Magnetic resonance images of the head are mostly used to detect brain tumors.

**HbO2**

d*I* = −*μ<sup>a</sup> I*d*z* (1)

(−*μa*) d*z* (2)

<sup>−</sup>*μaL* (3)

*μ<sup>a</sup>* = [*c*] *�* (4)

Fig. 1. The absorption spectrum, that is the molar extinction coefficient vs. the wavelength, of the oxygenated (red) and deoxygenated (blue) states of the hemoglobin. *I*<sup>1</sup> and *I*<sup>2</sup> are the

Functional Near Infrared Spectroscopy and Diffuse Optical Tomography in Neuroscience 53

When a radiation beam of known intensity *I*<sup>0</sup> crosses a layer of known width *L* of a certain medium (chromophore), part of the radiation is absorbed by the medium and the rest is transmitted on the other side of the layer (Fig. 2). Modeling this phenomenon supposing that each infinitesimal layer absorbs an amount of radiation d*I* and that this absorption is proportional to radiation intensity, brings to the following mathematical model. Calling *μ<sup>a</sup>* the absorption coefficient and z the radiation beam direction, the radiation intensity *I* measured

> *L* 0

The absorption coefficient is measured in m−<sup>1</sup> and has a straightforward interpretation:

*<sup>μ</sup><sup>a</sup>* is the mean free path a photon travels prior to being absorbed. Moreover, it can be decomposed into two distinct factors: the absorbent medium concentration and its molar

*I* = *I*0*e*

chosen wavelengths for Near InfraRed Spectroscopy measurements.

 *I I*0 d*I <sup>I</sup>* <sup>=</sup>

is modeled by the *Lambert-Beer law*:

*la* = <sup>1</sup>

extinction coefficient *�*:

2

**Hb** 

PET and SPECT are nuclear medicine techniques and they are used in cancer localization. A radionuclide is injected in the body via blood flow and then tracked with detectors placed around the patient.

Functional neuroimaging is a particular way to perform medical neuroimaging focused more on functionality rather than just resolve anatomical features. The goal of functional neuroimaging is to detect *spatial* and *temporal* changes of activated areas of the brain, by measuring related physiological or physical features. The most common techniques are *functional Magnetic Resonance Imaging* (fMRI), *Electroencephalography* (EEG), *Magnetoencephalography* (MEG) and, recently, *functional Near InfraRed Spectroscopy* (fNIRS) and *Diffuse Optical Tomography* (DOT). PET and SPECT can be considered functional images as they use a radionuclide concentration to distinguish between healthy tissue and tumoral tissue.

Usually functional neuroimaging techniques are divided into two main classes: those with high spatial resolution (∼ 1 cm or less) and those with high time resolution (∼ 1 s or less).


Having good time resolution and sufficient spatial resolution, DOT lays between the two classes and it is often considered the optimal compromise for functional neuroimaging studies. Spatial resolution in DOT is not an intrinsic feature, being the technique diffusion-limited. Nevertheless, in brain activation studies, integrating optical data and *a priori* anatomic data from MRI scans, brings to huge enhancements in spatial localization.

#### **2. Functional near infraRed spectroscopy (fNIRS)**

The simplest way to perform tissue oximetry using harmless electromagnetic radiation is to shine the tissue by means of a continuous beam of infrared light and to collect the re-emitted or transmitted light.

Hemoglobin is the key: oxygenated and deoxygenated states have considerably different absorption spectra and it turns out that, in first approximation, oxyhemoglobin and deoxyhemoglobin also are the two main chromophores in tissues into the InfraRed window. Because the underlining hypothesis states the presence of two chromophores, we need to inject radiation at two different wavelengths, possibly chosen where the spectra have the greatest differences (see Fig. 1). This spectroscopic technique uses constant light intensity and is known in literature as *Continuous Wave Near InfraRed Spectroscopy* (CW-NIRS) The most diffused wavelengths for CW-NIRS are 690 nm and 820 nm.

#### **2.1 Absorption spectroscopy**

Spectroscopy is the study of the interaction between radiation and matter, in particular absorption spectroscopy measures the absorption of radiation, in a selected range of frequency, in radiation-matter interaction. It is largely employed in analytical chemistry to check for the presence of elements or substances, being the absorption spectrum a sort of fingerprint and so characteristic of the substances. The simplest application is the detection of the amount of the substance present.

2 Will-be-set-by-IN-TECH

PET and SPECT are nuclear medicine techniques and they are used in cancer localization. A radionuclide is injected in the body via blood flow and then tracked with detectors placed

Functional neuroimaging is a particular way to perform medical neuroimaging focused more on functionality rather than just resolve anatomical features. The goal of functional neuroimaging is to detect *spatial* and *temporal* changes of activated areas of the brain, by measuring related physiological or physical features. The most common techniques are *functional Magnetic Resonance Imaging* (fMRI), *Electroencephalography* (EEG), *Magnetoencephalography* (MEG) and, recently, *functional Near InfraRed Spectroscopy* (fNIRS) and *Diffuse Optical Tomography* (DOT). PET and SPECT can be considered functional images as they use a radionuclide concentration to distinguish between healthy tissue and tumoral tissue. Usually functional neuroimaging techniques are divided into two main classes: those with high spatial resolution (∼ 1 cm or less) and those with high time resolution (∼ 1 s or less).

• High spatial resolution and poor time resolution techniques (fMRI, PET and SPECT).

Having good time resolution and sufficient spatial resolution, DOT lays between the two classes and it is often considered the optimal compromise for functional neuroimaging studies. Spatial resolution in DOT is not an intrinsic feature, being the technique diffusion-limited. Nevertheless, in brain activation studies, integrating optical data and *a priori* anatomic data from MRI scans, brings to huge enhancements in spatial localization.

The simplest way to perform tissue oximetry using harmless electromagnetic radiation is to shine the tissue by means of a continuous beam of infrared light and to collect the re-emitted

Hemoglobin is the key: oxygenated and deoxygenated states have considerably different absorption spectra and it turns out that, in first approximation, oxyhemoglobin and deoxyhemoglobin also are the two main chromophores in tissues into the InfraRed window. Because the underlining hypothesis states the presence of two chromophores, we need to inject radiation at two different wavelengths, possibly chosen where the spectra have the greatest differences (see Fig. 1). This spectroscopic technique uses constant light intensity and is known in literature as *Continuous Wave Near InfraRed Spectroscopy* (CW-NIRS) The most

Spectroscopy is the study of the interaction between radiation and matter, in particular absorption spectroscopy measures the absorption of radiation, in a selected range of frequency, in radiation-matter interaction. It is largely employed in analytical chemistry to check for the presence of elements or substances, being the absorption spectrum a sort of fingerprint and so characteristic of the substances. The simplest application is the detection of

• High time resolution and poor spatial resolution techniques (EEG and MEG).

**2. Functional near infraRed spectroscopy (fNIRS)**

diffused wavelengths for CW-NIRS are 690 nm and 820 nm.

around the patient.

or transmitted light.

**2.1 Absorption spectroscopy**

the amount of the substance present.

Fig. 1. The absorption spectrum, that is the molar extinction coefficient vs. the wavelength, of the oxygenated (red) and deoxygenated (blue) states of the hemoglobin. *I*<sup>1</sup> and *I*<sup>2</sup> are the chosen wavelengths for Near InfraRed Spectroscopy measurements.

When a radiation beam of known intensity *I*<sup>0</sup> crosses a layer of known width *L* of a certain medium (chromophore), part of the radiation is absorbed by the medium and the rest is transmitted on the other side of the layer (Fig. 2). Modeling this phenomenon supposing that each infinitesimal layer absorbs an amount of radiation d*I* and that this absorption is proportional to radiation intensity, brings to the following mathematical model. Calling *μ<sup>a</sup>* the absorption coefficient and z the radiation beam direction, the radiation intensity *I* measured is modeled by the *Lambert-Beer law*:

$$\mathbf{d}I = -\mu\_a I \mathbf{d}z \tag{1}$$

$$\int\_{I\_0}^{I} \frac{\mathbf{d}I}{I} = \int\_0^L (-\mu\_a) \,\mathrm{d}z \tag{2}$$

$$I = I\_0 e^{-\mu\_d L} \tag{3}$$

The absorption coefficient is measured in m−<sup>1</sup> and has a straightforward interpretation: *la* = <sup>1</sup> *<sup>μ</sup><sup>a</sup>* is the mean free path a photon travels prior to being absorbed. Moreover, it can be decomposed into two distinct factors: the absorbent medium concentration and its molar extinction coefficient *�*:

$$
\mu\_a = [c]\,\epsilon\tag{4}
$$

*coefficient*, is then introduced. The Lambert-Beer law it is written as follows:

*Turbid Medium (CW Light)*

(a) Absorption Spectroscopy using a turbid medium does not allow to compute the chromophore concentration within the sample because of scattering contributes

to intensity attenuation.

reflection geometry (Fig. 3(b)).

explained.

*ls* = <sup>1</sup>

(Fig. 3(b)).

*I* = *I*0*e*

The scattering coefficient here defined has a similar interpretation as the absorption coefficient:

Functional Near Infrared Spectroscopy and Diffuse Optical Tomography in Neuroscience 55

*<sup>μ</sup><sup>s</sup>* is the mean free path between two scattering events. Within this class of substances, absolute concentration measurements are possible only for thin samples (for thicker samples more complicated approaches are necessary, see section 2.3.1). Incidentally it can be noticed that this disturbing factor provides a helpful effect. Radiation scatters statistically in every direction, even backwards. It is then possible to collect radiation from the same side of the sample used for the radiation injection and to design experiments in reflectance geometry

<sup>−</sup>(*μa*+*μs*)*<sup>L</sup>* (6)

*Turbid Medium (CW Light)*

(b) When using turbid media experiments in reflection geometry are possible. Scattering provides a positive probability to have radiation traveling backwards after encountering scattering centers in

the medium.

Fig. 3. Absorption spectrocopy in turbid medium in trasmission geometry (Fig. 3(a)) and

A different approach is to irradiate the sample with a non-constant light intensity, but it becomes necessary to develop a more general model for light traveling inside turbid media. This will be done in section 2.3. First a simpler way to separate scattering contributes is

In eq. 6 it has been seen that both absorption and scattering take part in the same way to radiation attenuation and that is not trivial to discriminate between the two distinct contributions. There's at least a case in which a simple way the separation is easily done. Suppose to have a medium in which the chromophore concentration [*c* (*t*)] slightly varies in time and suppose to sample spectroscopic data in reflection geometry for a certain time interval. This case is not that far from reality as it may seem at first sight, for example in

Fig. 2. Absorption Spectroscopy in the most simple geometry. Radiation is injected in the sample from one side and radiation intensity is measure on the other side. The Lambert-Beer law then allows to compute the chromophore concentration.

The extinction coefficient depends on the radiation wavelength: *�* = *�* (*λ*) and so the absorption coefficient *μ<sup>a</sup>* = *μ<sup>a</sup>* (*λ*). Such a functional dependence constitutes the *absorption spectrum* of the substance considered1. When more than a chromophore is present, each contribute is summed:

$$\mu\_{\mathfrak{a}} = \left[c\_1\right]\mathfrak{e}\_1 + \left[c\_2\right]\mathfrak{e}\_2 + \left[c\_3\right]\mathfrak{e}\_3 + \dots = \sum\_{i} \left[c\_i\right]\mathfrak{e}\_i \tag{5}$$

Media for which this model holds are non-diffusing media, often called *clear media* and concentration measurements are easily conducted knowing the extinction coefficient of the chromophore and the length of the path traveled by the radiation. In the many chromophore experiment, performing multiple sessions, changing the radiation wavelength, allows to write a linear system and to obtain the unknown concentration of each chromophore.

Performing the same experiment described before using different media, it is found that for a vast class of substances this model is not accurate, as it doesn't take into account for the diffusion of radiation within the medium (Fig. 3(a)). Scattering events are due to discontinuities in dielectric properties of the medium and can be modeled as collisions between the photons of the radiation and particles in the medium. Many scattering models have been developed in physics: in classical electromagnetism by Rutherford, Thomson and Mie and in quantum physics by Compton and Raman [Feynman et al. (1964); Mie (1908); Rutherford (1911)]. A second important class of media then have to be described, media often called *turbid media*. Living tissues belong to this second class.

The Lambert-Beer model is still valid for this media, but the absorption coefficient alone is no more enough to describe the intensity drop found. A second coefficient *μs*, called *scattering*

<sup>1</sup> In many cases the absorption spectrum is given with respect to the frequency *ν* = *<sup>c</sup> <sup>n</sup><sup>λ</sup>* , being *c* the speed of light and *n* the refraction index of the medium.

*coefficient*, is then introduced. The Lambert-Beer law it is written as follows:

$$I = I\_0 \mathfrak{e}^{-(\mu\_s + \mu\_s)L} \tag{6}$$

The scattering coefficient here defined has a similar interpretation as the absorption coefficient: *ls* = <sup>1</sup> *<sup>μ</sup><sup>s</sup>* is the mean free path between two scattering events. Within this class of substances, absolute concentration measurements are possible only for thin samples (for thicker samples more complicated approaches are necessary, see section 2.3.1). Incidentally it can be noticed that this disturbing factor provides a helpful effect. Radiation scatters statistically in every direction, even backwards. It is then possible to collect radiation from the same side of the sample used for the radiation injection and to design experiments in reflectance geometry (Fig. 3(b)).

4 Will-be-set-by-IN-TECH

*Clear Medium (CW Light)*

Fig. 2. Absorption Spectroscopy in the most simple geometry. Radiation is injected in the sample from one side and radiation intensity is measure on the other side. The Lambert-Beer

*<sup>μ</sup><sup>a</sup>* = [*c*1] *�*<sup>1</sup> + [*c*2] *�*<sup>2</sup> + [*c*3] *�*<sup>3</sup> + ··· = ∑

a linear system and to obtain the unknown concentration of each chromophore.

<sup>1</sup> In many cases the absorption spectrum is given with respect to the frequency *ν* = *<sup>c</sup>*

Media for which this model holds are non-diffusing media, often called *clear media* and concentration measurements are easily conducted knowing the extinction coefficient of the chromophore and the length of the path traveled by the radiation. In the many chromophore experiment, performing multiple sessions, changing the radiation wavelength, allows to write

Performing the same experiment described before using different media, it is found that for a vast class of substances this model is not accurate, as it doesn't take into account for the diffusion of radiation within the medium (Fig. 3(a)). Scattering events are due to discontinuities in dielectric properties of the medium and can be modeled as collisions between the photons of the radiation and particles in the medium. Many scattering models have been developed in physics: in classical electromagnetism by Rutherford, Thomson and Mie and in quantum physics by Compton and Raman [Feynman et al. (1964); Mie (1908); Rutherford (1911)]. A second important class of media then have to be described, media often

The Lambert-Beer model is still valid for this media, but the absorption coefficient alone is no more enough to describe the intensity drop found. A second coefficient *μs*, called *scattering*

*i*

[*ci*] *�<sup>i</sup>* (5)

*<sup>n</sup><sup>λ</sup>* , being *c* the speed

The extinction coefficient depends on the radiation wavelength: *�* = *�* (*λ*) and so the absorption coefficient *μ<sup>a</sup>* = *μ<sup>a</sup>* (*λ*). Such a functional dependence constitutes the *absorption spectrum* of the substance considered1. When more than a chromophore is present, each

law then allows to compute the chromophore concentration.

called *turbid media*. Living tissues belong to this second class.

of light and *n* the refraction index of the medium.

contribute is summed:

(a) Absorption Spectroscopy using a turbid medium does not allow to compute the chromophore concentration within the sample because of scattering contributes to intensity attenuation.

(b) When using turbid media experiments in reflection geometry are possible. Scattering provides a positive probability to have radiation traveling backwards after encountering scattering centers in the medium.

Fig. 3. Absorption spectrocopy in turbid medium in trasmission geometry (Fig. 3(a)) and reflection geometry (Fig. 3(b)).

A different approach is to irradiate the sample with a non-constant light intensity, but it becomes necessary to develop a more general model for light traveling inside turbid media. This will be done in section 2.3. First a simpler way to separate scattering contributes is explained.

In eq. 6 it has been seen that both absorption and scattering take part in the same way to radiation attenuation and that is not trivial to discriminate between the two distinct contributions. There's at least a case in which a simple way the separation is easily done. Suppose to have a medium in which the chromophore concentration [*c* (*t*)] slightly varies in time and suppose to sample spectroscopic data in reflection geometry for a certain time interval. This case is not that far from reality as it may seem at first sight, for example in

this technique would be to measure what happens within the blood flow, neglecting the

Functional Near Infrared Spectroscopy and Diffuse Optical Tomography in Neuroscience 57

Collecting data at different times and then computing the absorbance variation values Δ*Aλ*<sup>1</sup>

*�Hb*,*λ*<sup>1</sup> *�*HbO2,*λ*<sup>1</sup> *�Hb*,*λ*<sup>2</sup> *�*HbO2,*λ*<sup>1</sup>

> *�Hb*,*λ*<sup>1</sup> *�*HbO2,*λ*<sup>1</sup> *�Hb*,*λ*<sup>2</sup> *�*HbO2,*λ*<sup>1</sup>

Being this a two equations systems with the two unknowns Δ [Hb] and Δ [HbO2], a solution to the linear system is possible and the oxyhemoglobin and deoxyhemoglobin concentration variations with respect to a reference value can be computed (Eq. 18). If data are sampled for an interval of time, solving the linear system for each time sample, the time course of the

Derived quantities such as total hemoglobin [tHb] are often plotted vs. time in order to

An important limitation of CW-NIRS is the possibility to obtain just concentration variations. Having absolute concentration values, hemoglobin saturation sO2 could be calculated:

sO2 <sup>=</sup> [HbO2]

In principle there is no depth sensitivity with a single source-detector pair in CW-NIRS. Information on the composition of the media crossed by the CW light can be measured as an average. To increase light penetration large source-detector separation have been used (4-5 cm) with the obvious consequence of increasing the sampling volume and thus worsening the spatial resolution. Increasing the complexity of the system, the use of multi-source and multi-detector tomographic or topographic arrangements could overcome this limitation and

Laboratory prototypes and commercial instrumentations have been developed through the years to measure oxygenation of the muscles and of the cortical regions of the brain. A

better depth sensitivity [Boas et al. (2004); Strangman et al. (2002)].

comprehensive list can be found in Wolf et al. (2007).

*�Hb*,*λ*1Δ [Hb] + *�*HbO2,*λ*1Δ [HbO2]

*�Hb*,*λ*2Δ [Hb] + *�*HbO2,*λ*2Δ [HbO2]

� *Dλ*<sup>1</sup> *L*

� *Dλ*<sup>2</sup> *L*

�

� � Δ [Hb]

�−<sup>1</sup> ⎡ ⎢ ⎣

Δ [HbO2]

Δ*Aλ*<sup>1</sup> *Dλ*<sup>1</sup> *L* Δ*Aλ*<sup>2</sup> *Dλ*<sup>2</sup> *L*

Δ [tHb] = Δ [Hb] + Δ [HbO2] (19)

[tHb] (20)

⎤ ⎥

<sup>⎦</sup> (18)

(16)

(17)

contributes from the surrounding environment.

In matrix form it becomes:

concentrations is given.

improve the data visualization.

and finally:

� <sup>Δ</sup>*Aλ*<sup>1</sup> <sup>=</sup> �

<sup>Δ</sup>*Aλ*<sup>2</sup> <sup>=</sup> �

⎡ ⎢ ⎣

Δ*Aλ*<sup>1</sup> *Dλ*<sup>1</sup> *L* Δ*Aλ*<sup>2</sup> *Dλ*<sup>2</sup> *L*

� Δ [Hb] Δ [HbO2]

⎤ ⎥ <sup>⎦</sup> <sup>=</sup>

�

� = �

and Δ*Aλ*<sup>2</sup> , it is possible to write the following system of equations:

biological tissues chromophores concentrations change in time because of blood flow. Then for each instant *tn* a Lambert-Beer equation can be written as follow:

$$I\left(t\_{\rm tr}\right) = I\_0 e^{-\left(\left[c\left(t\_{\rm tr}\right)\right]\varepsilon DL + G\right)}\tag{7}$$

where scattering contribution *μsL* has been grouped into the term *G* and the *differential pathlength factor D* has been introduced to take into account that photon travel paths longer than the source-detector separation because they penetrate inside the medium, scatters multiple times and then a little part reaches the detector traveling backwards. The differential pathlength factor depends only on the mean refractive index of the crossed medium. Hypothesizing that *G* doesn't vary in time and is constant for small variations of chromophore concentrations and introducing the physical quantity *absorbance A* = ln *<sup>I</sup> <sup>I</sup>*<sup>0</sup> = −([*c*] *�DL* + *G*), we can compute the variation in chromophore concentration between measurement collected at time *tn* and time *tm* in the following way:

$$
\Delta A\_{\rm mm} = A\_{\rm m} - A\_{\rm n} \tag{8}
$$

$$\mathbf{f} = -\mathbf{c}\_{\text{m}}\boldsymbol{\varepsilon}\mathbf{D}L - \mathbf{G} + \mathbf{c}\_{\text{n}}\boldsymbol{\varepsilon}\mathbf{D}L + \mathbf{G} \tag{9}$$

$$=-(\mathcal{c}\_{m} - \mathcal{c}\_{n})\varepsilon DL\tag{10}$$

$$=-\Delta c\_{mn}\varepsilon DL\tag{11}$$

The scattering contributes vanishes and the variation in chromophore concentration depends only on the two intensity measurements collected:

$$
\Delta \mathcal{L}\_{mn} = -\frac{\Delta A\_{mn}}{\epsilon D L} \tag{12}
$$

$$\epsilon = \frac{1}{\varepsilon DL} \left( A\_{\text{fl}} - A\_{\text{m}} \right) \tag{13}$$

$$=\frac{1}{\epsilon DL} \left( \ln \frac{I\_n}{I\_0} - \ln \frac{I\_m}{I\_0} \right) \tag{14}$$

$$=\frac{1}{\epsilon DL}\ln\left(\frac{I\_{\text{fl}}}{I\_{\text{m}}}\right)\tag{15}$$

Most of times a reference measurement is conducted and all the concentration values are then computed as variation from the reference value.

#### **2.2 Continuous wave NIRS (CW-NIRS)**

Hemodynamic activity of the brain is related through metabolic processes to the electrical activity of the neurons, then, performing hemoglobin concentration measurements on the cortical tissues, it is possible to obtain information about spatial localization and temporal behavior of neuronal activity. The spectroscopic tools described in section 2.1 provide a straightforward way to do this. It is possible to inject radiation into the tissues and collect the photons coming out, for example using a reflectance geometry. The usefulness of this technique would be to measure what happens within the blood flow, neglecting the contributes from the surrounding environment.

Collecting data at different times and then computing the absorbance variation values Δ*Aλ*<sup>1</sup> and Δ*Aλ*<sup>2</sup> , it is possible to write the following system of equations:

$$\begin{cases} \Delta A\_{\lambda\_1} = \left( \epsilon\_{Hb, \lambda\_1} \Delta \left[ \text{Hb} \right] + \epsilon\_{\text{HbO}\_2, \lambda\_1} \Delta \left[ \text{HbO}\_2 \right] \right) D\_{\lambda\_1} L \\ \Delta A\_{\lambda\_2} = \left( \epsilon\_{Hb, \lambda\_2} \Delta \left[ \text{Hb} \right] + \epsilon\_{\text{HbO}\_2, \lambda\_2} \Delta \left[ \text{HbO}\_2 \right] \right) D\_{\lambda\_2} L \end{cases} \tag{16}$$

In matrix form it becomes:

$$
\begin{bmatrix}
\frac{\Delta A\_{\lambda\_1}}{D\_{\lambda\_1}L} \\
\frac{\Delta A\_{\lambda\_2}}{D\_{\lambda\_2}L}
\end{bmatrix} = \begin{bmatrix}
\mathfrak{e}\_{Hb, \lambda\_1} \ \mathfrak{e}\_{HbO\_2, \lambda\_1} \\
\mathfrak{e}\_{Hb, \lambda\_2} \ \mathfrak{e}\_{HbO\_2, \lambda\_1}
\end{bmatrix} \begin{bmatrix}
\Delta \left[\text{Hb}\right] \\
\Delta \left[\text{HbO}\_2\right]
\end{bmatrix} \tag{17}
$$

and finally:

6 Will-be-set-by-IN-TECH

biological tissues chromophores concentrations change in time because of blood flow. Then

where scattering contribution *μsL* has been grouped into the term *G* and the *differential pathlength factor D* has been introduced to take into account that photon travel paths longer than the source-detector separation because they penetrate inside the medium, scatters multiple times and then a little part reaches the detector traveling backwards. The differential pathlength factor depends only on the mean refractive index of the crossed medium. Hypothesizing that *G* doesn't vary in time and is constant for small variations of chromophore

we can compute the variation in chromophore concentration between measurement collected

The scattering contributes vanishes and the variation in chromophore concentration depends

 ln *In I*0

Most of times a reference measurement is conducted and all the concentration values are then

Hemodynamic activity of the brain is related through metabolic processes to the electrical activity of the neurons, then, performing hemoglobin concentration measurements on the cortical tissues, it is possible to obtain information about spatial localization and temporal behavior of neuronal activity. The spectroscopic tools described in section 2.1 provide a straightforward way to do this. It is possible to inject radiation into the tissues and collect the photons coming out, for example using a reflectance geometry. The usefulness of

 *In Im* 

<sup>−</sup> ln *Im I*0 

<sup>Δ</sup>*cmn* <sup>=</sup> <sup>−</sup>Δ*Amn*

<sup>=</sup> <sup>1</sup>

<sup>=</sup> <sup>1</sup> *�DL*

<sup>=</sup> <sup>1</sup> *�DL* ln

<sup>−</sup>([*c*(*tn*)]*�DL*+*G*) (7)

Δ*Amn* = *Am* − *An* (8)

= −*cm�DL* − *G* + *cn�DL* + *G* (9)

= −(*cm* − *cn*)*�DL* (10)

= −Δ*cmn�DL* (11)

*�DL* (12)

*�DL* (*An* <sup>−</sup> *Am*) (13)

*<sup>I</sup>*<sup>0</sup> = −([*c*] *�DL* + *G*),

(14)

(15)

for each instant *tn* a Lambert-Beer equation can be written as follow:

*I* (*tn*) = *I*0*e*

concentrations and introducing the physical quantity *absorbance A* = ln *<sup>I</sup>*

at time *tn* and time *tm* in the following way:

only on the two intensity measurements collected:

computed as variation from the reference value.

**2.2 Continuous wave NIRS (CW-NIRS)**

$$
\begin{bmatrix}
\Delta \begin{bmatrix} \mathbf{H} \mathbf{b} \end{bmatrix} \\
\Delta \begin{bmatrix} \mathbf{H} \mathbf{b} \mathbf{O}\_{2} \end{bmatrix}
\end{bmatrix} = \begin{bmatrix}
\mathbf{c}\_{Hb\boldsymbol{\beta}\lambda\_{1}} & \mathbf{c}\_{\mathbf{H}\mathbf{b} \mathbf{O}\_{2}\lambda\_{1}} \\
\mathbf{c}\_{Hb\boldsymbol{\beta}\lambda\_{2}} & \mathbf{c}\_{\mathbf{H}\mathbf{b} \mathbf{O}\_{2}\lambda\_{1}}
\end{bmatrix}^{-1} \begin{bmatrix}
\frac{\Delta A\_{\lambda\_{1}}}{\mathbf{D}\_{\lambda\_{1}}L} \\
\Delta A\_{\lambda\_{2}} \\
\overline{\mathbf{D}\_{\lambda\_{2}}L}
\end{bmatrix} \tag{18}
$$

Being this a two equations systems with the two unknowns Δ [Hb] and Δ [HbO2], a solution to the linear system is possible and the oxyhemoglobin and deoxyhemoglobin concentration variations with respect to a reference value can be computed (Eq. 18). If data are sampled for an interval of time, solving the linear system for each time sample, the time course of the concentrations is given.

Derived quantities such as total hemoglobin [tHb] are often plotted vs. time in order to improve the data visualization.

$$
\Delta\left[\text{tHb}\right] = \Delta\left[\text{Hb}\right] + \Delta\left[\text{HbO}\_2\right] \tag{19}
$$

An important limitation of CW-NIRS is the possibility to obtain just concentration variations. Having absolute concentration values, hemoglobin saturation sO2 could be calculated:

$$\text{sO}\_2 = \frac{[\text{HbO}\_2]}{[\text{tHb}]} \tag{20}$$

In principle there is no depth sensitivity with a single source-detector pair in CW-NIRS. Information on the composition of the media crossed by the CW light can be measured as an average. To increase light penetration large source-detector separation have been used (4-5 cm) with the obvious consequence of increasing the sampling volume and thus worsening the spatial resolution. Increasing the complexity of the system, the use of multi-source and multi-detector tomographic or topographic arrangements could overcome this limitation and better depth sensitivity [Boas et al. (2004); Strangman et al. (2002)].

Laboratory prototypes and commercial instrumentations have been developed through the years to measure oxygenation of the muscles and of the cortical regions of the brain. A comprehensive list can be found in Wolf et al. (2007).

[Ishimaru (1978); Martelli et al. (2010)]:

Each term appearing has a precise significance:

*<sup>∂</sup><sup>t</sup>* ⇒ temporal change of energy. • −sˆ · ∇*I* ⇒ change due to energy flow.

) *I* (r, ˆs�

• *q* = *q* (r, ˆs, *t*) ⇒ Radiation source inside d*V*

The RTE has two immediate properties:

but scattered into the direction ˆs, being *p* (sˆ, ˆs�

*<sup>∂</sup><sup>t</sup>* <sup>=</sup> <sup>−</sup>s<sup>ˆ</sup> · ∇*<sup>I</sup>* <sup>−</sup> (*μ<sup>a</sup>* <sup>+</sup> *<sup>μ</sup>s*) *<sup>I</sup>* <sup>+</sup> *<sup>μ</sup><sup>s</sup>*

• − (*μ<sup>a</sup>* + *μs*) *I* ⇒ energy drop due to absorption and scattering.

that a photon flowing in the direction ˆs is scattered into the direction ˆs�

*q*<sup>0</sup> (r, ˆs) *δ* (*t*), then if the absorption coefficient is *μ<sup>a</sup>* the solution is:

� 4*π p* � sˆ, ˆs� � *I* � r, ˆs� , *t* �

Functional Near Infrared Spectroscopy and Diffuse Optical Tomography in Neuroscience 59

�

1. If *I*<sup>0</sup> is a solution for a non-absorbing medium with time impulsive source term *q* (r, ˆs, *t*) =

*I* = *I*0*e*

This property is a generalization of the Lambert-Beer law [Sassaroli & Fantini (2004);

*I* (r, ˆs, *t*)

Because of the high complexity of the RTE, no analytical solutions are available and approximation methods are then applied. Numerical methods like the *Finite Difference Method* or the *Finite Elements Method* (FEM) require a large amount of computational power, but are often applied [Arridge (1995); Arridge et al. (2000); Arridge & Hebden (1997); Arridge & Schweiger (1995); Arridge et al. (1993)]. Stochastic methods like the *Monte Carlo* are widely

Many ways to simplify the RTE have been developed through the years, but most methods are based on the *PN approximation*. Two derived quantities of interest are the *photon fluence* Φ

2. If I is the Green function for a medium having extinction coefficient *μ<sup>t</sup>* = *μ<sup>a</sup>* + *μ<sup>s</sup>* and

, then for a medium having extinction coefficient *μ*∗

�<sup>3</sup>

Wm−3sr−<sup>1</sup>

, *t*) dΩ� ⇒ energy gain from radiation coming from every direction,

� .

⎧ ⎨ ⎩ **r**∗ = � *μ<sup>t</sup> μ*∗ *t* � **r**

*t* ∗ = � *μ<sup>t</sup> μ*∗ *t* � *t*

) the *phase function*, that is the probability

<sup>−</sup>*μavt* (25)

*<sup>t</sup>* and the same albedo,

(26)

.

dΩ� + *q* (24)

1 *v ∂I*

• <sup>1</sup> *v ∂I*

• +*μ<sup>s</sup>*

�

<sup>4</sup>*<sup>π</sup> p* (sˆ, ˆs�

Tsuchiya (2001)].

*μt*

**2.3.2 The PN approximation**

the solution is scaled in this way:

*I* (r∗, ˆs, *t*

∗) =

� *μ*<sup>∗</sup> *t μt*

This scaling property is known as the *similarity principle* [Zege et al. (1991)].

used in many biological applications [Boas et al. (2002); Fang & Boas (2009)].

albedo *a* = *<sup>μ</sup><sup>s</sup>*

#### **2.3 Photon migration approach**

A different way to deal with diffusion and discriminate between scattering and absorption contributes is to use non-constant light intensities in the injected radiation. A possible way to describe the physics of the problem is to start from the electric and magnetic fields associated with the radiation and to develop a rigorous mathematical theory, taking into account scattering, diffraction and interference. This is quite complicated, especially because of multiple scattering.

A more practical physical model of photon propagation into diffuse media based on energy flow is the *Radiative Transfer Equation* (RTE). This model was first developed to describe the non-interacting neutrons diffusion in nuclear reactors [Sanchez & McCormick (1982)]. The hypothesis of non-interacting particles reasonably applies to photons in multiple elastic scattering regime. The RTE model neglects polarization effects. This is not a problem though, because even if light coherence exists, this is going to be completely lost after a few scattering events.

#### **2.3.1 The radiative transfer equation**

The average power that at position r and time *t* flows through the unit area oriented in the direction of the unit vector ˆs, due to photons within a unit frequency band centered at *ν*, that are moving within the unit solid angle around ˆs is the *spectral radiance Is* (r, ˆs, *t*, *ν*). Thus at time *t* through the unit area dΣ lying in r, oriented as ˆs, within the unit solid angle dΩ and in the frequency interval [*ν*, *ν* + d*ν*] flows a power d*P* given by:

$$\text{d}P = I\_{\text{s}} \left( r, \hat{\mathbf{s}}, t, \nu \right) \left| \hat{\mathbf{n}} \cdot \hat{\mathbf{s}} \right| \text{d}\Sigma \text{d}\Omega \text{d}\nu \tag{21}$$

Energy density <sup>d</sup>*<sup>E</sup>* <sup>d</sup>*<sup>V</sup>* is simply related to spectral radiance. The energy d*E* per unit frequency and per unit solid angle that crossed the area dΣ oriented as ˆs, in the time interval d*t* is:

$$\frac{\mathbf{d}E}{\mathbf{d}V} = \frac{I\_s \mathbf{d}\Sigma \mathbf{d}t}{\mathbf{d}\Sigma v \mathbf{d}t} = \frac{I\_s}{v} \tag{22}$$

Where *v* is the speed of radiation in the medium, *v* = *<sup>c</sup> <sup>n</sup>* . Energy density is the number of photons and so *Is* is proportional to the number of photons in the unit volume, with frequency *ν*, that at time *t* are moving in the direction ˆs. As we deal with media in which the radiation frequency doesn't change during propagation (elastic scattering), we can integrate *Is* over the frequency *ν* and obtain the *radiance I* (r, ˆs, *t*).

Integrating also over the solid angle and dividing by *v* gives the quantity *u* (r, *t*) = <sup>4</sup>*<sup>π</sup> I*(*r*, ˆ*s*,*t*)dΩ *<sup>v</sup>* , measured in Jm−<sup>3</sup> , that represents the energy density at position r and time *t*. Finally the photon density can be obtained dividing the energy density by the energy of a single photon:

$$m\left(r,t\right) = \frac{u\left(r,t\right)}{h\nu} \tag{23}$$

The RTE is an integro-differential equation stating the balance of the incoming and outgoing radiation along the direction ˆs, at the time *t*, inside the volume element d*V* at the position r

[Ishimaru (1978); Martelli et al. (2010)]:

8 Will-be-set-by-IN-TECH

A different way to deal with diffusion and discriminate between scattering and absorption contributes is to use non-constant light intensities in the injected radiation. A possible way to describe the physics of the problem is to start from the electric and magnetic fields associated with the radiation and to develop a rigorous mathematical theory, taking into account scattering, diffraction and interference. This is quite complicated, especially because

A more practical physical model of photon propagation into diffuse media based on energy flow is the *Radiative Transfer Equation* (RTE). This model was first developed to describe the non-interacting neutrons diffusion in nuclear reactors [Sanchez & McCormick (1982)]. The hypothesis of non-interacting particles reasonably applies to photons in multiple elastic scattering regime. The RTE model neglects polarization effects. This is not a problem though, because even if light coherence exists, this is going to be completely lost after a few scattering

The average power that at position r and time *t* flows through the unit area oriented in the direction of the unit vector ˆs, due to photons within a unit frequency band centered at *ν*, that are moving within the unit solid angle around ˆs is the *spectral radiance Is* (r, ˆs, *t*, *ν*). Thus at time *t* through the unit area dΣ lying in r, oriented as ˆs, within the unit solid angle dΩ and in

and per unit solid angle that crossed the area dΣ oriented as ˆs, in the time interval d*t* is:

<sup>d</sup>*<sup>V</sup>* <sup>=</sup> *Is*dΣd*<sup>t</sup>*

photons and so *Is* is proportional to the number of photons in the unit volume, with frequency *ν*, that at time *t* are moving in the direction ˆs. As we deal with media in which the radiation frequency doesn't change during propagation (elastic scattering), we can integrate *Is* over the

*t*. Finally the photon density can be obtained dividing the energy density by the energy of a

*<sup>n</sup>* (r, *<sup>t</sup>*) <sup>=</sup> *<sup>u</sup>* (r, *<sup>t</sup>*)

The RTE is an integro-differential equation stating the balance of the incoming and outgoing radiation along the direction ˆs, at the time *t*, inside the volume element d*V* at the position r

<sup>d</sup>Σ*v*d*<sup>t</sup>* <sup>=</sup> *Is*

d*E*

d*P* = *Is* (r, ˆs, *t*, *ν*)|nˆ · sˆ|dΣdΩd*ν* (21)

, that represents the energy density at position r and time

*<sup>v</sup>* (22)

*<sup>n</sup>* . Energy density is the number of

*v* gives the quantity *u* (r, *t*) =

*<sup>h</sup><sup>ν</sup>* (23)

<sup>d</sup>*<sup>V</sup>* is simply related to spectral radiance. The energy d*E* per unit frequency

**2.3 Photon migration approach**

**2.3.1 The radiative transfer equation**

the frequency interval [*ν*, *ν* + d*ν*] flows a power d*P* given by:

Where *v* is the speed of radiation in the medium, *v* = *<sup>c</sup>*

Integrating also over the solid angle and dividing by

 Jm−<sup>3</sup> 

frequency *ν* and obtain the *radiance I* (r, ˆs, *t*).

*<sup>v</sup>* , measured in

of multiple scattering.

events.

Energy density <sup>d</sup>*<sup>E</sup>*

<sup>4</sup>*<sup>π</sup> I*(*r*, ˆ*s*,*t*)dΩ

single photon:

$$\frac{1}{v}\frac{\partial I}{\partial t} = -\mathbf{\hat{s}} \cdot \nabla I - \left(\mu\_a + \mu\_s\right)I + \mu\_s \int\_{4\pi} p\left(\mathbf{\hat{s}}, \mathbf{\hat{s}}'\right) I\left(r, \mathbf{\hat{s}}', t\right) d\Omega' + q\tag{24}$$

Each term appearing has a precise significance:


The RTE has two immediate properties:

1. If *I*<sup>0</sup> is a solution for a non-absorbing medium with time impulsive source term *q* (r, ˆs, *t*) = *q*<sup>0</sup> (r, ˆs) *δ* (*t*), then if the absorption coefficient is *μ<sup>a</sup>* the solution is:

$$I = I\_0 \varepsilon^{-\mu\_d v t} \tag{25}$$

This property is a generalization of the Lambert-Beer law [Sassaroli & Fantini (2004); Tsuchiya (2001)].

2. If I is the Green function for a medium having extinction coefficient *μ<sup>t</sup>* = *μ<sup>a</sup>* + *μ<sup>s</sup>* and albedo *a* = *<sup>μ</sup><sup>s</sup> μt* , then for a medium having extinction coefficient *μ*∗ *<sup>t</sup>* and the same albedo, the solution is scaled in this way:

$$I\left(\mathbf{r}^\*, \mathbf{\dot{s}}, t^\*\right) = \left(\frac{\mu\_t^\*}{\mu\_t}\right)^3 I\left(\mathbf{r}, \mathbf{\dot{s}}, t\right) \qquad \begin{cases} \mathbf{r}^\* = \left(\frac{\mu\_t}{\mu\_t^\*}\right) \mathbf{r} \\ t^\* = \left(\frac{\mu\_t}{\mu\_t^\*}\right) t \end{cases} \tag{26}$$

This scaling property is known as the *similarity principle* [Zege et al. (1991)].

Because of the high complexity of the RTE, no analytical solutions are available and approximation methods are then applied. Numerical methods like the *Finite Difference Method* or the *Finite Elements Method* (FEM) require a large amount of computational power, but are often applied [Arridge (1995); Arridge et al. (2000); Arridge & Hebden (1997); Arridge & Schweiger (1995); Arridge et al. (1993)]. Stochastic methods like the *Monte Carlo* are widely used in many biological applications [Boas et al. (2002); Fang & Boas (2009)].

#### **2.3.2 The PN approximation**

Many ways to simplify the RTE have been developed through the years, but most methods are based on the *PN approximation*. Two derived quantities of interest are the *photon fluence* Φ

For the source term approximation the quantities monopole moment *q*<sup>0</sup> and dipole moment

Functional Near Infrared Spectroscopy and Diffuse Optical Tomography in Neuroscience 61

*<sup>g</sup>*<sup>1</sup> <sup>=</sup> <sup>&</sup>lt; <sup>s</sup><sup>ˆ</sup> · <sup>ˆ</sup>

<sup>J</sup> (r, *<sup>t</sup>*) <sup>−</sup> <sup>1</sup>

In a high-albedo (predominantly scattering) medium, the time for a substantial change in flux J is much longer than the time to traverse one transport mean free path. Thus, over one transport mean free path itself, the fractional change in flux is much less than unity. The

*∂*J (r, *t*)

<sup>J</sup> (r, *<sup>t</sup>*) <sup>=</sup> <sup>−</sup> <sup>1</sup>

Usually the factor multiplying the gradient is called *diffusion coefficient*, *D* = <sup>1</sup>

<sup>3</sup> This is consistent with the similarity principle previously introduced [Martelli et al. (2010)].

Dropping the dipole moment of the source is reasonable assuming an isotropic source. Eq. 43

3 (*μ<sup>a</sup>* + *μ*�

*s*)

*q*<sup>0</sup> (r, *t*) = *q*0,0 (r, *t*) (38)

*<sup>∂</sup><sup>t</sup>* <sup>=</sup> <sup>−</sup>*μa*<sup>Φ</sup> (r, *<sup>t</sup>*) −∇· <sup>J</sup> (r, *<sup>t</sup>*) <sup>+</sup> *<sup>q</sup>*<sup>0</sup> (r, *<sup>t</sup>*) (42)

3

⎞

*g*<sup>0</sup> = 1 (40)

s� > (41)

∇Φ (r, *t*) + q<sup>1</sup> (r, *t*) (43)

*s*:

*<sup>s</sup>* = (1 − *g*1) *μ<sup>s</sup>* (44)

*<sup>∂</sup><sup>t</sup>* <sup>=</sup> <sup>0</sup> (45) q<sup>1</sup> (r, *t*) = 0 (46)

*<sup>s</sup>* � *μa*) it doesn't depend on absorption and can be

∇Φ (r, *t*) (47)

3(*μa*+*μ*�

*s*). We

⎟⎠ (39)

<sup>2</sup> (*q*1,−<sup>1</sup> (r, *<sup>t</sup>*) <sup>−</sup> *<sup>q</sup>*1,1 (r, *<sup>t</sup>*))

<sup>√</sup><sup>2</sup> (*q*1,−<sup>1</sup> (r, *<sup>t</sup>*) <sup>+</sup> *<sup>q</sup>*1,1 (r, *<sup>t</sup>*)) *q*1,0 (r, *t*)

q<sup>1</sup> have been introduced:

1 *v*

> 1 *v*

**2.3.3 The diffusion approximation**

then becomes *Fick's law*:

calculated as *<sup>D</sup>* <sup>≈</sup> <sup>1</sup>

notice that in high-albedo medium (*μ*�

3*μ*� *s* 3.

*∂*Φ (r, *t*)

*∂*J (r, *t*)

*<sup>∂</sup><sup>t</sup>* <sup>=</sup> <sup>−</sup> �

q<sup>1</sup> (r, *t*) =

⎛

⎜⎝ √ 1

Normalizing *g*<sup>0</sup> to unity, *g*<sup>1</sup> is the average cosine of the scattering angle:

Finally, substituting these approximated expressions into 24, we obtain:

Where we made use of the definition of reduced scattering coefficient *μ*�

*diffusion approximation* results from making the following assumptions:

*μ<sup>a</sup>* + *μ*� *s* �

*μ*�

1 *i*

and the *photon flux* J:

$$\Phi\left(r,t\right) = \int\_{4\pi} I\left(r,\mathfrak{s},t\right)d\Omega\tag{27}$$

$$J\left(r,t\right) = \int\_{4\pi} I\left(r,\left\|t\right\|\right) \left\|d\Omega\right\|\tag{28}$$

Radiance and source term are expanded in spherical harmonics [Boas (1996); Jackson (1999)]:

$$I\left(r,\left.\dot{\mathbf{s}},t\right)=\sum\_{l=0}^{+\infty}\sum\_{m=-l}^{l}\sqrt{\frac{2l+1}{4\pi}}\phi\_{l,m}\left(r,t\right)\,\mathcal{Y}\_{l,m}\left(\dot{\mathbf{s}}\right)\tag{29}$$

$$q\left(r,\left.\dot{\mathbf{s}},t\right.\right) = \sum\_{l=0}^{+\infty} \sum\_{m=-l}^{l} \sqrt{\frac{2l+1}{4\pi}} q\_{l,m}\left(r,t\right) Y\_{l,m}\left(\dot{\mathbf{s}}\right) \tag{30}$$

Phase function is expanded in spherical harmonics too, but first the reasonable hypothesis that the scattering amplitude is only dependent on the change in direction of the photon is made. With this assumption the probability of scattering from a direction ˆs into a direction ˆs� depends only on the angle between ˆs and ˆs� :

$$p\left(\hat{\mathbf{s}}\cdot\hat{\mathbf{s}}'\right) = \sum\_{l=0}^{+\infty} \frac{2l+1}{4\pi} g\_l P\_l\left(\hat{\mathbf{s}}\cdot\hat{\mathbf{s}}'\right) \tag{31}$$

$$=\sum\_{l=0}^{+\infty}\sum\_{m=-l}^{l}g\_{l}Y\_{l,m}^{\*}\left(\mathring{\mathbf{s}}'\right)Y\_{l,m}\left(\mathring{\mathbf{s}}\right)\tag{32}$$

Where *Pl* is the Legendre polynomial of order *l* and the angular addition rule has been used [Jackson (1999)]. The normalization factors �2*l*+<sup>1</sup> <sup>4</sup>*<sup>π</sup>* and <sup>2</sup>*l*+<sup>1</sup> <sup>4</sup>*<sup>π</sup>* are introduced for convenience.

Truncating the series at the *l* = *N* term brings to the *PN* approximation of the RTE [Arridge (1999)]. In the *P*<sup>1</sup> approximation the radiance2, source and phase function reduce to:

$$I\left(r, \dot{\mathbf{s}}, t\right) = \frac{1}{4\pi} \Phi\left(r, t\right) + \frac{3}{4\pi} J\left(r, t\right) \cdot \dot{\mathbf{s}}\tag{33}$$

$$q\left(r, \dot{\mathbf{s}}, t\right) = \frac{1}{4\pi}q\_0\left(r, t\right) + \frac{3}{4\pi}q\_1\left(r, t\right) \cdot \dot{\mathbf{s}}\tag{34}$$

$$p\left(\hat{\mathbf{s}}\cdot\hat{\mathbf{s}}'\right) = \frac{1}{4\pi}\mathbf{g}\_0 + \frac{3}{4\pi}\mathbf{g}\_1\hat{\mathbf{s}}\cdot\hat{\mathbf{s}}'\tag{35}$$

Where we used the definition of fluence 27 and flux 28, obtaining:

$$
\Phi\left(r, t\right) = \phi\_{0,0}\left(r, t\right) \tag{36}
$$

$$J\left(r,t\right) = \begin{pmatrix} \frac{1}{\sqrt{2}}\left(\phi\_{1,-1}\left(r,t\right) - \phi\_{1,1}\left(r,t\right)\right) \\ \frac{1}{i\sqrt{2}}\left(\phi\_{1,-1}\left(r,t\right) + \phi\_{1,1}\left(r,t\right)\right) \\ \phi\_{1,0}\left(r,t\right) \end{pmatrix} \tag{37}$$

<sup>2</sup> As will be realized soon, this brings to assume the radiance to be nearly isotropic. This is a good approximation for biological tissues and in general for high-albedo media, because of the few absorption events relative to the scattering events.

For the source term approximation the quantities monopole moment *q*<sup>0</sup> and dipole moment q<sup>1</sup> have been introduced:

$$
\rho\_0\left(\mathbf{r},t\right) = \rho\_{0,0}\left(\mathbf{r},t\right) \tag{38}
$$

$$q\_1\left(r,t\right) = \begin{pmatrix} \frac{1}{\sqrt{2}}\left(q\_{1,-1}\left(r,t\right) - q\_{1,1}\left(r,t\right)\right) \\ \frac{1}{i\sqrt{2}}\left(q\_{1,-1}\left(r,t\right) + q\_{1,1}\left(r,t\right)\right) \\ q\_{1,0}\left(r,t\right) \end{pmatrix} \tag{39}$$

Normalizing *g*<sup>0</sup> to unity, *g*<sup>1</sup> is the average cosine of the scattering angle:

$$\mathbf{g\_0} = \mathbf{1} \tag{40}$$

$$
\mathbf{g}\_1 = <\mathbf{\hat{s}} \cdot \mathbf{\hat{s}'}>\tag{41}
$$

Finally, substituting these approximated expressions into 24, we obtain:

$$\frac{1}{v}\frac{\partial\Phi(\mathbf{r},t)}{\partial t} = -\mu\_d\Phi\left(\mathbf{r},t\right) - \nabla \cdot \mathbf{J}\left(\mathbf{r},t\right) + q\_0\left(\mathbf{r},t\right) \tag{42}$$

$$\frac{1}{v}\frac{\partial \mathbf{J}\left(r,t\right)}{\partial t} = -\left(\mu\_d + \mu\_s'\right)\mathbf{J}\left(r,t\right) - \frac{1}{3}\nabla \Phi\left(r,t\right) + q\_\mathbf{1}\left(r,t\right) \tag{43}$$

Where we made use of the definition of reduced scattering coefficient *μ*� *s*:

$$
\mu\_s' = \left(1 - g\_1\right)\mu\_s \tag{44}
$$

#### **2.3.3 The diffusion approximation**

10 Will-be-set-by-IN-TECH

*I* (r, ˆs, *t*) dΩ (27)

*I* (r, ˆs, *t*) sˆdΩ (28)

<sup>4</sup>*<sup>π</sup> <sup>φ</sup>l*,*<sup>m</sup>* (r, *<sup>t</sup>*) *Yl*,*<sup>m</sup>* (sˆ) (29)

<sup>4</sup>*<sup>π</sup> ql*,*<sup>m</sup>* (r, *<sup>t</sup>*) *Yl*,*<sup>m</sup>* (sˆ) (30)

� (31)

*Yl*,*<sup>m</sup>* (sˆ) (32)

<sup>4</sup>*<sup>π</sup>* are introduced for convenience.

J (r, *t*) · sˆ (33)

q<sup>1</sup> (r, *t*) · sˆ (34)

s� (35)

⎟⎠ (37)

� 4*π*

� 4*π*

Radiance and source term are expanded in spherical harmonics [Boas (1996); Jackson (1999)]:

�2*l* + 1

�2*l* + 1

Phase function is expanded in spherical harmonics too, but first the reasonable hypothesis that the scattering amplitude is only dependent on the change in direction of the photon is made. With this assumption the probability of scattering from a direction ˆs into a direction ˆs�

:

2*l* + 1 <sup>4</sup>*<sup>π</sup> glPl*

*l* ∑ *m*=−*l*

Where *Pl* is the Legendre polynomial of order *l* and the angular addition rule has been used

�2*l*+<sup>1</sup>

Truncating the series at the *l* = *N* term brings to the *PN* approximation of the RTE [Arridge

Φ (r, *t*) +

3 <sup>4</sup>*<sup>π</sup> <sup>g</sup>*1s<sup>ˆ</sup> · <sup>ˆ</sup>

<sup>4</sup>*<sup>π</sup> <sup>q</sup>*<sup>0</sup> (r, *<sup>t</sup>*) <sup>+</sup>

<sup>4</sup>*<sup>π</sup> <sup>g</sup>*<sup>0</sup> <sup>+</sup>

(1999)]. In the *P*<sup>1</sup> approximation the radiance2, source and phase function reduce to:

4*π*

� sˆ · sˆ �

<sup>4</sup>*<sup>π</sup>* and <sup>2</sup>*l*+<sup>1</sup>

3 4*π*

> 3 4*π*

<sup>2</sup> (*φ*1,−<sup>1</sup> (r, *<sup>t</sup>*) <sup>−</sup> *<sup>φ</sup>*1,1 (r, *<sup>t</sup>*))

<sup>√</sup><sup>2</sup> (*φ*1,−<sup>1</sup> (r, *<sup>t</sup>*) <sup>+</sup> *<sup>φ</sup>*1,1 (r, *<sup>t</sup>*)) *φ*1,0 (r, *t*)

<sup>2</sup> As will be realized soon, this brings to assume the radiance to be nearly isotropic. This is a good approximation for biological tissues and in general for high-albedo media, because of the few

Φ (r, *t*) = *φ*0,0 (r, *t*) (36)

⎞

*glY*<sup>∗</sup> *l*,*m* � sˆ � �

Φ (r, *t*) =

J (r, *t*) =

+∞ ∑ *l*=0

> +∞ ∑ *l*=0

> > = +∞ ∑ *l*=0

*<sup>I</sup>* (r, ˆs, *<sup>t</sup>*) <sup>=</sup> <sup>1</sup>

*<sup>q</sup>* (r, ˆs, *<sup>t</sup>*) <sup>=</sup> <sup>1</sup>

Where we used the definition of fluence 27 and flux 28, obtaining:

⎛

⎜⎝ √ 1

1 *i*

J (r, *t*) =

absorption events relative to the scattering events.

*p* � sˆ · sˆ � � <sup>=</sup> <sup>1</sup>

*l* ∑ *m*=−*l*

> *l* ∑ *m*=−*l*

*I* (r, ˆs, *t*) =

*q* (r, ˆs, *t*) =

*p* � sˆ · sˆ � � = +∞ ∑ *l*=0

depends only on the angle between ˆs and ˆs�

[Jackson (1999)]. The normalization factors

and the *photon flux* J:

In a high-albedo (predominantly scattering) medium, the time for a substantial change in flux J is much longer than the time to traverse one transport mean free path. Thus, over one transport mean free path itself, the fractional change in flux is much less than unity. The *diffusion approximation* results from making the following assumptions:

$$\frac{\partial \mathbf{J}(r,t)}{\partial t} = 0 \tag{45}$$

$$q\_1\left(r,t\right) = 0\tag{46}$$

Dropping the dipole moment of the source is reasonable assuming an isotropic source. Eq. 43 then becomes *Fick's law*:

$$J\left(r,t\right) = -\frac{1}{\Im\left(\mu\_a + \mu\_s'\right)} \nabla \Phi\left(r,t\right) \tag{47}$$

Usually the factor multiplying the gradient is called *diffusion coefficient*, *D* = <sup>1</sup> 3(*μa*+*μ*� *s*). We notice that in high-albedo medium (*μ*� *<sup>s</sup>* � *μa*) it doesn't depend on absorption and can be calculated as *<sup>D</sup>* <sup>≈</sup> <sup>1</sup> 3*μ*� *s* 3.

<sup>3</sup> This is consistent with the similarity principle previously introduced [Martelli et al. (2010)].

ρ

ρ

Fig. 5. Scheme of reflection geometry TD-NIRS measurement. Injected light travels in the medium and part of it can be collected by a detector at a source-detector distance *ρ*. Knowing the temporal shape of the injected light *I* (*t*), the temporal shape of the collected light *R* (*ρ*, *t*)

successive data to the baseline values, making the assumption that scattering is constant in

In section 2.3.1 the physical model for radiation propagation inside tissues has been introduced. It will be now shown how RTE gives a practically useful tool to perform a kind of

In many biological applications the semi-infinite geometry is assumed and light intensity measurement are conducted in reflectance geometry (see Fig. 5). In this case the Green function of the Diffusion Equation with the Extrapolated Boundary Condition, expressed as

time. This makes concentration variations the only possible obtainable data.

measurement in which absolute hemoglobin concentrations values can be sampled.

I (ω) R (ρ, ω)

Functional Near Infrared Spectroscopy and Diffuse Optical Tomography in Neuroscience 63

μa, μ s

R (ρ, t)

μa, μ

s

*s*.

z

z

O

I (t)

depends on *μ<sup>a</sup>* and *μ*

*s*.

Fig. 4. Scheme of reflection geometry FD-NIRS measurement. Intensity modulated light *I* (*ω*) travels in the medium and part of it can be collected by a detector at a source-detector

distance *ρ*. Amplitude and phase of the collected light *R* (*ρ*, *ω*) depend on *μ<sup>a</sup>* and *μ*

O

Substituting into Eq. 42 we finally have the *diffusion equation* (DE):

$$\frac{1}{v}\frac{\partial\Phi\left(r,t\right)}{\partial t} = -\mu\_d\Phi\left(r,t\right) + \nabla \cdot D\nabla\Phi\left(r,t\right) + q\_0\left(r,t\right) \tag{48}$$

The diffusion equation models many problems in Physics and Chemistry, such as heat propagation, ion diffusion and semiconductor doping. Analytical solution to the DE in photon migration have been carried out in symmetrical geometries such as the infinite medium, the semi-infinite medium or the semi-infinite medium with a high-absorbent inclusion.

In order to solve the DE, appropriate boundary conditions have to be declared. The most general boundary condition is the *Partial Current Boundary Condition* (PCBC):

$$
\Phi\left(r, t\right) + 2AD\nabla\Phi\left(r, t\right) \cdot \hat{n} = 0 \qquad r \in \partial V \tag{49}
$$

Where *A* is a coefficient depending on the refractive index and is due to Fresnel reflection [Contini et al. (1997); Martelli et al. (2010)]. In practical applications, though, the most commonly applied boundary conditions are the *Extrapolated Boundary Condition* (EBC) and the *Zero Boundary Condition* (ZBC):

$$
\Phi(\mathbf{r}, t) = \mathbf{0} \qquad \mathbf{r} \in \partial V\_{\text{ext}} \tag{50}
$$

$$
\Phi\left(r, t\right) = 0 \qquad r \in \partial V \tag{51}
$$

Where *∂Vext* is a surface distant 2*AD* from the boundary (extrapolated surface).

CW, TR and FD solutions for simple geometries have been calculated [Martelli et al. (2010)].

#### **2.3.4 FD-NIRS**

*Frequency Domain Near InfraRed Spectroscopy* (FD-NIRS) makes use of intensity modulated laser light (typically at radio frequencies), injecting it into the sample. The remitted wave is demodulated to obtain amplitude and phase as a function of the modulation frequency. The tissue acts like a low-pass filter: the amplitude is a decreasing function of the frequency and the phase typically increases. The analytical expressions for amplitude and phase can be obtained by a Fourier transform of the TR theoretical expression. Estimation of optical properties can thus be performed using the photon migration theory (see Fig. 4).

Some commercial device using frequency domain techniques is available and measurements and basic research often used this approach to study biological tissues oxygenation [D'Arceuil et al. (2005); Fantini et al. (1995)].

#### **2.3.5 TR-NIRS**

In section 2.2 a way to measure oxygenated and deoxygenated hemoglobin concentrations has been described. Using absorption spectroscopy techniques it is possible to investigate *in vivo* the oxygenation status of the superficial tissues and in particular of the cortical areas of the brain. Thus, CW-NIRS is a useful non-invasive technique for functional neuroimaging. As previously seen, the issue underlying Near InfraRed Spectroscopy is scattering. The way CW-NIRS deals with scattering is to perform a basal measurement and then to refer all 12 Will-be-set-by-IN-TECH

The diffusion equation models many problems in Physics and Chemistry, such as heat propagation, ion diffusion and semiconductor doping. Analytical solution to the DE in photon migration have been carried out in symmetrical geometries such as the infinite medium, the

In order to solve the DE, appropriate boundary conditions have to be declared. The most

Where *A* is a coefficient depending on the refractive index and is due to Fresnel reflection [Contini et al. (1997); Martelli et al. (2010)]. In practical applications, though, the most commonly applied boundary conditions are the *Extrapolated Boundary Condition* (EBC) and

CW, TR and FD solutions for simple geometries have been calculated [Martelli et al. (2010)].

*Frequency Domain Near InfraRed Spectroscopy* (FD-NIRS) makes use of intensity modulated laser light (typically at radio frequencies), injecting it into the sample. The remitted wave is demodulated to obtain amplitude and phase as a function of the modulation frequency. The tissue acts like a low-pass filter: the amplitude is a decreasing function of the frequency and the phase typically increases. The analytical expressions for amplitude and phase can be obtained by a Fourier transform of the TR theoretical expression. Estimation of optical

Some commercial device using frequency domain techniques is available and measurements and basic research often used this approach to study biological tissues oxygenation [D'Arceuil

In section 2.2 a way to measure oxygenated and deoxygenated hemoglobin concentrations has been described. Using absorption spectroscopy techniques it is possible to investigate *in vivo* the oxygenation status of the superficial tissues and in particular of the cortical areas of the brain. Thus, CW-NIRS is a useful non-invasive technique for functional neuroimaging. As previously seen, the issue underlying Near InfraRed Spectroscopy is scattering. The way CW-NIRS deals with scattering is to perform a basal measurement and then to refer all

semi-infinite medium or the semi-infinite medium with a high-absorbent inclusion.

general boundary condition is the *Partial Current Boundary Condition* (PCBC):

Where *∂Vext* is a surface distant 2*AD* from the boundary (extrapolated surface).

properties can thus be performed using the photon migration theory (see Fig. 4).

*<sup>∂</sup><sup>t</sup>* <sup>=</sup> <sup>−</sup>*μa*<sup>Φ</sup> (r, *<sup>t</sup>*) <sup>+</sup> ∇ · *<sup>D</sup>*∇<sup>Φ</sup> (r, *<sup>t</sup>*) <sup>+</sup> *<sup>q</sup>*<sup>0</sup> (r, *<sup>t</sup>*) (48)

Φ (r, *t*) + 2*AD*∇Φ (r, *t*) · nˆ = 0 r ∈ *∂V* (49)

Φ (r, *t*) = 0 r ∈ *∂Vext* (50) Φ (r, *t*) = 0 r ∈ *∂V* (51)

Substituting into Eq. 42 we finally have the *diffusion equation* (DE):

1 *v*

the *Zero Boundary Condition* (ZBC):

et al. (2005); Fantini et al. (1995)].

**2.3.4 FD-NIRS**

**2.3.5 TR-NIRS**

*∂*Φ (r, *t*)

Fig. 4. Scheme of reflection geometry FD-NIRS measurement. Intensity modulated light *I* (*ω*) travels in the medium and part of it can be collected by a detector at a source-detector distance *ρ*. Amplitude and phase of the collected light *R* (*ρ*, *ω*) depend on *μ<sup>a</sup>* and *μ s*.

Fig. 5. Scheme of reflection geometry TD-NIRS measurement. Injected light travels in the medium and part of it can be collected by a detector at a source-detector distance *ρ*. Knowing the temporal shape of the injected light *I* (*t*), the temporal shape of the collected light *R* (*ρ*, *t*) depends on *μ<sup>a</sup>* and *μ s*.

successive data to the baseline values, making the assumption that scattering is constant in time. This makes concentration variations the only possible obtainable data.

In section 2.3.1 the physical model for radiation propagation inside tissues has been introduced. It will be now shown how RTE gives a practically useful tool to perform a kind of measurement in which absolute hemoglobin concentrations values can be sampled.

In many biological applications the semi-infinite geometry is assumed and light intensity measurement are conducted in reflectance geometry (see Fig. 5). In this case the Green function of the Diffusion Equation with the Extrapolated Boundary Condition, expressed as

perturbation models for more complicated geometries, like multi-layered media, have been derived in the last years, but they would require the use of a priori information (e.g., the

Functional Near Infrared Spectroscopy and Diffuse Optical Tomography in Neuroscience 65

(a) Dependance from scattering (*μ<sup>a</sup>* = 0.09 cm−1, *μ*�



(b) Dependance from absorption (*μ<sup>a</sup>* = 0.04, 0.08, 0.16, 0.32 cm−1, *μ*�

Fig. 6. Normalized reflectance from the semi-infinite model, with source-detector distance

*ρ*=2 cm. Reflectance vs. time is plotted in semi-logarithmic scale.

cm−1).

 

 

*<sup>s</sup>* =4, 8, 16, 32 cm−1).

*<sup>s</sup>* = 20

the detected power per unit area (usually known as *reflectance R*), takes the form:

$$R\left(\rho, t\right) = \frac{1}{2\left(4\pi Dv\right)^{3/2}} e^{-\frac{\rho^2}{4Dvt}} e^{-\mu\_\mathrm{tr}\mathrm{tr}} \left[z\_0 e^{-\frac{z\_0^2}{4Dvt}} - z\_\mathrm{p} e^{-\frac{z\_p^2}{4Dvt}}\right] \tag{52}$$

Where *ρ* is the source-detector distance, *z*<sup>0</sup> = <sup>1</sup> *μ*� *s* , *zp* = *z*<sup>0</sup> + 2*ze* and *ze* = 2*AD*.

Looking at the the Green function of the diffusion equation it is possible to notice that absorption and scattering contributes are naturally separated in the model. Plotting these solutions for different values of the absorption coefficient and of the scattering coefficient evidence this property (see Fig. 6). This technique based on the injection of light with a known temporal shape (typical a pulse-like shape) and the estimation of the optical properties from the emitted light temporal shape is called *Time-Resolved Near InfraRed Spectroscopy* (TR-NIRS).

The light is collected by means of photomultipliers tubes (PMT) or avalanche photodiodes (APD) usually in Time Correlated Single Photon Counting (TCSPC) regime. Detailed descriptions of the way light is collected by light detectors is given in [Donati (1998)] and of TCSPC is given in [O'Connor & Phillips (1984)] and [Becker (2005)].

Estimation of the optical properties can be performed from all the analytical expressions for the TR response of a diffusive medium. Absorption and scattering coefficients can be computed by means of an inversion algorithm [Press et al. (1992)]. The instrumental response function (IRF) due to temporal dispersion in fiber, temporal jitter of detectors and finite pulse width of light sources has to be taken into account. Due to the linearity and time-invariance of the transport problem, the detected response is the convolution between the IRF and the theoretical response of the medium. Thus, a Levenberg-Marquardt non-linear minimization between experimental data and theoretical curve convoluted with the IRF is generally performed.

An estimation of the optical properties of a diffusive medium can be also performed by means of Monte Carlo simulations. A certain number of simulations are needed, then by using the scaling properties of the RTE it is possible to implement a fast search of the correct simulation [Pifferi et al. (1998)]. These two method enable an absolute estimation of optical properties.

A simple estimation of the variation of the absorption coefficient can be easily performed if only absorption changes are assumed (and usually are). If we collect at different times two reflectance curves *R*<sup>0</sup> (*ρ*, *t*) and *R*<sup>1</sup> (*ρ*, *t*), it is straightforward that:

$$
\Delta\mu\_a = -\frac{1}{vt}\ln\left(\frac{R\_1\left(\rho\_\prime t\right)}{R\_0\left(\rho\_\prime t\right)}\right) \tag{53}
$$

In principle TR techniques provides a richer insight than CW to the problem of non-invasively probing of a diffusive medium. These approaches can discriminate between absorption and scattering contributions and derive absolute values for the hemodynamic parameters [Patterson et al. (1989)]. This however can be obtained only in simple geometries like the infinite or the semi-infinite homogeneous models. In a real heterogeneous medium like the human tissues and in particular the human head it is easier to derive changes with respect to a baseline or effective average parameters rather than absolute values. Advanced time-resolved 14 Will-be-set-by-IN-TECH

*μ*� *s*

Looking at the the Green function of the diffusion equation it is possible to notice that absorption and scattering contributes are naturally separated in the model. Plotting these solutions for different values of the absorption coefficient and of the scattering coefficient evidence this property (see Fig. 6). This technique based on the injection of light with a known temporal shape (typical a pulse-like shape) and the estimation of the optical properties from the emitted light temporal shape is called *Time-Resolved Near InfraRed Spectroscopy* (TR-NIRS). The light is collected by means of photomultipliers tubes (PMT) or avalanche photodiodes (APD) usually in Time Correlated Single Photon Counting (TCSPC) regime. Detailed descriptions of the way light is collected by light detectors is given in [Donati (1998)] and

Estimation of the optical properties can be performed from all the analytical expressions for the TR response of a diffusive medium. Absorption and scattering coefficients can be computed by means of an inversion algorithm [Press et al. (1992)]. The instrumental response function (IRF) due to temporal dispersion in fiber, temporal jitter of detectors and finite pulse width of light sources has to be taken into account. Due to the linearity and time-invariance of the transport problem, the detected response is the convolution between the IRF and the theoretical response of the medium. Thus, a Levenberg-Marquardt non-linear minimization between experimental data and theoretical curve convoluted with the IRF is

An estimation of the optical properties of a diffusive medium can be also performed by means of Monte Carlo simulations. A certain number of simulations are needed, then by using the scaling properties of the RTE it is possible to implement a fast search of the correct simulation [Pifferi et al. (1998)]. These two method enable an absolute estimation of optical properties. A simple estimation of the variation of the absorption coefficient can be easily performed if only absorption changes are assumed (and usually are). If we collect at different times two

− *z*2 *p* 4*Dvt* 

, *zp* = *z*<sup>0</sup> + 2*ze* and *ze* = 2*AD*.

(52)

(53)

the detected power per unit area (usually known as *reflectance R*), takes the form:

3/2 *t*5/2

*e* <sup>−</sup> *<sup>ρ</sup>*<sup>2</sup> <sup>4</sup>*Dvt e* −*μavt z*0*e* <sup>−</sup> *<sup>z</sup>*<sup>2</sup> <sup>0</sup> <sup>4</sup>*Dvt* <sup>−</sup> *zpe*

*<sup>R</sup>* (*ρ*, *<sup>t</sup>*) <sup>=</sup> <sup>1</sup>

Where *ρ* is the source-detector distance, *z*<sup>0</sup> = <sup>1</sup>

generally performed.

2 (4*πDv*)

of TCSPC is given in [O'Connor & Phillips (1984)] and [Becker (2005)].

reflectance curves *R*<sup>0</sup> (*ρ*, *t*) and *R*<sup>1</sup> (*ρ*, *t*), it is straightforward that:

<sup>Δ</sup>*μ<sup>a</sup>* <sup>=</sup> <sup>−</sup> <sup>1</sup>

*vt* ln

In principle TR techniques provides a richer insight than CW to the problem of non-invasively probing of a diffusive medium. These approaches can discriminate between absorption and scattering contributions and derive absolute values for the hemodynamic parameters [Patterson et al. (1989)]. This however can be obtained only in simple geometries like the infinite or the semi-infinite homogeneous models. In a real heterogeneous medium like the human tissues and in particular the human head it is easier to derive changes with respect to a baseline or effective average parameters rather than absolute values. Advanced time-resolved

 *R*<sup>1</sup> (*ρ*, *t*) *R*<sup>0</sup> (*ρ*, *t*) perturbation models for more complicated geometries, like multi-layered media, have been derived in the last years, but they would require the use of a priori information (e.g., the

(a) Dependance from scattering (*μ<sup>a</sup>* = 0.09 cm−1, *μ*� *<sup>s</sup>* =4, 8, 16, 32 cm−1).

(b) Dependance from absorption (*μ<sup>a</sup>* = 0.04, 0.08, 0.16, 0.32 cm−1, *μ*� *<sup>s</sup>* = 20 cm−1).

Fig. 6. Normalized reflectance from the semi-infinite model, with source-detector distance *ρ*=2 cm. Reflectance vs. time is plotted in semi-logarithmic scale.

Fig. 8. Oxyhemoglobin (in red) and deoxyhemoglobin (in blue) concentrations for a hand grasping task. Single-subject, right hand movement, average of 10 repetitions, TR-NIRS

Functional Near Infrared Spectroscopy and Diffuse Optical Tomography in Neuroscience 67

To properly investigate a large area (larger than the common source-detector distance, that is 1-4 cm) a multi-channel approach is needed. The starting point is the arrangement of a number of sources and detectors to cover the area of interest and the management of the source-detector pairs. Depending on the physical dimensions of the area even a huge number could be arranged. No theoretical limitations on the number of optodes exists, however, actual technology features instrumentations with up to 32 sources and 32 detectors (CW) and 16 sources and 16 detection channels (TR). Shining NIR light on the skin and collecting the back-scattered light in a reflectance geometry from multiple points easily allows to draw a map of the hemoglobin concentrations around the area of interest. Spatial resolution is poor (depending on source-detector distances) and a little depth sensitivity (associated to the

Concentrations time series can be plotted vs time or a spatial map can be built (see Fig. 9). A movie-like map can be obtained if each frame corresponds to a time sample. More often, especially in books and publications, averaged concentrations data over time intervals are used to obtain maps relative to a time block of seconds or minutes, just as the example in Fig

device, contrast enhanced for the deep layers (see section 3.1).

overlapping measurements in CW or intrinsic in TR) is obtained.

9.

anatomy of the head as provided by a MRI scan) for their practical and effective use [Martelli et al. (2005)].

The actual potentiality of time-resolved techniques relies on an easier approach to the problem of depth sensitivity. As mentioned before, to enhance depth sensitivity, CW systems use large source-detector distance and multi-distance approaches. Conversely, depth sensitivity in TR-NIRS can be intrinsic in the measured reflectance curve having information on the time of flight of the photons in tissues.

Near InfraRed Spectroscopy applications to the brain spread from rehabilitation monitoring to neural plasticity studies, from localization of cortical activation in motor or cognitive tasks to stroke diagnosis, and more. Motor tasks such as finger tapping or hand grasping are easy to perform and hundreds of studies have been published [Torricelli et al. (2007)]. Cognitive tasks present more difficulties in activation interpretation, nevertheless a number of study is available in literature [Bandettini et al. (1997); Butti et al. (2009); Heekeren et al. (1997)]. Fig. 7 and Fig. 8 show two examples of TR-NIRS data.

Fig. 7. Δ [Hb] (blue) and Δ [HbO2] (red) concentrations (in arbitrary units) vs. time (in minutes) collected from the frontal area of the brain during a cognitive task of sustained attention. Vertical lines indicated the beginning and the end of the task. The average concentration value for the first two minutes block was used as baseline and then subtracted to all data.

#### **3. Diffuse Optical Tomography (DOT)**

Since the very beginning of the history of NIRS, many efforts have been made to improve space resolution and data accuracy. Near InfraRed Spectroscopy provides functional information about the oxygenation status of the explored tissues, but a single source-detector pair is only able to probe a small underlying area. In this section a way to perform functional neuroimaging starting form NIRS data will be introduced. Production of two dimensional maps, often called *topography*, by linear interpolation of NIRS data was the first application to be investigated. This imaging technique is usually called *Diffuse Optical Imaging* (DOI) and, adding depth sensitivity, is known as *Diffuse Optical Tomography* (DOT).

16 Will-be-set-by-IN-TECH

anatomy of the head as provided by a MRI scan) for their practical and effective use [Martelli

The actual potentiality of time-resolved techniques relies on an easier approach to the problem of depth sensitivity. As mentioned before, to enhance depth sensitivity, CW systems use large source-detector distance and multi-distance approaches. Conversely, depth sensitivity in TR-NIRS can be intrinsic in the measured reflectance curve having information on the time

Near InfraRed Spectroscopy applications to the brain spread from rehabilitation monitoring to neural plasticity studies, from localization of cortical activation in motor or cognitive tasks to stroke diagnosis, and more. Motor tasks such as finger tapping or hand grasping are easy to perform and hundreds of studies have been published [Torricelli et al. (2007)]. Cognitive tasks present more difficulties in activation interpretation, nevertheless a number of study is available in literature [Bandettini et al. (1997); Butti et al. (2009); Heekeren et al. (1997)]. Fig. 7

Fig. 7. Δ [Hb] (blue) and Δ [HbO2] (red) concentrations (in arbitrary units) vs. time (in minutes) collected from the frontal area of the brain during a cognitive task of sustained attention. Vertical lines indicated the beginning and the end of the task. The average

concentration value for the first two minutes block was used as baseline and then subtracted

Since the very beginning of the history of NIRS, many efforts have been made to improve space resolution and data accuracy. Near InfraRed Spectroscopy provides functional information about the oxygenation status of the explored tissues, but a single source-detector pair is only able to probe a small underlying area. In this section a way to perform functional neuroimaging starting form NIRS data will be introduced. Production of two dimensional maps, often called *topography*, by linear interpolation of NIRS data was the first application to be investigated. This imaging technique is usually called *Diffuse Optical Imaging* (DOI) and,

adding depth sensitivity, is known as *Diffuse Optical Tomography* (DOT).

et al. (2005)].

to all data.

of flight of the photons in tissues.

and Fig. 8 show two examples of TR-NIRS data.

**3. Diffuse Optical Tomography (DOT)**

Fig. 8. Oxyhemoglobin (in red) and deoxyhemoglobin (in blue) concentrations for a hand grasping task. Single-subject, right hand movement, average of 10 repetitions, TR-NIRS device, contrast enhanced for the deep layers (see section 3.1).

To properly investigate a large area (larger than the common source-detector distance, that is 1-4 cm) a multi-channel approach is needed. The starting point is the arrangement of a number of sources and detectors to cover the area of interest and the management of the source-detector pairs. Depending on the physical dimensions of the area even a huge number could be arranged. No theoretical limitations on the number of optodes exists, however, actual technology features instrumentations with up to 32 sources and 32 detectors (CW) and 16 sources and 16 detection channels (TR). Shining NIR light on the skin and collecting the back-scattered light in a reflectance geometry from multiple points easily allows to draw a map of the hemoglobin concentrations around the area of interest. Spatial resolution is poor (depending on source-detector distances) and a little depth sensitivity (associated to the overlapping measurements in CW or intrinsic in TR) is obtained.

Concentrations time series can be plotted vs time or a spatial map can be built (see Fig. 9). A movie-like map can be obtained if each frame corresponds to a time sample. More often, especially in books and publications, averaged concentrations data over time intervals are used to obtain maps relative to a time block of seconds or minutes, just as the example in Fig 9.

Time-resolved curves statistically contain information about tissues which light photons pass through. Photons reaching the detector at early times surely have been back-scattered from the superficial tissues and thus can bring information about the superficial layers only. Photons reaching the detector at late times have travelled inside the tissues for a longer time and statistically could have probed deeper layers and carry information about them. We usually

Functional Near Infrared Spectroscopy and Diffuse Optical Tomography in Neuroscience 69

A simple approach to the depth probing problem could be to develop a contrast function,

In order to develop models to discriminate between the variations of the absorption coefficient

time-dependent Path Length (MPL) is introduced [Steinbrink et al. (2001)]. As its own name suggests, MPL is the average length of the path travelled by photons in a specific layer and

The reflectance curve is divided into time intervals *τ*, called *time gates* (TG), and the total

Being *I*<sup>0</sup> (*t*) the intensity for the non-absorbent medium, *LUP* (*t*) and *LDOWN* (*t*) the mean time-dependent path lengths in the superficial layers and in the deep layers, respectively, the

In the simplest case, it is possible to think that early photons carry information only about the superficial layers and that late photons are not affected by a superficial inhomogeneity, and carry information only about the deep layers. Thus, we can write for the absorption coefficient

*LUP* (*τe*)

*LDOWN* (*τl*)

where *τ<sup>e</sup>* is the mean time of flight of photons in a early time gate and *τ<sup>l</sup>* is the mean time of flight for photons in a late time gate. Despite giving a general idea about how to obtain information about contributes of different layers, this model is pretty rough, indeed. The hypothesis of late photons not affected by superficial layers has to be rejected and more complicated expressions are often necessary to perform depth selection [Contini (2007)]. It can be proved that late photons carry both information about superficial and deep layers and that superficial layer absorption variations affect the deep layers absorption estimation. The subtraction of the superficial contribute is thus necessary in the estimation of Δ*μDOWN*

*<sup>a</sup> LUP*(*t*)+Δ*μDOWN*

ln *<sup>I</sup>* (*τe*) *I*<sup>0</sup> (*τe*)

> ln *<sup>I</sup>* (*τl*) *I*<sup>0</sup> (*τl*)

*tin R* (*t*) d*t*

expression for the intensity *I* (*t*) in a general absorbent medium can be expressed as:

<sup>−</sup>(Δ*μUP*

*<sup>a</sup>* ), the quantity time-dependent photon path length, usually called Mean

*<sup>a</sup>* ) and of the variation of the absorption coefficient in deep

*<sup>a</sup> <sup>L</sup>DOWN* (*t*)) (54)

(55)

(56)

*<sup>a</sup>* . Such

refer to them as *early photons* and *late photons*.

in the superficial layers (Δ*μUP*

layers (Δ*μDOWN*

changes:

considering just two domains: superficial layers and deep layers.

can be calculated from the mean time of flight of the photons.

*I* (*t*) = *I*<sup>0</sup> (*t*)*e*

Δ*μUP*

<sup>4</sup> The total counts are the time integral of the reflectance curve: *Igate* (*t*) <sup>=</sup> *tfin*

Δ*μDOWN*

*<sup>a</sup>* <sup>=</sup> <sup>−</sup> <sup>1</sup>

*<sup>a</sup>* <sup>=</sup> <sup>−</sup> <sup>1</sup>

counts (usually referred as *intensity*) are calculated for each time gate4.

Fig. 9. The concentration time series can be visualized as spatial map. Each time sample belongs to an average position between the associated source-detector pair. The map is built by linear interpolation with the neighboring points concentration values.

More complicated ways to image optical data have been developed. Using MRI anatomical information and using suitable values for the average optical properties of the tissues a inverse problem can be established. Localized alterations of the optical properties, corresponding to a increased absorption and related to the hemodynamic responses of neural activations, can be put in correspondence with cortical features and a better localization of brain activations can be performed.

DOT systems can consist of little more than a probe with fiber-optic sources and detectors, a piece of dedicated hardware about the size of a small suitcase and a laptop computer. Systems can be much larger, depending primarily on the type of laser source and detectors employed, but the approach generally offers a degree of portability unobtainable with many other modalities. For this reason, looking into the future, DOT may be ideally suited for clinical applications such as bedside monitoring of cerebral oxygenation.

#### **3.1 Time resolved DOT**

A fundamental point in NIRS measurements, whatever technique we are using, is the ability to separate systemic hemodynamic changes occurring in the superficial tissues, such as the skin, from functional hemodynamic changes related to brain activations. In order to reach this goal is fundamental to discuss NIRS depth sensitivity. Depth sensitivity is not intrinsic in CW- and TR-NIRS measurements, but can be achieved using proper optodes configurations or using time of flight information. Multi-distance source-detector approaches, both with CW [Saager & Berger (2005)] and with TR [Liebert et al. (2004)] techniques, have been proposed to improve depth selectivity and sensitivity. Single-distance approaches have also been discussed [Selb et al. (2005); Steinbrink et al. (2001)]. Contini et al. (2007) proposed a different approach to add depth information to the tissue probing problem, based on time-domain contrast functions. Wabnitz et al. (2008) discussed a method for depth selectivity based on time windows and moments of time-of-flight distributions for TD-NIRS.

18 Will-be-set-by-IN-TECH

Fig. 9. The concentration time series can be visualized as spatial map. Each time sample belongs to an average position between the associated source-detector pair. The map is built

More complicated ways to image optical data have been developed. Using MRI anatomical information and using suitable values for the average optical properties of the tissues a inverse problem can be established. Localized alterations of the optical properties, corresponding to a increased absorption and related to the hemodynamic responses of neural activations, can be put in correspondence with cortical features and a better localization of brain activations can

DOT systems can consist of little more than a probe with fiber-optic sources and detectors, a piece of dedicated hardware about the size of a small suitcase and a laptop computer. Systems can be much larger, depending primarily on the type of laser source and detectors employed, but the approach generally offers a degree of portability unobtainable with many other modalities. For this reason, looking into the future, DOT may be ideally suited for

A fundamental point in NIRS measurements, whatever technique we are using, is the ability to separate systemic hemodynamic changes occurring in the superficial tissues, such as the skin, from functional hemodynamic changes related to brain activations. In order to reach this goal is fundamental to discuss NIRS depth sensitivity. Depth sensitivity is not intrinsic in CW- and TR-NIRS measurements, but can be achieved using proper optodes configurations or using time of flight information. Multi-distance source-detector approaches, both with CW [Saager & Berger (2005)] and with TR [Liebert et al. (2004)] techniques, have been proposed to improve depth selectivity and sensitivity. Single-distance approaches have also been discussed [Selb et al. (2005); Steinbrink et al. (2001)]. Contini et al. (2007) proposed a different approach to add depth information to the tissue probing problem, based on time-domain contrast functions. Wabnitz et al. (2008) discussed a method for depth selectivity based on time windows and

by linear interpolation with the neighboring points concentration values.

clinical applications such as bedside monitoring of cerebral oxygenation.

moments of time-of-flight distributions for TD-NIRS.

be performed.

**3.1 Time resolved DOT**

Time-resolved curves statistically contain information about tissues which light photons pass through. Photons reaching the detector at early times surely have been back-scattered from the superficial tissues and thus can bring information about the superficial layers only. Photons reaching the detector at late times have travelled inside the tissues for a longer time and statistically could have probed deeper layers and carry information about them. We usually refer to them as *early photons* and *late photons*.

A simple approach to the depth probing problem could be to develop a contrast function, considering just two domains: superficial layers and deep layers.

In order to develop models to discriminate between the variations of the absorption coefficient in the superficial layers (Δ*μUP <sup>a</sup>* ) and of the variation of the absorption coefficient in deep layers (Δ*μDOWN <sup>a</sup>* ), the quantity time-dependent photon path length, usually called Mean time-dependent Path Length (MPL) is introduced [Steinbrink et al. (2001)]. As its own name suggests, MPL is the average length of the path travelled by photons in a specific layer and can be calculated from the mean time of flight of the photons.

The reflectance curve is divided into time intervals *τ*, called *time gates* (TG), and the total counts (usually referred as *intensity*) are calculated for each time gate4.

Being *I*<sup>0</sup> (*t*) the intensity for the non-absorbent medium, *LUP* (*t*) and *LDOWN* (*t*) the mean time-dependent path lengths in the superficial layers and in the deep layers, respectively, the expression for the intensity *I* (*t*) in a general absorbent medium can be expressed as:

$$I\left(t\right) = I\_0\left(t\right)e^{-\left(\Delta\mu\_d^{\rm IP}L^{\rm IP}(t) + \Delta\mu\_d^{\rm COWN}L^{\rm COWN}(t)\right)}\tag{54}$$

In the simplest case, it is possible to think that early photons carry information only about the superficial layers and that late photons are not affected by a superficial inhomogeneity, and carry information only about the deep layers. Thus, we can write for the absorption coefficient changes:

$$
\Delta \mu\_a^{IIP} = -\frac{1}{L^{LP}\left(\tau\_\ell\right)} \ln \left(\frac{I\left(\tau\_\ell\right)}{I\_0\left(\tau\_\ell\right)}\right) \tag{55}
$$

$$
\Delta\mu\_a^{DOWN} = -\frac{1}{L^{DOWN}\left(\tau\_l\right)}\ln\left(\frac{I\left(\tau\_l\right)}{I\_0\left(\tau\_l\right)}\right) \tag{56}
$$

where *τ<sup>e</sup>* is the mean time of flight of photons in a early time gate and *τ<sup>l</sup>* is the mean time of flight for photons in a late time gate. Despite giving a general idea about how to obtain information about contributes of different layers, this model is pretty rough, indeed. The hypothesis of late photons not affected by superficial layers has to be rejected and more complicated expressions are often necessary to perform depth selection [Contini (2007)]. It can be proved that late photons carry both information about superficial and deep layers and that superficial layer absorption variations affect the deep layers absorption estimation. The subtraction of the superficial contribute is thus necessary in the estimation of Δ*μDOWN <sup>a</sup>* . Such

<sup>4</sup> The total counts are the time integral of the reflectance curve: *Igate* (*t*) <sup>=</sup> *tfin tin R* (*t*) d*t*

Moreover, the comparison with standard clinical techniques such EEG or MEG can lead to a clinical standardization of the NIRS signal and push Near InfraRed Spectroscopy towards a

Functional Near Infrared Spectroscopy and Diffuse Optical Tomography in Neuroscience 71

Using anatomical magnetic resonance priors to perform MRI-guided optical reconstruction dramatically improves the spatial resolution of the diffuse optical tomography techniques [Boas & Dale (2005)]. Simulations of light propagation are run into the MRI head volume, using Monte Carlo statistical methods such as Fang & Boas (2009), and then the cortical activation profile is obtained solving an inverse problem. For a more detailed description

In the event of an unavailability of the subject-specific anatomy the efficient use of an MRI anatomical atlas has been demonstrated [Caffini et al. (2010)]. Fig. 11 reports an atlas reconstruction of a cortical activation during a visual protocol of pattern reversal

Fig. 11. Atlas map visualization of an MRI-guided optical reconstruction of the brain activation in the visual cortex (the brain is seen from the back), during a pattern reversal

through the last two decades and big steps have been done.

quantum efficiency in the near red region.

Near InfraRed Spectroscopy applied *in vivo* to cortical tissues has been widely investigated

From the technical point of view, the laser technology, and in particular the introduction of semiconductor laser diodes, has helped to fabricate compact and clinical instrumentations. Moreover, the improvements in light detectors has permitted to develop accurate devices, for example, in the time-resolved technology, the red-extended photocathodes well increased the

Nevertheless, a lot more has to be done. Pulsed lasers sources are getting better and better, especially concerning time stability and output power. Photomultiplier tubes is an efficient and well established technology, but, in NIRS application, suffers the need of high voltage

regular clinical use [*nEUROPt Project* (2008-2012)].

of it see Caffini et al. (2011).

checkerboard.

checkerboard protocol.

**4. Conclusions**

a subtraction can be performed as follow:

$$
\Delta \mu\_a^{\rm LP} = -\frac{1}{L^{\rm LP}\left(\tau\_\mathcal{e}\right)} \ln \left(\frac{I\left(\tau\_\mathcal{e}\right)}{I\_0\left(\tau\_\mathcal{e}\right)}\right) \tag{57}
$$

$$
\Delta\mu\_{\rm a}^{DOWN} = -\frac{1}{L^{DOWN}\left(\tau\_{l}\right)}\ln\left(1 + \frac{I\left(\tau\_{l}\right)}{I\_{0}\left(\tau\_{l}\right)} - \frac{I\left(\tau\_{\rm e}\right)}{I\_{0}\left(\tau\_{\rm e}\right)}\right) \tag{58}
$$

Oxy- and deoxy-hemoglobin concentrations for the upper layer and the lower layer can be easily obtained via Lambert-Beer law, from the estimated Δ*μa*. An example of application of this stratigraphic model is shown in Fig. 10, where data where collected during a Valsalva maneuver, that forces an increase of the systemic blood volume, resulting in an increase of both oxyhemoglobin and deoxyhemoglobin concentrations in the superficial layers of the skin.

Fig. 10. Estimated Δ*HbO*<sup>2</sup> and Δ*Hb* in the superficial layer (UP, continuous line) and in the brain (DOWN, dotted line) during the Valsalva maneuver without normalization. The two vertical dashed lines indicate the beginning and the end of the task period, respectively. During the Valsalva maneuver the sistemic blood volume dramatically increases, while the local blood volume in the brain is less affected by the maneuver effects. A proper depth selectivity remarks this differences.

#### **3.2 Multimodality approach**

Recently a great interest in multimodality approaches has grown. Merging the advantages of different imaging and functional techniques, such as co-registration of NIRS with MRI [Merritt et al. (2002)], blood flow monitors, fMRI [Torricelli et al. (2007)] and PET, gives the possibility to build anatomo-functional images and movies to largely improve information visualization. Moreover, the comparison with standard clinical techniques such EEG or MEG can lead to a clinical standardization of the NIRS signal and push Near InfraRed Spectroscopy towards a regular clinical use [*nEUROPt Project* (2008-2012)].

Using anatomical magnetic resonance priors to perform MRI-guided optical reconstruction dramatically improves the spatial resolution of the diffuse optical tomography techniques [Boas & Dale (2005)]. Simulations of light propagation are run into the MRI head volume, using Monte Carlo statistical methods such as Fang & Boas (2009), and then the cortical activation profile is obtained solving an inverse problem. For a more detailed description of it see Caffini et al. (2011).

In the event of an unavailability of the subject-specific anatomy the efficient use of an MRI anatomical atlas has been demonstrated [Caffini et al. (2010)]. Fig. 11 reports an atlas reconstruction of a cortical activation during a visual protocol of pattern reversal checkerboard.

Fig. 11. Atlas map visualization of an MRI-guided optical reconstruction of the brain activation in the visual cortex (the brain is seen from the back), during a pattern reversal checkerboard protocol.
